.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_examples_03_connectivity_plot_extract_regions_dictlearning_maps.py>`     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_03_connectivity_plot_extract_regions_dictlearning_maps.py:


Regions extraction using Dictionary Learning and functional connectomes
=======================================================================

This example shows how to use :class:`nilearn.regions.RegionExtractor`
to extract spatially constrained brain regions from whole brain maps decomposed
using dictionary learning and use them to build a functional connectome.

We used 20 movie-watching functional datasets from
:func:`nilearn.datasets.fetch_development_fmri` and
:class:`nilearn.decomposition.DictLearning` for set of brain atlas maps.

This example can also be inspired to apply the same steps to even regions extraction
using ICA maps. In that case, idea would be to replace dictionary learning to canonical
ICA decomposition using :class:`nilearn.decomposition.CanICA`

Please see the related documentation of :class:`nilearn.regions.RegionExtractor`
for more details.

.. note::

    The use of the attribute `components_img_` from dictionary learning
    estimator is implemented from version 0.4.1. For older versions,
    unmask the deprecated attribute `components_` to get the components
    image using attribute `masker_` embedded in estimator.
    See the :ref:`section Inverse transform: unmasking data <unmasking_step>`.

Fetch brain development functional datasets
------------------------------------------------------------

We use nilearn's datasets downloading utilities


.. code-block:: default

    from nilearn import datasets

    rest_dataset = datasets.fetch_development_fmri(n_subjects=20)
    func_filenames = rest_dataset.func
    confounds = rest_dataset.confounds








Extract functional networks with DictionaryLearning
-----------------------------------------------------------------------


.. code-block:: default


    # Import dictionary learning algorithm from decomposition module and call the
    # object and fit the model to the functional datasets
    from nilearn.decomposition import DictLearning

    # Initialize DictLearning object
    dict_learn = DictLearning(n_components=8, smoothing_fwhm=6.,
                              memory="nilearn_cache", memory_level=2,
                              random_state=0)
    # Fit to the data
    dict_learn.fit(func_filenames)
    # Resting state networks/maps in attribute `components_img_`
    # Note that this attribute is implemented from version 0.4.1.
    # For older versions, see the note section above for details.
    components_img = dict_learn.components_img_

    # Visualization of functional networks
    # Show networks using plotting utilities
    from nilearn import plotting

    plotting.plot_prob_atlas(components_img, view_type='filled_contours',
                             title='Dictionary Learning maps')




.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_001.png
    :alt: plot extract regions dictlearning maps
    :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    /home/nicolas/GitRepos/nilearn-fork/nilearn/plotting/displays.py:101: UserWarning: linewidths is ignored by contourf
      im = getattr(ax, type)(data_2d.copy(),
    /home/nicolas/GitRepos/nilearn-fork/nilearn/plotting/displays.py:101: UserWarning: No contour levels were found within the data range.
      im = getattr(ax, type)(data_2d.copy(),

    <nilearn.plotting.displays.OrthoSlicer object at 0x7f88205944c0>



Extract regions from networks
------------------------------


.. code-block:: default


    # Import Region Extractor algorithm from regions module
    # threshold=0.5 indicates that we keep nominal of amount nonzero voxels across all
    # maps, less the threshold means that more intense non-voxels will be survived.
    from nilearn.regions import RegionExtractor

    extractor = RegionExtractor(components_img, threshold=0.5,
                                thresholding_strategy='ratio_n_voxels',
                                extractor='local_regions',
                                standardize=True, min_region_size=1350)
    # Just call fit() to process for regions extraction
    extractor.fit()
    # Extracted regions are stored in regions_img_
    regions_extracted_img = extractor.regions_img_
    # Each region index is stored in index_
    regions_index = extractor.index_
    # Total number of regions extracted
    n_regions_extracted = regions_extracted_img.shape[-1]

    # Visualization of region extraction results
    title = ('%d regions are extracted from %d components.'
             '\nEach separate color of region indicates extracted region'
             % (n_regions_extracted, 8))
    plotting.plot_prob_atlas(regions_extracted_img, view_type='filled_contours',
                             title=title)




.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_002.png
    :alt: plot extract regions dictlearning maps
    :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    /home/nicolas/GitRepos/nilearn-fork/nilearn/plotting/displays.py:101: UserWarning: No contour levels were found within the data range.
      im = getattr(ax, type)(data_2d.copy(),
    /home/nicolas/GitRepos/nilearn-fork/nilearn/plotting/displays.py:101: UserWarning: linewidths is ignored by contourf
      im = getattr(ax, type)(data_2d.copy(),
    /home/nicolas/anaconda3/envs/nilearn/lib/python3.8/site-packages/numpy/ma/core.py:2825: UserWarning: Warning: converting a masked element to nan.
      _data = np.array(data, dtype=dtype, copy=copy,

    <nilearn.plotting.displays.OrthoSlicer object at 0x7f8820457b50>



Compute correlation coefficients
---------------------------------


.. code-block:: default


    # First we need to do subjects timeseries signals extraction and then estimating
    # correlation matrices on those signals.
    # To extract timeseries signals, we call transform() from RegionExtractor object
    # onto each subject functional data stored in func_filenames.
    # To estimate correlation matrices we import connectome utilities from nilearn
    from nilearn.connectome import ConnectivityMeasure

    correlations = []
    # Initializing ConnectivityMeasure object with kind='correlation'
    connectome_measure = ConnectivityMeasure(kind='correlation')
    for filename, confound in zip(func_filenames, confounds):
        # call transform from RegionExtractor object to extract timeseries signals
        timeseries_each_subject = extractor.transform(filename, confounds=confound)
        # call fit_transform from ConnectivityMeasure object
        correlation = connectome_measure.fit_transform([timeseries_each_subject])
        # saving each subject correlation to correlations
        correlations.append(correlation)

    # Mean of all correlations
    import numpy as np
    mean_correlations = np.mean(correlations, axis=0).reshape(n_regions_extracted,
                                                              n_regions_extracted)








Plot resulting connectomes
----------------------------


.. code-block:: default


    title = 'Correlation between %d regions' % n_regions_extracted

    # First plot the matrix
    display = plotting.plot_matrix(mean_correlations, vmax=1, vmin=-1,
                                   colorbar=True, title=title)

    # Then find the center of the regions and plot a connectome
    regions_img = regions_extracted_img
    coords_connectome = plotting.find_probabilistic_atlas_cut_coords(regions_img)

    plotting.plot_connectome(mean_correlations, coords_connectome,
                             edge_threshold='90%', title=title)




.. rst-class:: sphx-glr-horizontal


    *

      .. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_003.png
          :alt: plot extract regions dictlearning maps
          :class: sphx-glr-multi-img

    *

      .. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_004.png
          :alt: plot extract regions dictlearning maps
          :class: sphx-glr-multi-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none


    <nilearn.plotting.displays.OrthoProjector object at 0x7f88203030a0>



Plot regions extracted for only one specific network
----------------------------------------------------


.. code-block:: default


    # First, we plot a network of index=4 without region extraction (left plot)
    from nilearn import image

    img = image.index_img(components_img, 4)
    coords = plotting.find_xyz_cut_coords(img)
    display = plotting.plot_stat_map(img, cut_coords=coords, colorbar=False,
                                     title='Showing one specific network')




.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_005.png
    :alt: plot extract regions dictlearning maps
    :class: sphx-glr-single-img





Now, we plot (right side) same network after region extraction to show that
connected regions are nicely seperated.
Each brain extracted region is identified as separate color.


.. code-block:: default


    # For this, we take the indices of the all regions extracted related to original
    # network given as 4.
    regions_indices_of_map3 = np.where(np.array(regions_index) == 4)

    display = plotting.plot_anat(cut_coords=coords,
                                 title='Regions from this network')

    # Add as an overlay all the regions of index 4
    colors = 'rgbcmyk'
    for each_index_of_map3, color in zip(regions_indices_of_map3[0], colors):
        display.add_overlay(image.index_img(regions_extracted_img, each_index_of_map3),
                            cmap=plotting.cm.alpha_cmap(color))

    plotting.show()



.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_extract_regions_dictlearning_maps_006.png
    :alt: plot extract regions dictlearning maps
    :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 2 minutes  53.333 seconds)


.. _sphx_glr_download_auto_examples_03_connectivity_plot_extract_regions_dictlearning_maps.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example


  .. container:: binder-badge

    .. image:: images/binder_badge_logo.svg
      :target: https://mybinder.org/v2/gh/nilearn/nilearn.github.io/master?filepath=examples/auto_examples/03_connectivity/plot_extract_regions_dictlearning_maps.ipynb
      :alt: Launch binder
      :width: 150 px


  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: plot_extract_regions_dictlearning_maps.py <plot_extract_regions_dictlearning_maps.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_extract_regions_dictlearning_maps.ipynb <plot_extract_regions_dictlearning_maps.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_