[Insight-developers] Tensors - my thoughts so far

Wed Feb 2 16:14:50 EST 2005

So here's what I came up with. Note that what I am sending here does not 
claim to be a complete tensor implementation with bells and whistles. I 
am just trying to illustrate the fundamental desgin structure of a 
possible approach to the problem. Attached are the following files:

itkSymmetricTensor.{h,txx}
 - Class for symmetric tensor. This class should incorporate all 
functions that do _not_ benefit from knowledge of the tensor size, e.g., 
GetVnlMatrix()

itkSymmetricTensorTraits.{h,txx}
- Symmetric tensor traits. These specialized classes should implement 
all operations that _do_ depend on the tensor size, e.g., eigenvalue 
computation. Note that I copied and pasted 
vnl_symmetric_eigensystem_compute_eigenvals 
<http://paine.wiau.man.ac.uk/pub/doc_vxl/core/vnl/html/vnl__symmetric__eigensystem_8cxx.html#a0>, 
which does not seem to be in the VXL version shipped with ITK. 
Ultimately, the vnl function should simply be called here. Likewise, the 
eigenvalue computation for tensor sizes other than 3 need to be somehow 
implemented.

itkTensorToFractionalAnisotropyImageFilter.{h,txx}
- This is a simple example of an image filter that operates on tensor 
data. Note that if these was a different, more general tensor class 
"itkTensor" (ie., for non-symmetric tensors), this filter could also 
operate on images of that pixel type, as long as it implements 
itkTensor::ComputeEigenvalues with the same interface as itkSymmetricTensor.

DiffusionTensorToFractionalAnisotropy.cxx
- Example application that creates a pipeline to read a tensor image 
from a raw data file, computes the FA image, and writes the result. Note 
that the image file structure is hard coded into the source code by 
parameterization of itkRawImageIO. Sorry, but that's my kind of data ;))

So in summary, the points I am trying to make here are:

a) we can separate the operations which do depend on the tensor size (or 
benefit from knowing the size) from those that do not. The former go 
into tensor traits, the latter go into the basic tensor class.

b) there can be filter classes, which, by templating, operate on 
different tensor pixel type images. The tensor classes can be more or 
less specialized. In particular, the filters will be able to support 
more general tensor classes that may be implemented later, as long as 
they provide the necessary interface functions for any given filter.

Best,
  Torsten
<http://paine.wiau.man.ac.uk/pub/doc_vxl/core/vnl/html/vnl__symmetric__eigensystem_8cxx.html#a0> 

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

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

-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile$
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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.

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

#include "itkImage.h"
#include "itkSymmetricTensor.h"
#include "itkImageFileReader.h"
#include "itkRawImageIO.h"
#include "itkImageFileWriter.h"
#include "itkExceptionObject.h"
#include "itkCastImageFilter.h"
#include "itkTensorToFractionalAnisotropyImageFilter.h"

#include <iostream>

int
main( int argc, char * argv[] )
{
  typedef itk::SymmetricTensor<float,3> TensorType;
  typedef itk::Image<TensorType,3> TensorImageType;
  typedef itk::Image<float,3> OutputImageType;

  typedef itk::Vector<float,6> InputPixelType;
  typedef itk::Image<InputPixelType,3> InputImageType;

  typedef itk::RawImageIO<InputPixelType,3> TensorIOType;
  TensorIOType::Pointer io = TensorIOType::New();

  io->SetFileDimensionality( 3 );
  io->SetNumberOfComponents( 6 );
  io->SetPixelType( itk::ImageIOBase::VECTOR );
  io->SetComponentType( itk::ImageIOBase::FLOAT );

  io->SetHeaderSize( 0 );
  io->SetNumberOfDimensions( 3 );
  io->SetDimensions( 0, 256 );
  io->SetDimensions( 1, 256 );
  io->SetDimensions( 2, 62 );

  io->SetSpacing( 0, 0.9375 );
  io->SetSpacing( 1, 0.9375 );
  io->SetSpacing( 2, 2.5000 );

  io->SetByteOrderToLittleEndian();

  /// Read raw file into 6-d itkVector data.
  typedef itk::ImageFileReader<InputImageType> ImageReaderType;
  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName( argv[1] );
  reader->SetImageIO( io );

  /// Cast itkVector to itkSymmetricTensor
  typedef itk::CastImageFilter<InputImageType,TensorImageType> CasterType;
  CasterType::Pointer caster = CasterType::New();
  caster->SetInput( reader->GetOutput() );

  /// Filter tensor data to FA image.
  typedef itk::TensorToFractionalAnisotropyImageFilter<TensorImageType,float> FilterType;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput( caster->GetOutput() );

  typedef itk::ImageFileWriter<OutputImageType> ImageWriterType;
  ImageWriterType::Pointer writer = ImageWriterType::New();
  writer->SetFileName( argv[2] );
  writer->SetInput( filter->GetOutput() );

  try
    {
    writer->Update();
    }
  catch ( itk::ExceptionObject e )
    {
    e.Print( std::cerr );
    }

  return 0;
}

