[Insight-developers] Can someone clarify fractionalAnisotropy
for me ?
Karthik Krishnan
Karthik.Krishnan at kitware.com
Sun Jul 10 01:30:41 EDT 2005
On Sat, 2005-07-09 at 19:39 -0400, Gordon Kindlmann wrote:
> 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.
>
Thank you very much for your reply.
I must have been sleeping when I wrote that email. I was taking the
absolute value of the eigen values at some point. The two equations are
equivalent as the paper you've mentioned shows.
> I am interested in the filter you mention for estimating tensors from
> DWIs.
I've committed the filter in
Code/BasicFilters/itkDiffusionTensor3DReconstructionImageFilter.h , txx
An example for the same is at
Examples/Filtering/DiffusionTensor3DReconstructionImageFilter.cxx
You need tensor data to run the example.
> 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.
>
Jim Miller suggested using template specialization on itk::Image to
create a class that handles images of this nature. I'll try to chalk one
up.
> Can you describe more about the filter you currently have and the sort
> of input image it uses?
It takes a reference image (acquired in the absence of diffusion
sensitizing gradients) and n gradient images and gradient directions and
computes an image of tensors. (with DiffusionTensor3D as the pixel
type). Once that is done, you can apply filters on this tensor image to
compute FA, ADC, RGB weighted maps etc..
The example uses the DiffusionTensor3DReconstructionImageFilter and
the TensorFractionalAnisotropyImageFilter to do generate the FA map.
Thanks
Regards
Karthik
>
> 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
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
More information about the Insight-developers
mailing list