ENH: Add CIFTI surface plot (#663)

* ENH: Add CIFTI surface plot

* DEP: Add surfplot to installation requirements

* DEP: Bump matplotlib minimum and remove compatibility pin

* ENH: Plot average BOLD signal

* ENH: Add helper function / pytest fixture for data creation

* TST: Add test for CIFTI plot

* STY: Placate flake8

* FIX: Operate within tmp directory

* MAINT: Bump templateflow requirement to reflect new templates

* TST: Allow pytest to invoke xvfb

* MAINT: Add `tests` as installation extra

* DOC: Add intersphinx source for surfplot

* ENH: Add parameter to control data clipping

* Apply suggestions from code review

Co-authored-by: Oscar Esteban <[email protected]>

* FIX: Leverage existing numpy function; set default to clip at 0

* Update niworkflows/viz/plots.py

Co-authored-by: Chris Markiewicz <[email protected]>

* Update niworkflows/viz/plots.py

Co-authored-by: Chris Markiewicz <[email protected]>

Co-authored-by: Oscar Esteban <[email protected]>
Co-authored-by: Chris Markiewicz <[email protected]>

Author: 3 people

docs/conf.py

@@ -242,6 +242,7 @@
     "python": ("https://docs.python.org/3/", None),
     "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
     "smriprep": ("https://www.nipreps.org/smriprep/", None),
+    "surfplot": ("https://surfplot.readthedocs.io/en/latest/", None),
     "templateflow": ("https://www.templateflow.org/python-client", None),
 }
 

niworkflows/tests/generate_data.py

@@ -0,0 +1,70 @@
+from pathlib import Path
+
+import numpy as np
+
+
+def _create_dtseries_cifti(timepoints, models):
+    """Create a dense timeseries CIFTI-2 file"""
+    import nibabel.cifti2 as ci
+
+    def create_series_map():
+        return ci.Cifti2MatrixIndicesMap(
+            (0,),
+            'CIFTI_INDEX_TYPE_SERIES',
+            number_of_series_points=timepoints,
+            series_exponent=0,
+            series_start=0,
+            series_step=1,
+            series_unit='SECOND'
+        )
+
+    def create_geometry_map():
+        index_offset = 0
+        brain_models = []
+        timeseries = np.zeros((timepoints, 0))
+
+        for name, data in models:
+            if "CORTEX" in name:
+                model_type = "CIFTI_MODEL_TYPE_SURFACE"
+                attr = "vertex_indices"
+                indices = ci.Cifti2VertexIndices(np.arange(len(data)))
+            else:
+                model_type = "CIFTI_MODEL_TYPE_VOXELS"
+                attr = "voxel_indices_ijk"
+                indices = ci.Cifti2VoxelIndicesIJK(np.arange(len(data)))
+            bm = ci.Cifti2BrainModel(
+                index_offset=index_offset,
+                index_count=len(data),
+                model_type=model_type,
+                brain_structure=name,
+            )
+            setattr(bm, attr, indices)
+            index_offset += len(data)
+            brain_models.append(bm)
+            timeseries = np.column_stack((timeseries, data.T))
+
+        brain_models.append(
+            ci.Cifti2Volume(
+                (4, 4, 4),
+                ci.Cifti2TransformationMatrixVoxelIndicesIJKtoXYZ(-3, np.eye(4)),
+            )
+        )
+
+        return ci.Cifti2MatrixIndicesMap(
+            (1,),
+            "CIFTI_INDEX_TYPE_BRAIN_MODELS",
+            maps=brain_models,
+        ), timeseries
+
+    matrix = ci.Cifti2Matrix()
+    series_map = create_series_map()
+    geometry_map, ts = create_geometry_map()
+    matrix.append(series_map)
+    matrix.append(geometry_map)
+    hdr = ci.Cifti2Header(matrix)
+    img = ci.Cifti2Image(dataobj=ts, header=hdr)
+    img.nifti_header.set_intent("NIFTI_INTENT_CONNECTIVITY_DENSE_SERIES")
+
+    out_file = Path("test.dtseries.nii").absolute()
+    ci.save(img, out_file)
+    return out_file

niworkflows/tests/test_viz.py

@@ -22,11 +22,15 @@
 #
 """Test viz module"""
 import os
+from pathlib import Path
+
 import numpy as np
 import nibabel as nb
-from .. import viz
+import pytest
+
 from .conftest import datadir
-from pathlib import Path
+from .generate_data import _create_dtseries_cifti
+from .. import viz
 
 
 def test_carpetplot():
@@ -149,3 +153,25 @@ def test_compcor_variance_plot(tmp_path):
     out_file = str(out_dir / "variance_plot_short.svg")
     metadata_file = os.path.join(datadir, "confounds_metadata_short_test.tsv")
     viz.plots.compcor_variance_plot([metadata_file], output_file=out_file)
+
+
+@pytest.fixture
+def create_surface_dtseries():
+    """Create a dense timeseries CIFTI-2 file with only cortex structures"""
+    out_file = _create_dtseries_cifti(
+        timepoints=10,
+        models=[
+            ('CIFTI_STRUCTURE_CORTEX_LEFT', np.random.rand(29696, 10)),
+            ('CIFTI_STRUCTURE_CORTEX_RIGHT', np.random.rand(29716, 10)),
+        ],
+    )
+    yield str(out_file)
+    out_file.unlink()
+
+
+def test_cifti_surfaces_plot(tmp_path, create_surface_dtseries):
+    """Test plotting CIFTI-2 surfaces"""
+    os.chdir(tmp_path)
+    out_dir = Path(os.getenv("SAVE_CIRCLE_ARTIFACTS", str(tmp_path)))
+    out_file = str(out_dir / "cifti_surfaces_plot.svg")
+    viz.plots.cifti_surfaces_plot(create_surface_dtseries, output_file=out_file)

