From 593a29fa3094b3b9af70ac3933cea1ba412ccc17 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 12:26:02 -0400 Subject: [PATCH 01/16] Adds higher level carpot_plot function --- nireports/reportlets/modality/dwi.py | 95 ++++++++++++++++++++++++++++ nireports/tests/test_dwi.py | 21 +++++- 2 files changed, 115 insertions(+), 1 deletion(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 117b16f5..6f7b19d7 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -28,6 +28,7 @@ from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import art3d from nilearn.plotting import plot_anat +from nireports.reportlets.nuisance import plot_carpet as nw_plot_carpet def plot_dwi(dataobj, affine, gradient=None, **kwargs): @@ -405,3 +406,97 @@ 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 + ) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 6e43fea4..3b3bd154 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -27,7 +27,10 @@ import pytest from matplotlib import pyplot as plt -from nireports.reportlets.modality.dwi import plot_dwi, plot_gradients +import os.path as op +import dipy.core.gradients as dpg +from nireports.reportlets.modality.dwi \ + import plot_dwi, plot_gradients, plot_carpet def test_plot_dwi(tmp_path, testdata_path, outdir): @@ -69,3 +72,19 @@ def test_plot_gradients(tmp_path, testdata_path, dwi_btable, outdir): if outdir is not None: plt.savefig(outdir / f"{dwi_btable}.svg", bbox_inches="tight") + + +def test_plot_carpet(tmp_path, testdata_path, outdir): + """Check the carpet plot""" + + testdata_name = "ds000114_sub-01_ses-test_desc-trunc_dwi" + + nii_path = testdata_path / f'{testdata_name}.nii.gz' + bvec_path = testdata_path / f'{testdata_name}.bvec' + bval_path = testdata_path / f'{testdata_name}.bval' + + nii = nb.load(nii_path) + gtab = dpg.gradient_table(bval_path, bvec_path) + image_path = outdir / f'{testdata_name}_carpet.png' + + plot_carpet(nii, gtab, output_file=image_path) From 1c2aaff1e00e3c7d09fd207238a5b78f4093ccb2 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 13:49:39 -0400 Subject: [PATCH 02/16] Adds default if outdir None --- nireports/tests/test_dwi.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 3b3bd154..3623301f 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -85,6 +85,10 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): nii = nb.load(nii_path) gtab = dpg.gradient_table(bval_path, bvec_path) - image_path = outdir / f'{testdata_name}_carpet.png' + + image_path = None + + if outdir is not None: + image_path = outdir / f'{testdata_name}_carpet.svg' plot_carpet(nii, gtab, output_file=image_path) From 7783b9c5c7b8761635ecd6b32985023570c81dfc Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 14:10:57 -0400 Subject: [PATCH 03/16] Removes dipy dependency --- nireports/reportlets/modality/dwi.py | 14 +++++++------- nireports/tests/test_dwi.py | 11 +++-------- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 6f7b19d7..1a0c63bb 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -410,7 +410,7 @@ def plot_gradients( def plot_carpet( nii, - gtab, + bvals, segmentation=None, sort_by_bval=False, output_file=None, @@ -424,8 +424,8 @@ def plot_carpet( ---------- nii : Nifti1Image DW imaging data - gtab : :obj:`GradientTable` - DW imaging data gradient data + bvals : :obj:`numpy.ndarray` + Rounded bvals segmentation : Nifti1Image Boolean or segmentation mask of DW imaging data sort_by_bval : :obj:`bool` @@ -448,16 +448,16 @@ def plot_carpet( nii_data = nii.get_fdata() - b0_data = nii_data[..., gtab.b0s_mask] - dw_data = nii_data[..., ~gtab.b0s_mask] + b0_data = nii_data[..., bvals == 0] + dw_data = nii_data[..., bvals > 0] 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]))) + np.argsort(bvals[bvals > 0] if sort_by_bval + else np.arange(len(bvals[bvals > 0]))) ) nii_data_div_b0 = nii_data_div_b0[..., sort_inds] diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 3623301f..4cbebe55 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -28,7 +28,6 @@ from matplotlib import pyplot as plt import os.path as op -import dipy.core.gradients as dpg from nireports.reportlets.modality.dwi \ import plot_dwi, plot_gradients, plot_carpet @@ -79,16 +78,12 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): testdata_name = "ds000114_sub-01_ses-test_desc-trunc_dwi" - nii_path = testdata_path / f'{testdata_name}.nii.gz' - bvec_path = testdata_path / f'{testdata_name}.bvec' - bval_path = testdata_path / f'{testdata_name}.bval' - - nii = nb.load(nii_path) - gtab = dpg.gradient_table(bval_path, bvec_path) + nii = nb.load(testdata_path / f'{testdata_name}.nii.gz') + bvals = np.loadtxt(testdata_path / f'{testdata_name}.bval') image_path = None if outdir is not None: image_path = outdir / f'{testdata_name}_carpet.svg' - plot_carpet(nii, gtab, output_file=image_path) + plot_carpet(nii, bvals, output_file=image_path) From 51b7f80808a9135d0981352e664d401615073797 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 14:55:10 -0400 Subject: [PATCH 04/16] Adds segmentation to test and makes bval optional --- nireports/reportlets/modality/dwi.py | 36 +++++++++++++--------------- nireports/tests/test_dwi.py | 11 ++++++++- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 1a0c63bb..bc21a6dd 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -410,7 +410,7 @@ def plot_gradients( def plot_carpet( nii, - bvals, + bvals=None, segmentation=None, sort_by_bval=False, output_file=None, @@ -448,28 +448,29 @@ def plot_carpet( nii_data = nii.get_fdata() - b0_data = nii_data[..., bvals == 0] - dw_data = nii_data[..., bvals > 0] + if bvals is not None: + b0_data = nii_data[..., bvals == 0] + dw_data = nii_data[..., bvals > 0] - bzero = np.mean(b0_data, -1) + bzero = np.mean(b0_data, -1) - nii_data_div_b0 = dw_data / bzero[..., np.newaxis] + nii_data = dw_data / bzero[..., np.newaxis] - sort_inds = ( - np.argsort(bvals[bvals > 0] if sort_by_bval - else np.arange(len(bvals[bvals > 0]))) - ) - nii_data_div_b0 = nii_data_div_b0[..., sort_inds] + sort_inds = ( + np.argsort(bvals[bvals > 0] if sort_by_bval + else np.arange(len(bvals[bvals > 0]))) + ) + nii_data = nii_data[..., sort_inds] # Reshape - nii_data_reshaped = nii_data_div_b0.reshape(-1, nii_data_div_b0.shape[-1]) + nii_data = nii_data.reshape(-1, nii_data.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, :] + nii_data = nii_data[segmentation_reshaped > 0, :] segmentation_masked = segmentation_reshaped[segmentation_reshaped > 0] if segment_labels is not None: @@ -483,19 +484,16 @@ def plot_carpet( ) segments[label] = indices - else: - nii_data_masked = nii_data_reshaped - - bad_row_ind = np.where(~np.isfinite(nii_data_masked))[0] + bad_row_ind = np.where(~np.isfinite(nii_data))[0] - good_row_ind = np.ones(nii_data_masked.shape[0], dtype=bool) + good_row_ind = np.ones(nii_data.shape[0], dtype=bool) good_row_ind[bad_row_ind] = False - nii_data_masked = nii_data_masked[good_row_ind, :] + nii_data = nii_data[good_row_ind, :] # Plot return nw_plot_carpet( - nii_data_masked, + nii_data, detrend=detrend, segments=segments, output_file=output_file diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 4cbebe55..17e11b0e 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -81,9 +81,18 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): nii = nb.load(testdata_path / f'{testdata_name}.nii.gz') bvals = np.loadtxt(testdata_path / f'{testdata_name}.bval') + nii_data = nii.get_fdata() + segmentation = nii_data > 3000 + segment_labels = {"<3000": [0], ">3000": [1]} + image_path = None if outdir is not None: image_path = outdir / f'{testdata_name}_carpet.svg' - plot_carpet(nii, bvals, output_file=image_path) + plot_carpet(nii, + bvals=bvals, + segmentation=segmentation, + segment_labels=segment_labels, + output_file=image_path + ) From 2a135cc1b6a5a0c8b87bca11d90e1a687412564f Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 15:10:28 -0400 Subject: [PATCH 05/16] Updates test segmentation --- nireports/tests/test_dwi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 17e11b0e..fb14844d 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -82,8 +82,8 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): bvals = np.loadtxt(testdata_path / f'{testdata_name}.bval') nii_data = nii.get_fdata() - segmentation = nii_data > 3000 - segment_labels = {"<3000": [0], ">3000": [1]} + segmentation = round(3*np.rand(nii_data.shape)) + segment_labels = {"0": [0], "1": [1], "2": [2], "3": [3]} image_path = None From f22f6c15204987db79786e3c0d7bb60cbe8bdae0 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 15:26:00 -0400 Subject: [PATCH 06/16] Random segmentation --- nireports/tests/test_dwi.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index fb14844d..71a215ba 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -82,7 +82,10 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): bvals = np.loadtxt(testdata_path / f'{testdata_name}.bval') nii_data = nii.get_fdata() - segmentation = round(3*np.rand(nii_data.shape)) + segmentation = np.round(3*np.random.rand(nii_data.shape[0], + nii_data.shape[1], + nii_data.shape[2], + nii_data.shape[3])) segment_labels = {"0": [0], "1": [1], "2": [2], "3": [3]} image_path = None From 2a56cf9b27de51c66f2fe026631cc63f4d008dc7 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 15:35:05 -0400 Subject: [PATCH 07/16] Changes segmentation to np array --- nireports/reportlets/modality/dwi.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index bc21a6dd..1a1f5fbe 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -426,8 +426,9 @@ def plot_carpet( DW imaging data bvals : :obj:`numpy.ndarray` Rounded bvals - segmentation : Nifti1Image + segmentation : :obj:`numpy.ndarray` Boolean or segmentation mask of DW imaging data + e.g. np.asanyarray(Nifti1Image.dataobj, dtype=np.int16) sort_by_bval : :obj:`bool` Flag to reorder time points by bvalue output_file : :obj:`string` @@ -466,10 +467,8 @@ def plot_carpet( nii_data = nii_data.reshape(-1, nii_data.shape[-1]) if segmentation is not None: - segmentation_data = np.asanyarray(segmentation.dataobj, dtype=np.int16) - # Apply mask - segmentation_reshaped = segmentation_data.reshape(-1) + segmentation_reshaped = segmentation.reshape(-1) nii_data = nii_data[segmentation_reshaped > 0, :] segmentation_masked = segmentation_reshaped[segmentation_reshaped > 0] From 2faac2c90bfa0ca4bd9e313b8b32c263f341171f Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 15 May 2024 15:55:14 -0400 Subject: [PATCH 08/16] Segmentation fixes --- nireports/reportlets/modality/dwi.py | 1 + nireports/tests/test_dwi.py | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 1a1f5fbe..a2014133 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -428,6 +428,7 @@ def plot_carpet( Rounded bvals segmentation : :obj:`numpy.ndarray` Boolean or segmentation mask of DW imaging data + Data corresponding to 0 values will be dropped e.g. np.asanyarray(Nifti1Image.dataobj, dtype=np.int16) sort_by_bval : :obj:`bool` Flag to reorder time points by bvalue diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 71a215ba..9d8f8bb6 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -84,8 +84,7 @@ def test_plot_carpet(tmp_path, testdata_path, outdir): nii_data = nii.get_fdata() segmentation = np.round(3*np.random.rand(nii_data.shape[0], nii_data.shape[1], - nii_data.shape[2], - nii_data.shape[3])) + nii_data.shape[2])) segment_labels = {"0": [0], "1": [1], "2": [2], "3": [3]} image_path = None From 41f1daf0d06a64c2727da441eb1900081974626d Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 16 May 2024 10:54:29 -0400 Subject: [PATCH 09/16] Adds get_segment_labels() --- nireports/reportlets/modality/dwi.py | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index a2014133..25cda5c7 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -498,3 +498,52 @@ def plot_carpet( 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 = {} + + 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 From 1465cd5a905c2e7017976cb50043795c752c0bf2 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 16 May 2024 15:51:09 -0400 Subject: [PATCH 10/16] Breaks higher level function into preprocessing functions --- nireports/reportlets/modality/dwi.py | 160 +++++++++++------- .../aseg.auto_noCCseg.label_intensities.txt | 47 +++++ nireports/tests/test_dwi.py | 41 +++-- 3 files changed, 166 insertions(+), 82 deletions(-) create mode 100644 nireports/tests/data/aseg.auto_noCCseg.label_intensities.txt diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 25cda5c7..8315e101 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -28,7 +28,6 @@ from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import art3d from nilearn.plotting import plot_anat -from nireports.reportlets.nuisance import plot_carpet as nw_plot_carpet def plot_dwi(dataobj, affine, gradient=None, **kwargs): @@ -408,96 +407,128 @@ def plot_gradients( return ax -def plot_carpet( +def nii_to_carpetplot_data( nii, bvals=None, - segmentation=None, + divide_by_b0=True, + drop_b0=True, sort_by_bval=False, - output_file=None, - segment_labels=None, - detrend=False, + mask_nii=None, + segment_labels=None ): """ - Return carpet plot using niworkflows carpet_plot + Convert nii to data matrix for carpet plot Parameters ---------- nii : Nifti1Image DW imaging data - bvals : :obj:`numpy.ndarray` + bvals : :obj:`numpy.ndarray`, optional Rounded bvals - segmentation : :obj:`numpy.ndarray` - Boolean or segmentation mask of DW imaging data - Data corresponding to 0 values will be dropped - e.g. np.asanyarray(Nifti1Image.dataobj, dtype=np.int16) - sort_by_bval : :obj:`bool` + 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 - output_file : :obj:`string` - Path to save the plot - segment_labels : :obj:`dict` + 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]} - detrend : :obj:`bool` - niworkflows plot_carpet detrend flag Returns --------- - matplotlib GridSpec object + 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. """ - segments = None nii_data = nii.get_fdata() - if bvals is not None: - b0_data = nii_data[..., bvals == 0] - dw_data = nii_data[..., bvals > 0] + if mask_nii is None: + mask_data = np.ones(nii_data.shape[0:2]) + else: + mask_data = np.asanyarray(mask_nii.dataobj, dtype=np.int16) - bzero = np.mean(b0_data, -1) + 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] - nii_data = dw_data / bzero[..., np.newaxis] + if drop_b0: + nii_data = nii_data[..., bvals > 0] + bvals = bvals[bvals > 0] - sort_inds = ( - np.argsort(bvals[bvals > 0] if sort_by_bval - else np.arange(len(bvals[bvals > 0]))) - ) - nii_data = nii_data[..., sort_inds] + if sort_by_bval: + sort_inds = np.argsort(bvals) + nii_data = nii_data[..., sort_inds] # Reshape - nii_data = nii_data.reshape(-1, nii_data.shape[-1]) - - if segmentation is not None: - # Apply mask - segmentation_reshaped = segmentation.reshape(-1) - nii_data = nii_data[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 - - bad_row_ind = np.where(~np.isfinite(nii_data))[0] - - good_row_ind = np.ones(nii_data.shape[0], dtype=bool) + 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 - nii_data = nii_data[good_row_ind, :] + data = data[good_row_ind, :] + mask_data = mask_data[good_row_ind] - # Plot - return nw_plot_carpet( - nii_data, - detrend=detrend, - segments=segments, - output_file=output_file - ) + # 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 = dict() + + 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)] + ) + segments[label] = indices + + return segments def get_segment_labels( @@ -508,7 +539,8 @@ def get_segment_labels( label_position=1 ): """ - Return segment labels for plot_carpet function + Get segment labels from file by keyword for get_segments function + Parameters ---------- filepath : :obj:`string` @@ -517,13 +549,13 @@ def get_segment_labels( 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 : :obj:`string`, optional Delimiter between label index and label string in label file (' ' for freesurfer label file) - index_position : :obj:`int` + index_position : :obj:`int`, optional Position of label index in label file (0 for freesurfer label file) - label_position : :obj:`int` + label_position : :obj:`int`, optional Position of label string in label file (1 for freesurfer label file) Returns diff --git a/nireports/tests/data/aseg.auto_noCCseg.label_intensities.txt b/nireports/tests/data/aseg.auto_noCCseg.label_intensities.txt new file mode 100644 index 00000000..5c453366 --- /dev/null +++ b/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 diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 9d8f8bb6..3eb8bb7b 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -27,9 +27,10 @@ import pytest from matplotlib import pyplot as plt -import os.path as op from nireports.reportlets.modality.dwi \ - import plot_dwi, plot_gradients, plot_carpet + import plot_dwi, plot_gradients, nii_to_carpetplot_data, get_segment_labels + +from nireports.reportlets.nuisance import plot_carpet def test_plot_dwi(tmp_path, testdata_path, outdir): @@ -73,28 +74,32 @@ def test_plot_gradients(tmp_path, testdata_path, dwi_btable, outdir): plt.savefig(outdir / f"{dwi_btable}.svg", bbox_inches="tight") -def test_plot_carpet(tmp_path, testdata_path, outdir): - """Check the carpet plot""" +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') - nii_data = nii.get_fdata() - segmentation = np.round(3*np.random.rand(nii_data.shape[0], - nii_data.shape[1], - nii_data.shape[2])) - segment_labels = {"0": [0], "1": [1], "2": [2], "3": [3]} - image_path = None if outdir is not None: - image_path = outdir / f'{testdata_name}_carpet.svg' - - plot_carpet(nii, - bvals=bvals, - segmentation=segmentation, - segment_labels=segment_labels, - output_file=image_path - ) + image_path = outdir / f'{testdata_name}_nii_to_carpet.svg' + + data, segments = nii_to_carpetplot_data(nii, bvals=bvals) + + 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 From 51198f29e7756466d82d21ee78cdf640b2476eb7 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 16 May 2024 15:59:08 -0400 Subject: [PATCH 11/16] Updates default when no mask --- nireports/reportlets/modality/dwi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 8315e101..63258db5 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -452,7 +452,7 @@ def nii_to_carpetplot_data( nii_data = nii.get_fdata() if mask_nii is None: - mask_data = np.ones(nii_data.shape[0:2]) + mask_data = np.ones(nii_data.shape[0:3]) else: mask_data = np.asanyarray(mask_nii.dataobj, dtype=np.int16) From 3adb4cd8ab7c172f8d9b49642e379523bca18ac8 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 17 May 2024 10:23:17 -0400 Subject: [PATCH 12/16] Expands test --- nireports/tests/test_dwi.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 3eb8bb7b..df6b2a24 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -82,12 +82,27 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): 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], + nii.shape[3])) + + 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) + 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) From 901bd4deb0548c3cc9a493602c0e69b50e188e7e Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 17 May 2024 10:32:10 -0400 Subject: [PATCH 13/16] Updates test --- nireports/tests/test_dwi.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index df6b2a24..19af7560 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -101,6 +101,7 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): data, segments = nii_to_carpetplot_data(nii, bvals=bvals, + drop_b0=False, mask_nii=mask_nii, segment_labels=segment_labels) From aed7b041d98ec821ad24289f2dd892f90b186320 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 17 May 2024 10:37:28 -0400 Subject: [PATCH 14/16] Updates test segmentation --- nireports/tests/test_dwi.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 19af7560..88fa7097 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -84,8 +84,7 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): mask_data = np.round(82 * np.random.rand(nii.shape[0], nii.shape[1], - nii.shape[2], - nii.shape[3])) + nii.shape[2])) mask_nii = nb.Nifti1Image(mask_data, np.eye(4)) @@ -101,7 +100,6 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): data, segments = nii_to_carpetplot_data(nii, bvals=bvals, - drop_b0=False, mask_nii=mask_nii, segment_labels=segment_labels) From c3eee8d4d2748fdaac595b41eb8d839d04d0c238 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 17 May 2024 10:44:45 -0400 Subject: [PATCH 15/16] Updates get_segments() --- nireports/reportlets/modality/dwi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 63258db5..84465e82 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -524,7 +524,7 @@ def get_segments( indices = np.array([], dtype=int) for ii in idx: indices = np.concatenate( - [indices, np.nonzero(segment_mask == ii)] + [indices, np.nonzero(segment_mask == ii)[0]] ) segments[label] = indices From 75cc84eb35489984bac7a3ec868110c5471a2897 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 1 Jul 2024 15:09:19 +0200 Subject: [PATCH 16/16] sty: revise format --- nireports/reportlets/modality/dwi.py | 27 +++++++-------------------- nireports/tests/test_dwi.py | 26 +++++++++++++------------- 2 files changed, 20 insertions(+), 33 deletions(-) diff --git a/nireports/reportlets/modality/dwi.py b/nireports/reportlets/modality/dwi.py index 84465e82..1d45702b 100644 --- a/nireports/reportlets/modality/dwi.py +++ b/nireports/reportlets/modality/dwi.py @@ -414,7 +414,7 @@ def nii_to_carpetplot_data( drop_b0=True, sort_by_bval=False, mask_nii=None, - segment_labels=None + segment_labels=None, ): """ Convert nii to data matrix for carpet plot @@ -495,10 +495,7 @@ def nii_to_carpetplot_data( return data, segments -def get_segments( - segment_mask, - segment_labels -): +def get_segments(segment_mask, segment_labels): """ Return segments dict for plot_carpet function @@ -518,26 +515,18 @@ def get_segments( A mapping between segment labels (e.g., `"Left Cortex"`) and list of indexes in the data array. """ - segments = dict() + 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]] - ) + 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 -): +def get_segment_labels(filepath, keywords, delimiter=" ", index_position=0, label_position=1): """ Get segment labels from file by keyword for get_segments function @@ -570,12 +559,10 @@ def get_segment_labels( with open(filepath, "r") as f: labels = f.read() - labels_s = [label.split(delimiter) for label in labels.split("\n") - if label != ""] + 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]] + ind = [int(i[index_position]) for i in labels_s if keyword in i[label_position]] segment_labels[keyword] = ind return segment_labels diff --git a/nireports/tests/test_dwi.py b/nireports/tests/test_dwi.py index 88fa7097..c90c81ce 100644 --- a/nireports/tests/test_dwi.py +++ b/nireports/tests/test_dwi.py @@ -27,9 +27,12 @@ import pytest from matplotlib import pyplot as plt -from nireports.reportlets.modality.dwi \ - import plot_dwi, plot_gradients, nii_to_carpetplot_data, get_segment_labels - +from nireports.reportlets.modality.dwi import ( + get_segment_labels, + nii_to_carpetplot_data, + plot_dwi, + plot_gradients, +) from nireports.reportlets.nuisance import plot_carpet @@ -79,12 +82,10 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): 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') + 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_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)) @@ -96,12 +97,11 @@ def test_nii_to_carpetplot_data(tmp_path, testdata_path, outdir): image_path = None if outdir is not None: - image_path = outdir / f'{testdata_name}_nii_to_carpet.svg' + 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) + 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)