Merge pull request #1226 from nipreps/enh/ndc-iqm

ENH: Add new IQM for DWI -- NDC

mriqc/interfaces/diffusion.py

@@ -123,6 +123,12 @@ class _DiffusionQCInputSpec(_BaseInterfaceInputSpec):
         desc='FD threshold for orientation exclusion based on head motion'
     )
     in_fwhm = traits.List(traits.Float, desc='smoothness estimated with AFNI')
+    qspace_neighbors = traits.List(
+        traits.Tuple(traits.Int, traits.Int),
+        mandatory=True,
+        minlen=1,
+        desc='q-space nearest neighbor pairs',
+    )
 
 
 class _DiffusionQCOutputSpec(TraitedSpec):
@@ -132,8 +138,8 @@ class _DiffusionQCOutputSpec(TraitedSpec):
     fa_nans = traits.Float
     fber = traits.Dict
     fd = traits.Dict
+    ndc = traits.Float
     spikes_ppm = traits.Dict
-    # snr = traits.Float
     # gsr = traits.Dict
     # tsnr = traits.Float
     # fwhm = traits.Dict(desc='full width half-maximum measure')
@@ -172,7 +178,7 @@ def _run_interface(self, runtime):
         # Get brain mask data
         msknii = nb.load(self.inputs.brain_mask)
         mskdata = np.round(  # Protect the thresholding with a rounding for stability
-            np.asanyarray(msknii.dataobj),
+            msknii.get_fdata(),
             3,
         )
         if np.sum(mskdata) < 100:
@@ -256,6 +262,19 @@ def _run_interface(self, runtime):
             'perc': float(num_fd * 100 / (len(fd_data) + 1)),
         }
 
+        # NDC
+        dwidata = np.round(
+            np.nan_to_num(nb.load(self.inputs.in_file).get_fdata()),
+            3,
+        )
+        self._results['ndc'] = float(
+            dqc.neighboring_dwi_correlation(
+                dwidata,
+                neighbor_indices=self.inputs.qspace_neighbors,
+                mask=mskdata > 0.5,
+            )
+        )
+
         self._results['out_qc'] = _flatten_dict(self._results)
 
         return runtime
@@ -265,6 +284,10 @@ class _ReadDWIMetadataOutputSpec(_ReadSidecarJSONOutputSpec):
     out_bvec_file = File(desc='corresponding bvec file')
     out_bval_file = File(desc='corresponding bval file')
     out_bmatrix = traits.List(traits.List(traits.Float), desc='b-matrix')
+    qspace_neighbors = traits.List(
+        traits.Tuple(traits.Int, traits.Int),
+        desc='q-space nearest neighbor pairs',
+    )
 
 
 class ReadDWIMetadata(ReadSidecarJSON):
@@ -283,6 +306,7 @@ def _run_interface(self, runtime):
         bvecs = np.loadtxt(self._results['out_bvec_file']).T
         bvals = np.loadtxt(self._results['out_bval_file'])
 
+        self._results['qspace_neighbors'] = _find_qspace_neighbors(bvals, bvecs)
         self._results['out_bmatrix'] = np.hstack((bvecs, bvals[:, np.newaxis])).tolist()
 
         return runtime
