added uncertainty calculation

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)

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