niworkflows/viz/plots.py

@@ -1051,6 +1051,102 @@ def confounds_correlation_plot(
     return [ax0, ax1], gs
 
 
+def cifti_surfaces_plot(
+    in_cifti,
+    density="32k",
+    surface_type="inflated",
+    clip_range=(0, None),
+    output_file=None,
+    **splt_kwargs
+):
+    """
+    Plots a CIFTI-2 dense timeseries onto left/right mesh surfaces.
+
+    Parameters
+    ----------
+    in_cifti : str
+        CIFTI-2 dense timeseries (.dtseries.nii)
+    density : str
+        Surface density
+    surface_type : str
+        Inflation level of mesh surfaces. Supported: midthickness, inflated, veryinflated
+    clip_range : tuple or None
+        Range to clip `in_cifti` data prior to plotting.
+        If not None, two values must be provided as lower and upper bounds.
+        If values are None, no clipping is performed for that bound.
+    output_file: :obj:`str` or :obj:`None`
+        Path where the output figure should be saved. If this is not defined,
+        then the figure will be returned.
+    splt_kwargs : dict
+        Keyword arguments for :obj:`surfplot.Plot`
+
+    Outputs
+    -------
+    figure : matplotlib.pyplot.figure
+        Surface plot figure. Returned only if ``output_file`` is ``None``.
+    output_file: :obj:`str`
+        The file where the figure is saved.
+    """
+    import surfplot as splt
+    from surfplot.utils import add_fslr_medial_wall
+
+    def get_surface_meshes(density, surface_type):
+        import templateflow.api as tf
+        lh, rh = tf.get("fsLR", density=density, suffix=surface_type, extension=[".surf.gii"])
+        return str(lh), str(rh)
+
+    if density != "32k":
+        raise NotImplementedError("Only 32k density is currently supported.")
+
+    img = nb.cifti2.load(in_cifti)
+    if img.nifti_header.get_intent()[0] != "ConnDenseSeries":
+        raise TypeError(f"{in_cifti} is not a dense timeseries CIFTI file")
+
+    geo = img.header.get_index_map(1)
+    left_cortex, right_cortex = None, None
+    for bm in geo.brain_models:
+        if bm.brain_structure == "CIFTI_STRUCTURE_CORTEX_LEFT":
+            left_cortex = bm
+        elif bm.brain_structure == "CIFTI_STRUCTURE_CORTEX_RIGHT":
+            right_cortex = bm
+
+    if left_cortex is None or right_cortex is None:
+        raise RuntimeError("CIFTI is missing cortex information")
+
+    # calculate an average of the BOLD data, excluding the first 5 volumes
+    # as potential nonsteady states
+    data = img.dataobj[5:20].mean(axis=0)
+
+    cortex_data = _concat_brain_struct_data((left_cortex, right_cortex), data)
+    if density == "32k" and len(cortex_data) != 59412:
+        raise ValueError("Cortex data is not in fsLR space")
+    # medial wall needs to be added back in
+    cortex_data = add_fslr_medial_wall(cortex_data)
+    if clip_range:
+        cortex_data = np.clip(cortex_data, clip_range[0], clip_range[1], out=cortex_data)
+
+    lh_data, rh_data = np.array_split(cortex_data, 2)
+
+    # Build the figure
+    lh_mesh, rh_mesh = get_surface_meshes(density, surface_type)
+    p = splt.Plot(
+        surf_lh=lh_mesh,
+        surf_rh=rh_mesh,
+        layout=splt_kwargs.pop("layout", "row"),
+        **splt_kwargs
+    )
+    p.add_layer({'left': lh_data, 'right': rh_data}, cmap='YlOrRd_r')
+    figure = p.build()  # figsize - leave default?
+
+    if output_file is not None:
+        figure.savefig(output_file, bbox_inches="tight")
+        plt.close(figure)
+        figure = None
+        return output_file
+
+    return figure
+
+
 def _get_tr(img):
     """
     Attempt to extract repetition time from NIfTI/CIFTI header
@@ -1093,3 +1189,12 @@ def _decimate_data(data, seg, size):
     if t_dec:
         data = data[:, ::t_dec]
     return data, seg
+
+
+def _concat_brain_struct_data(structs, data):
+    concat_data = np.array([], dtype=data.dtype)
+    for struct in structs:
+        struct_upper_bound = struct.index_offset + struct.index_count
+        struct_data = data[struct.index_offset:struct_upper_bound]
+        concat_data = np.concatenate((concat_data, struct_data))
+    return concat_data

setup.cfg

@@ -27,7 +27,7 @@ python_requires = >= 3.7
 install_requires =
     attrs
     jinja2
-    matplotlib ~= 3.3, >= 3.3.1
+    matplotlib >= 3.4.2
     nibabel ~= 3.0
     nilearn >= 0.5.2
     nipype ~= 1.5
@@ -40,14 +40,16 @@ install_requires =
     scikit-image
     scipy
     seaborn
+    surfplot
     svgutils >= 0.3.4
     transforms3d
-    templateflow >= 0.6
+    templateflow >= 0.7.2
 test_requires =
     coverage >=5.2.1
     pytest >= 4.4
     pytest-cov
     pytest-xdist >= 1.28
+    pytest-xvfb
 packages = find:
 zip_safe = true
 
@@ -81,10 +83,17 @@ pointclouds =
     pyntcloud
 style =
     flake8 >= 3.7.0
+tests =
+    coverage >=5.2.1
+    pytest >= 4.4
+    pytest-cov
+    pytest-xdist >= 1.28
+    pytest-xvfb
 all =
     %(doc)s
     %(pointclouds)s
     %(style)s
+    %(tests)s
 
 [versioneer]
 VCS = git