[Insight-developers] Can someone clarify fractionalAnisotropy for me ?

[Karthik Krishnan]

Hi,

Coming back to a thread between Torsten and Simon Warfield a while ago. 

I am having some difficulty with the assumptions where the two equations are equivalent.

The definition of 
FA = sqrt(1.5 * ( \sum_i ( lambda_i - lambda_mean )^2 ) / \sum_i ( lambda_i^2 ) )
as in http://splweb.bwh.harvard.edu:8000/pages/papers/martha/DTI_Tech354.pdf
[lambda = eig(A)].

will be equivalent to  
FA = sqrt(1.5*sum(sum(N.*N))/sum((sum(D.*D))))
where N = D - ((1/3)*trace(D)*eye(3,3))
equation (28) in http://lmi.bwh.harvard.edu/papers/pdfs/2002/westinMEDIA02.pdf

only if the eigenvalues of the tensor D are positive.

I have a filter that reconstructs the tensor image from Diffusion weighted images and the FA computations on it
are done using the second equation simply because its much faster
(avoids the computational expense of the eigen value computation.)

The reconstructed tensor coefficients do not necessarily have positive eigen values. I do understand that the sign of 
the eigen values is merely a convention and can be reversed by reversing the sign of the eigen vectors. Nevertheless, 
the equivalence of the two equations requires it.

The results are then inconsistent. 

Am I missing something here ? I would appreciate any help.

Thanks
Regards
Karthik




>/ Torsten Rohlfing wrote:
/>/
/>>/
/>>/ Simon Warfield wrote:
/>>/
/>>>>/ The code (purged of numerical safeguards etc) that I don't get is 
/>>>>/ this:
/>>>>/
/>>>>/  const RealValueType anisotropy = 3.0 * isp - trace * trace;
/>>>>/  const RealValueType relativeAnisotropySquared =
/>>>>/        static_cast< RealValueType >( anisotropy / ( sqrt( 3.0 ) * 
/>>>>/ trace ) );
/>>>>/  const RealValueType relativeAnisotropy =
/>>>>/        static_cast< RealValueType >( sqrt( 
/>>>>/ relativeAnisotropySquared ) );
/>>>>/
/>>>>/ where "isp" is the inner scalar product and "trace" is the tensor 
/>>>>/ trace, both as returned by the respective class member functions.
/>>>>/
/>>>>/ My problem is with the definition of "anisotropy" above. I assume 
/>>>>/ that this is supposed to be the inner product of the anisotropic 
/>>>>/ part of the tensor, that is, the inner product of  (tensor - (trace 
/>>>>/ / 3) * identity). This is a double sum over the squares of 
/>>>>/ differences between elements of the original tensor and the 
/>>>>/ corresponding elements of (trace/3) times the identity tensor.
/>>>/
/>>>/
/>>>/
/>>>/
/>>>/ Check out Equation 27 here:
/>>>/ http://splweb.bwh.harvard.edu:8000/pages/papers/westin/media2002/Westin_Diffusion_MedAI_2002.pdf 
/>>>/
/>>>/ and the definition of the norm at the top of the column.
/>>>/
/>>/ Simon:
/>>/
/>>/ Thank you, but that's exactly my problem - how do you get from Eq. 27 
/>>/ in the Westin paper to the implementation in ITK? The Westin equation 
/>>/ is basically the same (except for the sqrt(2) factor) to Eq. 14 in 
/>>/ the 1996 JMR(B) paper by Basser and Pierpaoli, and their Eq. 10 is 
/>>/ the tensor norm as described in the Westin paper. But how do you get 
/>>/ from the numerator in the second line of Westin's Eq. 27 to 
/>>/ "3*isp-trace*trace" (or isp+1/3*trace*trace, which is basically the 
/>>/ same)?
/>/
/>/
/>/ I think we are allowed to consider the tensor after reorientation to 
/>/ the principal axes such that Dij = 0 for i != j.  Then trace(D) = D11 
/>/ + D22 + D33.
/>/ The Frobenius norm of
/>/ D - 1/3trace(D)I
/>/ is then the square root of
/>/ \sum_ij (D_ij - k)^2
/>/  for k = 1/3 trace(D) on the diagonal and zero elsewhere
/>/
/>/ Expanding this, we have
/>/ (D_11 - k)^2 + (D_22 - k)^2 + (D_33 - k)^2 + \sum_i != j D_ij^2
/>/  and this last term is 0.
/>/
/>/ Expanding the first three terms, I get:
/>/ D_11^2 + D_22^2 + D_33^2 + 3k^2 - 2k(D11 + D22 + D33)
/>/ = isp + 3k^2 - 2k.3k
/>/   for isp = D_11^2 + D_22^3 + D_33^2
/>/ = isp - 3k^2
/>/ = isp - 3.1/3 trace. 1/3 trace
/>/ = isp - 1/3 trace * trace
/>/
/>/ This is 1/3rd of what ITK is calling 'anisotropy'.  I didn't actually 
/>/ check the code for what isp is being computed as.
/>/ Does it make sense ?
/>/
/>/
/>>/
/>>/ I assume that the sum of the squared elements is broken up, that is, 
/>>/ the differences that make up each element are separated as (a-b)^2 -> 
/>>/ (a^2-2ab+b^2), and so sum_i sum_j (a-b)^2 becomes sum_ij a^2 - 
/>>/ 2sum_ij ab + sum_ij b^2.
/>>/
/>>/ Then with a_ij==D and b_ij==1/3 trace(D)I, that would give
/>>/
/>>/ sum_ij D_ij^2 - 2/3 trace(D) sum_ij I_ij + (1/3 trace(D))^2 sum_ij I_ij
/>>/
/>>/ and with the first term the norm of D and I the identity tensor we 
/>>/ would have
/>>/
/>>/ |D| - 2 trace(D) + 1/3 trace(D)^2
/>>/
/>>/ So where is the middle part going in ITK?
/>>/
/>>/ Sorry if I am suffering from a mental block here.
/>>/
/>>/ Again, thanks for your help!
/>>/  Torsten
/>>/
/>/
/>/
/