Adds WIP implementation of noise and SNR estimation.

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
+
+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)

mriqc/qc/tests/test_diffusion.py

@@ -21,13 +21,15 @@
 #     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
 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
+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)
+    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)