ENH: Higher-level carpetplot tooling for DWI

nireports/reportlets/modality/dwi.py

@@ -433,3 +433,164 @@
     """
 
     return plot_raincloud(data_file, group_name, feature, **kwargs)
+
+
+def nii_to_carpetplot_data(
+    nii,
+    bvals=None,
+    divide_by_b0=True,
+    drop_b0=True,
+    sort_by_bval=False,
+    mask_nii=None,
+    segment_labels=None,
+):
+    """
+    Convert nii to data matrix for carpet plot
+
+    Parameters
+    ----------
+    nii : Nifti1Image
+        DW imaging data
+    bvals : :obj:`numpy.ndarray`, optional
+        Rounded bvals
+    divide_by_b0 : :obj:`bool`, optional
+        Divide data by mean b0
+    drop_b0 : :obj:`bool`, optional
+        Flag to drop b0 data
+    sort_by_bval : :obj:`bool`, optional
+        Flag to reorder time points by bvalue
+    mask_nii : Nifti1Image, optional
+        Boolean or segmentation mask of DW imaging data
+    segment_labels : :obj:`dict`, optional
+        Dictionary of segment labels
+        e.g. {'Cerebral_White_Matter': [2, 41],
+              'Cerebral_Cortex': [3, 42],
+              'Ventricle': [4, 14, 15, 43, 72]}
+
+    Returns
+    ---------
+    data : N x T :obj:`numpy.array`
+        The functional data to be plotted
+        (*N* sampling locations by *T* timepoints).
+    segments: :obj:`dict`
+        A mapping between segment labels (e.g., `"Left Cortex"`)
+        and list of indexes in the data array.
+    """
+
+    nii_data = nii.get_fdata()
+
+    if mask_nii is None:
+        mask_data = np.ones(nii_data.shape[0:3])
+    else:
+        mask_data = np.asanyarray(mask_nii.dataobj, dtype=np.int16)
+
+    if bvals is not None:
+        if divide_by_b0:
+            b0_data = nii_data[..., bvals == 0]
+            bzero = np.mean(b0_data, -1)
+            nii_data = nii_data / bzero[..., np.newaxis]
+
+        if drop_b0:
+            nii_data = nii_data[..., bvals > 0]
+            bvals = bvals[bvals > 0]
+
+        if sort_by_bval:
+            sort_inds = np.argsort(bvals)
+            nii_data = nii_data[..., sort_inds]
+
+    # Reshape
+    data = nii_data.reshape(-1, nii_data.shape[-1])
+    mask_data = mask_data.reshape(-1)
+
+    # Apply mask
+    data = data[mask_data > 0, :]
+    mask_data = mask_data[mask_data > 0]
+
+    # Remove bad rows
+    bad_row_ind = np.where(~np.isfinite(data))[0]
+    good_row_ind = np.ones(data.shape[0], dtype=bool)
+    good_row_ind[bad_row_ind] = False
+
+    data = data[good_row_ind, :]
+    mask_data = mask_data[good_row_ind]
+
+    # Get segments dict
+    if mask_nii is None or segment_labels is None:
+        segments = None
+    else:
+        segments = get_segments(mask_data, segment_labels)
+
+    return data, segments
+
+
+def get_segments(segment_mask, segment_labels):
+    """
+    Return segments dict for plot_carpet function
+
+    Parameters
+    ----------
+    segment_mask : :obj:`numpy.ndarray`
+        Segmentation mask of DW imaging data
+    segment_labels : :obj:`dict`
+        Dictionary of segment labels
+        e.g. {'Cerebral_White_Matter': [2, 41],
+              'Cerebral_Cortex': [3, 42],
+              'Ventricle': [4, 14, 15, 43, 72]}
+
+    Returns
+    ---------
+    segments: :obj:`dict`
+        A mapping between segment labels (e.g., `"Left Cortex"`)
+        and list of indexes in the data array.
+    """
+    segments = {}
+
+    for label, idx in segment_labels.items():
+        indices = np.array([], dtype=int)
+        for ii in idx:
+            indices = np.concatenate([indices, np.nonzero(segment_mask == ii)[0]])
+        segments[label] = indices
+
+    return segments
+
+
+def get_segment_labels(filepath, keywords, delimiter=" ", index_position=0, label_position=1):
+    """
+    Get segment labels from file by keyword for get_segments function
+
+    Parameters
+    ----------
+    filepath : :obj:`string`
+        Path to segment label text file, such as freesurfer label file
+    keywords : list of :obj:`string`
+        List of label keywords.
+        All labels containing the keyword will be grouped together.
+        e.g. ["Cerebral_White_Matter", "Cerebral_Cortex", "Ventricle"]
+    delimiter : :obj:`string`, optional
+        Delimiter between label index and label string in label file
+        (' ' for freesurfer label file)
+    index_position : :obj:`int`, optional
+        Position of label index in label file
+        (0 for freesurfer label file)
+    label_position : :obj:`int`, optional
+        Position of label string in label file
+        (1 for freesurfer label file)
+    Returns
+    ---------
+    dict
+    e.g. {'Cerebral_White_Matter': [2, 41],
+          'Cerebral_Cortex': [3, 42],
+          'Ventricle': [4, 14, 15, 43, 72]}
+    """
+    segment_labels = {}
+
+    with open(filepath, "r") as f:
+        labels = f.read()
+
+    labels_s = [label.split(delimiter) for label in labels.split("\n") if label != ""]
+
+    for keyword in keywords:
+        ind = [int(i[index_position]) for i in labels_s if keyword in i[label_position]]
+        segment_labels[keyword] = ind
+
+    return segment_labels

nireports/tests/data/aseg.auto_noCCseg.label_intensities.txt

@@ -0,0 +1,47 @@
+1 Left_Cerebral_Exterior 1.64 0.0  0
+2 Left_Cerebral_White_Matter 1.00 0.0 105
+3 Left_Cerebral_Cortex 1.20 0.0 69
+4 Left_Lateral_Ventricle 1.30 0.0 24
+5 Left_Inf_Lat_Vent 1.27 0.0 41
+7 Left_Cerebellum_White_Matter 1.08 0.0 91
+8 Left_Cerebellum_Cortex 1.32 0.0 73
+9 Left_Thalamus 1.26 0.0 103
+10 Left_Thalamus_Proper 1.03 0.0 93
+11 Left_Caudate 1.10 0.0 75
+12 Left_Putamen 1.16 0.0 94
+13 Left_Pallidum 1.05 0.0 101
+14 Third_Ventricle 1.64 0.0 42
+15 Fourth_Ventricle 1.72 0.0 26
+16 Brain_Stem 1.07 0.0 83
+17 Left_Hippocampus 1.27 0.0 72
+18 Left_Amygdala 1.32 0.0 76
+24 CSF 1.64 0.0 60
+26 Left_Accumbens_area 1.10 0.0 68
+28 Left_VentralDC 1.07 0.0 91
+40 Right_Cerebral_Exterior 1.64 0.0  0
+41 Right_Cerebral_White_Matter 1.00 0.0 103
+42 Right_Cerebral_Cortex 1.20 0.0 69
+43 Right_Lateral_Ventricle 1.94 0.0 25
+44 Right_Inf_Lat_Vent 1.29 0.0 30
+46 Right_Cerebellum_White_Matter 1.05 0.0 89
+47 Right_Cerebellum_Cortex 1.28 0.0 70
+48 Right_Thalamus 1.26 0.0  0
+49 Right_Thalamus_Proper 1.05 0.0 85
+50 Right_Caudate 1.24 0.0 81
+51 Right_Putamen 1.16 0.0 88
+52 Right_Pallidum 1.03 0.0 99
+53 Right_Hippocampus 1.29 0.0 72
+54 Right_Amygdala 1.29 0.0 72
+58 Right_Accumbens_area 1.24 0.0 80
+60 Right_VentralDC 1.05 0.0 94
+72 Fifth_Ventricle 1.64 0.0 50
+75  1.64 0.0  0
+76  1.64 0.0  0
+77 WM_hypointensities 1.00 0.0 77
+78 Left_WM_hypointensities 1.00 0.0  0
+79 Right_WM_hypointensities 1.00 0.0  0
+80 non_WM_hypointensities 1.00 0.0 44
+81 Left_non_WM_hypointensities 1.00 0.0  0
+82 Right_non_WM_hypointensities 1.00 0.0  0
+186  1.00 0.0  0
+187  1.00 0.0  0

nireports/tests/test_dwi.py

@@ -27,7 +27,14 @@
 import pytest
 from matplotlib import pyplot as plt
 
-from nireports.reportlets.modality.dwi import plot_dwi, plot_gradients, plot_tissue_values
+from nireports.reportlets.modality.dwi import (
+    get_segment_labels,
+    nii_to_carpetplot_data,
+    plot_dwi,
+    plot_gradients,
+    plot_tissue_values,
+)
+from nireports.reportlets.nuisance import plot_carpet
 from nireports.tests.utils import _generate_raincloud_random_data
 
 
@@ -109,3 +116,45 @@ def test_plot_tissue_values(tmp_path):
         mark_nans=mark_nans,
         output_file=output_file,
     )
+
+
+def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir):
+    """Check the nii to carpet plot data function"""
+
+    testdata_name = "ds000114_sub-01_ses-test_desc-trunc_dwi"
+
+    nii = nb.load(testdata_path / f"{testdata_name}.nii.gz")
+    bvals = np.loadtxt(testdata_path / f"{testdata_name}.bval")
+
+    mask_data = np.round(82 * np.random.rand(nii.shape[0], nii.shape[1], nii.shape[2]))
+
+    mask_nii = nb.Nifti1Image(mask_data, np.eye(4))
+
+    filepath = testdata_path / "aseg.auto_noCCseg.label_intensities.txt"
+    keywords = ["Cerebral_White_Matter", "Cerebral_Cortex", "Ventricle"]
+
+    segment_labels = get_segment_labels(filepath, keywords)
+
+    image_path = None
+
+    if outdir is not None:
+        image_path = outdir / f"{testdata_name}_nii_to_carpet.svg"
+
+    data, segments = nii_to_carpetplot_data(
+        nii, bvals=bvals, mask_nii=mask_nii, segment_labels=segment_labels
+    )
+
+    plot_carpet(data, segments, output_file=image_path)
+
+
+def test_get_segment_labels(tmp_path, testdata_path):
+    """Check the segment label function"""
+
+    testdata_name = "aseg.auto_noCCseg.label_intensities.txt"
+
+    filepath = testdata_path / testdata_name
+    keywords = ["Cerebral_White_Matter", "Cerebral_Cortex", "Ventricle"]
+
+    segment_labels = get_segment_labels(filepath, keywords)
+
+    assert segment_labels is not None