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

[Karthik Krishnan]

Hi

To re-iterate what was just said. I am talking about the difference between
lines 172-187 vs 189-214. in
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/Common/itkDiffusionTensor3D.txx?annotate=1.6&root=Insight

Thanks
Regards
Karthik

Karthik Krishnan wrote:

> 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
> />>/
> />/
> />/
> /
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
>