.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/06_manipulating_images/plot_roi_extraction.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_06_manipulating_images_plot_roi_extraction.py: Computing a Region of Interest (ROI) mask manually ================================================== This example shows manual steps to create and further modify an ROI spatial mask. They represent a means for "data folding", i.e., extracting and then analyzing brain data from a subset of voxels rather than whole brain images. Masking can also help alleviate the curse of dimensionality (i.e., statistical problems that arise in the context of high-dimensional input variables). We demonstrate how to compute a ROI mask using a **T-test** and then how simple image operations can be used before and after computing the ROI to improve the quality of the computed mask. These chains of operations are easy to set up using Nilearn and Scipy Python libraries. Here we give clear guidelines about these steps, starting with pre-image operations to post-image operations. The main point is that visualization & results checking are possible at each step. .. seealso:: :doc:`plot_extract_rois_smith_atlas` for automatic ROI extraction of brain connected networks given in 4D image. .. GENERATED FROM PYTHON SOURCE LINES 29-31 Here are the coordinates of the slice we are interested in each direction. We will be using them for visualization. .. GENERATED FROM PYTHON SOURCE LINES 31-42 .. code-block:: Python # cut in x-direction sagittal = -25 # cut in y-direction coronal = -37 # cut in z-direction axial = -6 # coordinates displaying should be prepared as a list cut_coords = [sagittal, coronal, axial] .. GENERATED FROM PYTHON SOURCE LINES 43-48 Loading the data ---------------- We will use the Haxby dataset to demonstrate the complete list of operations. The data will then be automatically stored in our home directory under ``nilearn_data/``. .. GENERATED FROM PYTHON SOURCE LINES 48-76 .. code-block:: Python from nilearn import datasets # First, we fetch EPI images and masks images from the Haxby dataset. haxby_dataset = datasets.fetch_haxby() # Print basic information on the dataset. # Functional data fmri_filename = haxby_dataset.func[0] print( f"First subject functional nifti image (4D) is located at: {fmri_filename}" ) print( "Labels of the Haxby dataset (text file) is located " f"at: {haxby_dataset.session_target[0]}" ) # Second, load the labels stored in a text file into array using pandas. import pandas as pd run_target = pd.read_csv(haxby_dataset.session_target[0], sep=" ") # Now, we have the labels that will be useful while computing the # student's t-test. haxby_labels = run_target["labels"] .. rst-class:: sphx-glr-script-out .. code-block:: none [fetch_haxby] Dataset directory found: /home/runner/work/nilearn/nilearn/nilearn_data/haxby2001 First subject functional nifti image (4D) is located at: /home/runner/work/nilearn/nilearn/nilearn_data/haxby2001/subj2/bold.nii.gz Labels of the Haxby dataset (text file) is located at: /home/runner/work/nilearn/nilearn/nilearn_data/haxby2001/subj2/labels.txt .. GENERATED FROM PYTHON SOURCE LINES 77-80 We now have the paths to the images in this dataset. The next step is to do a simple pre-processing step called `image smoothing` on the functional images and then build a statistical test on smoothed images. .. GENERATED FROM PYTHON SOURCE LINES 82-96 Build a statistical test to find voxels of interest --------------------------------------------------- Smoothing ^^^^^^^^^ Functional MRI data have a low signal-to-noise ratio. When using methods that are not robust to noise, it is useful to apply a spatial filtering kernel on the data. Such data smoothing is usually applied using a Gaussian function with 4mm to 12mm :term:`full-width at half-maximum` (this is where the ``fwhm`` parameter below comes from). The function :func:`~nilearn.image.smooth_img` accounts for potential anisotropy in the image affine (i.e., non-identical :term:`voxel` size in all the three dimensions). Analogous to the majority of nilearn functions, :func:`~nilearn.image.smooth_img` can also use file names as input parameters. .. GENERATED FROM PYTHON SOURCE LINES 96-118 .. code-block:: Python # Smooth the data using image processing module from nilearn. # smoothing: first argument as functional data filename and smoothing value # (integer) in second argument. Output is a Nifti image. from nilearn.image import smooth_img fmri_img = smooth_img(fmri_filename, fwhm=6) # Visualize the mean of the smoothed EPI image using plotting function # `plot_epi`. # First, compute the voxel-wise mean of the smooth EPI image # (first argument) using the image processing module `image`. from nilearn.image import mean_img from nilearn.plotting import plot_epi, show mean_img = mean_img(fmri_img) # Second, we visualize the mean image with coordinates positioned manually. plot_epi(mean_img, title="Smoothed mean EPI", cut_coords=cut_coords) show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_001.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-124 Functional MRI data can be considered "high dimensional" given the p-versus-n ratio (e.g., p=~20,000-200,000 voxels for n=1000 samples or less). In this setting, machine-learning algorithms can perform poorly due to the so-called curse of dimensionality. However, simple means from classical statistics can help reduce the number of voxels. .. GENERATED FROM PYTHON SOURCE LINES 124-131 .. code-block:: Python from nilearn.image import get_data fmri_data = get_data(fmri_img) # number of voxels being x*y*z, samples in 4th dimension fmri_data.shape .. rst-class:: sphx-glr-script-out .. code-block:: none (40, 64, 64, 1452) .. GENERATED FROM PYTHON SOURCE LINES 132-147 Selecting features using a T-test ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Student's t-test (:func:`scipy.stats.ttest_ind`) is an established method to determine whether two distributions have a different mean value. It can be used to compare voxel time-series from two different experimental conditions (e.g., when houses or faces are shown to individuals during brain scanning). If the time-series distribution is similar in the two conditions, then the :term:`voxel` is not very interesting to discriminate the condition. This test returns p-values that represent probabilities that the two time-series were not drawn from the same distribution. The lower the p-value, the more discriminative is the voxel in distinguishing the two conditions (faces and houses). .. GENERATED FROM PYTHON SOURCE LINES 147-162 .. code-block:: Python import numpy as np from scipy import stats _, p_values = stats.ttest_ind( fmri_data[..., haxby_labels == "face"], fmri_data[..., haxby_labels == "house"], axis=-1, ) # Use a log scale for p-values log_p_values = -np.log10(p_values) # Set NAN values to zero log_p_values[np.isnan(log_p_values)] = 0.0 log_p_values[log_p_values > 10.0] = 10.0 .. GENERATED FROM PYTHON SOURCE LINES 163-171 Visualize statistical p-values .............................. Before visualizing, we transform the computed p-values to a Nifti-like image using function `new_img_like` from nilearn. First argument being a reference image and second argument should be p-values data to convert to a new image as output. This new image will have same header information as the reference image. .. GENERATED FROM PYTHON SOURCE LINES 171-175 .. code-block:: Python from nilearn.image import new_img_like log_p_values_img = new_img_like(fmri_img, log_p_values) .. GENERATED FROM PYTHON SOURCE LINES 176-179 Now, we visualize the log p-values image on the functional mean image as a background with coordinates given manually and a colorbar on the right side of the plot (by default, colorbar=True). .. GENERATED FROM PYTHON SOURCE LINES 179-191 .. code-block:: Python from nilearn.plotting import plot_stat_map plot_stat_map( log_p_values_img, mean_img, title="p-values", cut_coords=cut_coords, cmap="inferno", ) show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_002.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 192-197 Selecting features using f_classif ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ It is also possible to use the :func:`sklearn.feature_selection.f_classif` function, which works for feature selection in multi-class settings. .. GENERATED FROM PYTHON SOURCE LINES 199-211 Build a mask ------------ Thresholding ^^^^^^^^^^^^ We build the t-map to have better representation of voxels of interest, where voxels with lower p-values correspond to the most intense voxels. This can be done easily by applying a threshold to a t-map data in array. Note that we use log p-values data; we force values below 5 to 0 by thresholding. .. GENERATED FROM PYTHON SOURCE LINES 211-213 .. code-block:: Python log_p_values[log_p_values < 5] = 0 .. GENERATED FROM PYTHON SOURCE LINES 214-216 Visualize the reduced voxels of interest using statistical image plotting function. As shown above, we first transform data in array to Nifti image. .. GENERATED FROM PYTHON SOURCE LINES 216-229 .. code-block:: Python log_p_values_img = new_img_like(fmri_img, log_p_values) # Now, visualizing the created log p-values to image. plot_stat_map( log_p_values_img, mean_img, title="Thresholded p-values", cut_coords=cut_coords, cmap="inferno", ) show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_003.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 230-234 We can post-process the results obtained with simple operations such as mask intersection and :term:`dilation` to regularize the mask definition. The idea of using these operations are to have more compact or sparser blobs. .. GENERATED FROM PYTHON SOURCE LINES 236-241 Binarization and Intersection with Ventral Temporal (VT) mask ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We now want to restrict our investigation to the VT area. The corresponding spatial mask is provided in ``haxby_dataset.mask_vt``. We want to compute the intersection of this provided mask with our self-computed mask. .. GENERATED FROM PYTHON SOURCE LINES 241-253 .. code-block:: Python # self-computed mask bin_p_values = log_p_values != 0 # VT mask mask_vt_filename = haxby_dataset.mask_vt[0] # The first step is to load VT mask and at the same time to convert the # datatype from "number" to "boolean". from nilearn.image import load_img vt = get_data(load_img(mask_vt_filename)).astype(bool) .. GENERATED FROM PYTHON SOURCE LINES 254-257 We can then use a logical "and" operation - `numpy.logical_and` - to keep only voxels that have been selected in both masks. In neuroimaging jargon, this is called an "AND conjunction". .. GENERATED FROM PYTHON SOURCE LINES 257-259 .. code-block:: Python bin_p_values_and_vt = np.logical_and(bin_p_values, vt) .. GENERATED FROM PYTHON SOURCE LINES 260-268 Visualizing the mask intersection results using plotting function `plot_roi`, a function which can be used for visualizing target specific voxels. First, we create new image type of binarized and intersected mask (second argument) and use this created Nifti image type in visualization. Binarized values in data type boolean should be converted to int data type at the same time. Otherwise, an error will be raised. .. GENERATED FROM PYTHON SOURCE LINES 268-284 .. code-block:: Python bin_p_values_and_vt_img = new_img_like( fmri_img, bin_p_values_and_vt.astype(np.int32) ) # We visualize the mask using the computed mean of functional images as # background. from nilearn.plotting import plot_roi plot_roi( bin_p_values_and_vt_img, mean_img, cut_coords=cut_coords, title="Intersection with ventral temporal mask", ) show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_004.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 285-294 Dilation ^^^^^^^^ Thresholded functional brain images often contain scattered voxels across the brain. To consolidate such brain images towards more compact shapes, we use a `morphological dilation `_. This is a common step to be sure not to forget voxels located on the edge of a ROI. In other words, such operations can fill "holes" in masked :term:`voxel` representations. .. GENERATED FROM PYTHON SOURCE LINES 294-316 .. code-block:: Python # We use ndimage function from scipy Python library for mask dilation. from scipy.ndimage import binary_dilation # Input here is a binarized and intersected mask data # from the previous section. dil_bin_p_values_and_vt = binary_dilation(bin_p_values_and_vt) # Now, we visualize the same using `plot_roi` with the data being converted # to Nifti image. In all `new_img_like` calls, we use the same reference image. dil_bin_p_values_and_vt_img = new_img_like( fmri_img, dil_bin_p_values_and_vt.astype(np.int32) ) plot_roi( dil_bin_p_values_and_vt_img, mean_img, title="Dilated mask", cut_coords=cut_coords, ) show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_005.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 317-320 Finally, we end with splitting the connected ROIs to two hemispheres into two separate regions (ROIs). We use the function :func:`scipy.ndimage.label` from the scipy Python library. .. GENERATED FROM PYTHON SOURCE LINES 322-328 Identification of connected components ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The function :func:`scipy.ndimage.label` from the scipy Python library identifies immediately neighboring voxels in our voxels mask. It assigns a separate integer label to each one of them. .. GENERATED FROM PYTHON SOURCE LINES 328-338 .. code-block:: Python from scipy.ndimage import label labels, _ = label(dil_bin_p_values_and_vt) # we take first roi data with labels assigned as integer 1 first_roi_data = (labels == 5).astype(np.int32) # Similarly, second roi data is assigned as integer 2 second_roi_data = (labels == 3).astype(np.int32) .. GENERATED FROM PYTHON SOURCE LINES 339-342 Visualizing the connected components .................................... First, we create a Nifti image type from first roi data in a array. .. GENERATED FROM PYTHON SOURCE LINES 342-344 .. code-block:: Python first_roi_img = new_img_like(fmri_img, first_roi_data) .. GENERATED FROM PYTHON SOURCE LINES 345-349 Then, we visualize the same created Nifti image (first argument) with the mean of functional images as background (second argument). The cut_coords are the default now: coordinates are selected automatically and will be pointed exactly on the roi data. .. GENERATED FROM PYTHON SOURCE LINES 349-352 .. code-block:: Python plot_roi(first_roi_img, mean_img, title="Connected components: first ROI") .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_006.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 353-354 We do the same for the second roi data. .. GENERATED FROM PYTHON SOURCE LINES 354-359 .. code-block:: Python second_roi_img = new_img_like(fmri_img, second_roi_data) plot_roi(second_roi_img, mean_img, title="Connected components: second ROI") show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_007.png :alt: plot roi extraction :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 360-364 Use the new ROIs to extract data maps in both ROIs -------------------------------------------------- We extract data from ROIs using Nilearn's :class:`~nilearn.maskers.NiftiLabelsMasker`. .. GENERATED FROM PYTHON SOURCE LINES 364-366 .. code-block:: Python from nilearn.maskers import NiftiLabelsMasker .. GENERATED FROM PYTHON SOURCE LINES 367-371 Before data extraction, we convert array labels to a Nifti like image. All inputs to ``NiftiLabelsMasker`` must be Nifti-like images or filenames to Nifti images. We use the same reference image as used above in previous sections. .. GENERATED FROM PYTHON SOURCE LINES 371-373 .. code-block:: Python labels_img = new_img_like(fmri_img, labels) .. GENERATED FROM PYTHON SOURCE LINES 374-378 First, we initialize a masker with parameters suited for data extraction: labels as input image, ``resampling_target`` is None as the affine and shape/size are the same for all the data used here, time series signal processing parameters ``standardize`` and ``detrend`` are set to ``False``. .. GENERATED FROM PYTHON SOURCE LINES 378-386 .. code-block:: Python masker = NiftiLabelsMasker( labels_img, resampling_target=None, standardize=None, detrend=False, verbose=1, ) .. GENERATED FROM PYTHON SOURCE LINES 387-389 Preparing for data extraction: setting number of conditions, size, etc. from the Haxby dataset. .. GENERATED FROM PYTHON SOURCE LINES 389-395 .. code-block:: Python condition_names = haxby_labels.unique() n_cond_img = fmri_data[..., haxby_labels == "house"].shape[-1] n_conds = len(condition_names) X1, X2 = np.zeros((n_cond_img, n_conds)), np.zeros((n_cond_img, n_conds)) .. GENERATED FROM PYTHON SOURCE LINES 396-400 Gathering data for each condition and then use :meth:`~nilearn.maskers.NiftiLabelsMasker.fit_transform` on each data. The transformer extracts data in condition maps where the target regions are specified by labels images. .. GENERATED FROM PYTHON SOURCE LINES 400-408 .. code-block:: Python for i, cond in enumerate(condition_names): cond_maps = new_img_like( fmri_img, fmri_data[..., haxby_labels == cond][..., :n_cond_img] ) mask_data = masker.fit_transform(cond_maps) X1[:, i], X2[:, i] = mask_data[:, 0], mask_data[:, 1] condition_names[np.where(condition_names == "scrambledpix")] = "scrambled" .. rst-class:: sphx-glr-script-out .. code-block:: none [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Loading regions from [NiftiLabelsMasker.wrapped] Finished fit [NiftiLabelsMasker.wrapped] Loading data from [NiftiLabelsMasker.wrapped] Extracting region signals [NiftiLabelsMasker.wrapped] Cleaning extracted signals .. GENERATED FROM PYTHON SOURCE LINES 409-410 Save the ROI 'atlas' to a Nifti file. .. GENERATED FROM PYTHON SOURCE LINES 410-418 .. code-block:: Python from pathlib import Path output_dir = Path.cwd() / "results" / "plot_roi_extraction" output_dir.mkdir(exist_ok=True, parents=True) print(f"Output will be saved to: {output_dir}") new_img_like(fmri_img, labels).to_filename(output_dir / "mask_atlas.nii.gz") .. rst-class:: sphx-glr-script-out .. code-block:: none Output will be saved to: /home/runner/work/nilearn/nilearn/examples/06_manipulating_images/results/plot_roi_extraction .. GENERATED FROM PYTHON SOURCE LINES 419-420 Plot the average in the different condition names. .. GENERATED FROM PYTHON SOURCE LINES 420-432 .. code-block:: Python import matplotlib.pyplot as plt plt.figure(figsize=(15, 7)) for i in np.arange(2): plt.subplot(1, 2, i + 1) plt.boxplot(X1 if i == 0 else X2) plt.xticks( np.arange(len(condition_names)) + 1, condition_names, rotation=25 ) plt.title(f"Boxplots of data in ROI{int(i + 1)} per condition") show() .. image-sg:: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_008.png :alt: Boxplots of data in ROI1 per condition, Boxplots of data in ROI2 per condition :srcset: /auto_examples/06_manipulating_images/images/sphx_glr_plot_roi_extraction_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 24.583 seconds) **Estimated memory usage:** 2394 MB .. _sphx_glr_download_auto_examples_06_manipulating_images_plot_roi_extraction.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/nilearn/nilearn/0.14.1?urlpath=lab/tree/notebooks/auto_examples/06_manipulating_images/plot_roi_extraction.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_roi_extraction.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_roi_extraction.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_roi_extraction.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_