From afa0f5fe1e916a36c8e18094cf85cefa420e70c8 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 26 Jun 2023 17:51:01 -0500
Subject: [PATCH 01/13] added cone sdf implementation

---
 fury/actors/tensor.py          | 213 +++++++++++++++++++++++++++++++++
 fury/shaders/sdf/sd_cone.frag  |   5 +
 fury/shaders/sdf/sd_union.frag |   4 +
 3 files changed, 222 insertions(+)
 create mode 100644 fury/actors/tensor.py
 create mode 100644 fury/shaders/sdf/sd_cone.frag
 create mode 100644 fury/shaders/sdf/sd_union.frag

diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
new file mode 100644
index 000000000..53032be6e
--- /dev/null
+++ b/fury/actors/tensor.py
@@ -0,0 +1,213 @@
+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 double_cone(centers, axes, lengths, angles, colors, scales, opacity):
+    """
+    Visualize one or many Double Cones with different features.
+
+    Parameters
+    ----------
+    centers : ndarray(N, 3)
+    axes : ndarray (3, 3) or (N, 3, 3)
+    lengths : ndarray (3, ) or (N, 3)
+    angles : float or ndarray (N, )
+    colors : ndarray (N,3) or tuple (3,), optional
+    scales : float or ndarray (N, ), optional
+    opacity : float, optional
+
+    Returns
+    -------
+    box_actor: Actor
+
+    """
+
+    box_actor = actor.box(centers, colors=colors, scales=scales)
+    box_actor.GetMapper().SetVBOShiftScaleMethod(False)
+    box_actor.GetProperty().SetOpacity(opacity)
+
+    # Number of vertices that make up the box
+    n_verts = 8
+
+    big_centers = np.repeat(centers, n_verts, axis=0)
+    attribute_to_actor(box_actor, big_centers, 'center')
+
+    big_scales = np.repeat(scales, n_verts, axis=0)
+    attribute_to_actor(box_actor, big_scales, 'scale')
+
+    big_values = np.repeat(np.array(lengths, dtype=float), n_verts, axis=0)
+    attribute_to_actor(box_actor, big_values, 'evals')
+
+    evec1 = np.array([item[0] for item in axes])
+    evec2 = np.array([item[1] for item in axes])
+    evec3 = np.array([item[2] for item in axes])
+
+    big_vectors_1 = np.repeat(evec1, n_verts, axis=0)
+    attribute_to_actor(box_actor, big_vectors_1, 'evec1')
+    big_vectors_2 = np.repeat(evec2, n_verts, axis=0)
+    attribute_to_actor(box_actor, big_vectors_2, 'evec2')
+    big_vectors_3 = np.repeat(evec3, n_verts, axis=0)
+    attribute_to_actor(box_actor, big_vectors_3, 'evec3')
+
+    big_angles = np.repeat(np.array(angles, dtype=float), n_verts, axis=0)
+    attribute_to_actor(box_actor, big_angles, 'angle')
+
+    # Start of shader implementation
+
+    vs_dec = \
+        """
+        in vec3 center;
+        in float scale;
+        in vec3 evals;
+        in vec3 evec1;
+        in vec3 evec2;
+        in vec3 evec3;
+        in float angle;
+
+        out vec4 vertexMCVSOutput;
+        out vec3 centerMCVSOutput;
+        out float scaleVSOutput;
+        out vec3 evalsVSOutput;
+        out mat3 rotationMatrix;
+        out vec2 angleVSOutput;
+        """
+
+    # Variables assignment
+    v_assign = \
+        """
+        vertexMCVSOutput = vertexMC;
+        centerMCVSOutput = center;
+        scaleVSOutput = scale;
+        """
+
+    # Rotation matrix
+    rot_matrix = \
+        """
+        mat3 R = mat3(normalize(evec1), normalize(evec2), normalize(evec3));
+        float a = radians(90);
+        mat3 rot = mat3(cos(a),-sin(a),0,
+                        sin(a),cos(a), 0, 
+                        0,     0,      1);
+        rotationMatrix = transpose(R) * rot;
+        float ang = clamp(angle, 0, 6.283);
+        angleVSOutput = vec2(sin(ang), cos(ang));
+        """
+
+    vs_impl = compose_shader([v_assign, rot_matrix])
+
+    # Adding shader implementation to actor
+    shader_to_actor(box_actor, 'vertex', decl_code=vs_dec, impl_code=vs_impl)
+
+    fs_vars_dec = \
+        """
+        in vec4 vertexMCVSOutput;
+        in vec3 centerMCVSOutput;
+        in float scaleVSOutput;
+        in vec3 evalsVSOutput;
+        in mat3 rotationMatrix;
+        in vec2 angleVSOutput;
+
+        uniform mat4 MCVCMatrix;
+        """
+
+    # Importing the cone SDF
+    sd_cone = import_fury_shader(os.path.join('sdf', 'sd_cone.frag'))
+
+    # Importing the union SDF operation
+    sd_union = import_fury_shader(os.path.join('sdf', 'sd_union.frag'))
+
+    # SDF definition
+    sdf_map = \
+        """
+            float sdDoubleCone( vec3 p, vec2 c, float h )
+            {
+                return opUnion(sdCone(p,c,h),sdCone(-p,c,h));
+            }
+
+            float map(in vec3 position)
+            {
+                vec3 p = (position - centerMCVSOutput) /
+                                     scaleVSOutput * rotationMatrix;
+                vec2 a = angleVSOutput;
+                float h = .5 * a.y;
+                //return opUnion(sdCone(p,a,h),sdCone(-p,a,h)) * scaleVSOutput;
+                return sdDoubleCone((position - centerMCVSOutput)/scaleVSOutput
+                    *rotationMatrix, a, h) * scaleVSOutput;
+            }
+        """
+
+    # Importing central differences function for computing surface normals
+    central_diffs_normal = import_fury_shader(os.path.join(
+        'sdf', 'central_diffs.frag'))
+
+    # Importing raymarching function
+    cast_ray = import_fury_shader(os.path.join(
+        'ray_marching', 'cast_ray.frag'))
+
+    # Importing Blinn-Phong model for lighting
+    blinn_phong_model = import_fury_shader(os.path.join(
+        'lighting', 'blinn_phong_model.frag'))
+
+    # Full fragment shader declaration
+    fs_dec = compose_shader([fs_vars_dec, sd_cone, sd_union, sdf_map,
+                             central_diffs_normal, cast_ray,
+                             blinn_phong_model])
+
+    shader_to_actor(box_actor, 'fragment', decl_code=fs_dec)
+
+    # Vertex in Model Coordinates.
+    point = "vec3 point = vertexMCVSOutput.xyz;"
+
+    # Camera position in world space
+    ray_origin = "vec3 ro = (-MCVCMatrix[3] * MCVCMatrix).xyz;"
+
+    ray_direction = "vec3 rd = normalize(point - ro);"
+
+    light_direction = "vec3 ld = normalize(ro - point);"
+
+    ray_origin_update = "ro += point - ro;"
+
+    # Total distance traversed along the ray
+    distance = "float t = castRay(ro, rd);"
+
+    # Fragment shader output definition
+    # If surface is detected, color is assigned, otherwise, nothing is painted
+    frag_output_def = \
+        """
+        if(t < 20)
+        {
+            vec3 pos = ro + t * rd;
+            vec3 normal = centralDiffsNormals(pos, .0001);
+            // Light Attenuation
+            float la = dot(ld, normal);
+            vec3 color = blinnPhongIllumModel(la, lightColor0, 
+                diffuseColor, specularPower, specularColor, ambientColor);
+            fragOutput0 = vec4(color, opacity);
+        }
+        else
+        {
+            discard;
+        }
+        """
+
+    # Full fragment shader implementation
+    sdf_frag_impl = compose_shader([point, ray_origin, ray_direction,
+                                    light_direction, ray_origin_update,
+                                    distance, frag_output_def])
+
+    shader_to_actor(box_actor, 'fragment', impl_code=sdf_frag_impl,
+                    block='light')
+
+    return box_actor
+
+
diff --git a/fury/shaders/sdf/sd_cone.frag b/fury/shaders/sdf/sd_cone.frag
new file mode 100644
index 000000000..fe511612e
--- /dev/null
+++ b/fury/shaders/sdf/sd_cone.frag
@@ -0,0 +1,5 @@
+float sdCone( vec3 p, vec2 c, float h )
+{
+  float q = length(p.xz);
+  return max(dot(c.xy,vec2(q,p.y)),-h-p.y);
+}
\ No newline at end of file
diff --git a/fury/shaders/sdf/sd_union.frag b/fury/shaders/sdf/sd_union.frag
new file mode 100644
index 000000000..958c0e403
--- /dev/null
+++ b/fury/shaders/sdf/sd_union.frag
@@ -0,0 +1,4 @@
+float opUnion( float d1, float d2 )
+{
+    return min(d1,d2);
+}
\ No newline at end of file

