itkDiffusionTensor3D.txx

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkDiffusionTensor3D.txx,v $
  Language:  C++
  Date:      $Date: 2010-03-17 16:34:47 $
  Version:   $Revision: 1.9 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkDiffusionTensor3D_txx
#define __itkDiffusionTensor3D_txx

#include "itkDiffusionTensor3D.h"
#include "itkNumericTraits.h"

namespace itk
{


template<class T>
DiffusionTensor3D<T>
::DiffusionTensor3D()
{
}


template<class T>
DiffusionTensor3D<T>
::DiffusionTensor3D( const Superclass & r ):SymmetricSecondRankTensor<T,3>(r)
{
}


template<class T>
DiffusionTensor3D<T>
::DiffusionTensor3D( const ComponentType & r ):SymmetricSecondRankTensor<T,3>(r)
{
}


template<class T>
DiffusionTensor3D<T>
::DiffusionTensor3D( const ComponentArrayType r ):SymmetricSecondRankTensor<T,3>(r)
{
}


template<class T>
DiffusionTensor3D<T>&
DiffusionTensor3D<T>
::operator= (const ComponentType & r)
{
  Superclass::operator=(r);
  return *this;
}


template<class T>
DiffusionTensor3D<T>&
DiffusionTensor3D<T>
::operator= (const ComponentArrayType r)
{
  Superclass::operator=(r);
  return *this;
}


template<class T>
DiffusionTensor3D<T>&
DiffusionTensor3D<T>
::operator= (const Superclass & r)
{
  Superclass::operator=(r);
  return *this;
}


template<class T>
typename DiffusionTensor3D<T>::AccumulateValueType
DiffusionTensor3D<T>
::GetTrace() const
{
  AccumulateValueType trace = (*this)[0];
  trace += (*this)[3];
  trace += (*this)[5];
  return trace;
}


template<class T>
typename DiffusionTensor3D<T>::RealValueType
DiffusionTensor3D<T>
::GetFractionalAnisotropy() const
{
  // Computed as 
  // FA = vcl_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
  const RealValueType isp   = this->GetInnerScalarProduct();

  if( isp > 0.0 )
    {
    const RealValueType trace = this->GetTrace();
    const RealValueType anisotropy = 3.0 * isp - trace * trace;
    const RealValueType fractionalAnisotropy =
        static_cast< RealValueType >( vcl_sqrt(anisotropy / ( 2.0 * isp ) ) );
    return fractionalAnisotropy;
    }

  return 0.0;
}


template<class T>
typename DiffusionTensor3D<T>::RealValueType
DiffusionTensor3D<T>
::GetRelativeAnisotropy() const
{
  const RealValueType trace = this->GetTrace();
  const RealValueType isp   = this->GetInnerScalarProduct();

  // Avoid negative trace and traces small enough to look like a division by
  // zero.
  if( trace < NumericTraits< RealValueType >::min() )
    {
    return NumericTraits<RealValueType>::Zero;
    }

  const RealValueType anisotropy = 3.0 * isp - trace * trace;

  if( anisotropy  < NumericTraits< RealValueType >::Zero )
    {
    return NumericTraits<RealValueType>::Zero;
    }

  const RealValueType relativeAnisotropySquared =
        static_cast< RealValueType >( anisotropy / ( vcl_sqrt(3.0 ) * trace ) ); 

  const RealValueType relativeAnisotropy =
        static_cast< RealValueType >( vcl_sqrt(relativeAnisotropySquared ) );

  return relativeAnisotropy;
}


template<class T>
typename DiffusionTensor3D<T>::RealValueType
DiffusionTensor3D<T>
::GetInnerScalarProduct() const
{

  const RealValueType xx = (*this)[0];
  const RealValueType xy = (*this)[1];
  const RealValueType xz = (*this)[2];
  const RealValueType yy = (*this)[3];
  const RealValueType yz = (*this)[4];
  const RealValueType zz = (*this)[5];

  return ( xx*xx + yy*yy + zz*zz + 2.0*(xy*xy + xz*xz + yz*yz) );

}


} // end namespace itk

#endif