Allow not scaling background maps for surface plots

doc/changes/latest.rst

@@ -16,6 +16,8 @@ Fixes
 Enhancements
 ------------
 
+- Surface plotting methods no longer automatically rescale background maps, which, among other things, allows to use curvature sign as a background map (:gh:`3173` by `Alexis Thual`_).
+
 Changes
 -------
 - The behavior of :func:`~nilearn.datasets.fetch_atlas_craddock_2012`, :func:`~nilearn.datasets.fetch_atlas_smith_2009` and :func:`~nilearn.datasets.fetch_atlas_basc_multiscale_2015` is updated with their new parameters to return one map along with a deprecation cycle (:gh:`3353` by `Ahmad Chamma`_).

examples/01_plotting/plot_3d_map_to_surface_projection.py

@@ -27,11 +27,23 @@
 fsaverage = datasets.fetch_surf_fsaverage()
 
 ##############################################################################
-# Sample the 3D data around each node of the mesh
-# -----------------------------------------------
+# Use mesh curvature to display useful anatomical information
+# on inflated meshes
+#
+# Here, we load the curvature map of the hemisphere under study,
+# and define a surface map whose value for a given vertex is
+# 1 if the curvature is positive, -1 if the curvature is negative.
 
+import numpy as np
 from nilearn import surface
 
+curv_right = surface.load_surf_data(fsaverage.curv_right)
+curv_right_sign = np.sign(curv_right)
+
+##############################################################################
+# Sample the 3D data around each node of the mesh
+# -----------------------------------------------
+
 texture = surface.vol_to_surf(stat_img, fsaverage.pial_right)
 
 ##############################################################################
@@ -40,14 +52,14 @@
 #
 # You can visualize the texture on the surface using the function
 # :func:`~nilearn.plotting.plot_surf_stat_map` which uses ``matplotlib``
-# as the default plotting engine:
+# as the default plotting engine.
 
 from nilearn import plotting
 
 fig = plotting.plot_surf_stat_map(
     fsaverage.infl_right, texture, hemi='right',
     title='Surface right hemisphere', colorbar=True,
-    threshold=1., bg_map=fsaverage.sulc_right
+    threshold=1., bg_map=curv_right_sign,
 )
 fig.show()
 
@@ -72,7 +84,7 @@
 fig = plotting.plot_surf_stat_map(
     fsaverage.infl_right, texture, hemi='right',
     title='Surface right hemisphere', colorbar=True,
-    threshold=1., bg_map=fsaverage.sulc_right,
+    threshold=1., bg_map=curv_right_sign, bg_on_data=True,
     engine=engine  # Specify the plotting engine here
 )
 fig.show()  # Display the figure as with matplotlib figures
@@ -104,8 +116,6 @@
 # Use an atlas and choose regions to outline
 # ------------------------------------------
 
-import numpy as np
-
 destrieux_atlas = datasets.fetch_atlas_surf_destrieux()
 parcellation = destrieux_atlas['map_right']
 

nilearn/_utils/docs.py

