<jupyter_start><jupyter_text>Integrating Functional DataSo far most of our work has been examining anatomical images - the reason being is that it provides a nice visual way of exploring the effects of data manipulation and visualization is easy. In practice, you will most likely not analyze anatomical data using nilearn since there are other tools that are better suited for that kind of analysis (freesurfer, connectome-workbench, mindboggle, etc...). In this notebook we'll finally start working with functional MR data - the modality of interest in this workshop. First we'll cover some basics about how the data is organized (similar to T1s but slightly more complex), and then how we can integrate our anatomical and functional data together using tools provided by nilearn Functional data consists of full 3D brain volumes that are *sampled* at multiple time points. Therefore you have a sequence of 3D brain volumes, stepping through sequences is stepping through time and therefore time is our 4th dimension! Here's a visualization to make this concept more clear:  Each index along the 4th dimensions (called TR for "Repetition Time", or Sample) is a full 3D scan of the brain. Pulling out volumes from 4-dimensional images is similar to that of 3-dimensional images except you're now dealing with: img.slicer[x,y,z,time] !Let's try a couple of examples to familiarize ourselves with dealing with 4D images. But first, let's pull some functional data using PyBIDS!<jupyter_code>fmriprep_dir = '../data/ds000030/derivatives/fmriprep/'
layout=BIDSLayout(fmriprep_dir, validate=False)
T1w_files = layout.get(subject='10788', datatype='anat', suffix='preproc')
brainmask_files = layout.get(subject='10788', datatype='anat', suffix='brainmask')
func_files = layout.get(subject='10788', datatype='func', suffix='preproc')
func_mask_files = layout.get(subject='10788', datatype='func', suffix='brainmask')<jupyter_output><empty_output><jupyter_text>We'll be using functional files in MNI space rather than T1w space. Recall, that MNI space data is data that was been warped into standard space. These are the files you would typically use for a group-level functional imaging analysis! First, take a look at the shape of the functional image: Notice that the Functional MR scan contains *4 dimensions*. This is in the form of $(x,y,z,t)$, where $t$ is time. We can use slicer as usual where instead of using 3 dimensions we use 4. For example: func.slicer[x,y,z]  vs. func.slicer[x,y,z,t]   ExerciseTry pulling out the 5th TR and visualizing it using plot.plot_epi. plot_epi is exactly the same as plot_anat except it displays using colors that make more sense for functional images...<jupyter_code>#Pull the 5th TR
func_vol5 = func_mni_img.slicer[??,??,??,??]
plot.plot_epi(??)<jupyter_output><empty_output><jupyter_text>What fMRI actually represents We've represented fMRI as a snapshot of MR signal over multiple timepoints. This is a useful way of understanding the organization of fMRI, however it isn't typically how we think about the data when we analyze fMRI data. fMRI is typically thought of as **time-series** data. We can think of each voxel (x,y,z coordinate) as having a time-series of length T. The length T represents the number of volumes/timepoints in the data. Let's pick an example voxel and examine its time-series using func_mni_img.slicer:<jupyter_code>#Pick one voxel at coordinate (60,45,88)<jupyter_output><empty_output><jupyter_text>As you can see here, we pulled one voxel that contains 152 timepoints. For plotting purposes as 4-dimensional array is difficult to deal with so we'll flatten it to 1 dimension (time) for convenience: Here we've pulled out a voxel at a specific coordinate at every single time-point. This voxel has a single value for each timepoint and therefore is a time-series. We can visualize this time-series signal by using a standard python plotting library. We won't go into too much detail about python plotting, the intuition about what the data looks like is what is most important:<jupyter_code>import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>ResamplingRecall from our introductory exploration of neuroimaging data:- T1 images are typically composed of voxels that are 1x1x1 in dimension- Functional images are typically composed of voxels that are 4x4x4 in dimensionIf we'd like to overlay our functional on top of our T1 (for visualization purposes, or analyses), then we need to match the size of the voxels! Think of this like trying to overlay a 10x10 JPEG and a 20x20 JPEG on top of each other. To get perfect overlay we need to resize (or more accurately *resample*) our JPEGs to match!**Note**: Resampling is a method of interpolating in between data-points. When we stretch an image we need to figure out what goes in the spaces that are created via stretching - resampling does just that. In fact, resizing any type of image is actually just resampling to new dimensions.  Let's resampling some MRI data using nilearn. **Goal**: Match the dimensions of the structural image to that of the functional image<jupyter_code>#Files we'll be using (Notice that we're using _space-MNI..._ which means they are normalized brains)<jupyter_output><empty_output><jupyter_text>Let's take a look at the sizes of both our functional and structural files: Resampling in nilearn is as easy as telling it which image you want to sample and what the target image is.Structure of function:img.resample_to_img(source_img,target_img,interpolation) - source_img = the image you want to sample- target_img = the image you wish to *resample to* - interpolation = the method of interpolationA note on **interpolation**nilearn supports 3 types of interpolation, the one you'll use depends on the type of data you're resampling!1. **continuous** - Interpolate but maintain some edge features.  Ideal for structural images where edges are well-defined. Uses $3^\text{rd}$-order spline interpolation.2. **linear (default)** - Interpolate uses a combination of neighbouring voxels - will blur. Uses trilinear interpolation.3. **nearest** - matches value of closest voxel (majority vote from neighbours). This is ideal for masks which are binary since it will preserve the 0's and 1's and will not produce in-between values (ex: 0.342). Also ideal for numeric labels where values are 0,1,2,3... (parcellations). Uses nearest-neighbours interpolation with majority vote.<jupyter_code>#Try playing around with methods of interpolation
#options: 'linear','continuous','nearest'

