From 1c32a01f43f671ea904e6a1dcd3b0ce76120bfc5 Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 11:59:41 -0400 Subject: [PATCH 01/33] Adds files for diffusion-related IQMs. --- mriqc/qc/diffusion.py | 31 +++++++++++++++++++++++++++++++ mriqc/qc/tests/test_diffusion.py | 25 +++++++++++++++++++++++++ 2 files changed, 56 insertions(+) create mode 100644 mriqc/qc/diffusion.py create mode 100644 mriqc/qc/tests/test_diffusion.py diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py new file mode 100644 index 000000000..c31af8c58 --- /dev/null +++ b/mriqc/qc/diffusion.py @@ -0,0 +1,31 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# 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/ +# + +""" +Image quality metrics for diffusion MRI data +============================================ +""" +# XXX It's ok to require a mask (brain, wm, cc_mask) + +def noise_func(): + pass \ No newline at end of file diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py new file mode 100644 index 000000000..f31583d42 --- /dev/null +++ b/mriqc/qc/tests/test_diffusion.py @@ -0,0 +1,25 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# 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/ +# + +def test_noise_function(): + pass From e10224668093c6c57783423d170787f2353962e7 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:10:44 -0400 Subject: [PATCH 02/33] Adds initial spike functions --- mriqc/qc/diffusion.py | 73 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 2 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index c31af8c58..4eb7a3800 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -21,11 +21,80 @@ # https://www.nipreps.org/community/licensing/ # +import numpy as np + + """ Image quality metrics for diffusion MRI data ============================================ """ # XXX It's ok to require a mask (brain, wm, cc_mask) -def noise_func(): - pass \ No newline at end of file +def get_spike_mask(data, z_threshold=3): + """ + Return binary mask of spike/no spike + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + numpy array + """ + threshold = (z_threshold*np.std(np.ravel(data))) + np.mean(np.ravel(data)) + spike_mask = data > threshold + + return spike_mask + + +def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): + """ + Return percentage of slices spiking along each dimension + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + slice_threshold : :obj:`float` + Percentage of slice elements that need to be above spike threshold for slice to be considered spiking + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + + ndim = data.ndim + slice_spike_percentage = np.zeros(ndim) + + for ii in range(ndim): + slice_spike_percentage[ii] = np.mean(np.mean(spike_mask, ii) > slice_threshold) + + return slice_spike_percentage + + +def get_global_spike_percentage(data, z_threshold=3): + """ + Return percentage of array elements spiking + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + global_spike_percentage = np.mean(np.ravel(spike_mask)) + + return global_spike_percentage \ No newline at end of file From 82927b04f851ba7826a31d9d93484a5192a33bc2 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:12:22 -0400 Subject: [PATCH 03/33] Corrects docustring --- mriqc/qc/diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index 4eb7a3800..cd98df082 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -66,7 +66,7 @@ def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): Returns --------- - float + array """ spike_mask = get_spike_mask(data, z_threshold) From 9fb37170077e2565ac881f0d56ae3a209972e4ad Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:46:28 -0400 Subject: [PATCH 04/33] Adds initial test functions --- mriqc/qc/tests/test_diffusion.py | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index f31583d42..127be5e5e 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,6 +20,34 @@ # # https://www.nipreps.org/community/licensing/ # +import numpy as np +from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage -def test_noise_function(): - pass + +slice = np.array([[1, 1, 10, 1], [1, 10, 1, 1], [1, 10, 1, 1], [1, 1, 1, 1]]) +test_data = np.array([slice, slice, slice, slice, slice]) + +z_threshold = 2 +slice_threshold = .2 + +def test_get_spike_mask(): + spike_mask = get_spike_mask(test_data, z_threshold) + + assert np.min(np.ravel(spike_mask)) == 0 + assert np.max(np.ravel(spike_mask)) == 1 + assert spike_mask.shape == test_data.shape + + +def test_get_slice_spike_percentage(): + slice_spike_percentage = get_slice_spike_percentage(test_data, z_threshold, slice_threshold) + + assert np.min(np.ravel(slice_spike_percentage)) >= 0 + assert np.max(np.ravel(slice_spike_percentage)) <= 1 + assert len(slice_spike_percentage) == test_data.ndim + + +def test_get_global_spike_percentage(): + global_spike_percentage = get_global_spike_percentage(test_data, z_threshold) + + assert np.min(np.ravel(global_spike_percentage)) >= 0 + assert np.max(np.ravel(global_spike_percentage)) <= 1 From 9a1d62b43ead483f78941929295c58d9d74496ba Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 15:04:19 -0400 Subject: [PATCH 05/33] Adds test fixture to fetch and read dMRI data. --- mriqc/qc/diffusion.py | 6 +++--- mriqc/qc/tests/test_diffusion.py | 29 +++++++++++++++++++++++++++-- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index c31af8c58..e10871e32 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -25,7 +25,7 @@ Image quality metrics for diffusion MRI data ============================================ """ -# XXX It's ok to require a mask (brain, wm, cc_mask) -def noise_func(): - pass \ No newline at end of file + +def noise_func(img, gtab): + pass diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index f31583d42..1e39bc5de 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,6 +20,31 @@ # # https://www.nipreps.org/community/licensing/ # +import pytest +import nibabel as nib +from dipy.core.gradients import gradient_table +from dipy.data.fetcher import fetch_sherbrooke_3shell +import os.path as op +from ..diffusion import noise_func -def test_noise_function(): - pass +class DiffusionData(object): + def get_data(self): + """ + Generate test data + """ + _, path = fetch_sherbrooke_3shell() + fnifti = op.join(path, 'HARDI193.nii.gz') + fnifti, bval, bvec = [op.join(path, f'HARDI193.{ext}') for + ext in ["nii.gz", "bval", "bvec"]] + img = nib.load(fnifti) + gtab = gradient_table(bval, bvec) + return img, gtab + + +@pytest.fixture +def ddata(): + return DiffusionData() + +def test_noise_function(ddata): + img, gtab = ddata.get_fdata() + noise_func(img, gtab) \ No newline at end of file From 2e916ccd649b016d0f9b102604475b1408fd97f2 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:10:44 -0400 Subject: [PATCH 06/33] Adds initial spike functions --- mriqc/qc/diffusion.py | 73 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index e10871e32..8c008da2c 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -21,11 +21,82 @@ # https://www.nipreps.org/community/licensing/ # +import numpy as np + + """ Image quality metrics for diffusion MRI data ============================================ """ - def noise_func(img, gtab): pass + +def get_spike_mask(data, z_threshold=3): + """ + Return binary mask of spike/no spike + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + numpy array + """ + threshold = (z_threshold*np.std(np.ravel(data))) + np.mean(np.ravel(data)) + spike_mask = data > threshold + + return spike_mask + + +def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): + """ + Return percentage of slices spiking along each dimension + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + slice_threshold : :obj:`float` + Percentage of slice elements that need to be above spike threshold for slice to be considered spiking + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + + ndim = data.ndim + slice_spike_percentage = np.zeros(ndim) + + for ii in range(ndim): + slice_spike_percentage[ii] = np.mean(np.mean(spike_mask, ii) > slice_threshold) + + return slice_spike_percentage + + +def get_global_spike_percentage(data, z_threshold=3): + """ + Return percentage of array elements spiking + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + global_spike_percentage = np.mean(np.ravel(spike_mask)) + + return global_spike_percentage From 085b0c0d016906b5cdbe23e6dae42e78a04c5762 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:12:22 -0400 Subject: [PATCH 07/33] Corrects docustring --- mriqc/qc/diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index 8c008da2c..929e2efdb 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -68,7 +68,7 @@ def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): Returns --------- - float + array """ spike_mask = get_spike_mask(data, z_threshold) From 1b98d20b1490f33bb57ce25782bce3077f61f5e3 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:46:28 -0400 Subject: [PATCH 08/33] Adds initial test functions --- mriqc/qc/tests/test_diffusion.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 1e39bc5de..6355a6294 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,6 +20,7 @@ # # https://www.nipreps.org/community/licensing/ # + import pytest import nibabel as nib from dipy.core.gradients import gradient_table @@ -27,6 +28,10 @@ import os.path as op from ..diffusion import noise_func +import numpy as np +from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage + + class DiffusionData(object): def get_data(self): """ @@ -47,4 +52,26 @@ def ddata(): def test_noise_function(ddata): img, gtab = ddata.get_fdata() - noise_func(img, gtab) \ No newline at end of file + noise_func(img, gtab) + +def test_get_spike_mask(): + spike_mask = get_spike_mask(ddata(), 2) + + assert np.min(np.ravel(spike_mask)) == 0 + assert np.max(np.ravel(spike_mask)) == 1 + assert spike_mask.shape == test_data.shape + + +def test_get_slice_spike_percentage(): + slice_spike_percentage = get_slice_spike_percentage(ddata(), 2, .2) + + assert np.min(np.ravel(slice_spike_percentage)) >= 0 + assert np.max(np.ravel(slice_spike_percentage)) <= 1 + assert len(slice_spike_percentage) == ddata().ndim + + +def test_get_global_spike_percentage(): + global_spike_percentage = get_global_spike_percentage(ddata(), 2) + + assert np.min(np.ravel(global_spike_percentage)) >= 0 + assert np.max(np.ravel(global_spike_percentage)) <= 1 From faf7cccb8557b9226de407b3428798de97b08a9d Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 15:15:19 -0400 Subject: [PATCH 09/33] Merge with start_on_diffusion --- mriqc/qc/tests/test_diffusion.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 958edaff4..c2c223c7b 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,18 +20,13 @@ # # https://www.nipreps.org/community/licensing/ # -import numpy as np -from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage - import pytest import nibabel as nib from dipy.core.gradients import gradient_table from dipy.data.fetcher import fetch_sherbrooke_3shell import os.path as op -from ..diffusion import noise_func - +from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage import numpy as np -from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage class DiffusionData(object): From afe177d809a1adc680f08f5c6b8b55d363b502d2 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 15:22:14 -0400 Subject: [PATCH 10/33] Merge with start_on_diffusion --- mriqc/qc/tests/test_diffusion.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index c2c223c7b..6a5950d0d 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -51,24 +51,28 @@ def test_noise_function(ddata): img, gtab = ddata.get_fdata() noise_func(img, gtab) -def test_get_spike_mask(): - spike_mask = get_spike_mask(ddata(), 2) + +def test_get_spike_mask(ddata): + img, gtab = ddata.get_fdata() + spike_mask = get_spike_mask(img, 2) assert np.min(np.ravel(spike_mask)) == 0 assert np.max(np.ravel(spike_mask)) == 1 - assert spike_mask.shape == test_data.shape + assert spike_mask.shape == img.shape -def test_get_slice_spike_percentage(): - slice_spike_percentage = get_slice_spike_percentage(ddata(), 2, .2) +def test_get_slice_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) assert np.min(np.ravel(slice_spike_percentage)) >= 0 assert np.max(np.ravel(slice_spike_percentage)) <= 1 - assert len(slice_spike_percentage) == ddata().ndim + assert len(slice_spike_percentage) == img.ndim -def test_get_global_spike_percentage(): - global_spike_percentage = get_global_spike_percentage(ddata(), 2) +def test_get_global_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + global_spike_percentage = get_global_spike_percentage(img, 2) assert np.min(np.ravel(global_spike_percentage)) >= 0 assert np.max(np.ravel(global_spike_percentage)) <= 1 From cf3c0e5a2bf64a20e71f7fefcb2f138796d792bf Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 11:40:49 -0400 Subject: [PATCH 11/33] Corrects for ndim --- mriqc/qc/tests/test_diffusion.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 6a5950d0d..27baa6135 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -65,8 +65,8 @@ def test_get_slice_spike_percentage(ddata): img, gtab = ddata.get_fdata() slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) - assert np.min(np.ravel(slice_spike_percentage)) >= 0 - assert np.max(np.ravel(slice_spike_percentage)) <= 1 + assert np.min(slice_spike_percentage) >= 0 + assert np.max(slice_spike_percentage) <= 1 assert len(slice_spike_percentage) == img.ndim @@ -74,5 +74,5 @@ def test_get_global_spike_percentage(ddata): img, gtab = ddata.get_fdata() global_spike_percentage = get_global_spike_percentage(img, 2) - assert np.min(np.ravel(global_spike_percentage)) >= 0 - assert np.max(np.ravel(global_spike_percentage)) <= 1 + assert global_spike_percentage >= 0 + assert global_spike_percentage <= 1 From e89bcd91e3239969343ae8c5aed8ac000dfd8f0c Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Fri, 21 Jul 2023 11:48:16 -0400 Subject: [PATCH 12/33] Add function to get shelled data. --- mriqc/qc/diffusion.py | 4 ++++ mriqc/qc/tests/test_diffusion.py | 27 ++++++++++++++++++++++++--- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index e10871e32..0ae7d2b27 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -29,3 +29,7 @@ def noise_func(img, gtab): pass + + +def noise_func_for_shelled_data(shelled_data, gtab): + pass \ No newline at end of file diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 1e39bc5de..28e6ef6d3 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -21,12 +21,15 @@ # https://www.nipreps.org/community/licensing/ # import pytest +import numpy as np import nibabel as nib from dipy.core.gradients import gradient_table from dipy.data.fetcher import fetch_sherbrooke_3shell +from dipy.core.gradients import unique_bvals_magnitude, round_bvals import os.path as op from ..diffusion import noise_func + class DiffusionData(object): def get_data(self): """ @@ -37,14 +40,32 @@ def get_data(self): fnifti, bval, bvec = [op.join(path, f'HARDI193.{ext}') for ext in ["nii.gz", "bval", "bvec"]] img = nib.load(fnifti) + data = img.get_fdata() gtab = gradient_table(bval, bvec) - return img, gtab + return data, gtab + + def shelled_data(self): + data, gtab = self.get_data() + rounded_bvals = round_bvals(gtab.bvals) + unique_rounded_bvals = np.unique(rounded_bvals) + + out_data = [] + for u_bv in unique_rounded_bvals: + this = data[..., np.where(rounded_bvals == u_bv)] + out_data.append(this) + return out_data, gtab @pytest.fixture def ddata(): return DiffusionData() + def test_noise_function(ddata): - img, gtab = ddata.get_fdata() - noise_func(img, gtab) \ No newline at end of file + data, gtab = ddata.get_data() + noise_func(data, gtab) + + +def test_with_shelled_data(ddata): + shelled_data, gtab = ddata.shelled_data() + noise_func_for_shelled_data(shelled_data, gtab) \ No newline at end of file From e17645bfead54a026f5abdf31245fcbde4539c50 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 11:57:45 -0400 Subject: [PATCH 13/33] Removes redundant numpy import --- mriqc/qc/tests/test_diffusion.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 88f7151e9..0caf471b5 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -28,8 +28,6 @@ from dipy.core.gradients import unique_bvals_magnitude, round_bvals import os.path as op from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage -import numpy as np - class DiffusionData(object): From 5a4109c202799806392dc21a292569f664332738 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 12:00:40 -0400 Subject: [PATCH 14/33] Removes redundant lines --- mriqc/qc/tests/test_diffusion.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 0caf471b5..2c2711f6c 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -90,8 +90,6 @@ def test_get_global_spike_percentage(ddata): assert global_spike_percentage >= 0 assert global_spike_percentage <= 1 - data, gtab = ddata.get_data() - noise_func(data, gtab) def test_with_shelled_data(ddata): From fe1f29b05be3528cc1ec68a2f211d1b0325fffd7 Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 11:59:41 -0400 Subject: [PATCH 15/33] Adds files for diffusion-related IQMs. --- mriqc/qc/diffusion.py | 31 +++++++++++++++++++++++++++++++ mriqc/qc/tests/test_diffusion.py | 25 +++++++++++++++++++++++++ 2 files changed, 56 insertions(+) create mode 100644 mriqc/qc/diffusion.py create mode 100644 mriqc/qc/tests/test_diffusion.py diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py new file mode 100644 index 000000000..c31af8c58 --- /dev/null +++ b/mriqc/qc/diffusion.py @@ -0,0 +1,31 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# 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/ +# + +""" +Image quality metrics for diffusion MRI data +============================================ +""" +# XXX It's ok to require a mask (brain, wm, cc_mask) + +def noise_func(): + pass \ No newline at end of file diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py new file mode 100644 index 000000000..f31583d42 --- /dev/null +++ b/mriqc/qc/tests/test_diffusion.py @@ -0,0 +1,25 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# 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/ +# + +def test_noise_function(): + pass From 4e2c4e689ff60d7f645964e43ed125dce13a4b4f Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 15:04:19 -0400 Subject: [PATCH 16/33] Adds test fixture to fetch and read dMRI data. --- mriqc/qc/diffusion.py | 6 +++--- mriqc/qc/tests/test_diffusion.py | 29 +++++++++++++++++++++++++++-- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index c31af8c58..e10871e32 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -25,7 +25,7 @@ Image quality metrics for diffusion MRI data ============================================ """ -# XXX It's ok to require a mask (brain, wm, cc_mask) -def noise_func(): - pass \ No newline at end of file + +def noise_func(img, gtab): + pass diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index f31583d42..1e39bc5de 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,6 +20,31 @@ # # https://www.nipreps.org/community/licensing/ # +import pytest +import nibabel as nib +from dipy.core.gradients import gradient_table +from dipy.data.fetcher import fetch_sherbrooke_3shell +import os.path as op +from ..diffusion import noise_func -def test_noise_function(): - pass +class DiffusionData(object): + def get_data(self): + """ + Generate test data + """ + _, path = fetch_sherbrooke_3shell() + fnifti = op.join(path, 'HARDI193.nii.gz') + fnifti, bval, bvec = [op.join(path, f'HARDI193.{ext}') for + ext in ["nii.gz", "bval", "bvec"]] + img = nib.load(fnifti) + gtab = gradient_table(bval, bvec) + return img, gtab + + +@pytest.fixture +def ddata(): + return DiffusionData() + +def test_noise_function(ddata): + img, gtab = ddata.get_fdata() + noise_func(img, gtab) \ No newline at end of file From fd6670b892980a18e4dc94dd6f0cb3fc428efad5 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:10:44 -0400 Subject: [PATCH 17/33] Adds initial spike functions --- mriqc/qc/diffusion.py | 73 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index e10871e32..8c008da2c 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -21,11 +21,82 @@ # https://www.nipreps.org/community/licensing/ # +import numpy as np + + """ Image quality metrics for diffusion MRI data ============================================ """ - def noise_func(img, gtab): pass + +def get_spike_mask(data, z_threshold=3): + """ + Return binary mask of spike/no spike + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + numpy array + """ + threshold = (z_threshold*np.std(np.ravel(data))) + np.mean(np.ravel(data)) + spike_mask = data > threshold + + return spike_mask + + +def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): + """ + Return percentage of slices spiking along each dimension + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + slice_threshold : :obj:`float` + Percentage of slice elements that need to be above spike threshold for slice to be considered spiking + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + + ndim = data.ndim + slice_spike_percentage = np.zeros(ndim) + + for ii in range(ndim): + slice_spike_percentage[ii] = np.mean(np.mean(spike_mask, ii) > slice_threshold) + + return slice_spike_percentage + + +def get_global_spike_percentage(data, z_threshold=3): + """ + Return percentage of array elements spiking + + Parameters + ---------- + data : numpy array + Data to be thresholded + z_threshold : :obj:`float` + Number of standard deviations above the mean to use as spike threshold + + Returns + --------- + float + """ + spike_mask = get_spike_mask(data, z_threshold) + global_spike_percentage = np.mean(np.ravel(spike_mask)) + + return global_spike_percentage From 401ad7d99be6dc7c864044d6e8be3fe0573eb4bf Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:12:22 -0400 Subject: [PATCH 18/33] Corrects docustring --- mriqc/qc/diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index 8c008da2c..929e2efdb 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -68,7 +68,7 @@ def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): Returns --------- - float + array """ spike_mask = get_spike_mask(data, z_threshold) From 0a0e6de6289f08d08f1bfe60dac63ed6f2bbf972 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:46:28 -0400 Subject: [PATCH 19/33] Adds initial test functions --- mriqc/qc/tests/test_diffusion.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 1e39bc5de..6355a6294 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,6 +20,7 @@ # # https://www.nipreps.org/community/licensing/ # + import pytest import nibabel as nib from dipy.core.gradients import gradient_table @@ -27,6 +28,10 @@ import os.path as op from ..diffusion import noise_func +import numpy as np +from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage + + class DiffusionData(object): def get_data(self): """ @@ -47,4 +52,26 @@ def ddata(): def test_noise_function(ddata): img, gtab = ddata.get_fdata() - noise_func(img, gtab) \ No newline at end of file + noise_func(img, gtab) + +def test_get_spike_mask(): + spike_mask = get_spike_mask(ddata(), 2) + + assert np.min(np.ravel(spike_mask)) == 0 + assert np.max(np.ravel(spike_mask)) == 1 + assert spike_mask.shape == test_data.shape + + +def test_get_slice_spike_percentage(): + slice_spike_percentage = get_slice_spike_percentage(ddata(), 2, .2) + + assert np.min(np.ravel(slice_spike_percentage)) >= 0 + assert np.max(np.ravel(slice_spike_percentage)) <= 1 + assert len(slice_spike_percentage) == ddata().ndim + + +def test_get_global_spike_percentage(): + global_spike_percentage = get_global_spike_percentage(ddata(), 2) + + assert np.min(np.ravel(global_spike_percentage)) >= 0 + assert np.max(np.ravel(global_spike_percentage)) <= 1 From 67788092d4170cdd29f75393b4791ce6c06d1c1a Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:10:44 -0400 Subject: [PATCH 20/33] Adds initial spike functions --- mriqc/qc/diffusion.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index 929e2efdb..d29c20851 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -29,9 +29,6 @@ ============================================ """ -def noise_func(img, gtab): - pass - def get_spike_mask(data, z_threshold=3): """ Return binary mask of spike/no spike @@ -68,7 +65,7 @@ def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): Returns --------- - array + float """ spike_mask = get_spike_mask(data, z_threshold) @@ -99,4 +96,4 @@ def get_global_spike_percentage(data, z_threshold=3): spike_mask = get_spike_mask(data, z_threshold) global_spike_percentage = np.mean(np.ravel(spike_mask)) - return global_spike_percentage + return global_spike_percentage \ No newline at end of file From 3fd165c514766e5c998eb3932819a8914bc401f8 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:12:22 -0400 Subject: [PATCH 21/33] Corrects docustring --- mriqc/qc/diffusion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index d29c20851..dbdbc9eeb 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -65,7 +65,7 @@ def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): Returns --------- - float + array """ spike_mask = get_spike_mask(data, z_threshold) From a317e60ec30bb0fad80657335db4436d979065c8 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 14:46:28 -0400 Subject: [PATCH 22/33] Adds initial test functions --- mriqc/qc/tests/test_diffusion.py | 40 +++++++------------------------- 1 file changed, 8 insertions(+), 32 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 6355a6294..127be5e5e 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,42 +20,18 @@ # # https://www.nipreps.org/community/licensing/ # - -import pytest -import nibabel as nib -from dipy.core.gradients import gradient_table -from dipy.data.fetcher import fetch_sherbrooke_3shell -import os.path as op -from ..diffusion import noise_func - import numpy as np from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage -class DiffusionData(object): - def get_data(self): - """ - Generate test data - """ - _, path = fetch_sherbrooke_3shell() - fnifti = op.join(path, 'HARDI193.nii.gz') - fnifti, bval, bvec = [op.join(path, f'HARDI193.{ext}') for - ext in ["nii.gz", "bval", "bvec"]] - img = nib.load(fnifti) - gtab = gradient_table(bval, bvec) - return img, gtab - - -@pytest.fixture -def ddata(): - return DiffusionData() +slice = np.array([[1, 1, 10, 1], [1, 10, 1, 1], [1, 10, 1, 1], [1, 1, 1, 1]]) +test_data = np.array([slice, slice, slice, slice, slice]) -def test_noise_function(ddata): - img, gtab = ddata.get_fdata() - noise_func(img, gtab) +z_threshold = 2 +slice_threshold = .2 def test_get_spike_mask(): - spike_mask = get_spike_mask(ddata(), 2) + spike_mask = get_spike_mask(test_data, z_threshold) assert np.min(np.ravel(spike_mask)) == 0 assert np.max(np.ravel(spike_mask)) == 1 @@ -63,15 +39,15 @@ def test_get_spike_mask(): def test_get_slice_spike_percentage(): - slice_spike_percentage = get_slice_spike_percentage(ddata(), 2, .2) + slice_spike_percentage = get_slice_spike_percentage(test_data, z_threshold, slice_threshold) assert np.min(np.ravel(slice_spike_percentage)) >= 0 assert np.max(np.ravel(slice_spike_percentage)) <= 1 - assert len(slice_spike_percentage) == ddata().ndim + assert len(slice_spike_percentage) == test_data.ndim def test_get_global_spike_percentage(): - global_spike_percentage = get_global_spike_percentage(ddata(), 2) + global_spike_percentage = get_global_spike_percentage(test_data, z_threshold) assert np.min(np.ravel(global_spike_percentage)) >= 0 assert np.max(np.ravel(global_spike_percentage)) <= 1 From 395daecae22534c540b66e7fa2241a18beae7971 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 15:15:19 -0400 Subject: [PATCH 23/33] Merge with start_on_diffusion --- mriqc/qc/tests/test_diffusion.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 127be5e5e..0e08479b6 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -20,12 +20,13 @@ # # https://www.nipreps.org/community/licensing/ # +import pytest +import nibabel as nib +from dipy.core.gradients import gradient_table +from dipy.data.fetcher import fetch_sherbrooke_3shell +import os.path as op +from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage import numpy as np -from mriqc.mriqc.qc import get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage - - -slice = np.array([[1, 1, 10, 1], [1, 10, 1, 1], [1, 10, 1, 1], [1, 1, 1, 1]]) -test_data = np.array([slice, slice, slice, slice, slice]) z_threshold = 2 slice_threshold = .2 From ae9f1209dc201be35bc35731b514d007724e1a12 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Thu, 20 Jul 2023 15:22:14 -0400 Subject: [PATCH 24/33] Merge with start_on_diffusion --- mriqc/qc/tests/test_diffusion.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 0e08479b6..ec5cc257c 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -31,24 +31,28 @@ z_threshold = 2 slice_threshold = .2 -def test_get_spike_mask(): - spike_mask = get_spike_mask(test_data, z_threshold) + +def test_get_spike_mask(ddata): + img, gtab = ddata.get_fdata() + spike_mask = get_spike_mask(img, 2) assert np.min(np.ravel(spike_mask)) == 0 assert np.max(np.ravel(spike_mask)) == 1 - assert spike_mask.shape == test_data.shape + assert spike_mask.shape == img.shape -def test_get_slice_spike_percentage(): - slice_spike_percentage = get_slice_spike_percentage(test_data, z_threshold, slice_threshold) +def test_get_slice_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) assert np.min(np.ravel(slice_spike_percentage)) >= 0 assert np.max(np.ravel(slice_spike_percentage)) <= 1 - assert len(slice_spike_percentage) == test_data.ndim + assert len(slice_spike_percentage) == img.ndim -def test_get_global_spike_percentage(): - global_spike_percentage = get_global_spike_percentage(test_data, z_threshold) +def test_get_global_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + global_spike_percentage = get_global_spike_percentage(img, 2) assert np.min(np.ravel(global_spike_percentage)) >= 0 assert np.max(np.ravel(global_spike_percentage)) <= 1 From 25a8a06359223d8751ea58f16a270a2af00f6860 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 11:40:49 -0400 Subject: [PATCH 25/33] Corrects for ndim --- mriqc/qc/tests/test_diffusion.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index ec5cc257c..8b31c17e3 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -45,8 +45,8 @@ def test_get_slice_spike_percentage(ddata): img, gtab = ddata.get_fdata() slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) - assert np.min(np.ravel(slice_spike_percentage)) >= 0 - assert np.max(np.ravel(slice_spike_percentage)) <= 1 + assert np.min(slice_spike_percentage) >= 0 + assert np.max(slice_spike_percentage) <= 1 assert len(slice_spike_percentage) == img.ndim @@ -54,5 +54,5 @@ def test_get_global_spike_percentage(ddata): img, gtab = ddata.get_fdata() global_spike_percentage = get_global_spike_percentage(img, 2) - assert np.min(np.ravel(global_spike_percentage)) >= 0 - assert np.max(np.ravel(global_spike_percentage)) <= 1 + assert global_spike_percentage >= 0 + assert global_spike_percentage <= 1 From dac938be9e5c6315ba47176351b935af69d703c7 Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Fri, 21 Jul 2023 11:48:16 -0400 Subject: [PATCH 26/33] Add function to get shelled data. --- mriqc/qc/diffusion.py | 66 ++------------------------------ mriqc/qc/tests/test_diffusion.py | 52 ++++++++++++++++--------- 2 files changed, 37 insertions(+), 81 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index dbdbc9eeb..e09c49b23 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -33,67 +33,9 @@ def get_spike_mask(data, z_threshold=3): """ Return binary mask of spike/no spike - Parameters - ---------- - data : numpy array - Data to be thresholded - z_threshold : :obj:`float` - Number of standard deviations above the mean to use as spike threshold +def noise_func(img, gtab): + pass - Returns - --------- - numpy array - """ - threshold = (z_threshold*np.std(np.ravel(data))) + np.mean(np.ravel(data)) - spike_mask = data > threshold - - return spike_mask - - -def get_slice_spike_percentage(data, z_threshold=3, slice_threshold=.05): - """ - Return percentage of slices spiking along each dimension - - Parameters - ---------- - data : numpy array - Data to be thresholded - z_threshold : :obj:`float` - Number of standard deviations above the mean to use as spike threshold - slice_threshold : :obj:`float` - Percentage of slice elements that need to be above spike threshold for slice to be considered spiking - - Returns - --------- - array - """ - spike_mask = get_spike_mask(data, z_threshold) - - ndim = data.ndim - slice_spike_percentage = np.zeros(ndim) - - for ii in range(ndim): - slice_spike_percentage[ii] = np.mean(np.mean(spike_mask, ii) > slice_threshold) - - return slice_spike_percentage - - -def get_global_spike_percentage(data, z_threshold=3): - """ - Return percentage of array elements spiking - - Parameters - ---------- - data : numpy array - Data to be thresholded - z_threshold : :obj:`float` - Number of standard deviations above the mean to use as spike threshold - - Returns - --------- - float - """ - spike_mask = get_spike_mask(data, z_threshold) - global_spike_percentage = np.mean(np.ravel(spike_mask)) - return global_spike_percentage \ No newline at end of file +def noise_func_for_shelled_data(shelled_data, gtab): + pass \ No newline at end of file diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 8b31c17e3..d3d664a8b 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -21,38 +21,52 @@ # https://www.nipreps.org/community/licensing/ # import pytest +import numpy as np import nibabel as nib from dipy.core.gradients import gradient_table from dipy.data.fetcher import fetch_sherbrooke_3shell +from dipy.core.gradients import unique_bvals_magnitude, round_bvals import os.path as op from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage import numpy as np -z_threshold = 2 -slice_threshold = .2 +class DiffusionData(object): + def get_data(self): + """ + Generate test data + """ + _, path = fetch_sherbrooke_3shell() + fnifti = op.join(path, 'HARDI193.nii.gz') + fnifti, bval, bvec = [op.join(path, f'HARDI193.{ext}') for + ext in ["nii.gz", "bval", "bvec"]] + img = nib.load(fnifti) + data = img.get_fdata() + gtab = gradient_table(bval, bvec) + return data, gtab -def test_get_spike_mask(ddata): - img, gtab = ddata.get_fdata() - spike_mask = get_spike_mask(img, 2) + def shelled_data(self): + data, gtab = self.get_data() + rounded_bvals = round_bvals(gtab.bvals) + unique_rounded_bvals = np.unique(rounded_bvals) - assert np.min(np.ravel(spike_mask)) == 0 - assert np.max(np.ravel(spike_mask)) == 1 - assert spike_mask.shape == img.shape + out_data = [] + for u_bv in unique_rounded_bvals: + this = data[..., np.where(rounded_bvals == u_bv)] + out_data.append(this) + return out_data, gtab -def test_get_slice_spike_percentage(ddata): - img, gtab = ddata.get_fdata() - slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) +@pytest.fixture +def ddata(): + return DiffusionData() - assert np.min(slice_spike_percentage) >= 0 - assert np.max(slice_spike_percentage) <= 1 - assert len(slice_spike_percentage) == img.ndim +def test_noise_function(ddata): + data, gtab = ddata.get_data() + noise_func(data, gtab) -def test_get_global_spike_percentage(ddata): - img, gtab = ddata.get_fdata() - global_spike_percentage = get_global_spike_percentage(img, 2) - assert global_spike_percentage >= 0 - assert global_spike_percentage <= 1 +def test_with_shelled_data(ddata): + shelled_data, gtab = ddata.shelled_data() + noise_func_for_shelled_data(shelled_data, gtab) \ No newline at end of file From 278fc990c08e5248858bbfe3ede15c028b809696 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 11:57:45 -0400 Subject: [PATCH 27/33] Removes redundant numpy import --- mriqc/qc/tests/test_diffusion.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index d3d664a8b..145d09e68 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -28,8 +28,6 @@ from dipy.core.gradients import unique_bvals_magnitude, round_bvals import os.path as op from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage -import numpy as np - class DiffusionData(object): def get_data(self): From 7a5505847c1843d7fbe97069f9c84d106f54421e Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Fri, 21 Jul 2023 12:00:40 -0400 Subject: [PATCH 28/33] Removes redundant lines --- mriqc/qc/tests/test_diffusion.py | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 145d09e68..ed68b05b9 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -61,8 +61,34 @@ def ddata(): def test_noise_function(ddata): - data, gtab = ddata.get_data() - noise_func(data, gtab) + img, gtab = ddata.get_fdata() + noise_func(img, gtab) + + +def test_get_spike_mask(ddata): + img, gtab = ddata.get_fdata() + spike_mask = get_spike_mask(img, 2) + + assert np.min(np.ravel(spike_mask)) == 0 + assert np.max(np.ravel(spike_mask)) == 1 + assert spike_mask.shape == img.shape + + +def test_get_slice_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + slice_spike_percentage = get_slice_spike_percentage(img, 2, .2) + + assert np.min(slice_spike_percentage) >= 0 + assert np.max(slice_spike_percentage) <= 1 + assert len(slice_spike_percentage) == img.ndim + + +def test_get_global_spike_percentage(ddata): + img, gtab = ddata.get_fdata() + global_spike_percentage = get_global_spike_percentage(img, 2) + + assert global_spike_percentage >= 0 + assert global_spike_percentage <= 1 def test_with_shelled_data(ddata): From b05e226efdd42b794538561754b07019890e10b5 Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 14:14:39 -0400 Subject: [PATCH 29/33] Adds WIP implementation of noise and SNR estimation. --- mriqc/qc/diffusion.py | 84 +++++++++++++++++++++++++++++++- mriqc/qc/tests/test_diffusion.py | 21 +++++--- 2 files changed, 95 insertions(+), 10 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index e09c49b23..5d6f3cde6 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -29,6 +29,12 @@ ============================================ """ +import numpy as np +from dipy.core.gradients import gradient_table +from dipy.core.gradients import GradientTable +from dipy.reconst.dti import TensorModel +from dipy.denoise.noise_estimate import piesno + def get_spike_mask(data, z_threshold=3): """ Return binary mask of spike/no spike @@ -37,5 +43,79 @@ def noise_func(img, gtab): pass -def noise_func_for_shelled_data(shelled_data, gtab): - pass \ No newline at end of file + +def noise_b0(data, gtab, mask=None): + """ + Estimate noise in raw dMRI based on b0 variance. + + Parameters + ---------- + """ + if mask is None: + mask = np.ones(data.shape[:3], dtype=bool) + b0 = data[..., ~gtab.b0s_mask] + return np.percentile(np.var(b0[mask], -1), (25, 50, 75)) + + +def noise_piesno(data, n_channels=4): + """ + Estimate noise in raw dMRI data using the PIESNO [1]_ algorithm. + + + Parameters + ---------- + + Returns + ------- + + + Notes + ----- + + .. [1] Koay C.G., E. Ozarslan, C. Pierpaoli. Probabilistic Identification + and Estimation of Noise (PIESNO): A self-consistent approach and + its applications in MRI. JMR, 199(1):94-103, 2009. + """ + sigma, mask = piesno(data, N=n_channels, return_mask=True) + return sigma, mask + + +def cc_snr(data, gtab): + """ + Calculate worse-/best-case signal-to-noise ratio in the corpus callosum + + Parameters + ---------- + data : ndarray + + gtab : GradientTable class instance or tuple + + """ + if isinstance(gtab, GradientTable): + pass + + # XXX Per-shell calculation + tenmodel = TensorModel(gtab) + tensorfit = tenmodel.fit(data, mask=mask) + + from dipy.segment.mask import segment_from_cfa + from dipy.segment.mask import bounding_box + + threshold = (0.6, 1, 0, 0.1, 0, 0.1) + CC_box = np.zeros_like(data[..., 0]) + + mins, maxs = bounding_box(mask) + mins = np.array(mins) + maxs = np.array(maxs) + diff = (maxs - mins) // 4 + bounds_min = mins + diff + bounds_max = maxs - diff + + CC_box[bounds_min[0]:bounds_max[0], + bounds_min[1]:bounds_max[1], + bounds_min[2]:bounds_max[2]] = 1 + + mask_cc_part, cfa = segment_from_cfa(tensorfit, CC_box, threshold, + return_cfa=True) + + mean_signal = np.mean(data[mask_cc_part], axis=0) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index ed68b05b9..0008f79b0 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -21,6 +21,7 @@ # https://www.nipreps.org/community/licensing/ # import pytest +import os.path as op import numpy as np import nibabel as nib from dipy.core.gradients import gradient_table @@ -28,6 +29,7 @@ from dipy.core.gradients import unique_bvals_magnitude, round_bvals import os.path as op from ..diffusion import noise_func, get_spike_mask, get_slice_spike_percentage, get_global_spike_percentage +from ..diffusion import noise_b0, noise_piesno class DiffusionData(object): def get_data(self): @@ -54,17 +56,10 @@ def shelled_data(self): out_data.append(this) return out_data, gtab - @pytest.fixture def ddata(): return DiffusionData() - -def test_noise_function(ddata): - img, gtab = ddata.get_fdata() - noise_func(img, gtab) - - def test_get_spike_mask(ddata): img, gtab = ddata.get_fdata() spike_mask = get_spike_mask(img, 2) @@ -93,4 +88,14 @@ def test_get_global_spike_percentage(ddata): def test_with_shelled_data(ddata): shelled_data, gtab = ddata.shelled_data() - noise_func_for_shelled_data(shelled_data, gtab) \ No newline at end of file + noise_func_for_shelled_data(shelled_data, gtab) + + +def test_noise_b0(ddata): + data, gtab = ddata.get_data() + noise_b0(data, gtab) + + +def test_noise_piesno(ddata): + data, gtab = ddata.get_data() + noise_piesno(data) From a243f505435c3fb44e74fb89cbe6fa16292f0e0b Mon Sep 17 00:00:00 2001 From: Ariel Rokem Date: Thu, 20 Jul 2023 15:12:46 -0400 Subject: [PATCH 30/33] Adds whitespace. --- mriqc/qc/tests/test_diffusion.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mriqc/qc/tests/test_diffusion.py b/mriqc/qc/tests/test_diffusion.py index 0008f79b0..35af7f023 100644 --- a/mriqc/qc/tests/test_diffusion.py +++ b/mriqc/qc/tests/test_diffusion.py @@ -56,6 +56,7 @@ def shelled_data(self): out_data.append(this) return out_data, gtab + @pytest.fixture def ddata(): return DiffusionData() From 67f2aed12ae2fe5ae915671feba06c179bb7e1f0 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Mon, 5 Feb 2024 09:46:37 -0800 Subject: [PATCH 31/33] In progress, will update --- mriqc/qc/diffusion.py | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index d89ef0454..dd10c1aeb 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -32,7 +32,7 @@ def noise_func(img, gtab): pass -def get_spike_mask(data, z_threshold=3): +def get_spike_mask(data, grouping_vals, z_threshold=3): """ Return binary mask of spike/no spike @@ -40,6 +40,8 @@ def get_spike_mask(data, z_threshold=3): ---------- data : numpy array Data to be thresholded + grouping_vals : numpy array + Values by which to group data for thresholding (bvals or full mask) z_threshold : :obj:`float` Number of standard deviations above the mean to use as spike threshold @@ -47,8 +49,23 @@ def get_spike_mask(data, z_threshold=3): --------- numpy array """ - threshold = (z_threshold*np.std(np.ravel(data))) + np.mean(np.ravel(data)) - spike_mask = data > threshold + + threshold_mask = np.zeros(data.shape) + gvals = np.unique(grouping_vals) # Not the right function + gval_inds = np.where(grouping_vals == gval) # Not the right functions + + if grouping_vals.shape == data.shape: + for gval in gvals: + gval_data = data[gval_inds] + gval_threshold = (z_threshold*np.std(gval_data)) + np.mean(gval_data) + threshold_mask[gval_inds] = gval_threshold*np.ones(gval_data.shape) + else: + for gval in gvals: + gval_data = data[...,gval_inds] + gval_threshold = (z_threshold*np.std(gval_data)) + np.mean(gval_data) + threshold_mask[...,gval_inds] = gval_threshold*np.ones(gval_data.shape) + + spike_mask = data > threshold_mask return spike_mask From a88cfa128fca0cf4fca9636e0770e1cdce274b91 Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Mon, 5 Feb 2024 11:56:22 -0800 Subject: [PATCH 32/33] Adds start_on_diffusion functions --- mriqc/qc/diffusion.py | 119 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 107 insertions(+), 12 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index dd10c1aeb..c823148d7 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -1,3 +1,4 @@ + # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: # @@ -21,18 +22,100 @@ # https://www.nipreps.org/community/licensing/ # -import numpy as np - - """ Image quality metrics for diffusion MRI data ============================================ """ +import numpy as np +from dipy.core.gradients import gradient_table +from dipy.core.gradients import GradientTable +from dipy.reconst.dti import TensorModel +from dipy.denoise.noise_estimate import piesno +from dipy.core.gradients import unique_bvals_magnitude +from dipy.core.gradients import round_bvals + def noise_func(img, gtab): pass -def get_spike_mask(data, grouping_vals, z_threshold=3): +def noise_b0(data, gtab, mask=None): + """ + Estimate noise in raw dMRI based on b0 variance. + + Parameters + ---------- + """ + if mask is None: + mask = np.ones(data.shape[:3], dtype=bool) + b0 = data[..., ~gtab.b0s_mask] + return np.percentile(np.var(b0[mask], -1), (25, 50, 75)) + + +def noise_piesno(data, n_channels=4): + """ + Estimate noise in raw dMRI data using the PIESNO [1]_ algorithm. + + + Parameters + ---------- + + Returns + ------- + + + Notes + ----- + + .. [1] Koay C.G., E. Ozarslan, C. Pierpaoli. Probabilistic Identification + and Estimation of Noise (PIESNO): A self-consistent approach and + its applications in MRI. JMR, 199(1):94-103, 2009. + """ + sigma, mask = piesno(data, N=n_channels, return_mask=True) + return sigma, mask + + +def cc_snr(data, gtab): + """ + Calculate worse-/best-case signal-to-noise ratio in the corpus callosum + + Parameters + ---------- + data : ndarray + + gtab : GradientTable class instance or tuple + + """ + if isinstance(gtab, GradientTable): + pass + + # XXX Per-shell calculation + tenmodel = TensorModel(gtab) + tensorfit = tenmodel.fit(data, mask=mask) + + from dipy.segment.mask import segment_from_cfa + from dipy.segment.mask import bounding_box + + threshold = (0.6, 1, 0, 0.1, 0, 0.1) + CC_box = np.zeros_like(data[..., 0]) + + mins, maxs = bounding_box(mask) + mins = np.array(mins) + maxs = np.array(maxs) + diff = (maxs - mins) // 4 + bounds_min = mins + diff + bounds_max = maxs - diff + + CC_box[bounds_min[0]:bounds_max[0], + bounds_min[1]:bounds_max[1], + bounds_min[2]:bounds_max[2]] = 1 + + mask_cc_part, cfa = segment_from_cfa(tensorfit, CC_box, threshold, + return_cfa=True) + + mean_signal = np.mean(data[mask_cc_part], axis=0) + + +def get_spike_mask(data, z_threshold=3, grouping_vals=None, bmag=None): """ Return binary mask of spike/no spike @@ -40,30 +123,42 @@ def get_spike_mask(data, grouping_vals, z_threshold=3): ---------- data : numpy array Data to be thresholded - grouping_vals : numpy array - Values by which to group data for thresholding (bvals or full mask) z_threshold : :obj:`float` Number of standard deviations above the mean to use as spike threshold + grouping_vals : numpy array + Values by which to group data for thresholding (bvals or full mask) + bmag : int + From dipy.core.gradients: + The order of magnitude that the bvalues have to differ to be + considered an unique b-value. B-values are also rounded up to + this order of magnitude. Default: derive this value from the + maximal b-value provided: $bmag=log_{10}(max(bvals)) - 1$. Returns --------- numpy array """ + if grouping_vals is None: + threshold = (z_threshold*np.std(data)) + np.mean(data) + spike_mask = data > threshold + return spike_mask + threshold_mask = np.zeros(data.shape) - gvals = np.unique(grouping_vals) # Not the right function - gval_inds = np.where(grouping_vals == gval) # Not the right functions + + rounded_grouping_vals = round_bvals(grouping_vals, bmag) + gvals = unique_bvals_magnitude(grouping_vals, bmag) if grouping_vals.shape == data.shape: for gval in gvals: - gval_data = data[gval_inds] + gval_data = data[rounded_grouping_vals == gval] gval_threshold = (z_threshold*np.std(gval_data)) + np.mean(gval_data) - threshold_mask[gval_inds] = gval_threshold*np.ones(gval_data.shape) + threshold_mask[rounded_grouping_vals == gval] = gval_threshold*np.ones(gval_data.shape) else: for gval in gvals: - gval_data = data[...,gval_inds] + gval_data = data[..., rounded_grouping_vals == gval] gval_threshold = (z_threshold*np.std(gval_data)) + np.mean(gval_data) - threshold_mask[...,gval_inds] = gval_threshold*np.ones(gval_data.shape) + threshold_mask[..., rounded_grouping_vals == gval] = gval_threshold*np.ones(gval_data.shape) spike_mask = data > threshold_mask From 64f2513d40ff4265282a2f94703923969f30b54e Mon Sep 17 00:00:00 2001 From: Teresa Gomez <46339554+teresamg@users.noreply.github.com> Date: Wed, 7 Feb 2024 15:01:16 -0800 Subject: [PATCH 33/33] Updates cc_snr function with per-shell calculations --- mriqc/qc/diffusion.py | 61 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 54 insertions(+), 7 deletions(-) diff --git a/mriqc/qc/diffusion.py b/mriqc/qc/diffusion.py index c823148d7..c294659f9 100644 --- a/mriqc/qc/diffusion.py +++ b/mriqc/qc/diffusion.py @@ -34,6 +34,8 @@ from dipy.denoise.noise_estimate import piesno from dipy.core.gradients import unique_bvals_magnitude from dipy.core.gradients import round_bvals +from dipy.segment.mask import segment_from_cfa +from dipy.segment.mask import bounding_box def noise_func(img, gtab): pass @@ -74,7 +76,7 @@ def noise_piesno(data, n_channels=4): return sigma, mask -def cc_snr(data, gtab): +def cc_snr(data, gtab, bmag=None, mask=None): """ Calculate worse-/best-case signal-to-noise ratio in the corpus callosum @@ -84,21 +86,31 @@ def cc_snr(data, gtab): gtab : GradientTable class instance or tuple + bmag : int + From dipy.core.gradients: + The order of magnitude that the bvalues have to differ to be + considered an unique b-value. B-values are also rounded up to + this order of magnitude. Default: derive this value from the + maximal b-value provided: $bmag=log_{10}(max(bvals)) - 1$. + + mask : numpy array + Boolean brain mask + + """ if isinstance(gtab, GradientTable): pass - # XXX Per-shell calculation + if mask is None: + mask = np.ones(data.shape[:3]) + tenmodel = TensorModel(gtab) tensorfit = tenmodel.fit(data, mask=mask) - from dipy.segment.mask import segment_from_cfa - from dipy.segment.mask import bounding_box - threshold = (0.6, 1, 0, 0.1, 0, 0.1) CC_box = np.zeros_like(data[..., 0]) - mins, maxs = bounding_box(mask) + mins, maxs = bounding_box(mask) #mask needs to be volume mins = np.array(mins) maxs = np.array(maxs) diff = (maxs - mins) // 4 @@ -112,7 +124,42 @@ def cc_snr(data, gtab): mask_cc_part, cfa = segment_from_cfa(tensorfit, CC_box, threshold, return_cfa=True) - mean_signal = np.mean(data[mask_cc_part], axis=0) + b0_data = data[..., gtab.b0s_mask] + std_signal = np.std(b0_data[mask_cc_part], axis=-1) + + # Per-shell calculation + rounded_bvals = round_bvals(gtab.bvals, bmag) + bvals = unique_bvals_magnitude(gtab.bvals, bmag) + + cc_snr_best = np.zeros(gtab.bvals.shape) + cc_snr_worst = np.zeros(gtab.bvals.shape) + + for ind, bval in enumerate(bvals): + if bval == 0: + mean_signal = np.mean(data[..., rounded_bvals == 0], axis=-1) + cc_snr_worst[ind] = np.mean(mean_signal/std_signal) + cc_snr_best[ind] = np.mean(mean_signal/std_signal) + continue + + bval_data = data[..., rounded_bvals == bval] + bval_bvecs = gtab.bvecs[rounded_bvals == bval] + + axis_X = np.argmin(np.sum((bval_bvecs-np.array([1, 0, 0]))**2, axis=-1)) + axis_Y = np.argmin(np.sum((bval_bvecs-np.array([0, 1, 0]))**2, axis=-1)) + axis_Z = np.argmin(np.sum((bval_bvecs-np.array([0, 0, 1]))**2, axis=-1)) + + data_X = bval_data[..., axis_X] + data_Y = bval_data[..., axis_Y] + data_Z = bval_data[..., axis_Z] + + mean_signal_X = np.mean(data_X[mask_cc_part]) + mean_signal_Y = np.mean(data_Y[mask_cc_part]) + mean_signal_Z = np.mean(data_Z[mask_cc_part]) + + cc_snr_worst[ind] = np.mean(mean_signal_X/std_signal) + cc_snr_best[ind] = np.mean(np.mean(mean_signal_Y, mean_signal_Z)/std_signal) + + return cc_snr_worst, cc_snr_best def get_spike_mask(data, z_threshold=3, grouping_vals=None, bmag=None):