Merge pull request #687 from nipreps/fix/carpetplot-minimize-processing

ENH: Refactor carpetplot reasigning responsibilities

Author: oesteban

.circleci/config.yml

@@ -225,6 +225,7 @@ jobs:
             mkdir -p $PWD/artifacts $PWD/summaries
             docker run -u $( id -u ):$( id -g ) -it --rm -w /src/niworkflows \
               -e COVERAGE_FILE=/tmp/summaries/.pytest.coverage \
+              -e SAVE_CIRCLE_ARTIFACTS="/tmp/artifacts/" \
               -e TEST_DATA_HOME=/data -v /tmp/data:/data \
               -e FS_LICENSE=/etc/fslicense.txt \
               -v /tmp/fslicense/license.txt:/etc/fslicense.txt:ro \

niworkflows/interfaces/plotting.py

@@ -21,7 +21,9 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Visualization tools."""
+from collections import defaultdict
 import numpy as np
+import nibabel as nb
 
 from nipype.utils.filemanip import fname_presuffix
 from nipype.interfaces.base import (
@@ -37,14 +39,14 @@
 
 class _FMRISummaryInputSpec(BaseInterfaceInputSpec):
     in_func = File(exists=True, mandatory=True, desc="")
-    in_mask = File(exists=True, desc="")
+    in_spikes_bg = File(exists=True, desc="")
+    fd = File(exists=True, desc="")
+    dvars = File(exists=True, desc="")
+    outliers = File(exists=True, desc="")
     in_segm = File(exists=True, desc="")
-    in_spikes_bg = File(exists=True, mandatory=True, desc="")
-    fd = File(exists=True, mandatory=True, desc="")
-    fd_thres = traits.Float(0.2, usedefault=True, desc="")
-    dvars = File(exists=True, mandatory=True, desc="")
-    outliers = File(exists=True, mandatory=True, desc="")
     tr = traits.Either(None, traits.Float, usedefault=True, desc="the TR")
+    fd_thres = traits.Float(0.2, usedefault=True, desc="")
+    drop_trs = traits.Int(0, usedefault=True, desc="dummy scans")
 
 
 class _FMRISummaryOutputSpec(TraitedSpec):
@@ -67,28 +69,44 @@ def _run_interface(self, runtime):
             newpath=runtime.cwd,
         )
 
-        dataframe = pd.DataFrame(
-            {
-                "outliers": np.loadtxt(self.inputs.outliers, usecols=[0]).tolist(),
-                # Pick non-standardize dvars (col 1)
-                # First timepoint is NaN (difference)
-                "DVARS": [np.nan]
-                + np.loadtxt(self.inputs.dvars, skiprows=1, usecols=[1]).tolist(),
-                # First timepoint is zero (reference volume)
-                "FD": [0.0]
-                + np.loadtxt(self.inputs.fd, skiprows=1, usecols=[0]).tolist(),
-            }
+        dataframe = pd.DataFrame({
+            "outliers": np.loadtxt(self.inputs.outliers, usecols=[0]).tolist(),
+            # Pick non-standardize dvars (col 1)
+            # First timepoint is NaN (difference)
+            "DVARS": [np.nan]
+            + np.loadtxt(self.inputs.dvars, skiprows=1, usecols=[1]).tolist(),
+            # First timepoint is zero (reference volume)
+            "FD": [0.0]
+            + np.loadtxt(self.inputs.fd, skiprows=1, usecols=[0]).tolist(),
+        }) if (
+            isdefined(self.inputs.outliers)
+            and isdefined(self.inputs.dvars)
+            and isdefined(self.inputs.fd)
+        ) else None
+
+        input_data = nb.load(self.inputs.in_func)
+        seg_file = self.inputs.in_segm if isdefined(self.inputs.in_segm) else None
+        dataset, segments = (
+            _cifti_timeseries(input_data)
+            if isinstance(input_data, nb.Cifti2Image) else
+            _nifti_timeseries(input_data, seg_file)
         )
 
         fig = fMRIPlot(
-            self.inputs.in_func,
-            mask_file=self.inputs.in_mask if isdefined(self.inputs.in_mask) else None,
-            seg_file=self.inputs.in_segm if isdefined(self.inputs.in_segm) else None,
-            spikes_files=[self.inputs.in_spikes_bg],
-            tr=self.inputs.tr,
-            data=dataframe[["outliers", "DVARS", "FD"]],
+            dataset,
+            segments=segments,
+            spikes_files=(
+                [self.inputs.in_spikes_bg]
+                if isdefined(self.inputs.in_spikes_bg) else None
+            ),
+            tr=(
+                self.inputs.tr if isdefined(self.inputs.tr) else
+                _get_tr(self.inputs.in_func)
+            ),
+            confounds=dataframe,
             units={"outliers": "%", "FD": "mm"},
             vlines={"FD": [self.inputs.fd_thres]},
+            nskip=self.inputs.drop_trs,
         ).plot()
         fig.savefig(self._results["out_file"], bbox_inches="tight")
         return runtime
@@ -203,3 +221,78 @@ def _run_interface(self, runtime):
             reference=self.inputs.reference_column,
         )
         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
+
+    Examples
+    --------
+    >>> _get_tr(nb.load(Path(test_data) /
+    ...    'sub-ds205s03_task-functionallocalizer_run-01_bold_volreg.nii.gz'))
+    2.2
+    >>> _get_tr(nb.load(Path(test_data) /
+    ...    'sub-01_task-mixedgamblestask_run-02_space-fsLR_den-91k_bold.dtseries.nii'))
+    2.0
+
+    """
+
+    try:
+        return img.header.matrix.get_index_map(0).series_step
+    except AttributeError:
+        return img.header.get_zooms()[-1]
+    raise RuntimeError("Could not extract TR - unknown data structure type")