import matplotlib.animation
from IPython.display import HTML
%%capture
%matplotlib inline
#Resample the T1 to the size of the functional image!
resamp_t1 = img.resample_to_img(source_img=T1_mni_img, target_img=func_mni_img, interpolation='continuous')
fig, ax = plt.subplots()

def animate(image):
    plot.plot_anat(image, figure=fig, cut_coords=(0,0,0))
    ax.set_facecolor('black')
    
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=[resamp_t1, T1_mni_img]) 
#change the frames to look at the functional mask over the resampled T1
# ani = matplotlib.animation.FuncAnimation(fig, animate, frames=[resamp_t1, func]) 
# Display animation
HTML(ani.to_jshtml())<jupyter_output><empty_output><jupyter_text>**Exercise**Using **Native** T1 and **T1w** resting state functional do the following:1. Resample the native T1 image to resting state size2. Replace the brain in the T1 image with the first frame of the resting state brain<jupyter_code>#Files we'll need


####STRUCTURAL FILES

#T1 image
ex_t1 = img.load_img(T1w_files[0].path)

#mask file
ex_t1_bm = img.load_img(brainmask_files[0].path)


####FUNCTIONAL FILES

#This is the pre-processed resting state data that hasn't been standardized
ex_func = img.load_img(func_files[1].path)

#This is the associated mask for the resting state image.
ex_func_bm = img.load_img(func_mask_files[1].path)<jupyter_output><empty_output><jupyter_text>The first step we need to do is to make sure the dimensions for our T1 image and resting state image match each other:<jupyter_code>#Resample the T1 to the size of the functional image!
resamp_t1 = img.resample_to_img(source_img=??, target_img=??, interpolation='continuous')
plot.plot_anat(??)
print(resamp_t1.shape)<jupyter_output><empty_output><jupyter_text>Next we want to make sure that the brain mask for the T1 is also the same dimensions as the functional image. This is exactly the same as above, except we use the brain mask as the source.What kind of interpolation should we use for masks?<jupyter_code>resamp_bm = img.??(??)

#Plot the image
??

print(resamp_bm.shape)<jupyter_output><empty_output><jupyter_text>Once we've resampled both our T1 and our brain mask. We now want to remove the brain from the T1 image so that we can replace it with the funtional image instead. Remember to do this we need to:1. Invert the T1 mask2. Apply the inverted mask to the brain<jupyter_code>inverted_bm_t1 = img.math_img(??,a=resamp_bm)
plot.plot_anat(inverted_bm_t1)<jupyter_output><empty_output><jupyter_text>Now apply the mask:<jupyter_code>resamp_t1_nobrain = img.??(??)
plot.plot_anat(resamp_t1_nobrain)<jupyter_output><empty_output><jupyter_text>We now have a skull missing the structural T1 brain. The final steps is to stick in the brain from the functional image into the now brainless head. First we need to remove the surrounding signal from the functional image.Since a functional image is 4-Dimensional, we'll need to pull the first volume to work with. This is because the structural image is 3-dimensional and operations will fail if we try to mix 3D and 4D data.<jupyter_code>#Let's visualize the first volume of the functional image:
first_vol = ex_func.slicer[??,??,??,??]
plot.plot_epi(first_vol)<jupyter_output><empty_output><jupyter_text>As shown in the figure above, the image has some "signal" outside of the brain. In order to place this within the now brainless head we made earlier, we need to mask out the functional MR data as well!<jupyter_code>#Mask first_vol using ex_func_bm
masked_func = img.math_img('??', a=??, b=??)
plot.plot_epi(masked_func)<jupyter_output><empty_output><jupyter_text>The final step is to stick this data into the head of the T1 data. Since the hole in the T1 data is represented as $0$'s. We can add the two images together to place the functional data into the void:<jupyter_code>#Now overlay the functional image on top of the anatomical
combined_img = img.math_img(??)
plot.plot_anat(combined_img)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07184889901290813, 'fraction_numerical': 0.012433561123766135, 'mean_word_length': 5.618718592964824, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '10329', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3169, 'n_tokens_neox': 2961, 'n_words': 1365, 'path': 'code/04-integrating_functional_data.ipynb', 'repo_name': 'samwalkow/SDC-BIDS-fMRI', 'size': '10536'}