@@ -576,10 +576,12 @@ def custom_function(vertices):
 # bg_on_data
 docdict['bg_on_data'] = """
 bg_on_data : :obj:`bool`, optional
-    If ``True``, and a ``bg_map`` is specified,
+    If ``True`` and a ``bg_map`` is specified,
     the ``surf_data`` data is multiplied by the background
-    image, so that e.g. sulcal depth is visible beneath
-    the ``surf_data``.
+    image, so that e.g. sulcal depth is jointly visible with ``surf_data``.
+    Otherwise, the background image will only be visible where there
+    is no surface data (either because ``surf_data`` contains ``nan``\s
+    or because is was thresholded).
 
         .. note::
             This non-uniformly changes the surf_data values according

nilearn/plotting/html_surface.py

@@ -1,5 +1,6 @@
-import json
 import collections.abc
+import json
+import warnings
 
 import numpy as np
 import matplotlib as mpl
@@ -20,26 +21,86 @@ class SurfaceView(HTMLDocument):
     pass
 
 
+def _mix_colormaps(fg, bg):
+    """Mixes foreground and background arrays of RGBA colors.
+
+    Parameters
+    ----------
+    fg : numpy.ndarray
+        Array of shape (n, 4), foreground RGBA colors
+        represented as floats in [0, 1]
+    bg : numpy.ndarray
+        Array of shape (n, 4), background RGBA colors
+        represented as floats in [0, 1]
+
+    Returns
+    -------
+    mix : numpy.ndarray
+        Array of shape (n, 4), mixed colors
+        represented as floats in [0, 1]
+    """
+    # Adapted from https://stackoverflow.com/questions/726549/algorithm-for-additive-color-mixing-for-rgb-values/727339#727339 # noqa: E501
+    if fg.shape != bg.shape:
+        raise ValueError(
+            "Trying to mix colormaps with different shapes: "
+            f"{fg.shape}, {bg.shape}"
+        )
+
+    mix = np.empty_like(fg)
+
+    mix[:, 3] = 1 - (1 - fg[:, 3]) * (1 - bg[:, 3])
+
+    for color_index in range(0, 3):
+        mix[:, color_index] = (
+            fg[:, color_index] * fg[:, 3]
+            + bg[:, color_index] * bg[:, 3] * (1 - fg[:, 3])
+        ) / mix[:, 3]
+
+    return mix
+
+
 def _get_vertexcolor(surf_map, cmap, norm,
-                     absolute_threshold=None, bg_map=None):
-    vertexcolor = cmap(norm(surf_map).data)
-    if absolute_threshold is None:
-        return to_color_strings(vertexcolor)
+                     absolute_threshold=None, bg_map=None,
+                     bg_on_data=None, darkness=None):
     if bg_map is None:
-        bg_map = np.ones(len(surf_map)) * .5
+        bg_data = np.ones(len(surf_map)) * .5
         bg_vmin, bg_vmax = 0, 1
     else:
-        bg_map = surface.load_surf_data(bg_map)
-        bg_vmin, bg_vmax = np.min(bg_map), np.max(bg_map)
-    bg_norm = mpl.colors.Normalize(vmin=bg_vmin, vmax=bg_vmax)
-    bg_color = plt.get_cmap('Greys')(bg_norm(bg_map))
-    vertexcolor[np.abs(surf_map) < absolute_threshold] = bg_color[
-        np.abs(surf_map) < absolute_threshold]
-    return to_color_strings(vertexcolor)
+        bg_data = np.copy(surface.load_surf_data(bg_map))
+
+    # scale background map if need be
+    bg_vmin, bg_vmax = np.min(bg_data), np.max(bg_data)
+    if (bg_vmin < 0 or bg_vmax > 1):
+        bg_norm = mpl.colors.Normalize(vmin=bg_vmin, vmax=bg_vmax)
+        bg_data = bg_norm(bg_data)
+
+    if darkness is not None:
+        bg_data *= darkness
+
+    bg_colors = plt.get_cmap('Greys')(bg_data)
+
+    # select vertices which are filtered out by the threshold
+    if absolute_threshold is None:
+        under_threshold = np.zeros_like(surf_map, dtype=bool)
+    else:
+        under_threshold = np.abs(surf_map) < absolute_threshold
+
+    surf_colors = cmap(norm(surf_map).data)
+    # set transparency of voxels under threshold to 0
+    surf_colors[under_threshold, 3] = 0
+    if bg_on_data:
+        # if need be, set transparency of voxels above threshold to 0.7
+        # so that background map becomes visible
+        surf_colors[~under_threshold, 3] = 0.7
+
+    vertex_colors = _mix_colormaps(surf_colors, bg_colors)
+
+    return to_color_strings(vertex_colors)
 
 
 def one_mesh_info(surf_map, surf_mesh, threshold=None, cmap=cm.cold_hot,
                   black_bg=False, bg_map=None, symmetric_cmap=True,
+                  bg_on_data=False, darkness=.7,
                   vmax=None, vmin=None):
     """Prepare info for plotting one surface map on a single mesh.
 
@@ -55,7 +116,9 @@ def one_mesh_info(surf_map, surf_mesh, threshold=None, cmap=cm.cold_hot,
     info['inflated_left'] = mesh_to_plotly(surf_mesh)
     info['vertexcolor_left'] = _get_vertexcolor(
         surf_map, colors['cmap'], colors['norm'],
-        colors['abs_threshold'], bg_map)
+        absolute_threshold=colors['abs_threshold'], bg_map=bg_map,
+        bg_on_data=bg_on_data, darkness=darkness,
+    )
     info["cmin"], info["cmax"] = float(colors['vmin']), float(colors['vmax'])
     info['black_bg'] = black_bg
     info['full_brain_mesh'] = False