From 36072de6594d2c4bcf5be6b22a0eae4fecf1a69a Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 26 Jun 2023 18:35:44 -0500
Subject: [PATCH 02/13] added uncertainty calculation

---
 fury/actor.py         |  27 ++++++++
 fury/actors/tensor.py | 149 +++++++++++++++++++++++++++++++-----------
 2 files changed, 139 insertions(+), 37 deletions(-)

diff --git a/fury/actor.py b/fury/actor.py
index e802b4368..a0107cea6 100644
--- a/fury/actor.py
+++ b/fury/actor.py
@@ -10,6 +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.colormap import colormap_lookup_table
 from fury.deprecator import deprecate_with_version, deprecated_params
 from fury.io import load_image
@@ -3795,3 +3796,29 @@ def callback(
     shader_to_actor(sq_actor, 'fragment', impl_code=fs_impl_code, block='light')
 
     return sq_actor
+
+
+def dti_uncertainty(
+    data,
+    bvals,
+    bvecs,
+    scales=.6,
+    opacity=1.0
+):
+    """
+    VTK actor for visualizing tensor ellipsoids.
+
+    Parameters
+    ----------
+    scales : float or ndarray (N, ), optional
+        Tensor size, default(1).
+    opacity : float, optional
+        Takes values from 0 (fully transparent) to 1 (opaque). Default is 1.
+
+    Returns
+    -------
+    uncertainty_cone: Actor
+
+    """
+
+    return uncertainty_cone(data, bvals, bvecs, scales, opacity)
\ No newline at end of file
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index 53032be6e..ed04d0837 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -12,6 +12,12 @@
                           shader_to_actor, compose_shader)
 
 
+def uncertainty_cone(data, bvals, bvecs, scales, opacity):
+    angles, centers, axes, lengths = main_dir_uncertainty(data, bvals, bvecs)
+    colors = np.array([107, 107, 107])
+    return double_cone(centers, axes, lengths, angles, colors, scales, opacity)
+
+
 def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
     """
     Visualize one or many Double Cones with different features.
@@ -79,34 +85,25 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         out float scaleVSOutput;
         out vec3 evalsVSOutput;
         out mat3 rotationMatrix;
-        out vec2 angleVSOutput;
+        out float angleVSOutput;
         """
 
-    # Variables assignment
-    v_assign = \
+    vs_impl = \
         """
         vertexMCVSOutput = vertexMC;
         centerMCVSOutput = center;
         scaleVSOutput = scale;
-        """
-
-    # Rotation matrix
-    rot_matrix = \
-        """
         mat3 R = mat3(normalize(evec1), normalize(evec2), normalize(evec3));
         float a = radians(90);
         mat3 rot = mat3(cos(a),-sin(a),0,
                         sin(a),cos(a), 0, 
                         0,     0,      1);
         rotationMatrix = transpose(R) * rot;
-        float ang = clamp(angle, 0, 6.283);
-        angleVSOutput = vec2(sin(ang), cos(ang));
+        angleVSOutput = angle;
         """
 
-    vs_impl = compose_shader([v_assign, rot_matrix])
-
-    # Adding shader implementation to actor
-    shader_to_actor(box_actor, 'vertex', decl_code=vs_dec, impl_code=vs_impl)
+    shader_to_actor(box_actor, 'vertex', decl_code=vs_dec,
+                    impl_code=vs_impl)
 
     fs_vars_dec = \
         """
@@ -115,35 +112,47 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         in float scaleVSOutput;
         in vec3 evalsVSOutput;
         in mat3 rotationMatrix;
-        in vec2 angleVSOutput;
+        in float angleVSOutput;
 
         uniform mat4 MCVCMatrix;
         """
 
-    # Importing the cone SDF
-    sd_cone = import_fury_shader(os.path.join('sdf', 'sd_cone.frag'))
-
-    # Importing the union SDF operation
-    sd_union = import_fury_shader(os.path.join('sdf', 'sd_union.frag'))
+    sd_sphere = import_fury_shader(os.path.join('sdf', 'sd_sphere.frag'))
 
-    # SDF definition
     sdf_map = \
         """
-            float sdDoubleCone( vec3 p, vec2 c, float h )
-            {
-                return opUnion(sdCone(p,c,h),sdCone(-p,c,h));
-            }
-
-            float map(in vec3 position)
-            {
-                vec3 p = (position - centerMCVSOutput) /
-                                     scaleVSOutput * rotationMatrix;
-                vec2 a = angleVSOutput;
-                float h = .5 * a.y;
-                //return opUnion(sdCone(p,a,h),sdCone(-p,a,h)) * scaleVSOutput;
-                return sdDoubleCone((position - centerMCVSOutput)/scaleVSOutput
-                    *rotationMatrix, a, h) * scaleVSOutput;
-            }
+        float opUnion( float d1, float d2 ) { return min(d1,d2); }
+
+        float sdCone( vec3 p, vec2 c, float h )
+        {
+            // c is the sin/cos of the angle, h is height
+            // Alternatively pass q instead of (c,h),
+            // which is the point at the base in 2D
+            vec2 q = h*vec2(c.x/c.y,-1.0);
+
+            vec2 w = vec2( length(p.xz), p.y );
+            vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
+            vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
+            float k = sign( q.y );
+            float d = min(dot( a, a ),dot(b, b));
+            float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
+            return sqrt(d)*sign(s);
+        }
+
+        float sdDoubleCone( vec3 p, vec2 c, float h )
+        {
+            return opUnion(sdCone(p,c,h),sdCone(-p,c,h));
+        }
+
+        float map(in vec3 position)
+        {
+            float a = clamp(angleVSOutput, 0, 6.283);
+            //float a = angleVSOutput;
+            vec2 angle = vec2(sin(a), cos(a));
+            return sdDoubleCone((position - centerMCVSOutput)/scaleVSOutput
+                *rotationMatrix, angle, .5*angle.y) * scaleVSOutput;
+
+        }
         """
 
     # Importing central differences function for computing surface normals
@@ -159,7 +168,7 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         'lighting', 'blinn_phong_model.frag'))
 
     # Full fragment shader declaration
-    fs_dec = compose_shader([fs_vars_dec, sd_cone, sd_union, sdf_map,
+    fs_dec = compose_shader([fs_vars_dec, sdf_map,
                              central_diffs_normal, cast_ray,
                              blinn_phong_model])
 
@@ -211,3 +220,69 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
     return box_actor
 
 
+def main_dir_uncertainty(data, bvals, bvecs):
+    gtab = gradient_table(bvals, bvecs)
+
+    tenmodel = dti.TensorModel(gtab)
+    tenfit = tenmodel.fit(data)
+    evals = tenfit.evals
+    evecs = tenfit.evecs
+
+    tensor_vals = dti.lower_triangular(tenfit.quadratic_form)
+    dti_params = dti.eig_from_lo_tri(tensor_vals)
+
+    signalOI = dti.tensor_prediction(dti_params, gtab, 1.0)
+    b_matrix = dti.design_matrix(gtab)
+
+    valid_mask = np.abs(evecs).max(axis=(-2, -1)) > 0
+    indices = np.nonzero(valid_mask)
+
+    dofs_vecs = evecs[indices]
+    dofs_vals = evals[indices]
+    signal = signalOI[indices]
+
+    # Uncertainty calculations ------------------------------------------------
+
+    sigma = ne.estimate_sigma(data)  # standard deviation of the noise
+
+    # Angles for cone of uncertainty ------------------------------------------
+
+    angles = np.ones(dofs_vecs.shape[0])
+    for i in range(dofs_vecs.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]]))
+
+        try:
+            delta_D = np.linalg.inv(delta_DD)
+        except:
+            delta_D = delta_DD
+
+        D_ = dofs_vecs
+        eigen_vals = dofs_vals[i]
+
+        e1, e2, e3 = np.array(D_[i, :, 0]), np.array(D_[i, :, 1]), \
+                     np.array(D_[i, :, 2])
+        lambda1, lambda2, lambda3 = eigen_vals[0], eigen_vals[1], eigen_vals[2]
+
+        if (lambda1 > lambda2 and lambda1 > lambda3):
+            # The perturbation of the eigenvector associated with the largest
+            # eigenvalue is given by
+            a = np.dot(np.outer(np.dot(e1, delta_D), np.transpose(e2)) /
+                       (lambda1 - lambda2), e2)
+            b = np.dot(np.outer(np.dot(e1, delta_D), np.transpose(e3)) /
+                       (lambda1 - lambda3), e3)
+            delta_e1 = a + b
+
+            # The angle \theta between the perturbed principal eigenvector of D
+            theta = np.arctan(np.linalg.norm(delta_e1))
+            angles[i] = theta
+        else:
+            theta = 0.0872665
+
+    centers = np.asarray(indices).T
+    return angles, centers, dofs_vecs, dofs_vals

