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

[Torsten Rohlfing]

Simon:

Yes, that makes sense now. I don't think you even need to assume D has 
been diagonalized, but that's probably not worth pondering now. I 
suspect I also figured out where I went wrong in my original thinking - 
when I split the difference that is squared inside a double summation, I 
failed to account for the second nested summation level, so I was 
missing some terms.

I appreciate your help with this!

Best,
  Torsten

Simon Warfield wrote:

> 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
>>
>
>


-- 
Torsten Rohlfing, PhD          SRI International, Neuroscience Program
 Research Scientist             333 Ravenswood Ave, Menlo Park, CA 94025
  Phone: ++1 (650) 859-3379      Fax: ++1 (650) 859-2743
   torsten at synapse.sri.com        http://www.stanford.edu/~rohlfing/

     "Though this be madness, yet there is a method in't"