[Insight-developers] Can someone explain DiffusionTensor3D::GetRelativeAnisotropy() to me?

[Simon Warfield]

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
>


-- 
Simon K. Warfield, Ph.D. warfield at bwh.harvard.edu Phone:617-732-7090
http://www.spl.harvard.edu/~warfield           FAX:  617-582-6033
Associate Professor of Radiology,          Harvard Medical School
Director, Computational Radiology Laboratory
Thorn 329, Dept Radiology,  Brigham and Women's Hospital 
75 Francis St, Boston, MA, 02115
MA 280, Dept Radiology, Children's Hospital Phone: 617-355-4566