From 29649035cf9a76433904933be372c36e94fa1438 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 26 Jun 2023 21:22:02 -0500
Subject: [PATCH 03/13] added initial test

---
 fury/actors/tensor.py                 | 70 ++++++++++++++-------------
 fury/actors/tests/test_uncertainty.py | 36 ++++++++++++++
 fury/shaders/sdf/sd_cone.frag         | 14 +++++-
 3 files changed, 85 insertions(+), 35 deletions(-)
 create mode 100644 fury/actors/tests/test_uncertainty.py

diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index ed04d0837..d19dc1dcb 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -15,6 +15,20 @@
 def uncertainty_cone(data, bvals, bvecs, scales, opacity):
     angles, centers, axes, lengths = main_dir_uncertainty(data, bvals, bvecs)
     colors = np.array([107, 107, 107])
+
+    if centers.ndim != 2:
+        centers = np.array([centers])
+        axes = np.array([axes])
+        lengths = np.array([lengths])
+        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, lengths, angles, colors, scales, opacity)
 
 
@@ -88,20 +102,28 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         out float angleVSOutput;
         """
 
-    vs_impl = \
+    # Variables assignment
+    v_assign = \
         """
         vertexMCVSOutput = vertexMC;
         centerMCVSOutput = center;
         scaleVSOutput = scale;
+        angleVSOutput = angle;
+        """
+
+    # Rotation matrix
+    rot_matrix = \
+        """
         mat3 R = mat3(normalize(evec1), normalize(evec2), normalize(evec3));
         float a = radians(90);
         mat3 rot = mat3(cos(a),-sin(a),0,
                         sin(a),cos(a), 0, 
                         0,     0,      1);
         rotationMatrix = transpose(R) * rot;
-        angleVSOutput = angle;
         """
 
+    vs_impl = compose_shader([v_assign, rot_matrix])
+
     shader_to_actor(box_actor, 'vertex', decl_code=vs_dec,
                     impl_code=vs_impl)
 
@@ -117,41 +139,23 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         uniform mat4 MCVCMatrix;
         """
 
-    sd_sphere = import_fury_shader(os.path.join('sdf', 'sd_sphere.frag'))
+    # Importing the cone SDF
+    sd_cone = import_fury_shader(os.path.join('sdf', 'sd_cone.frag'))
 
+    # Importing the union operation SDF
+    sd_union = import_fury_shader(os.path.join('sdf', 'sd_union.frag'))
+
+    # SDF definition
     sdf_map = \
         """
-        float opUnion( float d1, float d2 ) { return min(d1,d2); }
-
-        float sdCone( vec3 p, vec2 c, float h )
-        {
-            // c is the sin/cos of the angle, h is height
-            // Alternatively pass q instead of (c,h),
-            // which is the point at the base in 2D
-            vec2 q = h*vec2(c.x/c.y,-1.0);
-
-            vec2 w = vec2( length(p.xz), p.y );
-            vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
-            vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
-            float k = sign( q.y );
-            float d = min(dot( a, a ),dot(b, b));
-            float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
-            return sqrt(d)*sign(s);
-        }
-
-        float sdDoubleCone( vec3 p, vec2 c, float h )
-        {
-            return opUnion(sdCone(p,c,h),sdCone(-p,c,h));
-        }
-
         float map(in vec3 position)
         {
-            float a = clamp(angleVSOutput, 0, 6.283);
-            //float a = angleVSOutput;
-            vec2 angle = vec2(sin(a), cos(a));
-            return sdDoubleCone((position - centerMCVSOutput)/scaleVSOutput
-                *rotationMatrix, angle, .5*angle.y) * scaleVSOutput;
-
+            vec3 p = (position - centerMCVSOutput)/scaleVSOutput
+                *rotationMatrix;
+            float angle = clamp(angleVSOutput, 0, 6.283);
+            vec2 a = vec2(sin(angle), cos(angle));
+            float h = .5 * a.y;
+            return opUnion(sdCone(p,a,h), sdCone(-p,a,h)) * scaleVSOutput;
         }
         """
 
@@ -168,7 +172,7 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         'lighting', 'blinn_phong_model.frag'))
 
     # Full fragment shader declaration
-    fs_dec = compose_shader([fs_vars_dec, sdf_map,
+    fs_dec = compose_shader([fs_vars_dec, sd_cone, sd_union, sdf_map,
                              central_diffs_normal, cast_ray,
                              blinn_phong_model])
 
diff --git a/fury/actors/tests/test_uncertainty.py b/fury/actors/tests/test_uncertainty.py
new file mode 100644
index 000000000..191c925f7
--- /dev/null
+++ b/fury/actors/tests/test_uncertainty.py
@@ -0,0 +1,36 @@
+"""
+This spript includes the implementation of cone of uncertainty using matrix
+perturbation analysis
+"""
+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
+
+
+def test_uncertainty():
+    hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
+
+    data, affine = load_nifti(hardi_fname)
+
+    bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
+
+    from dipy.segment.mask import median_otsu
+
+    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)
+
+    scene = window.Scene()
+    scene.background([255, 255, 255])
+
+    scene.add(uncertainty_cones)
+
+    scene.reset_camera()
+    scene.reset_clipping_range()
+
+    window.show(scene, reset_camera=False)
diff --git a/fury/shaders/sdf/sd_cone.frag b/fury/shaders/sdf/sd_cone.frag
index fe511612e..6ff3fd40f 100644
--- a/fury/shaders/sdf/sd_cone.frag
+++ b/fury/shaders/sdf/sd_cone.frag
@@ -1,5 +1,15 @@
 float sdCone( vec3 p, vec2 c, float h )
 {
-  float q = length(p.xz);
-  return max(dot(c.xy,vec2(q,p.y)),-h-p.y);
+    // c is the sin/cos of the angle, h is height
+    // Alternatively pass q instead of (c,h),
+    // which is the point at the base in 2D
+    vec2 q = h*vec2(c.x/c.y,-1.0);
+
+    vec2 w = vec2( length(p.xz), p.y );
+    vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 );
+    vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 );
+    float k = sign( q.y );
+    float d = min(dot( a, a ),dot(b, b));
+    float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
+    return sqrt(d)*sign(s);
 }
\ No newline at end of file

From 9445de65fd3c26417977b50e98dea0061cacfb1c Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Tue, 4 Jul 2023 17:08:27 -0500
Subject: [PATCH 04/13] made adjustments on documentation

---
 fury/actor.py                         |  31 +++++-
 fury/actors/tensor.py                 | 133 +++++++++++++++++++-------
 fury/actors/tests/test_uncertainty.py |  14 ++-
 fury/tests/test_actors.py             |  18 ++++
 4 files changed, 155 insertions(+), 41 deletions(-)

diff --git a/fury/actor.py b/fury/actor.py
index a0107cea6..39b4e95c3 100644
--- a/fury/actor.py
+++ b/fury/actor.py
@@ -3806,19 +3806,42 @@ def dti_uncertainty(
     opacity=1.0
 ):
     """
-    VTK actor for visualizing tensor ellipsoids.
+    VTK actor for visualizing the cone of uncertainty representing the
+    variance of the main direction of diffusion.
 
     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, ), optional
-        Tensor size, default(1).
+        Cones of uncertainty size, default(1.0).
     opacity : float, optional
-        Takes values from 0 (fully transparent) to 1 (opaque). Default is 1.
+        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)
+
     """
 