@@ -1152,3 +1176,79 @@ def get_spike_mask(
         spike_mask[..., b_mask] = shelldata > a_thres
 
     return spike_mask
+
+
+def _find_qspace_neighbors(
+    bvals: np.ndarray, bvecs: np.ndarray
+) -> list[tuple[int, int]]:
+    """
+    Create a mapping of dwi volume index to its nearest neighbor in q-space.
+
+    This function implements an approximate nearest neighbor search in q-space
+    (excluding delta encoding). It calculates the Cartesian distance between
+    q-space representations of each diffusion-weighted imaging (DWI) volume
+    (represented by b-values and b-vectors) and identifies the closest one
+    (excluding the volume itself and b=0 volumes).
+
+    Parameters
+    ----------
+    bvals : :obj:`~numpy.ndarray`
+        List of b-values.
+    bvecs : :obj:`~numpy.ndarray`
+        Table of b-vectors.
+
+    Returns
+    -------
+    :obj:`list` of :obj:`tuple`
+        A list of 2-tuples indicating the nearest q-space neighbor
+        of each dwi volume.
+
+    Examples
+    --------
+    >>> _find_qspace_neighbors(
+    ...     np.array([0, 1000, 1000, 2000]),
+    ...     np.array([
+    ...         [1, 0, 0],
+    ...         [1, 0, 0],
+    ...         [0.99, 0.0001, 0.0001],
+    ...         [1, 0, 0]
+    ...     ]),
+    ... )
+    [(1, 2), (2, 1), (3, 1)]
+
+    Notes
+    -----
+    This is a copy of DIPY's code to be removed (and just imported) as soon as
+    a new release of DIPY is cut including
+    `dipy/dipy#3156 <https://github.com/dipy/dipy/pull/3156>`__.
+
+    """
+    from dipy.core.geometry import cart_distance
+    from dipy.core.gradients import gradient_table
+
+    gtab = gradient_table(bvals, bvecs)
+
+    dwi_neighbors: list[tuple[int, int]] = []
+
+    # Only correlate the b>0 images
+    dwi_indices = np.flatnonzero(~gtab.b0s_mask)
+
+    # Get a pseudo-qspace value for b>0s
+    qvecs = np.sqrt(gtab.bvals)[:, np.newaxis] * gtab.bvecs
+
+    for dwi_index in dwi_indices:
+        qvec = qvecs[dwi_index]
+
+        # Calculate distance in q-space, accounting for symmetry
+        pos_dist = cart_distance(qvec[np.newaxis, :], qvecs)
+        neg_dist = cart_distance(qvec[np.newaxis, :], -qvecs)
+        distances = np.min(np.column_stack([pos_dist, neg_dist]), axis=1)
+
+        # Be sure we don't select the image as its own neighbor
+        distances[dwi_index] = np.inf
+        # Or a b=0
+        distances[gtab.b0s_mask] = np.inf
+        neighbor_index = np.argmin(distances)
+        dwi_neighbors.append((dwi_index, neighbor_index))
+
+    return dwi_neighbors

mriqc/qc/diffusion.py

@@ -316,3 +316,60 @@ def spike_ppm(
     ]
 
     return {'global': spike_perc_global, 'slice': spike_perc_slice}
+
+
+def neighboring_dwi_correlation(
+    dwi_data: np.ndarray,
+    neighbor_indices: list[tuple[int, int]],
+    mask: np.ndarray | None = None,
+) -> float:
+    """
+    Calculates the Neighboring DWI Correlation (NDC) from diffusion MRI (dMRI) data.
+
+    The NDC is a measure of the correlation between signal intensities in neighboring
+    diffusion-weighted images (DWIs) within a mask. A low NDC (typically below 0.4)
+    can indicate poor image quality, according to Yeh et al. [Yeh2019]_.
+
+    Parameters
+    ----------
+    dwi_data : 4D :obj:`~numpy.ndarray`
+        DWI data on which to calculate NDC
+    neighbor_indices : :obj:`list` of :obj:`tuple`
+        List of (from, to) index neighbors.
+    mask : 3D :obj:`~numpy.ndarray`, optional
+        optional mask of voxels to include in the NDC calculation
+
+    Returns
+    -------
+    :obj:`float`
+        The NDC value.
+
+    References
+    ----------
+    .. [Yeh2019] Yeh, Fang-Cheng, et al. "Differential tractography as a
+                 track-based biomarker for neuronal injury."
+                 NeuroImage 202 (2019): 116131.
+
+    Notes
+    -----
+    This is a copy of DIPY's code to be removed (and just imported) as soon as
+    a new release of DIPY is cut including
+    `dipy/dipy#3156 <https://github.com/dipy/dipy/pull/3156>`__.
+
+    """
+
+    neighbor_correlations = []
+
+    mask = np.ones_like(dwi_data[..., 0], dtype=bool) if mask is None else mask
+
+    dwi_data = dwi_data[mask]
+
+    for from_index, to_index in neighbor_indices:
+        flat_from_image = dwi_data[from_index]
+        flat_to_image = dwi_data[to_index]
+
+        neighbor_correlations.append(
+            np.corrcoef(flat_from_image, flat_to_image)[0, 1]
+        )
+
+    return np.round(np.mean(neighbor_correlations), 4)

mriqc/workflows/diffusion/base.py

@@ -249,6 +249,7 @@ def dmri_qc_workflow(name='dwiMRIQC'):
         (hmcwf, get_hmc_shells, [('outputnode.out_file', 'in_file')]),
         (dti, cc_mask, [('out_fa', 'in_fa'),
                         ('out_cfa', 'in_cfa')]),
+        (meta, iqms_wf, [('qspace_neighbors', 'inputnode.qspace_neighbors')]),
         (averages, iqms_wf, [(('out_file', _first), 'inputnode.in_b0')]),
         (sp_mask, iqms_wf, [('out_mask', 'inputnode.spikes_mask')]),
         (hmcwf, iqms_wf, [('outputnode.out_fd', 'inputnode.framewise_displacement')]),
@@ -322,6 +323,7 @@ def compute_iqms(name='ComputeIQMs'):
                 'cc_mask',
                 'spikes_mask',
                 'framewise_displacement',
+                'qspace_neighbors',
             ]
         ),
         name='inputnode',
@@ -377,7 +379,8 @@ def compute_iqms(name='ComputeIQMs'):
                                ('in_cfa', 'in_cfa'),
                                ('in_fa_nans', 'in_fa_nans'),
                                ('in_fa_degenerate', 'in_fa_degenerate'),
-                               ('framewise_displacement', 'in_fd')]),
+                               ('framewise_displacement', 'in_fd'),
+                               ('qspace_neighbors', 'qspace_neighbors')]),
         (inputnode, addprov, [('in_file', 'in_file')]),
         (addprov, datasink, [('out_prov', 'provenance')]),
         (meta, datasink, [('subject', 'subject_id'),