Merge pull request #1516 from jeromedockes/cortex_surface_projections

[MRG] Cortex surface projections

doc/modules/reference.rst

@@ -284,6 +284,7 @@ uses.
 
 .. _plotting_ref:
 
+
 :mod:`nilearn.plotting`: Plotting brain data
 ================================================
 
@@ -349,3 +350,25 @@ uses.
 
    clean
    high_variance_confounds
+
+
+:mod:`nilearn.surface`: Manipulating surface data
+===================================================
+
+.. automodule:: nilearn.surface
+   :no-members:
+   :no-inherited-members:
+
+.. No relevant user manual section yet.
+
+**Functions**:
+
+.. currentmodule:: nilearn.surface
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   load_surf_data
+   load_surf_mesh
+   vol_to_surf

doc/plotting/index.rst

@@ -346,7 +346,6 @@ NOTE: These functions works for only with matplotlib higher than 1.3.1.
      :target: ../auto_examples/01_plotting/plot_surf_stat_map.html
      :scale: 50
 
-
 =====================   ===================================================================
 =====================   ===================================================================
 |plot_surf_roi|          :func:`plot_surf_roi`
@@ -362,5 +361,5 @@ NOTE: These functions works for only with matplotlib higher than 1.3.1.
                          |hack|
                          **Example:**
                          :ref:`sphx_glr_auto_examples_01_plotting_plot_surf_stat_map.py`
-=====================   ===================================================================
 
+=====================   ===================================================================

doc/whats_new.rst

@@ -61,6 +61,9 @@ Enhancements
       :func:`nilearn.plotting.plot_surf_stat_map` and
       :func:`nilearn.plotting.plot_surf_roi`
 
+    - :func:`nilearn.surface.vol_to_surf` projects a 3d or 4d brain
+      volume on the cortical surface.
+
     - :class:`nilearn.decoding.SearchLight` has new parameter "groups" to
       do LeaveOneGroupOut type cv with new scikit-learn module model selection.
 

examples/01_plotting/plot_3d_map_to_surface_projection.py

@@ -0,0 +1,57 @@
+"""
+Making a surface plot of a 3D statistical map
+=============================================
+
+project a 3D statistical map onto a cortical mesh using
+:func:`nilearn.surface.vol_to_surf`. Display a surface plot of the projected
+map using :func:`nilearn.plotting.plot_surf_stat_map`.
+
+NOTE: Example needs matplotlib version higher than 1.3.1.
+
+"""
+
+##############################################################################
+# Get a statistical map
+# ---------------------
+
+from nilearn import datasets
+
+localizer_dataset = datasets.fetch_localizer_button_task()
+localizer_tmap = localizer_dataset.tmaps[0]
+
+##############################################################################
+# Get a cortical mesh
+# -------------------
+
+fsaverage = datasets.fetch_surf_fsaverage5()
+
+##############################################################################
+# Sample the 3D data around each node of the mesh
+# -----------------------------------------------
+
+from nilearn import surface
+
+texture = surface.vol_to_surf(localizer_tmap, fsaverage.pial_right)
+
+##############################################################################
+# Plot the result
+# ---------------
+
+from nilearn import plotting
+
+plotting.plot_surf_stat_map(fsaverage.infl_right, texture, hemi='right',
+                            title='Surface right hemisphere',
+                            threshold=1., bg_map=fsaverage.sulc_right,
+                            cmap='cold_hot')
+
+##############################################################################
+# Plot 3D image for comparison
+# ----------------------------
+
+plotting.plot_glass_brain(localizer_tmap, display_mode='r', plot_abs=False,
+                          title='Glass brain', threshold=2.)
+
+plotting.plot_stat_map(localizer_tmap, display_mode='x', threshold=1.,
+                       cut_coords=range(0, 51, 10), title='Slices')
+
+plotting.show()

examples/01_plotting/plot_surf_stat_map.py

@@ -85,9 +85,9 @@
 # --------------------------------
 
 # Load resting state time series from nilearn
-from nilearn import plotting
+from nilearn import surface
 
-timeseries = plotting.surf_plotting.load_surf_data(nki_dataset['func_left'][0])
+timeseries = surface.load_surf_data(nki_dataset['func_left'][0])
 
 # Extract seed region via label
 pcc_region = b'G_cingul-Post-dorsal'
@@ -115,6 +115,9 @@
 
 ###############################################################################
 # Display ROI on surface