-    return uncertainty_cone(data, bvals, bvecs, scales, opacity)
\ No newline at end of file
+    return uncertainty_cone(data, bvals, bvecs, scales, opacity)
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index d19dc1dcb..ed5bdf93e 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -13,13 +13,33 @@
 
 
 def uncertainty_cone(data, bvals, bvecs, scales, opacity):
-    angles, centers, axes, lengths = main_dir_uncertainty(data, bvals, bvecs)
+    """
+    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])
-        lengths = np.array([lengths])
         colors = np.array([colors])
 
     x, y, z = axes.shape
@@ -29,22 +49,27 @@ def uncertainty_cone(data, bvals, bvecs, scales, opacity):
     if scales.size == 1:
         scales = np.repeat(scales, x)
 
-    return double_cone(centers, axes, lengths, angles, colors, scales, opacity)
+    return double_cone(centers, axes, angles, colors, scales, opacity)
 
 
-def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
+def double_cone(centers, axes, angles, colors, scales, opacity):
     """
     Visualize one or many Double Cones with different features.
 
     Parameters
     ----------
     centers : ndarray(N, 3)
+        Cone positions.
     axes : ndarray (3, 3) or (N, 3, 3)
-    lengths : ndarray (3, ) or (N, 3)
+        Axes of the cone.
     angles : float or ndarray (N, )
-    colors : ndarray (N,3) or tuple (3,), optional
-    scales : float or ndarray (N, ), optional
-    opacity : float, optional
+        Angles of the cone.
+    colors : ndarray (N, 3) or tuple (3,)
+        R, G and B should be in the range [0, 1].
+    scales : float or ndarray (N, )
+        Cone size.
+    opacity : float
+        Takes values from 0 (fully transparent) to 1 (opaque).
 
     Returns
     -------
@@ -65,9 +90,6 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
     big_scales = np.repeat(scales, n_verts, axis=0)
     attribute_to_actor(box_actor, big_scales, 'scale')
 
-    big_values = np.repeat(np.array(lengths, dtype=float), n_verts, axis=0)
-    attribute_to_actor(box_actor, big_values, 'evals')
-
     evec1 = np.array([item[0] for item in axes])
     evec2 = np.array([item[1] for item in axes])
     evec3 = np.array([item[2] for item in axes])
@@ -88,7 +110,6 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         """
         in vec3 center;
         in float scale;
-        in vec3 evals;
         in vec3 evec1;
         in vec3 evec2;
         in vec3 evec3;
@@ -97,7 +118,6 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         out vec4 vertexMCVSOutput;
         out vec3 centerMCVSOutput;
         out float scaleVSOutput;
-        out vec3 evalsVSOutput;
         out mat3 rotationMatrix;
         out float angleVSOutput;
         """
@@ -116,9 +136,9 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         """
         mat3 R = mat3(normalize(evec1), normalize(evec2), normalize(evec3));
         float a = radians(90);
-        mat3 rot = mat3(cos(a),-sin(a),0,
-                        sin(a),cos(a), 0, 
-                        0,     0,      1);
+        mat3 rot = mat3(cos(a),-sin(a), 0,
+                        sin(a), cos(a), 0, 
+                            0 ,      0, 1);
         rotationMatrix = transpose(R) * rot;
         """
 
@@ -132,7 +152,6 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
         in vec4 vertexMCVSOutput;
         in vec3 centerMCVSOutput;
         in float scaleVSOutput;
-        in vec3 evalsVSOutput;
         in mat3 rotationMatrix;
         in float angleVSOutput;
 
@@ -225,34 +244,81 @@ def double_cone(centers, axes, lengths, angles, colors, scales, opacity):
 
 
 def main_dir_uncertainty(data, bvals, bvecs):
+    """
+    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.
+
+    Returns
+    -------
+    angles, centers, evecs
+
+    Notes
+    -----
+    The uncertainty calculation is based on first-order matrix perturbation
+    analysis described in [1]_. The idea is to estimate the variance of the
+    main eigenvector which corresponds to the main direction of diffusion,
+    directly from estimated D and its estimated covariance matrix \Delta D (see
+    [2]_, equation 4). The subtended angle, $\Delta\theta_i$, between the ith
+    perturbed eigenvector of D, $\varepsilon_i+\Delta\varepsilon_i$, and the
+    estimated eigenvector $\varepsilon_i$, measures the angular deviation of
+    the fiber direction, $\Delta\theta_i$:
+
+    .. math::
+        \theta=\tan^{-1}(\|\Delta\varepsilon_1\|)
+
+    Giving way to a graphical construct for displaying both the main
+    eigenvector of D and its associated uncertainty, with the so-called
+    "uncertainty cone".
+
+    References
+    ----------
+    .. [1] Basser, P. J. (1997). Quantifying errors in fiber direction and
+    diffusion tensor field maps resulting from MR noise. In 5th Scientific
+    Meeting of the ISMRM (Vol. 1740).
+
+    .. [2] Chang, L. C., Koay, C. G., Pierpaoli, C., & Basser, P. J. (2007).
+    Variance of estimated DTI‐derived parameters via first‐order perturbation
+    methods. Magnetic Resonance in Medicine: An Official Journal of the
+    International Society for Magnetic Resonance in Medicine, 57(1), 141-149.
+    """
+
     gtab = gradient_table(bvals, bvecs)
 
     tenmodel = dti.TensorModel(gtab)
     tenfit = tenmodel.fit(data)
-    evals = tenfit.evals
-    evecs = tenfit.evecs
+    fevals = tenfit.evals
+    fevecs = tenfit.evecs
 
     tensor_vals = dti.lower_triangular(tenfit.quadratic_form)
     dti_params = dti.eig_from_lo_tri(tensor_vals)
 
-    signalOI = dti.tensor_prediction(dti_params, gtab, 1.0)
+    fsignal = dti.tensor_prediction(dti_params, gtab, 1.0)
     b_matrix = dti.design_matrix(gtab)
 
-    valid_mask = np.abs(evecs).max(axis=(-2, -1)) > 0
+    valid_mask = np.abs(fevecs).max(axis=(-2, -1)) > 0
     indices = np.nonzero(valid_mask)
 
-    dofs_vecs = evecs[indices]
-    dofs_vals = evals[indices]
-    signal = signalOI[indices]
+    evecs = fevecs[indices]
+    evals = fevals[indices]
+    signal = fsignal[indices]
 
-    # Uncertainty calculations ------------------------------------------------
+    # Uncertainty calculations
 
     sigma = ne.estimate_sigma(data)  # standard deviation of the noise
 
-    # Angles for cone of uncertainty ------------------------------------------
-
-    angles = np.ones(dofs_vecs.shape[0])
-    for i in range(dofs_vecs.shape[0]):
+    # 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)
@@ -261,13 +327,14 @@ def main_dir_uncertainty(data, bvals, bvecs):
         delta_DD = dti.from_lower_triangular(
             np.array([dd[0], dd[3], dd[1], dd[4], dd[5], dd[2]]))
 
+        # perturbation matrix of tensor D
         try:
             delta_D = np.linalg.inv(delta_DD)
         except:
             delta_D = delta_DD
 
-        D_ = dofs_vecs
-        eigen_vals = dofs_vals[i]
+        D_ = evecs
+        eigen_vals = evals[i]
 
         e1, e2, e3 = np.array(D_[i, :, 0]), np.array(D_[i, :, 1]), \
                      np.array(D_[i, :, 2])
@@ -286,7 +353,7 @@ def main_dir_uncertainty(data, bvals, bvecs):
             theta = np.arctan(np.linalg.norm(delta_e1))
             angles[i] = theta
         else:
-            theta = 0.0872665
+            theta = 1.39626
 
     centers = np.asarray(indices).T
-    return angles, centers, dofs_vecs, dofs_vals
+    return angles, centers, evecs
diff --git a/fury/actors/tests/test_uncertainty.py b/fury/actors/tests/test_uncertainty.py
index 191c925f7..e99156d67 100644
--- a/fury/actors/tests/test_uncertainty.py
+++ b/fury/actors/tests/test_uncertainty.py
@@ -1,6 +1,7 @@
 """
-This spript includes the implementation of cone of uncertainty using matrix
-perturbation analysis
+This spript includes the implementation of dti_uncertainty actor for the
+visualization of the cones of uncertainty which uses matrix perturbation
+analysis for its calculation
 """
 from fury import actor, window
 
