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

[Gordon Kindlmann]

hello,

For reference, the paper which originally defined FA is:
http://dir2.nichd.nih.gov/nichd/stbb/Microstruc.pdf
(equation 15, page 213)
This is the same as equation 28 of the Westin paper, and it suggests  
that you should go with the version which does not rely on eigenvalue  
computation.  It is not uncommon for estimated tensors to have negative  
eigenvalues.  And the sign of eigenvalues is *not* merely convention,  
though the sign of eigenvectors is.

I am interested in the filter you mention for estimating tensors from  
DWIs.  During the recent NA-MIC Programming Week there was interest in  
having an ITK image type which could handle a variable number of DWI  
values (number of pixel components known only at run-time), while also  
having all DWI values contiguous in memory, and with the DWIs for each  
pixel along the "fastest" axis in memory (as opposed to being volume  
interleaved).  My understanding from the impromptu discussion with Luis  
et al was that this was slightly beyond the current representational  
abilities in ITK.

Can you describe more about the filter you currently have and the sort  
of input image it uses?

Gordon

On Jul 9, 2005, at 6:37 PM, 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