+
+from nilearn import plotting
+
 plotting.plot_surf_roi(fsaverage['pial_left'], roi_map=pcc_labels,
                        hemi='left', view='medial',
                        bg_map=fsaverage['sulc_left'], bg_on_data=True,

examples/01_plotting/plot_surface_projection_strategies.py

@@ -0,0 +1,69 @@
+"""
+Technical point: Illustration of the volume to surface sampling schemes
+=======================================================================
+
+In nilearn, :func:`nilearn.surface.vol_to_surf` allows us to measure values of
+a 3d volume at the nodes of a cortical mesh, transforming it into surface data.
+This data can then be plotted with :func:`nilearn.plotting.plot_surf_stat_map`
+for example.
+
+This script shows, on a toy example, where samples are drawn around each mesh
+vertex. Image values are interpolated at each sample location, then these
+samples are averaged to produce a value for the vertex.
+
+Two strategies are available to choose sample locations: they can be spread
+along the normal to the mesh, or inside a ball around the vertex. Don't worry
+too much about choosing one or the other: they take a similar amount of time
+and give almost identical results for most images.
+
+"""
+
+import numpy as np
+
+import matplotlib
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+from nilearn.surface import surface
+
+
+######################################################################
+# Build a mesh (of a cylinder)
+######################################################################
+
+N_Z = 5
+N_T = 10
+u, v = np.mgrid[:N_T, :N_Z]
+triangulation = matplotlib.tri.Triangulation(u.flatten(), v.flatten())
+angles = u.flatten() * 2 * np.pi / N_T
+x, y = np.cos(angles), np.sin(angles)
+z = v.flatten() * 2 / N_Z
+
+mesh = [np.asarray([x, y, z]).T, triangulation.triangles]
+
+
+#########################################################################
+# Get the locations from which vol_to_surf would draw its samples
+#########################################################################
+
+line_sample_points = surface._line_sample_locations(
+    mesh, np.eye(4), segment_half_width=.2, n_points=6)
+
+ball_sample_points = surface._ball_sample_locations(
+    mesh, np.eye(4), ball_radius=.15, n_points=20)
+
+
+######################################################################
+# Plot the mesh and the sample locations
+######################################################################
+
+for sample_points in [line_sample_points, ball_sample_points]:
+    fig = plt.figure()
+    ax = plt.subplot(projection='3d')
+    ax.set_aspect(1)
+
+    ax.plot_trisurf(x, y, z, triangles=triangulation.triangles)
+
+    ax.scatter(*sample_points.T, color='r')
+
+plt.show()

nilearn/__init__.py

@@ -58,4 +58,4 @@
 # list all submodules available in nilearn and version
 __all__ = ['datasets', 'decoding', 'decomposition', 'connectome',
            'image', 'input_data', 'masking', 'mass_univariate', 'plotting',
-           'region', 'signal', '__version__']
+           'region', 'signal', 'surface', '__version__']

nilearn/datasets/tests/test_func.py

@@ -390,7 +390,7 @@ def test__load_mixed_gambles():
         assert_equal(len(zmaps), len(gain))
 
 