@@ -11,11 +12,16 @@
 
 
 def test_uncertainty():
-    hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
+    hardi_fname, hardi_bval_fname, hardi_bvec_fname =\
+        get_fnames('stanford_hardi')
 
     data, affine = load_nifti(hardi_fname)
+    data_ = data[20:24, 68:72, 28:29]
+    print(data_)
+    print(data_.shape)
 
     bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
+    print(bvals, bvecs)
 
     from dipy.segment.mask import median_otsu
 
@@ -23,7 +29,7 @@ def test_uncertainty():
                                  numpass=1, autocrop=True, dilate=2)
 
     uncertainty_cones = actor.dti_uncertainty(
-        data=maskdata[13:43, 44:74, 28:29], bvals=bvals, bvecs=bvecs)
+        data=maskdata[20:24, 68:72, 28:29], bvals=bvals, bvecs=bvecs)
 
     scene = window.Scene()
     scene.background([255, 255, 255])
diff --git a/fury/tests/test_actors.py b/fury/tests/test_actors.py
index e21593bb7..bc3c502d4 100644
--- a/fury/tests/test_actors.py
+++ b/fury/tests/test_actors.py
@@ -1734,6 +1734,24 @@ def test_marker_actor(interactive=False):
     npt.assert_equal(report.objects, 12)
 
 
+def test_uncertainty_actor(interactive=False):
+    scene = window.Scene()
+    scene.background((0, 0, 0))
+
+    data = np.array([0, 0, 0])
+    bvals = np.array([0, 0, 0])
+    bvecs = np.array([0, 0, 0])
+
+    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs,
+                                         scales=1.0, opacity=1.0)
+    scene.add(uncert_cones)
+
+    if interactive:
+        window.show(scene)
+
+    report = window.analyze_scene(scene)
+    npt.assert_equal(report.actors, 1)
+
 def test_actors_primitives_count():
     centers = np.array([[1, 1, 1], [2, 2, 2]])
     directions = np.array([[1, 0, 0], [1, 0, 0]])

From 291b4d2dc589cea03f971c8fd599fd8bd5414e07 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 10 Jul 2023 01:06:22 -0500
Subject: [PATCH 05/13] updated tests

---
 fury/actors/tests/test_uncertainty.py | 53 +++++++++++++++++++++------
 fury/tests/test_actors.py             | 31 +++++++++++++---
 2 files changed, 67 insertions(+), 17 deletions(-)

diff --git a/fury/actors/tests/test_uncertainty.py b/fury/actors/tests/test_uncertainty.py
index e99156d67..f4ebdd5d6 100644
--- a/fury/actors/tests/test_uncertainty.py
+++ b/fury/actors/tests/test_uncertainty.py
@@ -1,14 +1,19 @@
 """
 This spript includes the implementation of dti_uncertainty actor for the
-visualization of the cones of uncertainty which uses matrix perturbation
-analysis for its calculation
+visualization of the cones of uncertainty along with the diffusion tensors for
+comparison
 """
+from dipy.reconst import dti
+from dipy.segment.mask import median_otsu
+
 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
+from dipy.data import get_fnames, read_stanford_hardi
+
+from fury.primitive import prim_sphere
 
 
 def test_uncertainty():
@@ -16,12 +21,9 @@ def test_uncertainty():
         get_fnames('stanford_hardi')
 
     data, affine = load_nifti(hardi_fname)
-    data_ = data[20:24, 68:72, 28:29]
-    print(data_)
-    print(data_.shape)
 
+    # load the b-values and b-vectors
     bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
-    print(bvals, bvecs)
 
     from dipy.segment.mask import median_otsu
 
@@ -29,14 +31,43 @@ def test_uncertainty():
                                  numpass=1, autocrop=True, dilate=2)
 
     uncertainty_cones = actor.dti_uncertainty(
-        data=maskdata[20:24, 68:72, 28:29], bvals=bvals, bvecs=bvecs)
+        data=maskdata[13:43, 44:74, 28:29], bvals=bvals, bvecs=bvecs)
 
     scene = window.Scene()
     scene.background([255, 255, 255])
 
+    scene.add(diffusion_tensors())
+    window.show(scene, reset_camera=False)
+    scene.clear()
+
     scene.add(uncertainty_cones)
+    window.show(scene, reset_camera=False)
 
-    scene.reset_camera()
-    scene.reset_clipping_range()
 
-    window.show(scene, reset_camera=False)
+class Sphere:
+
+    vertices = None
+    faces = None
+
+def diffusion_tensors():
+    # https://dipy.org/documentation/1.0.0./examples_built/reconst_dti/
+    img, gtab = read_stanford_hardi()
+    data = img.get_data()
+
+    maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
+                                 numpass=1, autocrop=True, dilate=2)
+    tenmodel = dti.TensorModel(gtab)
+    tenfit = tenmodel.fit(maskdata)
+
+    evals = tenfit.evals[13:43, 44:74, 28:29]
+    evecs = tenfit.evecs[13:43, 44:74, 28:29]
+
+    vertices, faces = prim_sphere('symmetric724', True)
+    sphere = Sphere()
+    sphere.vertices = vertices
+    sphere.faces = faces
+
+    from dipy.data import get_sphere
+    sphere = get_sphere('symmetric724')
+
+    return actor.tensor_slicer(evals, evecs, sphere=sphere, scale=0.3)
diff --git a/fury/tests/test_actors.py b/fury/tests/test_actors.py
index bc3c502d4..751a1c32e 100644
--- a/fury/tests/test_actors.py
+++ b/fury/tests/test_actors.py
@@ -1736,14 +1736,31 @@ def test_marker_actor(interactive=False):
 
 def test_uncertainty_actor(interactive=False):
     scene = window.Scene()
-    scene.background((0, 0, 0))
 
-    data = np.array([0, 0, 0])
-    bvals = np.array([0, 0, 0])
-    bvecs = np.array([0, 0, 0])
+    from dipy.data import get_fnames
+    from dipy.io import read_bvals_bvecs
+    _, fbvals, fbvecs = get_fnames('small_101D')
+    bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
 
-    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs,
-                                         scales=1.0, opacity=1.0)
+    n = bvals.shape[0]
+    data = np.ones((10, 10, 1, n)) * 0.1
+
+    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs)
+    scene.add(uncert_cones)
+
+    if interactive:
+        window.show(scene)
+
+    report = window.analyze_scene(scene)
+    npt.assert_equal(report.actors, 1)
+    arr = window.snapshot(scene, offscreen=True)
+    report = window.analyze_snapshot(arr)
+    npt.assert_equal(report.objects, 100)
+    scene.clear()
+
+    n = bvals.shape[0]
+    data = np.ones((5, 5, 5, n)) * 0.5
+    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs)
     scene.add(uncert_cones)
 
     if interactive:
@@ -1751,6 +1768,8 @@ def test_uncertainty_actor(interactive=False):
 
     report = window.analyze_scene(scene)
     npt.assert_equal(report.actors, 1)
+    scene.clear()
+
 
 def test_actors_primitives_count():
     centers = np.array([[1, 1, 1], [2, 2, 2]])

From 2a69b4a96c458e6a55c18baa6ddc614e2d423c02 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 10 Jul 2023 13:13:43 -0500
Subject: [PATCH 06/13] line commented in test

---
 fury/actors/tests/test_uncertainty.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fury/actors/tests/test_uncertainty.py b/fury/actors/tests/test_uncertainty.py
index f4ebdd5d6..9e40c1672 100644
--- a/fury/actors/tests/test_uncertainty.py
+++ b/fury/actors/tests/test_uncertainty.py
@@ -67,7 +67,7 @@ def diffusion_tensors():
     sphere.vertices = vertices
     sphere.faces = faces
 
-    from dipy.data import get_sphere
-    sphere = get_sphere('symmetric724')
+    #from dipy.data import get_sphere
+    #sphere = get_sphere('symmetric724')
 
     return actor.tensor_slicer(evals, evecs, sphere=sphere, scale=0.3)

From f445b3907e7e5f98d2629e51a2696b6973c1ba57 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Fri, 14 Jul 2023 18:33:23 -0500
Subject: [PATCH 07/13] updated uncertainty calculation with different
 parameters

---
 fury/actor.py                         | 66 ++++++++++--------
 fury/actors/tensor.py                 | 99 +++++----------------------
 fury/actors/tests/test_uncertainty.py | 29 +++++++-
 3 files changed, 80 insertions(+), 114 deletions(-)

