remove old dipy funcs

docs/conf.py

@@ -25,7 +25,7 @@
 # -- General configuration ------------------------------------------------
 # from recommonmark.parser import CommonMarkParser
 
-autodoc_mock_imports = ["scipy, spams, dipy, nibabel"]
+autodoc_mock_imports = ["scipy, spams, nibabel"]
 language = 'en'
 
 # Add any Sphinx extension module names here, as strings. They can be

nlsam/angular_tools.py

@@ -1,6 +1,5 @@
 import numpy as np
 
-
 def angular_neighbors(vec, n):
     """
     Returns the indices of the n closest neighbors (excluding the vector itself)
@@ -131,3 +130,41 @@ def split_per_shell(bvals, bvecs, angular_size, dwis, is_symmetric=False, bval_t
         current_shell = [(dwi,) + tuple(current_shell[pos]) for pos, dwi in enumerate(new_positions) if dwi in dwis]
         neighbors[shell] = current_shell
     return neighbors
+
+
+def read_bvals_bvecs(fbvals, fbvecs):
+    """Read bvalues and bvectors from disk.
+
+    Parameters
+    ----------
+    fbvals : str
+       Path of bvalues
+    fbvecs : str
+       Path of bvectors
+
+    Returns
+    -------
+    bvals : array, (N,)
+    bvecs : array, (N, 3)
+    """
+
+    bvals = np.loadtxt(fbvals)
+    bvecs = np.loadtxt(fbvecs)
+
+    if 3 not in bvecs.shape:
+        raise ValueError("bvec file should have three rows")
+
+    if bvecs.shape[1] != 3:
+        bvecs = bvecs.T
+
+    if len(bvals.shape) > 1:
+        raise ValueError("bval file should have one row")
+
+    if bvals.shape[0] != bvecs.shape[0]:
+        raise ValueError(f"bvalues shape {bvals.shape} and bvectors shape {bvecs.shape} do not correspond")
+
+    norm = np.linalg.norm(bvecs, axis=1, keepdims=True)
+    norm[norm == 0] = 1
+    bvecs /= norm
+
+    return bvals, bvecs

nlsam/smoothing.py

@@ -4,9 +4,6 @@
 from joblib import Parallel, delayed
 from tqdm.autonotebook import tqdm
 
-from dipy.core.geometry import cart2sphere
-from dipy.reconst.shm import sph_harm_ind_list, real_sph_harm, smooth_pinv
-
 from scipy.ndimage import convolve, gaussian_filter
 
 logger = logging.getLogger('nlsam')
@@ -49,7 +46,6 @@ def sh_smooth(data, bvals, bvecs, sh_order=4, b0_threshold=1.0, similarity_thres
     where_b0s = bvals <= b0_threshold
     pred_sig = np.zeros_like(data, dtype=np.float32)
 
-
     # Round similar bvals together for identifying similar shells
     rounded_bvals = np.zeros_like(bvals)
 
@@ -69,16 +65,16 @@ def sh_smooth(data, bvals, bvecs, sh_order=4, b0_threshold=1.0, similarity_thres
                 pred_sig[..., idx] = data[..., idx]
             continue
 
-
         x, y, z = bvecs[:, idx]
         _, theta, phi = cart2sphere(x, y, z)
 
-        # Find the sh coefficients to smooth the signal
-        B_dwi = real_sph_harm(m, n, theta[:, None], phi[:, None])
+        # Find the sh coefficients to predict the signal
+        B_dwi = real_sh_descoteaux_from_index(m, n, theta[:, None], phi[:, None])
         invB = smooth_pinv(B_dwi, np.sqrt(regul) * L)
-        sh_coeff = np.dot(data[..., idx], invB.T)
-        # Find the smoothed signal from the sh fit for the given gtab
-        pred_sig[..., idx] = np.dot(sh_coeff, B_dwi.T)
+        sh_coeff = data[..., idx] @ invB.T
+
+        # Find the predicted signal from the sh fit for the given bvecs
+        pred_sig[..., idx] = sh_coeff @ B_dwi.T
     return pred_sig
 
 
@@ -170,3 +166,183 @@ def local_standard_deviation(arr, n_cores=-1, verbose=False):
     blur = full_width_at_half_max / np.sqrt(8 * np.log(2))
 
     return gaussian_filter(sigma, blur, mode='reflect')
+
+
+# Copyright (c) 2008-2024, dipy developers
+# All rights reserved.
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+
+#     * Redistributions of source code must retain the above copyright
+#        notice, this list of conditions and the following disclaimer.
+
+#     * Redistributions in binary form must reproduce the above
+#        copyright notice, this list of conditions and the following
+#        disclaimer in the documentation and/or other materials provided
+#        with the distribution.
+
+#     * Neither the name of the dipy developers nor the names of any
+#        contributors may be used to endorse or promote products derived
+#        from this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+def cart2sphere(x, y, z):
+    r"""Return angles for Cartesian 3D coordinates `x`, `y`, and `z`
+
+    This is the standard physics convention where `theta` is the
+    inclination (polar) angle, and `phi` is the azimuth angle.
+
+    $0\le\theta\mathrm{(theta)}\le\pi$ and $-\pi\le\phi\mathrm{(phi)}\le\pi$
+
+    Parameters
+    ----------
+    x : array_like
+       x coordinate in Cartesian space
+    y : array_like
+       y coordinate in Cartesian space
+    z : array_like
+       z coordinate
+
+    Returns
+    -------
+    r : array
+       radius
+    theta : array
+       inclination (polar) angle
+    phi : array
+       azimuth angle
+
+    """
+    r = np.sqrt(x * x + y * y + z * z)
+    cos = np.divide(z, r, where=r > 0)
+    theta = np.arccos(cos, where=(cos >= -1) & (cos <= 1))
+    theta = np.where(r > 0, theta, 0.0)
+    phi = np.arctan2(y, x)
+    r, theta, phi = np.broadcast_arrays(r, theta, phi)
+    return r, theta, phi
+
+
+def sph_harm_ind_list(sh_order_max):
+    """
+    Returns the order (``l``) and phase_factor (``m``) of all the symmetric
+    spherical harmonics of order less then or equal to ``sh_order_max``.
+    The results, ``m_list`` and ``l_list`` are kx1 arrays, where k depends on
+    ``sh_order_max``.
+
+    Parameters
+    ----------
+    sh_order_max : int
+        The maximum order ($l$) of the spherical harmonic basis.
+        Even int > 0, max order to return
+
+    Returns
+    -------
+    m_list : array of int
+        phase factors ($m$) of even spherical harmonics
+    l_list : array of int
+        orders ($l$) of even spherical harmonics
+    """
+    if sh_order_max % 2 != 0:
+        raise ValueError("sh_order_max must be an even integer >= 0")
+    l_range = np.arange(0, sh_order_max + 1, 2, dtype=int)
+    ncoef = int((sh_order_max + 2) * (sh_order_max + 1) // 2)
+
+    l_list = np.repeat(l_range, l_range * 2 + 1)
+    offset = 0
+    m_list = np.empty(ncoef, "int")
+    for ii in l_range:
+        m_list[offset : offset + 2 * ii + 1] = np.arange(-ii, ii + 1)
+        offset = offset + 2 * ii + 1
+
+    # makes the arrays ncoef by 1, allows for easy broadcasting later in code
+    return m_list, l_list
+
+
+def real_sh_descoteaux_from_index(m_values, l_values, theta, phi):
+    r"""Compute real spherical harmonics.
+
+    The definition adopted here follows Descoteaux2007 where the
+    real harmonic $Y_l^m$ is defined to be:
+
+    .. math::
+       :nowrap:
+
+        Y_l^m =
+        \begin{cases}
+            \sqrt{2} * \Im(Y_l^m) \; if m > 0 \\
+            Y^0_l \; if m = 0 \\
+            \sqrt{2} * \Re(Y_l^m)  \; if m < 0 \\
+        \end{cases}
+
+    This may take scalar or array arguments. The inputs will be broadcast
+    against each other.
+
+    Parameters
+    ----------
+    m_values : array of int ``|m| <= l``
+        The phase factors ($m$) of the harmonics.
+    l_values : array of int ``l >= 0``
+        The orders ($l$) of the harmonics.
+    theta : float [0, pi]
+        The polar (colatitudinal) coordinate.
+    phi : float [0, 2*pi]
+        The azimuthal (longitudinal) coordinate.
+
+    Returns
+    -------
+    real_sh : real float
+        The real harmonic $Y_l^m$ sampled at ``theta`` and ``phi``.
+    """
+
+    def spherical_harmonics(m_values, l_values, theta, phi):
+         from scipy.special import sph_harm
+         return sph_harm(m_values, l_values, theta, phi, dtype=complex)
+
+    sh = spherical_harmonics(np.abs(m_values), l_values, phi, theta)
+
+    real_sh = np.where(m_values > 0, sh.imag, sh.real)
+    real_sh *= np.where(m_values == 0, 1.0, np.sqrt(2))
+
+    return real_sh
+
+
+def smooth_pinv(B, L):
+    """Regularized pseudo-inverse
+
+    Computes a regularized least square inverse of B
+
+    Parameters
+    ----------
+    B : array_like (n, m)
+        Matrix to be inverted
+    L : array_like (m,)
+
+    Returns
+    -------
+    inv : ndarray (m, n)
+        regularized least square inverse of B
+
+    Notes
+    -----
+    In the literature this inverse is often written $(B^{T}B+L^{2})^{-1}B^{T}$.
+    However here this inverse is implemented using the pseudo-inverse because
+    it is more numerically stable than the direct implementation of the matrix
+    product.
+
+    """
+    L = np.diag(L)
+    inv = np.linalg.pinv(np.concatenate((B, L)))
+    return inv[:, : len(B)]

nlsam/tests/test_angular_tools.py

@@ -1,8 +1,8 @@
 import numpy as np
-from numpy.testing import assert_almost_equal, assert_equal
-
-from nlsam.angular_tools import angular_neighbors, _angle, greedy_set_finder, split_per_shell
+from numpy.testing import assert_allclose, assert_equal
 
+from nlsam.angular_tools import angular_neighbors, _angle, greedy_set_finder, split_per_shell, read_bvals_bvecs
+from os.path import join, dirname
 
 def test_angular_neighbors():
 
@@ -29,7 +29,7 @@ def test_angle():
                [np.pi / 2,     np.pi / 2, np.pi / 2],
                [np.pi / 2,     np.pi / 2, 0]]
 
-    assert_almost_equal(angles, correct)
+    assert_allclose(angles, correct)
 
 
 def test_greedy_set_finder():
@@ -87,3 +87,17 @@ def test_split_per_shell():
 
     for n, ii in enumerate(idx):
         assert_equal(np.unique(greedy_set_finder(ii)), true_idx[n])
+
+
+def test_read_bvals_bvecs():
+    path = dirname(__file__)
+    fbvals = join(path, '../../example', 'bvals')
+    fbvecs = join(path, '../../example', 'bvecs')
+
+    bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
+    zeros = bvals == 0
+    nonzeros = np.logical_not(zeros)
+    norm = np.linalg.norm(bvecs, axis=1)
+
+    assert_allclose(norm[nonzeros], 1.0)
+    assert_allclose(norm[zeros], 0)