# Source code for dipy.stats.qc

    import numpy as np

    from dipy.core.geometry import cart_distance
    from dipy.testing.decorators import warning_for_keywords

    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]),
        ...     bvecs=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

    @warning_for_keywords()
    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 :footcite:t:`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
            Mask of voxels to include in the NDC calculation

        Returns
        -------
        ndc : float
            The neighboring DWI correlation

        References
        ----------
        .. footbibliography::

        """

        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)