diff --git a/fury/actor.py b/fury/actor.py
index 39b4e95c3..e340cbab9 100644
--- a/fury/actor.py
+++ b/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,12 +3813,16 @@ 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
@@ -3824,24 +3830,26 @@ def dti_uncertainty(
 
     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)
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index ed5bdf93e..059c0deef 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -1,10 +1,5 @@
 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
@@ -12,46 +7,6 @@
                           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,31 +250,6 @@ 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)
@@ -324,8 +257,9 @@ def main_dir_uncertainty(data, bvals, bvecs):
         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
diff --git a/fury/actors/tests/test_uncertainty.py b/fury/actors/tests/test_uncertainty.py
index 9e40c1672..d0914c971 100644
--- a/fury/actors/tests/test_uncertainty.py
+++ b/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])

From 057c22d237ad4794440c90854931ee9f3bf4454b Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Mon, 17 Jul 2023 08:57:48 -0500
Subject: [PATCH 08/13] updated test

---
 fury/tests/test_actors.py | 50 ++++++++++++++++++++++++++++++---------
 1 file changed, 39 insertions(+), 11 deletions(-)

diff --git a/fury/tests/test_actors.py b/fury/tests/test_actors.py
index 751a1c32e..b1796de04 100644
--- a/fury/tests/test_actors.py
+++ b/fury/tests/test_actors.py
@@ -1734,18 +1734,36 @@ def test_marker_actor(interactive=False):
     npt.assert_equal(report.objects, 12)
 
 
-def test_uncertainty_actor(interactive=False):
+def test_uncertainty_cone_actor(interactive=False):
     scene = window.Scene()
 
-    from dipy.data import get_fnames
-    from dipy.io import read_bvals_bvecs
-    _, fbvals, fbvecs = get_fnames('small_101D')
-    bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
+    evals = np.array([1.4, 0.5, 0.35])
+    evecs = np.eye(3)
 
-    n = bvals.shape[0]
-    data = np.ones((10, 10, 1, n)) * 0.1
+    mevals = np.zeros((10, 10, 1, 3))
+    mevecs = np.zeros((10, 10, 1, 3, 3))
 
-    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs)
+    mevals[..., :] = evals
+    mevecs[..., :, :] = evecs
+
+    signal = np.ones((10, 10, 1, 10))
+    sigma = np.array([14.791911, 14.999622, 14.880976, 14.933881, 14.392784,
+                      14.132468, 14.334953, 14.409375, 14.514647, 14.409275])
+
+    b_matrix = np.array([[-1.8, -1.9, -4.8, -4.4, -2.3, -1.2, -1.0],
+                         [-5.4, -1.8, -1.6, -1.7, -6.1, -1.3, -1.0],
+                         [-6.2, -5.1, -1.0, -1.9, -9.3, -2.2, -1.0],
+                         [-2.8, -1.9, -4.8, -1.4, -2.1, -3.6, -1.0],
+                         [-5.6, -1.3, -7.8, -2.4, -5.2, -4.2, -1.0],
+                         [-1.8, -2.5, -1.8, -1.2, -2.3, -4.8, -1.0],
+                         [-2.3, -1.9, -6.8, -4.4, -6.4, -1.9, -1.0],
+                         [-1.8, -2.6, -4.8, -6.5, -7.7, -3.1, -1.0],
+                         [-6.2, -1.9, -5.6, -4.6, -1.5, -2.0, -1.0],
+                         [-2.4, -1.9, -4.5, -3.6, -2.5, -1.2, -1.0]])
+
+    uncert_cones = actor.uncertainty_cone(evecs=mevecs, evals=mevals,
+                                          signal=signal, sigma=sigma,
+                                          b_matrix=b_matrix)
     scene.add(uncert_cones)
 
     if interactive:
@@ -1758,9 +1776,19 @@ def test_uncertainty_actor(interactive=False):
     npt.assert_equal(report.objects, 100)
     scene.clear()
 
-    n = bvals.shape[0]
-    data = np.ones((5, 5, 5, n)) * 0.5
-    uncert_cones = actor.dti_uncertainty(data=data, bvals=bvals, bvecs=bvecs)
+    evals = np.array([1.4, 0.5, 0.35])
+    evecs = np.eye(3)
+
+    mevals = np.zeros((4, 4, 4, 3))
+    mevecs = np.zeros((4, 4, 4, 3, 3))
+
+    mevals[..., :] = evals
+    mevecs[..., :, :] = evecs
+
+    signal = np.ones((4, 4, 4, 10))
+    uncert_cones = actor.uncertainty_cone(evecs=mevecs, evals=mevals,
+                                          signal=signal, sigma=sigma,
+                                          b_matrix=b_matrix)
     scene.add(uncert_cones)
 
     if interactive:

From 942584c031af14f1c9a1ebb6089d47d7b11366cf Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Wed, 19 Jul 2023 16:33:09 -0500
Subject: [PATCH 09/13] moved uncertainty script example and fixed formula
 display

---
 .../examples/viz_uncertainty.py               | 63 +++++++++----------
 fury/actors/tensor.py                         | 12 ++--
 2 files changed, 36 insertions(+), 39 deletions(-)
 rename fury/actors/tests/test_uncertainty.py => docs/examples/viz_uncertainty.py (94%)

diff --git a/fury/actors/tests/test_uncertainty.py b/docs/examples/viz_uncertainty.py
similarity index 94%
rename from fury/actors/tests/test_uncertainty.py
rename to docs/examples/viz_uncertainty.py
index d0914c971..33ac15159 100644
--- a/fury/actors/tests/test_uncertainty.py
+++ b/docs/examples/viz_uncertainty.py
@@ -1,5 +1,5 @@
 """
-This spript includes the implementation of dti_uncertainty actor for the
+This script includes the implementation of dti_uncertainty actor for the
 visualization of the cones of uncertainty along with the diffusion tensors for
 comparison
 """
@@ -19,7 +19,35 @@
 from fury.primitive import prim_sphere
 
 
-def test_uncertainty():
+class Sphere:
+
+    vertices = None
+    faces = None
+
+def diffusion_tensors():
+    # https://dipy.org/documentation/1.0.0./examples_built/reconst_dti/
+    img, gtab = read_stanford_hardi()
+    data = img.get_fdata()
+
+    maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
+                                 numpass=1, autocrop=True, dilate=2)
+    tenmodel = dti.TensorModel(gtab)
+    tenfit = tenmodel.fit(maskdata)
+
+    evals = tenfit.evals[13:43, 44:74, 28:29]
+    evecs = tenfit.evecs[13:43, 44:74, 28:29]
+
+    vertices, faces = prim_sphere('symmetric724', True)
+    sphere = Sphere()
+    sphere.vertices = vertices
+    sphere.faces = faces
+
+    #from dipy.data import get_sphere
+    #sphere = get_sphere('symmetric724')
+
+    return actor.tensor_slicer(evals, evecs, sphere=sphere, scale=0.3)
+
+if __name__ == '__main__':
     hardi_fname, hardi_bval_fname, hardi_bvec_fname =\
         get_fnames('stanford_hardi')
 
@@ -28,8 +56,6 @@ def test_uncertainty():
     # load the b-values and b-vectors
     bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
 
-    from dipy.segment.mask import median_otsu
-
     maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
                                  numpass=1, autocrop=True, dilate=2)
 
@@ -67,32 +93,3 @@ def test_uncertainty():
 
     scene.add(uncertainty_cones)
     window.show(scene, reset_camera=False)
-
-
-class Sphere:
-
-    vertices = None
-    faces = None
-
-def diffusion_tensors():
-    # https://dipy.org/documentation/1.0.0./examples_built/reconst_dti/
-    img, gtab = read_stanford_hardi()
-    data = img.get_data()
-
-    maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
-                                 numpass=1, autocrop=True, dilate=2)
-    tenmodel = dti.TensorModel(gtab)
-    tenfit = tenmodel.fit(maskdata)
-
-    evals = tenfit.evals[13:43, 44:74, 28:29]
-    evecs = tenfit.evecs[13:43, 44:74, 28:29]
-
-    vertices, faces = prim_sphere('symmetric724', True)
-    sphere = Sphere()
-    sphere.vertices = vertices
-    sphere.faces = faces
-
-    #from dipy.data import get_sphere
-    #sphere = get_sphere('symmetric724')
-
-    return actor.tensor_slicer(evals, evecs, sphere=sphere, scale=0.3)
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index 1059066a7..e83e0222f 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -431,14 +431,14 @@ def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
     The uncertainty calculation is based on first-order matrix perturbation
     analysis described in [1]_. The idea is to estimate the variance of the
     main eigenvector which corresponds to the main direction of diffusion,