@@ -80,6 +143,7 @@ def _check_mesh(mesh):
 
 def full_brain_info(volume_img, mesh='fsaverage5', threshold=None,
                     cmap=cm.cold_hot, black_bg=False, symmetric_cmap=True,
+                    bg_on_data=False, darkness=.7,
                     vmax=None, vmin=None, vol_to_surf_kwargs={}):
     """Project 3D map on cortex; prepare info to plot both hemispheres.
 
@@ -101,15 +165,19 @@ def full_brain_info(volume_img, mesh='fsaverage5', threshold=None,
         symmetric_cmap=symmetric_cmap, vmax=vmax, vmin=vmin)
 
     for hemi, surf_map in surface_maps.items():
-        bg_map = surface.load_surf_data(mesh['sulc_{}'.format(hemi)])
+        curv_map = surface.load_surf_data(mesh["curv_{}".format(hemi)])
+        bg_map = np.sign(curv_map)
+
         info['pial_{}'.format(hemi)] = mesh_to_plotly(
             mesh['pial_{}'.format(hemi)])
         info['inflated_{}'.format(hemi)] = mesh_to_plotly(
             mesh['infl_{}'.format(hemi)])
 
         info['vertexcolor_{}'.format(hemi)] = _get_vertexcolor(
             surf_map, colors['cmap'], colors['norm'],
-            colors['abs_threshold'], bg_map)
+            absolute_threshold=colors['abs_threshold'], bg_map=bg_map,
+            bg_on_data=bg_on_data, darkness=darkness,
+        )
     info["cmin"], info["cmax"] = float(colors['vmin']), float(colors['vmax'])
     info['black_bg'] = black_bg
     info['full_brain_mesh'] = True
@@ -130,6 +198,7 @@ def _fill_html_template(info, embed_js=True):
 def view_img_on_surf(stat_map_img, surf_mesh='fsaverage5',
                      threshold=None, cmap=cm.cold_hot,
                      black_bg=False, vmax=None, vmin=None, symmetric_cmap=True,
+                     bg_on_data=False, darkness=.7,
                      colorbar=True, colorbar_height=.5, colorbar_fontsize=25,
                      title=None, title_fontsize=25, vol_to_surf_kwargs={}):
     """Insert a surface plot of a statistical map into an HTML page.
@@ -164,6 +233,12 @@ def view_img_on_surf(stat_map_img, surf_mesh='fsaverage5',
         If True, image is plotted on a black background. Otherwise on a
         white background. Default=False.
 
+    %(bg_on_data)s
+        Default=False.
+
+    %(darkness)s
+        Default=1.
+
     vmax : float or None, optional
         upper bound for the colorbar. if None, use the absolute max of the
         brain map.
@@ -221,7 +296,9 @@ def view_img_on_surf(stat_map_img, surf_mesh='fsaverage5',
     info = full_brain_info(
         volume_img=stat_map_img, mesh=surf_mesh, threshold=threshold,
         cmap=cmap, black_bg=black_bg, vmax=vmax, vmin=vmin,
-        symmetric_cmap=symmetric_cmap, vol_to_surf_kwargs=vol_to_surf_kwargs)
+        bg_on_data=bg_on_data, darkness=darkness,
+        symmetric_cmap=symmetric_cmap, vol_to_surf_kwargs=vol_to_surf_kwargs
+    )
     info['colorbar'] = colorbar
     info['cbar_height'] = colorbar_height
     info['cbar_fontsize'] = colorbar_fontsize
@@ -232,8 +309,9 @@ def view_img_on_surf(stat_map_img, surf_mesh='fsaverage5',
 
 def view_surf(surf_mesh, surf_map=None, bg_map=None, threshold=None,
               cmap=cm.cold_hot, black_bg=False, vmax=None, vmin=None,
-              symmetric_cmap=True, colorbar=True, colorbar_height=.5,
-              colorbar_fontsize=25, title=None, title_fontsize=25):
+              bg_on_data=False, darkness=.7, symmetric_cmap=True,
+              colorbar=True, colorbar_height=.5, colorbar_fontsize=25,
+              title=None, title_fontsize=25):
     """Insert a surface plot of a surface map into an HTML page.
 
     Parameters
@@ -252,10 +330,20 @@ def view_surf(surf_mesh, surf_map=None, bg_map=None, threshold=None,
         .thickness, .area, .curv, .sulc, .annot, .label) or
         a Numpy array
 
-    bg_map : Surface data, optional
+    bg_map : str or numpy.ndarray, optional
         Background image to be plotted on the mesh underneath the
         surf_data in greyscale, most likely a sulcal depth map for
         realistic shading.
+        If the map contains values outside [0, 1], it will be
+        rescaled such that all values are in [0, 1]. Otherwise,
+        it will not be modified.
+        Default=None.
+
+    %(bg_on_data)s
+        Default=False.
+
+    %(darkness)s
+        Default=1.
 
     threshold : str, number or None, optional
         If None, no thresholding.
@@ -329,6 +417,7 @@ def view_surf(surf_mesh, surf_map=None, bg_map=None, threshold=None,
     info = one_mesh_info(
         surf_map=surf_map, surf_mesh=surf_mesh, threshold=threshold,
         cmap=cmap, black_bg=black_bg, bg_map=bg_map,
+        bg_on_data=bg_on_data, darkness=darkness,
         symmetric_cmap=symmetric_cmap, vmax=vmax, vmin=vmin)
     info['colorbar'] = colorbar
     info['cbar_height'] = colorbar_height