Introduction to BrainIAK tutorials

Congratulations, if you are viewing this Jupyter notebook, you have already acquired many of the skills necessary to excel in this course and you are well on your way to learning cutting-edge methods for cognitive neuroscience!

For users on NeuroLibre, you are seeing a ready to run version of these tutorials that needed no installation or configuration on your part. If you would like to install and use these tutorials on your own machines, follow the instructions on brainiak tutorials.

In this course we will use a variety of tools, many of which will likely be new to you. Don't worry if you are having trouble wrapping your head around them now: by the end of this course you will be proficient in not only these useful skills but also the exciting analyses that use them.

Goal of this notebook

1. Familiarize yourself with the tools that will be used in these notebooks. 

Table of Contents

Exercises

Exercise 1

Contributions

Resources

Here are some resources (Python, fMRI and machine learning, etc.): BrainIAK tutorials resource page

Import necessary packages

While importing packages, you may see warning messages. It is safe to ignore these warnings as they will not impact your execution of the tutorials.

# suppress warnings
import warnings
import sys 
if not sys.warnoptions:
    warnings.simplefilter("ignore")

# The plotting tool we will be using in this course
import matplotlib.pyplot as plt

# Module essential for data organization and manipulation
import numpy as np #numpy's "nickname" is np

# Import a function from BrainIAK to simulate fMRI data
import brainiak.utils.fmrisim as sim  

# display the plots inline 
%matplotlib inline 
# autosave for every 5 secs
%autosave 5

Autosaving every 5 seconds

Brain template

We are now going to use some of the tools we just loaded. First we'll call a function from brainiak to load a gray matter mask from the MNI152 standard brain. Here's an article talking about different anatomical standards, including MNI152: Structural Brain Atlases: Design, Rationale, and Applications in Normal and Pathological Cohorts

# Set the size (in terms of X, Y, Z) of the volume we want to create
dimensions = np.asarray([64, 64, 64])

# Generate an anatomical image with the size above of brain voxels in gray matter
# This outputs variables for two versions of the image, binary (mask) and probabilistic (template)
mask, template = sim.mask_brain(dimensions, mask_self=False)

Congrats, you just ran a command from BrainIAK!!

We are now going to take a slice from that template and display it.

# Get an axial (a.k.a. transverse or horizontal) slice halfway through the brain
mid_idx = dimensions[2] // 2
axial_slice = template[:, :, mid_idx]

# imshow can visualize a 2d array 
plt.imshow(axial_slice)
plt.title('An axial brain slice');

There are great tools in Python for exploring brains in notebooks. One such tool is niwidgets. Below we use that tool to look at the brain interactively. If you cannot install it, there are other options to consider, like nibabel.viewers.OrthoSlicer3D

try:
    # Load in the new variable
    from niwidgets import NiftiWidget

    template_nii = nib.Nifti1Image(template, np.eye(4))
    viewer = NiftiWidget(template_nii)
    viewer.nifti_plotter();

except:
    print('niwidgets cannot run, try installing it or some other viewing tool')

niwidgets cannot run, try installing it or some other viewing tool

"help()"

help is a very useful function in Python. If you type help(function_name) in Python, you will get some basic infomation about how to use this function. If you run the following line, you will see that sim.mask_brain takes the dimension of x, y, and z, and then output a MNI152 template with the specified dimensions. Note, you can also do this by typing [SHIFT] + [TAB] while the cursor is hovering over a function name.

Note: The [SHIFT] + [TAB] works in Jupyter environments, but you will see small differences in this functionality when these notebooks are are used in other environments such as NeuroLibre that use Binder.

help(sim.mask_brain)

Help on function mask_brain in module brainiak.utils.fmrisim:

mask_brain(volume, template_name=None, mask_threshold=None, mask_self=True)
    Mask the simulated volume
    This creates a mask specifying the approximate likelihood that a voxel is
    part of the brain. All values are bounded to the range of 0 to 1. An
    appropriate threshold to isolate brain voxels is >0.2. Critically,
    the data that should be used to create a template shouldn't already be
    masked/skull stripped. If it is then it will give in accurate estimates
    of non-brain noise and corrupt estimations of SNR.
    
    Parameters
    ----------
    
    volume : multidimensional array
        Either numpy array of a volume or a tuple describing the dimensions
        of the mask to be created
    
    template_name : str
        What is the path to the template to be loaded? If empty then it
        defaults to an MNI152 grey matter mask. This is ignored if mask_self
        is True.
    
    mask_threshold : float
        What is the threshold (0 -> 1) for including a voxel in the mask? If
        None then the program will try and identify the last wide peak in a
        histogram of the template (assumed to be the brain voxels) and takes
        the minima before that peak as the threshold. Won't work when the
        data is not bimodal.
    
    mask_self : bool or None
        If set to true then it makes a mask from the volume supplied (by
        averaging across time points and changing the range). If it is set
        to false then it will use the template_name as an input.
    
    Returns
    ----------
    
    mask : 3 dimensional array, binary
        The masked brain, thresholded to distinguish brain and non-brain
    
    template : 3 dimensional array, float
        A continuous (0 -> 1) volume describing the likelihood a voxel is in
        the brain. This can be used to contrast the brain and non brain.

Look at the source code

If you want to see the source code, you can use the getsource function from the inspect package.