-@with_setup(setup_mock)
+@with_setup(setup_mock, teardown_mock)
 @with_setup(tst.setup_tmpdata, tst.teardown_tmpdata)
 def test_fetch_mixed_gambles():
     local_url = "file://" + os.path.join(tst.datadir,

nilearn/plotting/__init__.py

@@ -41,7 +41,6 @@ def _set_mpl_backend():
 from .img_plotting import plot_img, plot_anat, plot_epi, \
     plot_roi, plot_stat_map, plot_glass_brain, plot_connectome, \
     plot_prob_atlas, show
-from .surf_plotting import plot_surf, plot_surf_stat_map, plot_surf_roi
 from .find_cuts import find_xyz_cut_coords, find_cut_slices
 from .matrix_plotting import plot_matrix
 __all__ = ['cm', 'plot_img', 'plot_anat', 'plot_epi',

nilearn/plotting/surf_plotting.py

@@ -2,144 +2,17 @@
 Functions for surface visualization.
 Only matplotlib is required.
 """
-
-# Import libraries
-import nibabel
 import numpy as np
+
 import matplotlib.pyplot as plt
 
 from mpl_toolkits.mplot3d import Axes3D
-from nibabel import gifti
 
+from ..surface import load_surf_data, load_surf_mesh
 from .._utils.compat import _basestring
 from .img_plotting import _get_colorbar_and_data_ranges
 
 
-# function to figure out datatype and load data
-def load_surf_data(surf_data):
-    """Loading data to be represented on a surface mesh.
-
-    Parameters
-    ----------
-    surf_data : str or numpy.ndarray
-        Either a file containing surface data (valid format are .gii,
-        .mgz, .nii, .nii.gz, or Freesurfer specific files such as
-        .thickness, .curv, .sulc, .annot, .label) or
-        a Numpy array containing surface data.
-    Returns
-    -------
-    data : numpy.ndarray
-        An array containing surface data
-    """
-    # if the input is a filename, load it
-    if isinstance(surf_data, _basestring):
-        if (surf_data.endswith('nii') or surf_data.endswith('nii.gz') or
-                surf_data.endswith('mgz')):
-            data = np.squeeze(nibabel.load(surf_data).get_data())
-        elif (surf_data.endswith('curv') or surf_data.endswith('sulc') or
-                surf_data.endswith('thickness')):
-            data = nibabel.freesurfer.io.read_morph_data(surf_data)
-        elif surf_data.endswith('annot'):
-            data = nibabel.freesurfer.io.read_annot(surf_data)[0]
-        elif surf_data.endswith('label'):
-            data = nibabel.freesurfer.io.read_label(surf_data)
-        elif surf_data.endswith('gii'):
-            gii = gifti.read(surf_data)
-            try:
-                data = np.zeros((len(gii.darrays[0].data), len(gii.darrays)))
-                for arr in range(len(gii.darrays)):
-                    data[:, arr] = gii.darrays[arr].data
-                data = np.squeeze(data)
-            except IndexError:
-                raise ValueError('Gifti must contain at least one data array')
-        else:
-            raise ValueError(('The input type is not recognized. %r was given '
-                              'while valid inputs are a Numpy array or one of '
-                              'the following file formats: .gii, .mgz, .nii, '
-                              '.nii.gz, Freesurfer specific files such as '
-                              '.curv, .sulc, .thickness, .annot, '
-                              '.label') % surf_data)
-    # if the input is a numpy array
-    elif isinstance(surf_data, np.ndarray):
-        data = np.squeeze(surf_data)
-    else:
-        raise ValueError('The input type is not recognized. '
-                         'Valid inputs are a Numpy array or one of the '
-                         'following file formats: .gii, .mgz, .nii, .nii.gz, '
-                         'Freesurfer specific files such as .curv,  .sulc, '
-                         '.thickness, .annot, .label')
-    return data
-
-
-# function to figure out datatype and load data
-def load_surf_mesh(surf_mesh):
-    """Loading a surface mesh geometry
-
-    Parameters
-    ----------
-    surf_mesh : str or numpy.ndarray
-        Either a file containing surface mesh geometry (valid formats
-        are .gii or Freesurfer specific files such as .orig, .pial,
-        .sphere, .white, .inflated) or a list of two Numpy arrays,
-        the first containing the x-y-z coordinates of the mesh
-        vertices, the second containing the indices (into coords)
-        of the mesh faces.
-
-    Returns
-    --------
-    [coords, faces] : List of two numpy.ndarray
-        The first containing the x-y-z coordinates of the mesh vertices,
-        the second containing the indices (into coords) of the mesh faces.
-    """
-    # if input is a filename, try to load it
-    if isinstance(surf_mesh, _basestring):
-        if (surf_mesh.endswith('orig') or surf_mesh.endswith('pial') or
-                surf_mesh.endswith('white') or surf_mesh.endswith('sphere') or
-                surf_mesh.endswith('inflated')):
-            coords, faces = nibabel.freesurfer.io.read_geometry(surf_mesh)
-        elif surf_mesh.endswith('gii'):
-            try:
-                coords = gifti.read(surf_mesh).getArraysFromIntent(
-                    nibabel.nifti1.intent_codes['NIFTI_INTENT_POINTSET'])[0].data
-            except IndexError:
-                raise ValueError('Gifti file needs to contain a data array '
-                                 'with intent NIFTI_INTENT_POINTSET')
-            try:
-                faces = gifti.read(surf_mesh).getArraysFromIntent(
-                    nibabel.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE'])[0].data
-            except IndexError:
-                raise ValueError('Gifti file needs to contain a data array '
-                                 'with intent NIFTI_INTENT_TRIANGLE')
-        else:
-            raise ValueError(('The input type is not recognized. %r was given '
-                              'while valid inputs are one of the following '
-                              'file formats: .gii, Freesurfer specific files '
-                              'such as .orig, .pial, .sphere, .white, '
-                              '.inflated or a list containing two Numpy '
-                              'arrays [vertex coordinates, face indices]'
-                              ) % surf_mesh)
-    elif isinstance(surf_mesh, list):
-        if len(surf_mesh) == 2:
-            coords, faces = surf_mesh[0], surf_mesh[1]
-        else:
-            raise ValueError(('If a list is given as input, it must have '
-                              'two elements, the first is a Numpy array '
-                              'containing the x-y-z coordinates of the mesh '
-                              'vertices, the second is a Numpy array '
-                              'containing  the indices (into coords) of the '
-                              'mesh faces. The input was a list with '
-                              '%r elements.') % len(surf_mesh))
-    else:
-        raise ValueError('The input type is not recognized. '
-                         'Valid inputs are one of the following file '
-                         'formats: .gii, Freesurfer specific files such as '
-                         '.orig, .pial, .sphere, .white, .inflated '
-                         'or a list containing two Numpy arrays '
-                         '[vertex coordinates, face indices]')
-
-    return [coords, faces]
-
-
 def plot_surf(surf_mesh, surf_map=None, bg_map=None,
               hemi='left', view='lateral', cmap=None,
               avg_method='mean', threshold=None, alpha='auto',