[RNH post #13] Start wrapping up - Test singularities of kurtosis statistics.

As we are reaching the end of the GSoC coding period, I am starting to wrap up the code that I developed this summer.

When reviewing the code implementing the kurtosis standard statistics, I detected some problems on the performance of the analytical solution of the mean kurtosis function.

In this post, I am reporting how I overcome these issues! This post is extremely relevant for whom is interested in knowing the full details of the implementation of kurtosis standard measures.

Problematic performance near to function singularities

As I mention in previous posts, the function to compute the mean kurtosis was implemented according to an analytical solution proposed by Tabesh et al., 2011. The mathematical formulas of this analytical solution have some singularities, particularly for cases that the diffusion tensor has equal eigenvalues.

To illustrate theses singularities, I am plotting below the diffusion and kurtosis tensors of crossing fiber simulates with different intersection angles. Simulates were preformed based on the modules implemented during the GSoC coding period.

Figure 1 - Diffusion tensor (upper panels) and kurtosis tensors (lower panels) for crossing fibers intersecting at different angles (the ground truth fiber direction are shown in red).

The values of the eigenvalues of the diffusion tensor as function of the intersection angle are shown below.

Figure 2 - Diffusion eigenvalues as function of crossing fibers intersection angle. First eigenvalues is plotted in red while the second and third are plotted in green and blue.

From the figure above, we can detect two problematic cases for the MK analytical solution:

1. When intersection angle is zero (i.e. when the two fibers are aligned), the second diffusion eigenvalue is equal to the third eigenvalue.
2. When the intersection angle is 90 degrees, the first diffusion eigenvalue is equal to the second eigenvalue.

Based on the work done by Tabesh et al., 2011, these MK estimation singularities can be mathematically revolved by detecting the problematic cases and using specific formulas for each detected situation. In the previous version of my codes, I was detecting the cases that two or three eigenvalues were equal by analysis if their differences were three orders of magnitude larger than system's epslon. For example, to automatically check if the first eigenvalue equals the second eigenvalue, the following lines of code were used:

import numpy as np
er = np.finfo(L1.ravel()[0]).eps * 1e3
cond1 = (abs(L1 - L2) < er)

Although, my testing modules were showing that this procedure was successfully solving the singularities for eigenvalue differences three orders of magnitude smaller than the system's epslon, when plotting MK as function of the intersection angle, some unexpected underestimated were present on the regions near to the singularities (see the figures below).

Figure 3 - MK values as function of the crossing angle. The blue line shows the MK values estimated from the analytical solution while the red line show the MK values estimated from a numerical method described in previous posts.

Figure 4 - MK values as function of the crossing angle (range between 85 and 90 degrees). The blue line shows the MK values estimated from the analytical solution while the red line show the MK values estimated from a numerical method described in previous posts. This figure was produced for a better visualization of the underestimations still present near to the crossing angle of 90 degrees.

After some analysis, I noticed that MK underestimations were still present if eigenvalues were not 2% different to each other. Given this, I was able to solve this underestimation by adjusting the criteria of eigenvalue comparison. As example, to compare the first eigenvalue with the second, the following lines of code are now used:

er = 2.5e-2  # difference (in %) between two eigenvalues to be considered as different
cond1 = (abs((L1 - L2) / L1) > er)

Below, I am showing the new MK estimates as function of the crossing angle, where all underestimations seem to be corrected. Moreover, discontinuities on the limits between the problematic and the non-problematic eigenvalue regime are relatively small. The most significant differences are now between different MK estimation methods (for details on the difference between these methods revisit post #9).

Figure 5 - Corrected MK values as function of the crossing angle. The blue line shows the MK values estimated from the analytical solution while the red line show the MK values estimated from the numerical method described in previous posts.