Run the following code to see the source code of sim.mask_brain.

import inspect # this "inspect" package can let you peek what's inside a function
source_code = inspect.getsource(sim.mask_brain)
print(source_code)

def mask_brain(volume,
               template_name=None,
               mask_threshold=None,
               mask_self=True,
               ):
    """ Mask the simulated volume
    This creates a mask specifying the approximate likelihood that a voxel is
    part of the brain. All values are bounded to the range of 0 to 1. An
    appropriate threshold to isolate brain voxels is >0.2. Critically,
    the data that should be used to create a template shouldn't already be
    masked/skull stripped. If it is then it will give in accurate estimates
    of non-brain noise and corrupt estimations of SNR.

    Parameters
    ----------

    volume : multidimensional array
        Either numpy array of a volume or a tuple describing the dimensions
        of the mask to be created

    template_name : str
        What is the path to the template to be loaded? If empty then it
        defaults to an MNI152 grey matter mask. This is ignored if mask_self
        is True.

    mask_threshold : float
        What is the threshold (0 -> 1) for including a voxel in the mask? If
        None then the program will try and identify the last wide peak in a
        histogram of the template (assumed to be the brain voxels) and takes
        the minima before that peak as the threshold. Won't work when the
        data is not bimodal.

    mask_self : bool or None
        If set to true then it makes a mask from the volume supplied (by
        averaging across time points and changing the range). If it is set
        to false then it will use the template_name as an input.

    Returns
    ----------

    mask : 3 dimensional array, binary
        The masked brain, thresholded to distinguish brain and non-brain

    template : 3 dimensional array, float
        A continuous (0 -> 1) volume describing the likelihood a voxel is in
        the brain. This can be used to contrast the brain and non brain.

    """

    # If the volume supplied is a 1d array then output a volume of the
    # supplied dimensions
    if len(volume.shape) == 1:
        volume = np.ones(volume)

    # Load in the mask
    if mask_self is True:
        mask_raw = volume
    elif template_name is None:
        mask_raw = np.load(resource_stream(__name__, "grey_matter_mask.npy"))
    else:
        mask_raw = np.load(template_name)

    # Make the masks 3dremove_baseline
    if len(mask_raw.shape) == 3:
        mask_raw = np.array(mask_raw)
    elif len(mask_raw.shape) == 4 and mask_raw.shape[3] == 1:
        mask_raw = np.array(mask_raw[:, :, :, 0])
    else:
        mask_raw = np.mean(mask_raw, 3)

    # Find the max value (so you can calulate these as proportions)
    mask_max = mask_raw.max()

    # Make sure the mask values range from 0 to 1 (make out of max of volume
    #  so that this is invertible later)
    mask_raw = mask_raw / mask_max

    # If there is only one brain volume then make this a forth dimension
    if len(volume.shape) == 3:
        temp = np.zeros([volume.shape[0], volume.shape[1], volume.shape[2], 1])
        temp[:, :, :, 0] = volume
        volume = temp

    # Reshape the mask to be the size as the brain
    brain_dim = volume.shape
    mask_dim = mask_raw.shape

    zoom_factor = (brain_dim[0] / mask_dim[0],
                   brain_dim[1] / mask_dim[1],
                   brain_dim[2] / mask_dim[2],
                   )

    # Scale the mask according to the input brain
    # You might get a warning if the zoom_factor is not an integer but you
    # can safely ignore that.
    template = ndimage.zoom(mask_raw, zoom_factor, order=2)
    template[template < 0] = 0

    # If the mask threshold is not supplied then guess it is a minima
    # between the two peaks of the bimodal distribution of voxel activity
    if mask_threshold is None:

        # How many bins on either side of a peak will be compared
        order = 5

        # Make the histogram
        template_vector = template.reshape(brain_dim[0] * brain_dim[1] *
                                           brain_dim[2])
        template_hist = np.histogram(template_vector, 100)

        # Zero pad the values
        binval = np.concatenate([np.zeros((order,)), template_hist[0]])
        bins = np.concatenate([np.zeros((order,)), template_hist[1]])

        # Identify the first two peaks
        peaks = signal.argrelmax(binval, order=order)[0][0:2]

        # What is the minima between peaks
        minima = binval[peaks[0]:peaks[1]].min()

        # What is the index of the last idx with this min value (since if
        # zero, there may be many)
        minima_idx = (np.where(binval[peaks[0]:peaks[1]] == minima) + peaks[
            0])[-1]

        # Convert the minima into a threshold
        mask_threshold = bins[minima_idx][0]

    # Mask the template based on the threshold
    mask = np.zeros(template.shape)
    mask[template > mask_threshold] = 1

    return mask, template

Creating a Python function

sim.mask_brain() is a Python "function". In general, a Python function has the following structure:

def function_name(input_1, input_2, ..., input_m):
    some code 
    some code
    ...
    some code
    return output1, output2, ... output_n

Exercise 1: Change the above script in at least 3 ways (examples: add a subplot of different slices, change the colors, show a histogram of values, etc.):

Contributions

M. Kumar, C. Ellis and N. Turk-Browne produced the initial notebook 01/2018
T. Meissner minor edits
Q. Lu: switch to matplotlib, fix dead links, add resources, encapsulate brainiak fmrisim
C. Ellis updated with comments from cmhn-s19