updated uncertainty calculation with different parameters

fury/actor.py

@@ -10,7 +10,7 @@
 from fury import layout
 from fury.actors.odf_slicer import OdfSlicerActor
 from fury.actors.peak import PeakActor
-from fury.actors.tensor import uncertainty_cone
+from fury.actors.tensor import double_cone, main_dir_uncertainty
 from fury.colormap import colormap_lookup_table
 from fury.deprecator import deprecate_with_version, deprecated_params
 from fury.io import load_image
@@ -3798,10 +3798,12 @@ def callback(
     return sq_actor
 
 
-def dti_uncertainty(
-    data,
-    bvals,
-    bvecs,
+def uncertainty_cone(
+    evals,
+    evecs,
+    signal,
+    sigma,
+    b_matrix,
     scales=.6,
     opacity=1.0
 ):
@@ -3811,37 +3813,43 @@ def dti_uncertainty(
 
     Parameters
     ----------
-    data : 3D or 4D ndarray
-        Diffusion data.
-    bvals : array, (N,) or None
-        Array containing the b-values.
-    bvecs : array, (N, 3) or None
-        Array containing the b-vectors.
+    evals : ndarray (3, ) or (N, 3)
+        Eigenvalues.
+    evecs : ndarray (3, 3) or (N, 3, 3)
+        Eigenvectors.
+    signal : 3D or 4D ndarray
+        Predicted signal.
+    sigma : ndarray
+        Standard deviation of the noise.
+    b_matrix : array (N, 7)
+        Design matrix for DTI.
     scales : float or ndarray (N, ), optional
         Cones of uncertainty size, default(1.0).
     opacity : float, optional
         Takes values from 0 (fully transparent) to 1 (opaque), default(1.0).
 
     Returns
     -------
-    uncertainty_cone: Actor
-
-    Examples
-    --------
-    >>> from fury import actor, window
-    >>> from dipy.io.image import load_nifti
-    >>> from dipy.io.gradients import read_bvals_bvecs
-    >>> from dipy.data import get_fnames
-    >>> scene = window.Scene()
-    >>> hardi_fname, hardi_bval_fname, hardi_bvec_fname =\
-    >>>     get_fnames("stanford_hardi")
-    >>> data, _ = load_nifti(hardi_fname)
-    >>> bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
-    >>> uncertainty_cones = actor.dti_uncertainty(
-    >>>     data=data[13:43, 44:74, 28:29], bvals=bvals, bvecs=bvecs)
-    >>> scene.add(uncertainty_cones)
-    >>> #window.show(scene)
+    double_cone: Actor
 
     """
 
-    return uncertainty_cone(data, bvals, bvecs, scales, opacity)
+    valid_mask = np.abs(evecs).max(axis=(-2, -1)) > 0
+    indices = np.nonzero(valid_mask)
+
+    evecs = evecs[indices]
+    evals = evals[indices]
+    signal = signal[indices]
+
+    centers = np.asarray(indices).T
+    colors = np.array([107, 107, 107])
+
+    x, y, z = evecs.shape
+    if not isinstance(scales, np.ndarray):
+        scales = np.array(scales)
+    if scales.size == 1:
+        scales = np.repeat(scales, x)
+
+    angles = main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix)
+
+    return double_cone(centers, evecs, angles, colors, scales, opacity)

fury/actors/tensor.py

@@ -1,57 +1,12 @@
 import os
 
-from dipy.core.gradients import gradient_table
-from dipy.reconst import dti
-
-import dipy.denoise.noise_estimate as ne
-
 import numpy as np
 
 from fury import actor
 from fury.shaders import (attribute_to_actor, import_fury_shader,
                           shader_to_actor, compose_shader)
 
 
-def uncertainty_cone(data, bvals, bvecs, scales, opacity):
-    """
-    Visualize the cones of uncertainty for DTI.
-
-    Parameters
-    ----------
-    data : 3D or 4D ndarray
-        Diffusion data.
-    bvals : array, (N,) or None
-        Array containing the b-values.
-    bvecs : array, (N, 3) or None
-        Array containing the b-vectors.
-    scales : float or ndarray (N, )
-        Cones of uncertainty size.
-    opacity : float, optional
-        Takes values from 0 (fully transparent) to 1 (opaque).
-
-    Returns
-    -------
-    box_actor: Actor
-
-    """
-    angles, centers, axes = main_dir_uncertainty(data, bvals, bvecs)
-    colors = np.array([107, 107, 107])
-
-    if centers.ndim != 2:
-        centers = np.array([centers])
-        axes = np.array([axes])
-        colors = np.array([colors])
-
-    x, y, z = axes.shape
-
-    if not isinstance(scales, np.ndarray):
-        scales = np.array(scales)
-    if scales.size == 1:
-        scales = np.repeat(scales, x)
-
-    return double_cone(centers, axes, angles, colors, scales, opacity)
-
-
 def double_cone(centers, axes, angles, colors, scales, opacity):
     """
     Visualize one or many Double Cones with different features.
@@ -243,24 +198,27 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
     return box_actor
 
 
-def main_dir_uncertainty(data, bvals, bvecs):
+def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
     """
     Calculates the angle of the cone of uncertainty that represents the
     perturbation of the main eigenvector of the diffusion tensor matrix.
-    Additionally, it gives information needed for the cone visualization.
 
     Parameters
     ----------
-    data : 3D or 4D ndarray
-        Diffusion data.
-    bvals : array, (N,) or None
-        Array containing the b-values.
-    bvecs : array, (N, 3) or None
-        Array containing the b-vectors.
+    evals : ndarray (3, ) or (N, 3)
+        Eigenvalues.
+    evecs : ndarray (3, 3) or (N, 3, 3)
+        Eigenvectors.
+    signal : 3D or 4D ndarray
+        Predicted signal.
+    sigma : ndarray
+        Standard deviation of the noise.
+    b_matrix : array (N, 7)
+        Design matrix for DTI.
 
     Returns
     -------
-    angles, centers, evecs
+    angles: array
 
     Notes
     -----
@@ -292,40 +250,16 @@ def main_dir_uncertainty(data, bvals, bvecs):
     International Society for Magnetic Resonance in Medicine, 57(1), 141-149.
     """
 
-    gtab = gradient_table(bvals, bvecs)
-
-    tenmodel = dti.TensorModel(gtab)
-    tenfit = tenmodel.fit(data)
-    fevals = tenfit.evals
-    fevecs = tenfit.evecs
-
-    tensor_vals = dti.lower_triangular(tenfit.quadratic_form)
-    dti_params = dti.eig_from_lo_tri(tensor_vals)
-
-    fsignal = dti.tensor_prediction(dti_params, gtab, 1.0)
-    b_matrix = dti.design_matrix(gtab)
-
-    valid_mask = np.abs(fevecs).max(axis=(-2, -1)) > 0
-    indices = np.nonzero(valid_mask)
-
-    evecs = fevecs[indices]
-    evals = fevals[indices]
-    signal = fsignal[indices]
-
-    # Uncertainty calculations
-
-    sigma = ne.estimate_sigma(data)  # standard deviation of the noise
-
-    # Angles for cone of uncertainty
     angles = np.ones(evecs.shape[0])
     for i in range(evecs.shape[0]):
         sigma_e = np.diag(signal[i] / sigma ** 2)
         k = np.dot(np.transpose(b_matrix), sigma_e)
         sigma_ = np.dot(k, b_matrix)
 
         dd = np.diag(sigma_)
-        delta_DD = dti.from_lower_triangular(
-            np.array([dd[0], dd[3], dd[1], dd[4], dd[5], dd[2]]))
+        delta_DD = np.array([[dd[0], dd[3], dd[4]],
+                             [dd[3], dd[1], dd[5]],
+                             [dd[4], dd[5], dd[2]]])
 
         # perturbation matrix of tensor D
         try:
@@ -355,5 +289,4 @@ def main_dir_uncertainty(data, bvals, bvecs):
         else:
             theta = 1.39626
 
-    centers = np.asarray(indices).T
-    return angles, centers, evecs
+    return angles

fury/actors/tests/test_uncertainty.py

@@ -3,9 +3,12 @@
 visualization of the cones of uncertainty along with the diffusion tensors for
 comparison
 """
+from dipy.core.gradients import gradient_table
 from dipy.reconst import dti
 from dipy.segment.mask import median_otsu
 
+import dipy.denoise.noise_estimate as ne
+
 from fury import actor, window
 
 from dipy.io.image import load_nifti
@@ -30,8 +33,30 @@ def test_uncertainty():
     maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
                                  numpass=1, autocrop=True, dilate=2)
 
-    uncertainty_cones = actor.dti_uncertainty(
-        data=maskdata[13:43, 44:74, 28:29], bvals=bvals, bvecs=bvecs)
+    gtab = gradient_table(bvals, bvecs)
+
+    tenmodel = dti.TensorModel(gtab)
+    tenfit = tenmodel.fit(maskdata[13:43, 44:74, 28:29])
+
+    # Eigenvalues and eigenvectors
+    fevals = tenfit.evals
+    fevecs = tenfit.evecs
+
+    tensor_vals = dti.lower_triangular(tenfit.quadratic_form)
+    dti_params = dti.eig_from_lo_tri(tensor_vals)
+
+    # Predicted signal given tensor parameters
+    fsignal = dti.tensor_prediction(dti_params, gtab, 1.0)
+
+    # Design matrix or B matrix
+    b_matrix = dti.design_matrix(gtab)
+
+    # Standard deviation of the noise
+    sigma = ne.estimate_sigma(maskdata[13:43, 44:74, 28:29])
+
+    uncertainty_cones = actor.uncertainty_cone(evecs=fevecs, evals=fevals,
+                                               signal=fsignal, sigma=sigma,
+                                               b_matrix=b_matrix)
 
     scene = window.Scene()
     scene.background([255, 255, 255])