itkDiffusionTensor3DNearestCorrection.h

Go to the documentation of this file.

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

  Program:   Diffusion Applications
  Module:    $HeadURL: http://svn.slicer.org/Slicer3/trunk/Applications/CLI/DiffusionApplications/ResampleDTI/itkDiffusionTensor3DNearestCorrection.h $
  Language:  C++
  Date:      $Date: 2010-04-29 11:58:49 -0400 (Thu, 29 Apr 2010) $
  Version:   $Revision: 13073 $

  Copyright (c) Brigham and Women's Hospital (BWH) All Rights Reserved.

  See License.txt or http://www.slicer.org/copyright/copyright.txt for details.

==========================================================================*/
#ifndef __itkDiffusionTensor3DNearestCorrectionFilter_h
#define __itkDiffusionTensor3DNearestCorrectionFilter_h

#include "itkUnaryFunctorImageFilter.h"
#include "vnl/vnl_math.h"
#include <itkMatrix.h>
#include "itkDiffusionTensor3DExtended.h"

namespace itk
{
  
namespace Functor {  
  
template< class TInput , class TOutput >
class DiffusionTensor3DNearest
{
public:
  DiffusionTensor3DNearest() {}
  ~DiffusionTensor3DNearest() {}
  bool operator!=( const DiffusionTensor3DNearest & other ) const
  {
    return ( *this != other ) ;
  }
  bool operator==( const DiffusionTensor3DNearest & other ) const
  {
    return !( *this != other ) ;
  }
  inline DiffusionTensor3D< TOutput > operator()
       ( const DiffusionTensor3D< TInput > & tensorA )
  {
    DiffusionTensor3D< TOutput > tensor;
    DiffusionTensor3DExtended< double > tensorDouble ;
    tensorDouble = ( DiffusionTensor3DExtended< TInput > ) tensorA ;
    Matrix< double , 3 , 3 > B ;
    Matrix< double , 3 , 3 > A ;
    Matrix< double , 3 , 3 > transpose ;
    Matrix< double , 3 , 3 > H ;
    Matrix< double , 3 , 3 > mat ;
    A=tensorDouble.GetTensor2Matrix();
    transpose=A.GetTranspose();
    B=(A+transpose)/2;
    transpose=B.GetTranspose();
    H=transpose*B;
    tensorDouble.SetTensorFromMatrix(H);
    
    typename DiffusionTensor3DExtended< double >::EigenValuesArrayType eigenValues ;
    typename DiffusionTensor3DExtended< double >::EigenVectorsMatrixType eigenVectors ;
    tensorDouble.ComputeEigenAnalysis( eigenValues , eigenVectors ) ;
    for( int i = 0 ; i < 3 ; i++ )
      {
      mat[ i ][ i ] = sqrt(eigenValues[ i ]) ;
      }
    eigenVectors = eigenVectors.GetTranspose();
    H = eigenVectors * mat * eigenVectors.GetInverse() ;
    mat=(B+H)/2;
    tensorDouble.SetTensorFromMatrix( mat ) ;
    tensorDouble.ComputeEigenAnalysis( eigenValues , eigenVectors ) ;//sometimes very small negative eigenvalues appear; we suppress them
    mat.Fill(0);
    for( int i = 0 ; i < 3 ; i++ )
      {
      mat[ i ][ i ] = ( eigenValues[ i ] <= 0 ? ZERO : eigenValues[ i ] ) ;
      }
    eigenVectors = eigenVectors.GetTranspose();
    tensorDouble.SetTensorFromMatrix< double >( eigenVectors * mat * eigenVectors.GetInverse() );
    for( int i = 0 ; i < 6 ; i++ )
      { tensor[ i ] = ( TOutput ) tensorDouble[ i ] ; }
    return tensor ;
  }
}; 
}//end of Functor namespace

template < class TInputImage, class TOutputImage >
class DiffusionTensor3DNearestCorrectionFilter :
    public
  UnaryFunctorImageFilter< TInputImage , TOutputImage , 
                        Functor::DiffusionTensor3DNearest< 
  typename TInputImage::PixelType::ComponentType , 
  typename TOutputImage::PixelType::ComponentType> >
{
public:
  typedef DiffusionTensor3DNearestCorrectionFilter Self ;
  typedef UnaryFunctorImageFilter< TInputImage , TOutputImage , 
           Functor::DiffusionTensor3DNearest< typename TInputImage::PixelType , 
                                typename TOutputImage::PixelType> >  Superclass ;
  typedef SmartPointer< Self > Pointer ;
  typedef SmartPointer< const Self >  ConstPointer ;

  itkNewMacro( Self ) ;

#ifdef ITK_USE_CONCEPT_CHECKING

  itkConceptMacro( InputTensorTypeCheck ,
    ( Concept::SameType< DiffusionTensor3D< typename TInputImage::PixelType::ComponentType > ,
                                            typename TInputImage::PixelType > ) ) ;
  itkConceptMacro( OutputTensorTypeCheck ,
    ( Concept::SameType< DiffusionTensor3D< typename TOutputImage::PixelType::ComponentType > ,
                                            typename TOutputImage::PixelType>) ) ;

#endif

protected:
  DiffusionTensor3DNearestCorrectionFilter() {}
  virtual ~DiffusionTensor3DNearestCorrectionFilter() {}

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

};

} // end namespace itk


#endif