<jupyter_start><jupyter_text>`fMRI-04`: Software, spaces, and maskingThis notebook is a mixed bag of basic concepts in fMRI research. First, we'll highlight the variety of software packages used in fMRI research. Next, we'll focus on the idea of a "coordinate space" with accompanying concepts like templates, atlases, and masking.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')
sns.set_context('notebook', font_scale=1.5)
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Part 1: The fMRI software ecosystemTraditionally, fMRI researchers have relied on one of several software libraries for performing fMRI analysis. Packages such as [SPM](https://www.fil.ion.ucl.ac.uk/spm/) (MATLAB, UCL), [FSL](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) (C/Python, Oxford), [AFNI](https://afni.nimh.nih.gov/) (C, NIMH), and [FreeSurfer](https://surfer.nmr.mgh.harvard.edu/) (C/MATLAB, MGH) are still the major workhouses for fMRI analysis. As the field of fMRI analysis has advanced (and as the wider fMRI community has grown more computationally savvy), there's been an increasing push to centralize and standardize tools and incorporate fMRI analysis into larger analytic workflows.Recently, there's been a collosal effort to move fMRI analysis into the Python ecosystem (and also R). The most prominent example of this is [NIPY](https://nipy.org/), a constellation of python packages for (f)MRI preprocessing and analysis including:- [nipype](https://nipype.readthedocs.io/en/latest/): a python pipeline framework that interfaces with existing neuroimaging software (the basis for [fMRIPrep](https://fmriprep.org/en/stable/))- [nibabel](https://nipy.org/nibabel/): provides read/write protocols to many common neuroimaging file formats (e.g. [NIfTI](https://nifti.nimh.nih.gov/))- [nilearn](https://nilearn.github.io/): comprehensive analytic toolkit for fMRI data- [BrainIAK](https://brainiak.org/): advanced computational analysis tools for fMRI- [nitime](http://nipy.org/nitime/): time-series analysis of data from neuroscience experiments- [PyMVPA](http://www.pymvpa.org/): multivariate pattern analysis for neuroimaging data- [CoSMoMVPA](https://www.cosmomvpa.org/): MATLAB port of PyMVPA- [mindboggle](https://mindboggle.info/): automated labeling and morphometric analysis of human brain image data - [mne](http://martinos.org/mne/stable/index.html): software for M/EEG data analysis- [lyman](https://github.com/mwaskom/lyman): python port of FSL- [FIRdeconvolution](https://github.com/tknapen/FIRDeconvolution): deconvolution algorithms for timeseries data- [neurodesign](https://github.com/neuropower/neurodesign): tools for fMRI experimental design and power analysisMany software tools provide some basic visualization as well. However, there are also a variety of dedicated visualization tools:- [MRIcron](https://people.cas.sc.edu/rorden/mricron/index.html): basic fMRI visualization- [pycortex](https://gallantlab.github.io/pycortex/): advanced cortical surface visualizaiton- [pysurfer](https://pysurfer.github.io/): cortical surface visualization in Python- [Surf-Ice](https://github.com/neurolabusc/surf-ice): standalone software for cortical surface visualization  The Brain Imaging Data Structure (BIDS)There has (finally) been widespread adoption of a standardized format for fMRI data as well as the accompanying metadata and derivatives: the [Brain Imaging Data Structure (BIDS)](https://bids.neuroimaging.io/). This has accelerated data sharing and allowed for the development of automated processing pipelines:- [fMRIPrep](https://fmriprep.org/en/stable/): adaptive, automated preprocessing for fMRI data- [MRIQC](https://mriqc.readthedocs.io/en/stable/): automated quality control for fMRI data- [fitlins](https://github.com/poldracklab/fitlins): automated first-level analysis via fMRIPrep  Where to begin...With so many options, the ecosystem for fMRI analysis can feel overwhelming. For this course, we'll be relying on only a handful of software packages with an emphasis on current best practices and the ubiquity of Python. We'll use [fMRIPrep](https://fmriprep.org/en/stable/) for automated preprocessing and [nilearn](https://nilearn.github.io/stable/index.html) for most downstream analyses. The [Princeton Handbook for Reproducible Neuroimaging](https://brainhack-princeton.github.io/handbook/) provides a good starting point.  Part 2: Coordinate space(s)Over the course of fMRI processing, the functional data are often spatially transformed into another "space"—typically a standard space. We'll need to keep the following distinctions in mind during analysis:* native space vs. standard space* volume space vs. surface space* standard space and templates, e.g. MNI152, MNI305, fsaverage* voxel space vs. world space Native space vs. standard spaceThe first distinction is between native and standard space. Native space refers to the subject-specific configuration under which the data were collected. Analyzing data in native space requires the minimal amount of data manipulation (e.g. aligning and resampling), but the lack of voxelwise correspondence across subjects prevents some types of group analysis. To achieve voxelwise correspondence across subjects, we typically perform spatial normalization to align each subject's data to a template image (e.g. *MNI152NLin2009cAsym*) in a standard space (e.g. *MNI space*). There are several standard spaces. Volume space vs. surface spaceThe second distinction is between volume- and surface-based analysis. MRI data are collected in a volumetric format where each data point corresponds to a voxel at a particular location in a three-dimensional grid; volumetric analysis retains this three-dimensional grid representation of the data. In contrast, surface-based analysis is argued to better capture the sheet-like structure of the cortical mantle. In surace-based analysis, the functional data are resampled onto two-dimensional graph-like triangular mesh where each data point corresponds to a vertex at a particular location in the mesh. Importantly, two voxels on either side of a sulcus may appear to be close according to their Euclidean distance in volume space, whereas surface space more accurately that they are farther apart in terms of geodesic distance on the folded cortical surface. The *fsaverage* template is typically used for spatial normalization of surface-based data. Standard space and templatesThe third distinction is between a standard space and the specific *template* in that space. A template is an actual brain image (typically the average of many subjects) that serves as a target for alignment. Although "MNI space" is very widely used, there are in fact many specific *templates* in MNI space that each differ in small ways (read more [about the MNI space(s)](https://www.lead-dbs.org/about-the-mni-spaces/) at Lead-DBS). Voxel space vs. world spaceThe fourth distinction is between voxel space and world space. A standard space (e.g. MNI space) is defined in real-world coordinates. However, any given brain image is rendered in a voxel space (or "image space") with a particular geometry and grid: i.e. voxel size, field of view, etc. The translation between the voxel space of a given brain image and real-world coordinates are stored in an matrix accompanying the data (called an *affine matrix*). When we align a brain image to a standard space, we align it in real-world coordinates, but then we resample (or slice it) onto a particular grid of voxels. Different brain images with different geometries (voxel size, field of view, etc) can all be *in* (i.e. *aligned to*) MNI space.  Part 3: Working with Volume Data Reading / Writing Volume DataReading in MRI data is easy with `nibabel`. For a more complete tutorial, see [here](https://nipy.org/nibabel/gettingstarted.html).<jupyter_code>import nibabel as nib

# Load anatomical data.
img = nib.load('T1.nii.gz')<jupyter_output><empty_output><jupyter_text>Nibabel objects store many useful pieces of information, including the data, the affine transformation, and the header.<jupyter_code># Show header metadata
print(img.header.keys()[:10])
# Show affine matrix
print(img.affine)<jupyter_output>[[  -2.    0.    0.   74.]
 [   0.    0.    2. -109.]
 [   0.   -2.    0.   85.]
 [   0.    0.    0.    1.]]<jupyter_text>Extracting the actual data is similarly easy.<jupyter_code># Extract data
data = img.get_fdata()
print(data.shape)

# Visualize one slice
sns.heatmap(data[50], cmap='binary_r', square=True, cbar=False, 
            xticklabels=[], yticklabels=[]);<jupyter_output>(76, 76, 93)<jupyter_text>Saving data is also easy. Whenever we save data, we must also store its associated affine transformation (otherwise we will not know how to interpret its spatial reference frame).<jupyter_code># Convert data back into a NIfTI image
new_img = nib.Nifti1Image(data, img.affine)

# Save data to disk
filename = 'example.nii.gz'
# nib.save(nifti, filename)<jupyter_output><empty_output><jupyter_text>Masking volumetric dataTypically, we don't actually want to analyize an entire MRI image. We want to extract some subset of voxels based on a particular "mask"—for example, the entire brain (excluding everything outside the brain), only gray matter (excluding white matter, CSF, etc), or a particular region of interest (ROI; e.g. ventral temporal cortex). As long as our data and mask are in the same space and on the same grid, we can mask our data using simple NumPy indexing. Here, we demonstrate masking with the brainmask and with an anatomical atlas.<jupyter_code># Load brain mask
mask = nib.load('brainmask.nii.gz').get_fdata()
print(f"Data shape: {data.shape}")
print(f"Mask shape: {mask.shape}")

# Visualize one slice of mask
plt.matshow(mask[50], cmap='binary_r')
plt.xticks([])
plt.yticks([])
plt.show()
# Make copy of data
copy = data.copy()

# Zero out non-brain voxels for visualization
copy[mask == 0] = 0

# Visualize one slice
plt.matshow(copy[50], cmap='binary_r')
plt.xticks([])
plt.yticks([])
plt.show()<jupyter_output><empty_output><jupyter_text>When we extract data based on a mask, the original 3-dimensional data array is collapsed into a single dimension. We can use `np.where` to keep track of where the masked voxels came from in the original image.<jupyter_code># Count the voxels in our mask
mask_count = np.sum(mask != 0)
print(f"{mask.size} total voxels in image")
print(f"{mask_count} nonzero voxels in mask")

# Masked data is no longer 3-dimensional
masked_data = data[mask != 0]
print(f"Masked data shape: {masked_data.shape}")
# Track the original mask coordinates
mask_coords = np.where(mask != 0)

# Create an empty new image with original shape
new_data = np.zeros(img.shape)

# Modify our masked data in some way
result = masked_data + np.random.randn(*masked_data.shape) * 20
new_data[mask_coords] = result

# Visualize one slice
plt.matshow(new_data[50])
plt.xticks([])
plt.yticks([])
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09095034945992557, 'fraction_numerical': 0.006081510393028955, 'mean_word_length': 5.13132999443517, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 29, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '664380', 'license': 'MIT', 'n_tokens_mistral': 3278, 'n_tokens_neox': 3008, 'n_words': 1312, 'path': 'fmri-04/fmri-04-spaces.ipynb', 'repo_name': '2022-NEU502b/demos', 'size': '11017'}