Merge pull request #3156 from PennLINC/add_ndc

dipy/stats/meson.build

@@ -2,6 +2,7 @@ python_sources = [
   '__init__.py',
   'analysis.py',
   'resampling.py',
+  'qc.py',
   ]
 
 

dipy/stats/qc.py

@@ -0,0 +1,116 @@
+import numpy as np
+from dipy.core.geometry import cart_distance
+
+
+def find_qspace_neighbors(gtab):
+    """Create a mapping of dwi volume index to its nearest neighbor.
+
+    An approximate q-space is used (the deltas are not included).
+    Note that neighborhood is not necessarily bijective. One neighbor
+    is found per dwi volume.
+
+    Parameters
+    ----------
+
+    gtab: dipy.core.gradients.GradientTable
+        Gradient table.
+
+    Returns
+    -------
+
+    neighbors: list of tuple
+        A list of 2-tuples indicating the nearest q-space neighbor
+        of each dwi volume.
+
+    Examples
+    --------
+    >>> from dipy.core.gradients import gradient_table
+    >>> import numpy as np
+    >>> gtab = gradient_table(
+    ...     np.array([0, 1000, 1000, 2000]),
+    ...     np.array([
+    ...         [1, 0, 0],
+    ...         [1, 0, 0],
+    ...         [0.99, 0.0001, 0.0001],
+    ...         [1, 0, 0]]))
+    >>> find_qspace_neighbors(gtab)
+    [(1, 2), (2, 1), (3, 1)]
+
+    """
+    dwi_neighbors = []
+
+    # Only correlate the b>0 images
+    dwi_mask = np.logical_not(gtab.b0s_mask)
+    dwi_indices = np.flatnonzero(dwi_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
+
+
+def neighboring_dwi_correlation(dwi_data, gtab, mask=None):
+    """Calculate the Neighboring DWI Correlation (NDC) from dMRI data.
+
+    Using a mask is highly recommended, otherwise the FOV will influence the
+    correlations. According to [Yeh2019], an NDC less than 0.4 indicates a
+    low quality image.
+
+    Parameters
+    ----------
+    dwi_data : 4D ndarray
+        dwi data on which to calculate NDC
+    gtab : dipy.core.gradients.GradientTable
+        Gradient table.
+    mask : 3D ndarray, optional
+        optional mask of voxels to include in the NDC calculation
+
+    Returns
+    -------
+    ndc : float
+        The neighboring DWI correlation
+
+    References
+    ----------
+
+    .. [Yeh2019] Yeh, Fang-Cheng, et al. "Differential tractography as a
+                 track-based biomarker for neuronal injury."
+                 NeuroImage 202 (2019): 116131.
+
+    """
+
+    neighbor_indices = find_qspace_neighbors(gtab)
+    neighbor_correlations = []
+
+    if mask is not None:
+        binary_mask = mask > 0
+
+    for from_index, to_index in neighbor_indices:
+
+        # Flatten the dwi images
+        if mask is not None:
+            flat_from_image = dwi_data[..., from_index][binary_mask]
+            flat_to_image = dwi_data[..., to_index][binary_mask]
+        else:
+            flat_from_image = dwi_data[..., from_index].flatten()
+            flat_to_image = dwi_data[..., to_index].flatten()
+
+        neighbor_correlations.append(
+            np.corrcoef(flat_from_image, flat_to_image)[0, 1])
+
+    return np.mean(neighbor_correlations)

dipy/stats/tests/meson.build

@@ -2,6 +2,7 @@ python_sources = [
   '__init__.py',
   'test_analysis.py',
   'test_resampling.py',
+  'test_qc.py',
   ]
 
 

dipy/stats/tests/test_qc.py

@@ -0,0 +1,135 @@
+import numpy as np
+from dipy.stats.qc import neighboring_dwi_correlation
+from dipy.core.gradients import gradient_table
+from dipy.core.geometry import normalized_vector
+
+rng = np.random.default_rng()
+
+
+def create_test_data(test_r, cube_size, mask_size, num_dwi_vols, num_b0s):
+    """Create testing data with a known neighbor structure and a known NDC.
+
+    The b>0 images have 2 images per shell, separated by a very small angle,
+    guaranteeing they will be neighbors. The within-mask data is filled with
+    random data with correlation value of approximately ``test_r``.
+
+    Parameters
+    ----------
+
+    test_r: float
+        The approximate NDC that the simulated data should have
+    cube_size: int
+        The simulated data will be a cube with this many voxels per dim
+    mask_size: int
+        A cubic "brain" is this size per side and filled with data. Must
+        be less than ``cube_size``
+    num_dwi_vols: int
+        The number of b>0 images to simulate. Must be even to ensure we
+        can make known neighbors
+    num_b0s: int
+        The number of b=0 images to prepend to the b>0 images
+
+    Returns
+    -------
+
+    real_r: float
+        The ground-truth neighbor correlation of the simulated data
+    dwi_data: np.ndarray
+        A 4D array containing simulated data
+    mask_data: np.ndarray
+        A 3D array indicating which voxels in ``dwi_data`` contain
+        brain data
+    gtab: dipy.core.gradients.GradientTable
+        Gradient table with known neighbors
+
+    """
+
+    if not num_dwi_vols % 2 == 0:
+        raise Exception(
+            "Needs an even number of dwi vols to ensure known neighbors")
+
+    # Create a volume mask
+    test_mask = np.zeros((cube_size, cube_size, cube_size))
+    test_mask[:mask_size, :mask_size, :mask_size] = 1
+    n_voxels_in_mask = mask_size ** 3
+
+    # 4D Data array
+    dwi_data = np.zeros(
+        (cube_size, cube_size, cube_size, num_b0s + num_dwi_vols))
+
+    # Create a sampling scheme where we know what volumes will be neighbors
+    n_known = num_dwi_vols // 2
+    dwi_bvals = np.column_stack(
+        [np.arange(n_known) + 1] * 2).flatten(order="C").tolist()
+    bvals = np.array([0] * num_b0s + dwi_bvals) * 1000
+
+    # The bvecs will be a straight line with a minor perturbance every other
+    ref_vec = np.array([1., 0., 0.])
+    nbr_vec = normalized_vector(ref_vec + 0.00001)
+    bvecs = np.row_stack(
+        [ref_vec] * num_b0s + [np.row_stack([ref_vec, nbr_vec])] * n_known)
+
+    cor = np.ones((2, 2)) * test_r
+    np.fill_diagonal(cor, 1)
+    L = np.linalg.cholesky(cor)
+
+    known_correlations = []
+    for starting_vol in np.arange(n_known) * 2 + num_b0s:
+        uncorrelated = rng.standard_normal((2, n_voxels_in_mask))
+        correlated = np.dot(L, uncorrelated)
+
+        dwi_data[:, :, :, starting_vol][test_mask > 0] = correlated[0]
+        dwi_data[:, :, :, starting_vol+1][test_mask > 0] = correlated[1]
+
+        known_correlations += [np.corrcoef(correlated)[0, 1]] * 2
+
+    gtab = gradient_table(bvals, bvecs, b0_threshold=50)
+
+    return np.mean(known_correlations), dwi_data, test_mask, gtab
+
+
+def test_neighboring_dwi_correlation():
+    """Test NDC under various conditions."""
+
+    # Test data with b=0s, low correlation, using mask
+    real_r, dwi_data, mask, gtab = create_test_data(
+        test_r=0.3,
+        cube_size=10,
+        mask_size=6,
+        num_dwi_vols=10,
+        num_b0s=2)
+    estimated_ndc = neighboring_dwi_correlation(dwi_data, gtab, mask)
+    assert np.allclose(real_r, estimated_ndc)
+
+    maskless_ndc = neighboring_dwi_correlation(dwi_data, gtab)
+    assert maskless_ndc != real_r
+
+    # Try with no b=0s
+    real_r, dwi_data, mask, gtab = create_test_data(
+        test_r=0.3,
+        cube_size=10,
+        mask_size=6,
+        num_dwi_vols=10,
+        num_b0s=0)
+    estimated_ndc = neighboring_dwi_correlation(dwi_data, gtab, mask)
+    assert np.allclose(real_r, estimated_ndc)
+
+    # Try with realistic correlation value
+    real_r, dwi_data, mask, gtab = create_test_data(
+        test_r=0.8,
+        cube_size=10,
+        mask_size=6,
+        num_dwi_vols=10,
+        num_b0s=2)
+    estimated_ndc = neighboring_dwi_correlation(dwi_data, gtab, mask)
+    assert np.allclose(real_r, estimated_ndc)
+
+    # Try with a bigger volume, lower correlation
+    real_r, dwi_data, mask, gtab = create_test_data(
+        test_r=0.5,
+        cube_size=100,
+        mask_size=49,
+        num_dwi_vols=160,
+        num_b0s=2)
+    estimated_ndc = neighboring_dwi_correlation(dwi_data, gtab, mask)
+    assert np.allclose(real_r, estimated_ndc)