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]])