-    directly from estimated D and its estimated covariance matrix \Delta D (see
-    [2]_, equation 4). The subtended angle, $\Delta\theta_i$, between the ith
-    perturbed eigenvector of D, $\varepsilon_i+\Delta\varepsilon_i$, and the
-    estimated eigenvector $\varepsilon_i$, measures the angular deviation of
-    the fiber direction, $\Delta\theta_i$:
+    directly from estimated D and its estimated covariance matrix :math:`\Delta
+    D` (see [2]_, equation 4). The angle :math:`\Theta` between the perturbed
+    principal eigenvector of D, :math:`\epsilon_1+\Delta\epsilon_1`, and the
+    estimated eigenvector :math:`\epsilon_1`, measures the angular deviation of
+    the main fiber direction and can be approximated by:
 
     .. math::
-        \theta=\tan^{-1}(\|\Delta\varepsilon_1\|)
+        \Theta=tan^{-1}(\|\Delta\epsilon_1\|)
 
     Giving way to a graphical construct for displaying both the main
     eigenvector of D and its associated uncertainty, with the so-called

From dd6320014640717e6a8f367c5eb303f34cb3600b Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Sat, 29 Jul 2023 01:17:23 -0500
Subject: [PATCH 10/13] minor changes on format

---
 fury/actor.py         |  2 +-
 fury/actors/tensor.py | 13 +++++++------
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/fury/actor.py b/fury/actor.py
index a0f0cf9f2..10648c124 100644
--- a/fury/actor.py
+++ b/fury/actor.py
@@ -3896,7 +3896,7 @@ def uncertainty_cone(
     b_matrix : array (N, 7)
         Design matrix for DTI.
     scales : float or ndarray (N, ), optional
-        Cones of uncertainty size, default(1.0).
+        Cones of uncertainty size.
     opacity : float, optional
         Takes values from 0 (fully transparent) to 1 (opaque), default(1.0).
 
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index e83e0222f..8442a997c 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -6,6 +6,7 @@
 from fury.shaders import (attribute_to_actor, compose_shader,
                           import_fury_shader, shader_to_actor)
 
+
 def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
     """
     Visualize one or many Tensor Ellipsoids with different features.
@@ -136,7 +137,7 @@ def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
             */
             float scFactor = min(evalsVSOutput.x, min(evalsVSOutput.y,
                                      evalsVSOutput.z));
-                                     
+
             /*
             The approximation of distance is calculated by stretching the
             space such that the ellipsoid becomes a sphere (multiplying by
@@ -192,7 +193,7 @@ def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
             vec3 normal = centralDiffsNormals(pos, .0001);
             // Light Attenuation
             float la = dot(ld, normal);
-            vec3 color = blinnPhongIllumModel(la, lightColor0, 
+            vec3 color = blinnPhongIllumModel(la, lightColor0,
                 diffuseColor, specularPower, specularColor, ambientColor);
             fragOutput0 = vec4(color, opacity);
         }
@@ -298,7 +299,7 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
         mat3 R = mat3(normalize(evec1), normalize(evec2), normalize(evec3));
         float a = radians(90);
         mat3 rot = mat3(cos(a),-sin(a), 0,
-                        sin(a), cos(a), 0, 
+                        sin(a), cos(a), 0,
                             0 ,      0, 1);
         rotationMatrix = transpose(R) * rot;
         """
@@ -383,7 +384,7 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
             vec3 normal = centralDiffsNormals(pos, .0001);
             // Light Attenuation
             float la = dot(ld, normal);
-            vec3 color = blinnPhongIllumModel(la, lightColor0, 
+            vec3 color = blinnPhongIllumModel(la, lightColor0,
                 diffuseColor, specularPower, specularColor, ambientColor);
             fragOutput0 = vec4(color, opacity);
         }
@@ -470,7 +471,7 @@ def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
         # perturbation matrix of tensor D
         try:
             delta_D = np.linalg.inv(delta_DD)
-        except:
+        except np.linalg.LinAlgError:
             delta_D = delta_DD
 
         D_ = evecs
@@ -480,7 +481,7 @@ def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
                      np.array(D_[i, :, 2])
         lambda1, lambda2, lambda3 = eigen_vals[0], eigen_vals[1], eigen_vals[2]
 
-        if (lambda1 > lambda2 and lambda1 > lambda3):
+        if lambda1 > lambda2 and lambda1 > lambda3:
             # The perturbation of the eigenvector associated with the largest
             # eigenvalue is given by
             a = np.dot(np.outer(np.dot(e1, delta_D), np.transpose(e2)) /

From baa18daf227b5cf2a4be47bcc8c8c2540c872807 Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Sun, 6 Aug 2023 22:32:19 -0500
Subject: [PATCH 11/13] updated text

---
 fury/tests/test_actors.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/fury/tests/test_actors.py b/fury/tests/test_actors.py
index f73c8585b..c09862b98 100644
--- a/fury/tests/test_actors.py
+++ b/fury/tests/test_actors.py
@@ -1810,9 +1810,6 @@ def test_uncertainty_cone_actor(interactive=False):
 
     report = window.analyze_scene(scene)
     npt.assert_equal(report.actors, 1)
-    arr = window.snapshot(scene, offscreen=True)
-    report = window.analyze_snapshot(arr)
-    npt.assert_equal(report.objects, 100)
     scene.clear()
 
     evals = np.array([1.4, 0.5, 0.35])

From d555495dad542d76c4c290dcc95d3d20d3bf350b Mon Sep 17 00:00:00 2001
From: Tania Castillo <31288525+tvcastillod@users.noreply.github.com>
Date: Sat, 12 Aug 2023 22:49:12 -0500
Subject: [PATCH 12/13] added function for ray generation

---
 docs/examples/viz_uncertainty.py       | 95 --------------------------
 fury/actors/tensor.py                  | 64 ++++++++---------
 fury/shaders/ray_marching/gen_ray.frag | 17 +++++
 3 files changed, 44 insertions(+), 132 deletions(-)
 delete mode 100644 docs/examples/viz_uncertainty.py
 create mode 100644 fury/shaders/ray_marching/gen_ray.frag

diff --git a/docs/examples/viz_uncertainty.py b/docs/examples/viz_uncertainty.py
deleted file mode 100644
index 33ac15159..000000000
--- a/docs/examples/viz_uncertainty.py
+++ /dev/null
@@ -1,95 +0,0 @@
-"""
-This script includes the implementation of dti_uncertainty actor for the
-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
-from dipy.io.gradients import read_bvals_bvecs
-
-from dipy.data import get_fnames, read_stanford_hardi
-
-from fury.primitive import prim_sphere
-
-
-class Sphere:
-
-    vertices = None
-    faces = None
-
-def diffusion_tensors():
-    # https://dipy.org/documentation/1.0.0./examples_built/reconst_dti/
-    img, gtab = read_stanford_hardi()
-    data = img.get_fdata()
-
-    maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
-                                 numpass=1, autocrop=True, dilate=2)
-    tenmodel = dti.TensorModel(gtab)
-    tenfit = tenmodel.fit(maskdata)
-
-    evals = tenfit.evals[13:43, 44:74, 28:29]
-    evecs = tenfit.evecs[13:43, 44:74, 28:29]
-
-    vertices, faces = prim_sphere('symmetric724', True)
-    sphere = Sphere()
-    sphere.vertices = vertices
-    sphere.faces = faces
-
-    #from dipy.data import get_sphere
-    #sphere = get_sphere('symmetric724')
-
-    return actor.tensor_slicer(evals, evecs, sphere=sphere, scale=0.3)
-
-if __name__ == '__main__':
-    hardi_fname, hardi_bval_fname, hardi_bvec_fname =\
-        get_fnames('stanford_hardi')
-
-    data, affine = load_nifti(hardi_fname)
-
-    # load the b-values and b-vectors
-    bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
-
-    maskdata, mask = median_otsu(data, vol_idx=range(10, 50), median_radius=3,
-                                 numpass=1, autocrop=True, dilate=2)
-
-    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])
-
-    scene.add(diffusion_tensors())
-    window.show(scene, reset_camera=False)
-    scene.clear()
-
-    scene.add(uncertainty_cones)
-    window.show(scene, reset_camera=False)
diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index 8442a997c..e72492f4d 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -157,31 +157,26 @@ def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
     cast_ray = import_fury_shader(os.path.join(
         'ray_marching', 'cast_ray.frag'))
 
+    # Importing the function that generates the ray components
+    ray_generation = import_fury_shader(os.path.join(
+        'ray_marching', 'gen_ray.frag'))
+
     # Importing Blinn-Phong model for lighting
     blinn_phong_model = import_fury_shader(os.path.join(
         'lighting', 'blinn_phong_model.frag'))
 
     # Full fragment shader declaration
     fs_dec = compose_shader([fs_vars_dec, sd_sphere, sdf_map,
-                             central_diffs_normal, cast_ray,
+                             central_diffs_normal, cast_ray, ray_generation,
                              blinn_phong_model])
 
     shader_to_actor(box_actor, 'fragment', decl_code=fs_dec)
 
-    # Vertex in Model Coordinates.
-    point = "vec3 point = vertexMCVSOutput.xyz;"
-
-    # Camera position in world space
-    ray_origin = "vec3 ro = (-MCVCMatrix[3] * MCVCMatrix).xyz;"
-
-    ray_direction = "vec3 rd = normalize(point - ro);"
-
-    light_direction = "vec3 ld = normalize(ro - point);"
-
-    ray_origin_update = "ro += point - ro;"
-
-    # Total distance traversed along the ray
-    distance = "float t = castRay(ro, rd);"
+    ray_components = \
+        """
+        vec3 ro; vec3 rd; float t;
+        gen_ray(ro, rd, t);
+        """
 
     # Fragment shader output definition
     # If surface is detected, color is assigned, otherwise, nothing is painted
@@ -191,6 +186,8 @@ def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
         {
             vec3 pos = ro + t * rd;
             vec3 normal = centralDiffsNormals(pos, .0001);
+            // Light Direction
+            vec3 ld = normalize(ro - pos);
             // Light Attenuation
             float la = dot(ld, normal);
             vec3 color = blinnPhongIllumModel(la, lightColor0,
@@ -204,9 +201,7 @@ def tensor_ellipsoid(centers, axes, lengths, colors, scales, opacity):
         """
 
     # Full fragment shader implementation
-    sdf_frag_impl = compose_shader([point, ray_origin, ray_direction,
-                                    light_direction, ray_origin_update,
-                                    distance, frag_output_def])
+    sdf_frag_impl = compose_shader([ray_components, frag_output_def])
 
     shader_to_actor(box_actor, 'fragment', impl_code=sdf_frag_impl,
                     block='light')
@@ -348,31 +343,26 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
     cast_ray = import_fury_shader(os.path.join(
         'ray_marching', 'cast_ray.frag'))
 
+    # Importing the function that generates the ray components
+    ray_generation = import_fury_shader(os.path.join(
+        'ray_marching', 'gen_ray.frag'))
+
     # Importing Blinn-Phong model for lighting
     blinn_phong_model = import_fury_shader(os.path.join(
         'lighting', 'blinn_phong_model.frag'))
 
     # Full fragment shader declaration
     fs_dec = compose_shader([fs_vars_dec, sd_cone, sd_union, sdf_map,
-                             central_diffs_normal, cast_ray,
+                             central_diffs_normal, cast_ray, ray_generation,
                              blinn_phong_model])
 
     shader_to_actor(box_actor, 'fragment', decl_code=fs_dec)
 
-    # Vertex in Model Coordinates.
-    point = "vec3 point = vertexMCVSOutput.xyz;"
-
-    # Camera position in world space
-    ray_origin = "vec3 ro = (-MCVCMatrix[3] * MCVCMatrix).xyz;"
-
-    ray_direction = "vec3 rd = normalize(point - ro);"
-
-    light_direction = "vec3 ld = normalize(ro - point);"
-
-    ray_origin_update = "ro += point - ro;"
-
-    # Total distance traversed along the ray
-    distance = "float t = castRay(ro, rd);"
+    ray_components = \
+        """
+        vec3 ro; vec3 rd; float t;
+        gen_ray(ro, rd, t);
+        """
 
     # Fragment shader output definition
     # If surface is detected, color is assigned, otherwise, nothing is painted
@@ -382,6 +372,8 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
         {
             vec3 pos = ro + t * rd;
             vec3 normal = centralDiffsNormals(pos, .0001);
+            // Light Direction
+            vec3 ld = normalize(ro - pos);
             // Light Attenuation
             float la = dot(ld, normal);
             vec3 color = blinnPhongIllumModel(la, lightColor0,
@@ -395,9 +387,7 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
         """
 
     # Full fragment shader implementation
-    sdf_frag_impl = compose_shader([point, ray_origin, ray_direction,
-                                    light_direction, ray_origin_update,
-                                    distance, frag_output_def])
+    sdf_frag_impl = compose_shader([ray_components, frag_output_def])
 
     shader_to_actor(box_actor, 'fragment', impl_code=sdf_frag_impl,
                     block='light')
@@ -407,7 +397,7 @@ def double_cone(centers, axes, angles, colors, scales, opacity):
 
 def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
     """
-    Calculates the angle of the cone of uncertainty that represents the
+    Calculate the angle of the cone of uncertainty that represents the
     perturbation of the main eigenvector of the diffusion tensor matrix.
 
     Parameters
diff --git a/fury/shaders/ray_marching/gen_ray.frag b/fury/shaders/ray_marching/gen_ray.frag
new file mode 100644
index 000000000..1ddd92801
--- /dev/null
+++ b/fury/shaders/ray_marching/gen_ray.frag
@@ -0,0 +1,17 @@
+void gen_ray(out vec3 ro, out vec3 rd, out float t)
+{
+    // Vertex in Model Coordinates
+    vec3 point = vertexMCVSOutput.xyz;
+
+    // Ray Origin
+    // Camera position in world space
+    ro = (-MCVCMatrix[3] * MCVCMatrix).xyz;
+
+    // Ray Direction
+    rd = normalize(point - ro);
+
+    ro += point - ro;
+
+    // Total distance traversed along the ray
+    t = castRay(ro, rd);
+}

From 022a73ba209f5614d688bf841e7d2e3bbc9fe4f8 Mon Sep 17 00:00:00 2001
From: tvcastillod <tvcastillod@unal.edu.co>
Date: Tue, 31 Oct 2023 15:17:25 -0500
Subject: [PATCH 13/13] added comment on main_dir_uncertainty function

---
 fury/actors/tensor.py          | 3 +++
 fury/shaders/sdf/sd_cone.frag  | 2 +-
 fury/shaders/sdf/sd_union.frag | 2 +-
 3 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/fury/actors/tensor.py b/fury/actors/tensor.py
index e72492f4d..ab6aa3c3a 100644
--- a/fury/actors/tensor.py
+++ b/fury/actors/tensor.py
@@ -484,6 +484,9 @@ def main_dir_uncertainty(evals, evecs, signal, sigma, b_matrix):
             theta = np.arctan(np.linalg.norm(delta_e1))
             angles[i] = theta
         else:
+            # If the condition is not satisfied it means that there is not a
+            # predominant diffusion direction, hence the uncertainty will be
+            # higher and a default value close to pi/2 is assigned
             theta = 1.39626
 
     return angles
diff --git a/fury/shaders/sdf/sd_cone.frag b/fury/shaders/sdf/sd_cone.frag
index 6ff3fd40f..678e5107d 100644
--- a/fury/shaders/sdf/sd_cone.frag
+++ b/fury/shaders/sdf/sd_cone.frag
@@ -12,4 +12,4 @@ float sdCone( vec3 p, vec2 c, float h )
     float d = min(dot( a, a ),dot(b, b));
     float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y)  );
     return sqrt(d)*sign(s);
-}
\ No newline at end of file
+}
diff --git a/fury/shaders/sdf/sd_union.frag b/fury/shaders/sdf/sd_union.frag
index 958c0e403..6bcfe8903 100644
--- a/fury/shaders/sdf/sd_union.frag
+++ b/fury/shaders/sdf/sd_union.frag
@@ -1,4 +1,4 @@
 float opUnion( float d1, float d2 )
 {
     return min(d1,d2);
-}
\ No newline at end of file
+}