-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSymmetricTensor.txx,v $
  Language:  C++
  Date:      $Date: 2004/04/15 22:37:33 $
  Version:   $Revision: 1.10 $

  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 __itkSymmetricTensor_h
#define __itkSymmetricTensor_h

#include "itkFixedArray.h"
#include "itkSymmetricTensorTraits.h"
#include "vnl/vnl_matrix.h"

namespace itk
{
/** \class itkSymmetricTensor
 *
 * \brief
 * 
 * \par Overview
 *
 * \par Template Parameters
 *
 * \par Filter Parameters 
 *
 * \par Constraints
 *
 * \par Performance
 *
 * \par References
 *
 * \ingroup DataRepresentation
 *
 */
template
 < typename TRealType = float, unsigned int VTensorDimension = 3,
  typename TTensorTraits = SymmetricTensorTraits<TRealType,VTensorDimension>
>
class ITK_EXPORT SymmetricTensor :
    public FixedArray< TRealType, VTensorDimension*(VTensorDimension+1)/2 >
{
 public:
  /// The compact linear vector dimension of the tensor.
  static const unsigned int NVectorDimension = 
    VTensorDimension*(VTensorDimension+1)/2;

  /** Standard class typedefs. */
  typedef SymmetricTensor Self;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self>  ConstPointer;
  typedef FixedArray< TRealType, NVectorDimension > Superclass;

  /** Run-time type information (and related methods) */
  itkTypeMacro(SymmetricTensor, FixedArray);

  /** Define the data type and the vector of data type used in calculations. 
  */
  typedef TRealType RealType;  

  typedef vnl_matrix_fixed<TRealType,VTensorDimension,VTensorDimension> VnlMatrixType;
  /* Return vnl_matrix object. We do this here, since it's always the same,
   * regardless of tensor size.
   */
  VnlMatrixType GetVnlMatrix() const;

  /** Compute eigenvalues. This is done in the traits class, since the best
   * implementation depends on the actual tensor size.
   */
  void ComputeEigenvalues( double (&lambda)[VTensorDimension] ) const 
    {
    TTensorTraits::ComputeEigenvalues( *this, lambda );
    }

  SymmetricTensor() : Superclass() {}
  SymmetricTensor(const Self& other) : Superclass( other ) {}
  void operator=(const Self& other) { this->Superclass::operator=( other ); }

  /** This is part of a hack that allows reading a raw vector field using
    * itkRawImageIO and subsequently casting to itkSymmetricTensor
    * using itkCastImageFilter.
    */    
  SymmetricTensor(const Vector<TRealType,NVectorDimension>& other) :
    Superclass( other ) {}

  /** This is part of a hack that allows reading a raw vector field using
    * itkRawImageIO and subsequently casting to itkSymmetricTensor
    * using itkCastImageFilter.
    */    
  typedef TRealType ComponentType;

  /** This is part of a hack that allows reading a raw vector field using
    * itkRawImageIO and subsequently casting to itkSymmetricTensor
    * using itkCastImageFilter.
    */    
  static int GetNumberOfComponents()
    { return VTensorDimension*(VTensorDimension+1)/2;}

  /** This is part of a hack that allows reading a raw vector field using
    * itkRawImageIO and subsequently casting to itkSymmetricTensor
    * using itkCastImageFilter.
    */    
  void SetNthComponent(int c, const ComponentType& v)  
  {  this->operator[](c) = v; }

protected:
  void PrintSelf(std::ostream& os, Indent indent) const;
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkSymmetricTensor.txx"
#endif

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSymmetricTensor.txx,v $
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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.

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

namespace itk
{

template
 < typename TRealType, unsigned int VTensorDimension, typename TTensorTraits >
vnl_matrix_fixed<TRealType,VTensorDimension,VTensorDimension>
SymmetricTensor<TRealType,VTensorDimension,TTensorTraits>
::GetVnlMatrix() const
{
  vnl_matrix_fixed<TRealType,VTensorDimension,VTensorDimension> matrix;

  unsigned int idx = 0;
  for ( unsigned int i = 0; i < VTensorDimension; ++i )
    {
    for ( unsigned int j = 0; j < VTensorDimension; ++j )
      {
      if ( i <= j )
	{
	matrix[i][j] = (*this)[idx++];
	}
      else
	{
	matrix[i][j] = matrix[j][i];
	}
      }
    }

  return matrix;
}

} // namespace itk
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSymmetricTensorTraits.h,v $
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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 __itkSymmetricTensorTraits_h
#define __itkSymmetricTensorTraits_h

#include "itkFixedArray.h"

namespace itk
{
/** \class itkSymmetricTensorTraits
 *
 * \brief
 * 
 * \par Overview
 *
 * \par Template Parameters (Input and Output)
 *
 * \par Filter Parameters 
 *
 * \par Constraints
 *
 * \par Performance
 *
 * \par References
 *
 * \ingroup DataRepresentation
 *
 */
template<typename TRealType, unsigned int VTensorDimension>
class ITK_EXPORT SymmetricTensorTraits
{
public:
  static const unsigned int NTensorDimension = VTensorDimension;
  static const unsigned int NVectorDimension = 
  VTensorDimension * (VTensorDimension+1) / 2;

  /** Standard class typedefs. */
  typedef SymmetricTensorTraits Self;
  typedef FixedArray<TRealType,NVectorDimension> SymmetricTensorBaseType;

  /// Perform eigenvalue computation.
  static void ComputeEigenvalues( const SymmetricTensorBaseType& tensor,
				  double (&lambda)[VTensorDimension] );
};

/// Specialization for 3D tensors.
template<typename TRealType>
class ITK_EXPORT SymmetricTensorTraits<TRealType,3>
{
public:
  static const unsigned int NTensorDimension = 3;
  static const unsigned int NVectorDimension = 6;

  /** Standard class typedefs. */
  typedef SymmetricTensorTraits Self;
  typedef FixedArray<TRealType,NVectorDimension> SymmetricTensorBaseType;

  /// Perform eigenvalue computation.
  static void ComputeEigenvalues( const SymmetricTensorBaseType& tensor,
				  double (&lambda)[3] );
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkSymmetricTensorTraits.txx"
#endif

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkSymmetricTensorTraits.txx,v $
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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.

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

namespace itk
{

template < typename TRealType, unsigned int VTensorDimension >
void
SymmetricTensorTraits<TRealType,VTensorDimension>
::ComputeEigenvalues( const SymmetricTensorBaseType& tensor,
		      double (&lambda)[VTensorDimension] )
{
}

template <typename TRealType>
void
SymmetricTensorTraits<TRealType,3>
::ComputeEigenvalues( const SymmetricTensorBaseType& tensor,
		      double (&lambda)[3] )
{
  const double M11 = tensor[0];
  const double M12 = tensor[1];
  const double M13 = tensor[2];
  const double M22 = tensor[3];
  const double M23 = tensor[4];
  const double M33 = tensor[5];

// <begin copyright notice>
// ---------------------------------------------------------------------------
//
//                Copyright (c) 2000-2003 TargetJr Consortium
//               GE Corporate Research and Development (GE CRD)
//                             1 Research Circle
//                            Niskayuna, NY 12309
//                            All Rights Reserved
//              Reproduction rights limited as described below.
//
//      Permission to use, copy, modify, distribute, and sell this software
//      and its documentation for any purpose is hereby granted without fee,
//      provided that (i) the above copyright notice and this permission
//      notice appear in all copies of the software and related documentation,
//      (ii) the name TargetJr Consortium (represented by GE CRD), may not be
//      used in any advertising or publicity relating to the software without
//      the specific, prior written permission of GE CRD, and (iii) any
//      modifications are clearly marked and summarized in a change history
//      log.
//
//      THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
//      EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
//      WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
//      IN NO EVENT SHALL THE TARGETJR CONSORTIUM BE LIABLE FOR ANY SPECIAL,
//      INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND OR ANY
//      DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
//      WHETHER OR NOT ADVISED OF THE POSSIBILITY OF SUCH DAMAGES, OR ON
//      ANY THEORY OF LIABILITY ARISING OUT OF OR IN CONNECTION WITH THE
//      USE OR PERFORMANCE OF THIS SOFTWARE.
//
// ---------------------------------------------------------------------------
// <end copyright notice>

  // Characteristic eqtn |M - xI| = 0
  // x^3 + b x^2 + c x + d = 0
  const double b = -M11-M22-M33;
  const double c =  M11*M22 +M11*M33 +M22*M33  -M12*M12 -M13*M13 -M23*M23;
  const double d = M11*M23*M23 +M12*M12*M33 +M13*M13*M22 -2.0*M12*M13*M23 -M11*M22*M33;
  // Using a numerically tweaked version of the real cubic solver http://www.1728.com/cubic2.htm
  const double b_3 = b/3.0;
  const double f = b_3*b_3 -  c/3.0 ;
  const double g = b*c/6.0 - b_3*b_3*b_3 - 0.5*d;

  if (f == 0.0 && g == 0.0)
    {
    lambda[0] = lambda[1] = lambda[2] = - b_3 ;
    return;
    }

  const double f3 = f*f*f;
  const double g2 = g*g;
  const double sqrt_f = -vcl_sqrt(f);

   // deal explicitly with repeated root and treat
   // complex conjugate roots as numerically inaccurate repeated roots.

   // first check we are not too numerically innacurate
//  assert((g2 - f3) / vnl_math_sqr(b*b*b) < 1e-8);  

  if (g2 >= f3)
    {
    if (g < 0.0)
      {
      lambda[0] = 2.0 * sqrt_f  - b_3;
      lambda[1] = lambda[2] = - sqrt_f - b_3;
      }
    else
      {
      lambda[0] = lambda[1] = sqrt_f  - b_3;
      lambda[2] = -2.0 * sqrt_f - b_3;
      }
    return;
    }

  const double sqrt_f3 = sqrt_f * sqrt_f * sqrt_f;
  const double k = vcl_acos(g / sqrt_f3) / 3.0;
  const double j = 2.0 * sqrt_f;
  lambda[0] = j * vcl_cos(k) - b_3;
  lambda[1] = j * vcl_cos(k + vnl_math::pi * 2.0 / 3.0) - b_3;
  lambda[2] = j * vcl_cos(k - vnl_math::pi * 2.0 / 3.0) - b_3;

  double tmp;
  if (lambda[1] < lambda[0]) 
    {
    tmp = lambda[1]; lambda[1] = lambda[0]; lambda[0] = tmp;
    }

  if (lambda[2] < lambda[1])
    {
    tmp = lambda[1]; lambda[1] = lambda[2]; lambda[2] = tmp;
    if (lambda[1] < lambda[0]) 
      {
      tmp = lambda[1]; lambda[1] = lambda[0]; lambda[0] = tmp;
      }
    }
}

}
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile$
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

     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 __itkTensorToFractionalAnisotropyImageFilter_h
#define __itkTensorToFractionalAnisotropyImageFilter_h

#include "itkImageRegionConstIterator.h"
#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkSymmetricTensor.h"

namespace itk
{
/** \class TensorToFractionalAnisotropyImageFilter
 *
 * \brief This filter computes the fractional anisotropy for each pixel in a
 * diffusion tensor image.
 * 
 * \par Overview
 * The fractional anisotropy coefficient is computed from the eigenvalues of 
 * a 3x3 symmetric tensor matrix. The eigenvalues are computed using VXL's
 * vnl_symmetric_eigensystem class.
 *
 * \par Template Parameters (Input and Output)
 *
 * \par Filter Parameters 
 *
 * \par Constraints
 *
 * \par Performance
 *
 * \par References
 * E. R. Melhem, S. Mori, G. Mukundan, M. A. Kraut, M. G. Pomper, and 
 * P. C. M. van Zijl, "Diffusion tensor MR imaging of the brain and white 
 * matter tractography," Am. J. Roentgenol., vol. 178, pp. 3-16, 2002.
 *
 * \ingroup TensorFilters
 *
 */
template <
  typename TInputImage = Image< SymmetricTensor< float, 3 >, 3 >,
  typename TRealType = typename TInputImage::RealType,
  typename TOutputImage = Image<TRealType,::itk::GetImageDimension<TInputImage>::ImageDimension > 
>
class ITK_EXPORT TensorToFractionalAnisotropyImageFilter :
    public ImageToImageFilter< TInputImage, TOutputImage>
{
public:
  /** Standard class typedefs. */
  typedef TInputImage InputImageType;
  typedef TOutputImage OutputImageType;

  typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
  typedef TensorToFractionalAnisotropyImageFilter Self;

  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  typedef typename InputImageType::PixelType TensorType;
  typedef TRealType RealType;

  typedef typename InputImageType::Pointer InputImagePointer;
  typedef typename OutputImageType::Pointer OutputImagePointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods) */
  itkTypeMacro(TensorToFractionalAnisotropyImageFilter, ImageToImageFilter);

  /** The dimensionality of the input and output images. */
  itkStaticConstMacro(ImageDimension, unsigned int,
                      OutputImageType::ImageDimension);

  /** Type of the iterator that will be used to move through the image.  Also
      the type which will be passed to the evaluate function */
  typedef ImageRegionConstIterator<InputImageType> ImageRegionConstIteratorType;

  /** Superclass typedefs. */
  typedef typename Superclass::OutputImageRegionType OutputImageRegionType;

  /** TensorToFractionalAnisotropyImageFilter needs a larger input requested
   * region than the output requested region (larger by the kernel
   * size to calculate derivatives).  As such,
   * TensorToFractionalAnisotropyImageFilter needs to provide an implementation
   * for GenerateInputRequestedRegion() in order to inform the
   * pipeline execution model.
   *
   * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
  virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);

protected:
  TensorToFractionalAnisotropyImageFilter();
  virtual ~TensorToFractionalAnisotropyImageFilter() {}

  /** Do any necessary casting/copying of the input data.  Input pixel types
     whose value types are not real number types must be cast to real number
     types.*/
  void BeforeThreadedGenerateData ();

  /** TensorToFractionalAnisotropyImageFilter can be implemented as a
   * multithreaded filter.  Therefore, this implementation provides a
   * ThreadedGenerateData() routine which is called for each
   * processing thread. The output image data is allocated
   * automatically by the superclass prior to calling
   * ThreadedGenerateData().  ThreadedGenerateData can only write to
   * the portion of the output image specified by the parameter
   * "outputRegionForThread"
   *
   * \sa ImageToImageFilter::ThreadedGenerateData(),
   *     ImageToImageFilter::GenerateData() */
  void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                            int threadId );

  void PrintSelf(std::ostream& os, Indent indent) const;

  /// Evaluate fractional anisotropy at iterator location.
  TRealType Evaluate(const ImageRegionConstIteratorType &it) const;

private:
  TensorToFractionalAnisotropyImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTensorToFractionalAnisotropyImageFilter.txx"
#endif

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkTensorToFractionalAnisotropyImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2004/04/15 22:37:33 $
  Version:   $Revision: 1.10 $

     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 _itkTensorToFractionalAnisotropyImageFilter_txx
#define _itkTensorToFractionalAnisotropyImageFilter_txx

#include "itkTensorToFractionalAnisotropyImageFilter.h"

#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkProgressReporter.h"
#include "itkVectorCastImageFilter.h"

#include "vnl/vnl_math.h"
#include "vnl/vnl_copy.h"

namespace itk
{

template 
< typename TInputImage, typename TRealType, typename TOutputImage >
void
TensorToFractionalAnisotropyImageFilter< 
TInputImage,TRealType,TOutputImage
>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
}

template 
< typename TInputImage, typename TRealType, typename TOutputImage >
TensorToFractionalAnisotropyImageFilter<
TInputImage,TRealType,TOutputImage
>
::TensorToFractionalAnisotropyImageFilter()
{
}

template 
< typename TInputImage, typename TRealType, typename TOutputImage >
void 
TensorToFractionalAnisotropyImageFilter<
TInputImage,TRealType,TOutputImage
>
::GenerateInputRequestedRegion() throw(InvalidRequestedRegionError)
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

  // get pointers to the input and output
  InputImagePointer  inputPtr = 
    const_cast< InputImageType * >( this->GetInput());
  OutputImagePointer outputPtr = this->GetOutput();

  if ( !inputPtr || !outputPtr )
    {
    return;
    }

  // get a copy of the input requested region (should equal the output
  // requested region)
  typename InputImageType::RegionType inputRequestedRegion;
  inputRequestedRegion = inputPtr->GetRequestedRegion();

  // crop the input requested region at the input's largest possible region
  if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
    {
    inputPtr->SetRequestedRegion( inputRequestedRegion );
    return;
    }
  else
    {
    // Couldn't crop the region (requested region is outside the largest
    // possible region).  Throw an exception.

    // store what we tried to request (prior to trying to crop)
    inputPtr->SetRequestedRegion( inputRequestedRegion );

    // build an exception
    InvalidRequestedRegionError e(__FILE__, __LINE__);
    OStringStream msg;
    msg << static_cast<const char *>(this->GetNameOfClass())
        << "::GenerateInputRequestedRegion()";
    e.SetLocation(msg.str().c_str());
    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
    e.SetDataObject(inputPtr);
    throw e;
    }
}

template
< typename TInputImage, typename TRealType, typename TOutputImage >
void
TensorToFractionalAnisotropyImageFilter<
TInputImage,TRealType,TOutputImage
>
::BeforeThreadedGenerateData()
{
  Superclass::BeforeThreadedGenerateData();
}

template
< typename TInputImage, typename TRealType, typename TOutputImage >
void
TensorToFractionalAnisotropyImageFilter<
TInputImage,TRealType,TOutputImage
>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                       int threadId)
{
  ImageRegionConstIterator<InputImageType> itIn;
  ImageRegionIterator<OutputImageType> itOut;

  itOut = ImageRegionIterator<OutputImageType>
    ( this->GetOutput(), outputRegionForThread );
  itIn = ImageRegionConstIterator<InputImageType>
    ( this->GetInput(), outputRegionForThread );

  itOut.GoToBegin();
  itIn.GoToBegin();

  while ( ! itIn.IsAtEnd() ) 
    {
    itOut.Set( this->Evaluate( itIn ) );
    ++itIn;
    ++itOut;
    }
}

template
< typename TInputImage, typename TRealType, typename TOutputImage >
TRealType
TensorToFractionalAnisotropyImageFilter<TInputImage,TRealType,TOutputImage>
::Evaluate(const ImageRegionConstIteratorType &it) const
{
  double lambda[3];
  const TensorType tensor = it.Get();
  tensor.ComputeEigenvalues( lambda );

  // trace of the diagonal eigenvalue matrix
  const TRealType lambdaMean = (lambda[0] + lambda[1] + lambda[2])/3.0;

  const TRealType tmp0 = lambda[0] - lambdaMean;
  const TRealType tmp1 = lambda[1] - lambdaMean;
  const TRealType tmp2 = lambda[2] - lambdaMean;

  const TRealType faSquare =
    (3 * (tmp0*tmp0 + tmp1*tmp1 + tmp2*tmp2)) / 
    (2 * (lambda[0]*lambda[0] + lambda[1]*lambda[1] + lambda[2]*lambda[2]));

  return sqrt( faSquare );
}

} // end namespace itk

#endif