Merge pull request #688 from nipreps/enh/carpetplot-improvements

ENH: Miscellaneous improvements to carpetplot

Author: oesteban

niworkflows/interfaces/plotting.py

@@ -21,7 +21,6 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Visualization tools."""
-from collections import defaultdict
 import numpy as np
 import nibabel as nb
 
@@ -34,7 +33,12 @@
     traits,
     isdefined,
 )
-from ..viz.plots import fMRIPlot, compcor_variance_plot, confounds_correlation_plot
+from niworkflows.utils.timeseries import _cifti_timeseries, _nifti_timeseries
+from niworkflows.viz.plots import (
+    fMRIPlot,
+    compcor_variance_plot,
+    confounds_correlation_plot,
+)
 
 
 class _FMRISummaryInputSpec(BaseInterfaceInputSpec):
@@ -101,7 +105,7 @@ def _run_interface(self, runtime):
             ),
             tr=(
                 self.inputs.tr if isdefined(self.inputs.tr) else
-                _get_tr(self.inputs.in_func)
+                _get_tr(input_data)
             ),
             confounds=dataframe,
             units={"outliers": "%", "FD": "mm"},
@@ -223,59 +227,6 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def _cifti_timeseries(dataset):
-    """Extract timeseries from CIFTI2 dataset."""
-    dataset = nb.load(dataset) if isinstance(dataset, str) else dataset
-
-    if dataset.nifti_header.get_intent()[0] != "ConnDenseSeries":
-        raise ValueError("Not a dense timeseries")
-
-    matrix = dataset.header.matrix
-    seg = defaultdict(list)
-    for bm in matrix.get_index_map(1).brain_models:
-        label = bm.brain_structure.replace("CIFTI_STRUCTURE_", "").replace("_", " ").title()
-        if "CORTEX" not in bm.brain_structure and "CEREBELLUM" not in bm.brain_structure:
-            label = "Other"
-
-        seg[label] += list(range(
-            bm.index_offset, bm.index_offset + bm.index_count
-        ))
-
-    return dataset.get_fdata(dtype="float32").T, seg
-
-
-def _nifti_timeseries(
-    dataset,
-    segmentation=None,
-    lut=None,
-    labels=("CSF", "WM", "Cerebellum", "Cortex")
-):
-    """Extract timeseries from NIfTI1/2 datasets."""
-    dataset = nb.load(dataset) if isinstance(dataset, str) else dataset
-    data = dataset.get_fdata(dtype="float32").reshape((-1, dataset.shape[-1]))
-
-    if segmentation is None:
-        return data, None
-
-    segmentation = nb.load(segmentation) if isinstance(segmentation, str) else segmentation
-    # Map segmentation
-    if lut is None:
-        lut = np.zeros((256,), dtype="int")
-        lut[1:11] = 4
-        lut[255] = 3
-        lut[30:99] = 2
-        lut[100:201] = 1
-    # Apply lookup table
-    seg = lut[np.asanyarray(segmentation.dataobj, dtype=int)].reshape(-1)
-    fgmask = seg > 0
-    seg = seg[fgmask]
-    seg_dict = {}
-    for i in np.unique(seg):
-        seg_dict[labels[i - 1]] = np.argwhere(seg == i).squeeze()
-
-    return data[fgmask], seg_dict
-
-
 def _get_tr(img):
     """
     Attempt to extract repetition time from NIfTI/CIFTI header

niworkflows/interfaces/tests/test_plotting.py

@@ -24,11 +24,8 @@
 import os
 import nibabel as nb
 from niworkflows import viz
-from niworkflows.interfaces.plotting import (
-    _get_tr,
-    _cifti_timeseries,
-    _nifti_timeseries,
-)
+from niworkflows.utils.timeseries import _cifti_timeseries, _nifti_timeseries
+from niworkflows.interfaces.plotting import _get_tr
 from niworkflows.tests.conftest import datadir
 
 

niworkflows/tests/test_viz.py

@@ -33,11 +33,8 @@
 from .generate_data import _create_dtseries_cifti
 from .. import viz
 from niworkflows.viz.plots import fMRIPlot
-from niworkflows.interfaces.plotting import (
-    _cifti_timeseries,
-    _nifti_timeseries,
-    _get_tr,
-)
+from niworkflows.utils.timeseries import _cifti_timeseries, _nifti_timeseries
+from niworkflows.interfaces.plotting import _get_tr
 
 
 @pytest.mark.parametrize("tr", (None, 0.7))

niworkflows/utils/timeseries.py