niworkflows/interfaces/tests/test_plotting.py

@@ -0,0 +1,81 @@
+# 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/
+#
+"""Tests plotting interfaces."""
+import os
+import nibabel as nb
+from niworkflows import viz
+from niworkflows.interfaces.plotting import (
+    _get_tr,
+    _cifti_timeseries,
+    _nifti_timeseries,
+)
+from niworkflows.tests.conftest import datadir
+
+
+def test_cifti_carpetplot():
+    """Exercise extraction of timeseries from CIFTI2."""
+    save_artifacts = os.getenv("SAVE_CIRCLE_ARTIFACTS", False)
+
+    cifti_file = os.path.join(
+        datadir,
+        "sub-01_task-mixedgamblestask_run-02_space-fsLR_den-91k_bold.dtseries.nii",
+    )
+    data, segments = _cifti_timeseries(cifti_file)
+    viz.plot_carpet(
+        data,
+        segments,
+        tr=_get_tr(nb.load(cifti_file)),
+        output_file=(
+            os.path.join(
+                save_artifacts, "carpetplot_cifti.svg"
+            ) if save_artifacts else None
+        ),
+        drop_trs=0,
+        cmap="paired",
+    )
+
+
+def test_nifti_carpetplot():
+    """Exercise extraction of timeseries from CIFTI2."""
+    save_artifacts = os.getenv("SAVE_CIRCLE_ARTIFACTS", False)
+
+    nifti_file = os.path.join(
+        datadir,
+        "sub-ds205s03_task-functionallocalizer_run-01_bold_volreg.nii.gz",
+    )
+    seg_file = os.path.join(
+        datadir,
+        "sub-ds205s03_task-functionallocalizer_run-01_bold_parc.nii.gz",
+    )
+    data, segments = _nifti_timeseries(nifti_file, seg_file)
+    viz.plot_carpet(
+        data,
+        segments,
+        tr=_get_tr(nb.load(nifti_file)),
+        output_file=(
+            os.path.join(
+                save_artifacts, "carpetplot_nifti.svg"
+            ) if save_artifacts else None
+        ),
+        drop_trs=0,
+    )