diff --git a/setup.cfg b/setup.cfg index f093ade8..c75129a8 100755 --- a/setup.cfg +++ b/setup.cfg @@ -30,6 +30,7 @@ install_requires = joblib nipype>= 1.5.1, < 2.0 nitransforms>=21.0.0 + nireports numpy>=1.17.3 nest-asyncio>=1.5.1 scikit-image>=0.14.2 diff --git a/src/eddymotion/viz.py b/src/eddymotion/viz.py index a46c777d..a298d092 100644 --- a/src/eddymotion/viz.py +++ b/src/eddymotion/viz.py @@ -23,6 +23,7 @@ """Visualization utilities.""" import nibabel as nb import numpy as np +from nireports.reportlets.nuisance import plot_carpet as nw_plot_carpet def plot_dwi(dataobj, affine, gradient=None, **kwargs): @@ -277,3 +278,139 @@ def plot_gradients( plt.suptitle(title) return ax + + +def plot_carpet( + nii, + gtab, + segmentation=None, + sort_by_bval=False, + output_file=None, + segment_labels=None, + detrend=False, +): + """ + Return carpet plot using niworkflows carpet_plot + + Parameters + ---------- + nii : Nifti1Image + DW imaging data + gtab : :obj:`GradientTable` + DW imaging data gradient data + segmentation : Nifti1Image + Boolean or segmentation mask of DW imaging data + sort_by_bval : :obj:`bool` + Flag to reorder time points by bvalue + output_file : :obj:`string` + Path to save the plot + 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]} + detrend : :obj:`bool` + niworkflows plot_carpet detrend flag + + Returns + --------- + matplotlib GridSpec object + """ + segments = None + + nii_data = nii.get_fdata() + + b0_data = nii_data[..., gtab.b0s_mask] + dw_data = nii_data[..., ~gtab.b0s_mask] + + bzero = np.mean(b0_data, -1) + + nii_data_div_b0 = dw_data / bzero[..., np.newaxis] + + sort_inds = np.argsort( + gtab.bvals[~gtab.b0s_mask] + if sort_by_bval + else np.arange(len(gtab.bvals[~gtab.b0s_mask])) + ) + nii_data_div_b0 = nii_data_div_b0[..., sort_inds] + + # Reshape + nii_data_reshaped = nii_data_div_b0.reshape(-1, nii_data_div_b0.shape[-1]) + + if segmentation is not None: + segmentation_data = np.asanyarray(segmentation.dataobj, dtype=np.int16) + + # Apply mask + segmentation_reshaped = segmentation_data.reshape(-1) + nii_data_masked = nii_data_reshaped[segmentation_reshaped > 0, :] + segmentation_masked = segmentation_reshaped[segmentation_reshaped > 0] + + if segment_labels is not None: + segments = dict() + labels = list(segment_labels.keys()) + for label in labels: + indices = np.array([], dtype=int) + for ii in segment_labels[label]: + indices = np.concatenate( + [indices, np.where(segmentation_masked == ii)[0]] + ) + segments[label] = indices + + else: + nii_data_masked = nii_data_reshaped + + bad_row_ind = np.where(~np.isfinite(nii_data_masked))[0] + + good_row_ind = np.ones(nii_data_masked.shape[0], dtype=bool) + good_row_ind[bad_row_ind] = False + + nii_data_masked = nii_data_masked[good_row_ind, :] + + # Plot + return nw_plot_carpet( + nii_data_masked, detrend=detrend, segments=segments, output_file=output_file + ) + + +def get_segment_labels( + filepath, keywords, delimiter=" ", index_position=0, label_position=1 +): + """ + Return segment labels for plot_carpet 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` + Delimiter between label index and label string in label file + (' ' for freesurfer label file) + index_position : :obj:`int` + Position of label index in label file + (0 for freesurfer label file) + label_position : :obj:`int` + 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 = dict() + + 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