@@ -0,0 +1,87 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2022 The NiPreps Developers <nipreps@gmail.com>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Extracting signals from NIfTI and CIFTI2 files."""
+import numpy as np
+import nibabel as nb
+
+
+def _cifti_timeseries(dataset):
+    """Extract timeseries from CIFTI2 dataset."""
+    dataset = nb.load(dataset) if isinstance(dataset, str) else dataset
+
+    if dataset.nifti_header.get_intent()[0] != "ConnDenseSeries":
+        raise ValueError("Not a dense timeseries")
+
+    matrix = dataset.header.matrix
+    labels = {
+        "CIFTI_STRUCTURE_CORTEX_LEFT": "CtxL",
+        "CIFTI_STRUCTURE_CORTEX_RIGHT": "CtxR",
+        "CIFTI_STRUCTURE_CEREBELLUM_LEFT": "CbL",
+        "CIFTI_STRUCTURE_CEREBELLUM_RIGHT": "CbR",
+    }
+    seg = {label: [] for label in list(labels.values()) + ["Other"]}
+    for bm in matrix.get_index_map(1).brain_models:
+        label = (
+            "Other" if bm.brain_structure not in labels else
+            labels[bm.brain_structure]
+        )
+        seg[label] += list(range(
+            bm.index_offset, bm.index_offset + bm.index_count
+        ))
+
+    return dataset.get_fdata(dtype="float32").T, seg
+
+
+def _nifti_timeseries(
+    dataset,
+    segmentation=None,
+    labels=("Ctx GM", "dGM", "WM+CSF", "Cb"),
+    remap_rois=True,
+    lut=None,
+):
+    """Extract timeseries from NIfTI1/2 datasets."""
+    dataset = nb.load(dataset) if isinstance(dataset, str) else dataset
+    data = dataset.get_fdata(dtype="float32").reshape((-1, dataset.shape[-1]))
+
+    if segmentation is None:
+        return data, None
+
+    segmentation = nb.load(segmentation) if isinstance(segmentation, str) else segmentation
+    # Map segmentation
+    if remap_rois or lut is not None:
+        if lut is None:
+            lut = np.zeros((256,), dtype="int")
+            lut[100:201] = 1  # Ctx GM
+            lut[30:99] = 2    # dGM
+            lut[1:11] = 3     # WM+CSF
+            lut[255] = 4      # Cerebellum
+        # Apply lookup table
+        segmentation = lut[np.asanyarray(segmentation.dataobj, dtype=int)].reshape(-1)
+
+    fgmask = segmentation > 0
+    segmentation = segmentation[fgmask]
+    seg_dict = {}
+    for i in np.unique(segmentation):
+        seg_dict[labels[i - 1]] = np.argwhere(segmentation == i).squeeze()
+
+    return data[fgmask], seg_dict

niworkflows/viz/plots.py

@@ -29,7 +29,7 @@
 import matplotlib.pyplot as plt
 from matplotlib import gridspec as mgs
 import matplotlib.cm as cm
-from matplotlib.colors import ListedColormap, Normalize
+from matplotlib.colors import Normalize
 from matplotlib.colorbar import ColorbarBase
 
 DINA4_LANDSCAPE = (11.69, 8.27)
@@ -192,13 +192,15 @@ def plot_carpet(
     """
     if segments is None:
         segments = {
-            "brain": list(range(data.shape[0]))
+            "whole brain (voxels)": list(range(data.shape[0]))
         }
 
     if cmap is None:
-        cmap = ListedColormap(cm.get_cmap("tab10").colors[:len(segments)])
+        colors = cm.get_cmap("tab10").colors
     elif cmap == "paired":
-        cmap = ListedColormap([cm.get_cmap("Paired").colors[i] for i in (1, 0, 7, 3)])
+        colors = list(cm.get_cmap("Paired").colors)
+        colors[0], colors[1] = colors[1], colors[0]
+        colors[2], colors[7] = colors[7], colors[2]
 
     vminmax = (None, None)
     if detrend:
@@ -239,76 +241,75 @@ def plot_carpet(
     # Length before decimation
     n_trs = data.shape[-1] - drop_trs
 
-    # Define nested GridSpec
-    gs = mgs.GridSpecFromSubplotSpec(
-        1, 2, subplot_spec=subplot, width_ratios=(1, 100), wspace=0.0,
-    )
-
-    # Segmentation colorbar
-    colors = np.hstack([[i + 1] * len(v) for i, v in enumerate(segments.values())])
-    ax0 = plt.subplot(gs[0])
-    ax0.set_yticks([])
-    ax0.set_xticks([])
-    ax0.imshow(
-        colors[:, np.newaxis],
-        interpolation="none",
-        aspect="auto",
-        cmap=cmap
-    )
-
-    ax0.grid(False)
-    ax0.spines["left"].set_visible(False)
-    ax0.spines["bottom"].set_color("none")
-    ax0.spines["bottom"].set_visible(False)
-
     # Calculate time decimation factor
     t_dec = max(int((1.8 * n_trs) // size[1]), 1)
-    data = data[np.hstack(list(segments.values())), drop_trs::t_dec]
-
-    # Carpet plot
-    ax1 = plt.subplot(gs[1])
-    ax1.imshow(
-        data,
-        interpolation="nearest",
-        aspect="auto",
-        cmap="gray",
-        vmin=vminmax[0],
-        vmax=vminmax[1],
-    )
+    data = data[:, drop_trs::t_dec]
 
-    ax1.grid(False)
-    ax1.set_yticks([])
-    ax1.set_yticklabels([])
-
-    xticks = np.linspace(0, data.shape[-1], endpoint=True, num=7)
-    xlabel = "time-points (index)"
-    xticklabels = (xticks * n_trs / data.shape[-1]).astype("uint32") + drop_trs
-    if tr is not None:
-        xlabel = "time (mm:ss)"
-        xticklabels = [
-            f"{int(t // 60):02d}:{(t % 60).round(0).astype(int):02d}"
-            for t in (tr * xticklabels)
-        ]
+    nsegments = len(segments)
+
+    # Define nested GridSpec
+    gs = mgs.GridSpecFromSubplotSpec(
+        nsegments,
+        1,
+        subplot_spec=subplot,
+        hspace=0.05,
+        height_ratios=[len(v) for v in segments.values()]
+    )
 
-    ax1.set_xticks(xticks)
-    ax1.set_xlabel(xlabel)
-    ax1.set_xticklabels(xticklabels)
+    for i, (label, indices) in enumerate(segments.items()):
+        # Carpet plot
+        ax = plt.subplot(gs[i])
+
+        ax.imshow(
+            data[indices, :],
+            interpolation="nearest",
+            aspect="auto",
+            cmap="gray",
+            vmin=vminmax[0],
+            vmax=vminmax[1],
+        )
 
-    # Remove and redefine spines
-    for side in ["top", "right"]:
         # Toggle the spine objects
-        ax0.spines[side].set_color("none")
-        ax0.spines[side].set_visible(False)
-        ax1.spines[side].set_color("none")
-        ax1.spines[side].set_visible(False)
-
-    ax1.yaxis.set_ticks_position("left")
-    ax1.xaxis.set_ticks_position("bottom")
-    ax1.spines["bottom"].set_visible(False)
-    ax1.spines["left"].set_color("none")
-    ax1.spines["left"].set_visible(False)
-    if title:
-        ax1.set_title(title)
+        ax.spines["top"].set_color("none")
+        ax.spines["top"].set_visible(False)
+        ax.spines["right"].set_color("none")
+        ax.spines["right"].set_visible(False)
+
+        # Make colored left axis
+        ax.spines["left"].set_linewidth(3)
+        ax.spines["left"].set_color(colors[i])
+        ax.spines["left"].set_capstyle("butt")
+        ax.spines["left"].set_position(("outward", 2))
+
+        # Make all subplots have same xticks
+        xticks = np.linspace(0, data.shape[-1], endpoint=True, num=7)
+        ax.set_xticks(xticks)
+        ax.set_yticks([])
+        ax.set_ylabel(label)
+        ax.grid(False)
+
+        if i < (nsegments - 1):
+            ax.spines["bottom"].set_color("none")
+            ax.spines["bottom"].set_visible(False)
+            ax.set_xticklabels([])
+        else:
+            xlabel = "time-points (index)"
+            xticklabels = (xticks * n_trs / data.shape[-1]).astype("uint32") + drop_trs
+            if tr is not None:
+                xlabel = "time (mm:ss)"
+                xticklabels = [
+                    f"{int(t // 60):02d}:{(t % 60).round(0).astype(int):02d}"
+                    for t in (tr * xticklabels)
+                ]
+
+            ax.set_xlabel(xlabel)
+            ax.set_xticklabels(xticklabels)
+            ax.spines["bottom"].set_position(("outward", 5))
+            ax.spines["bottom"].set_color("k")
+            ax.spines["bottom"].set_linewidth(.8)
+
+        if title and i == 0:
+            ax.set_title(title)
 
     if output_file is not None:
         figure = plt.gcf()
@@ -317,7 +318,7 @@ def plot_carpet(
         figure = None
         return output_file
 
-    return (ax0, ax1), gs
+    return gs
 
 
 def spikesplot(