Slicer  4.10
Slicer is a multi-platform, free and open source software package for visualization and medical image computing
itkDiffusionTensor3DNearestCorrection.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Diffusion Applications
4  Module: $HeadURL$
5  Language: C++
6  Date: $Date$
7  Version: $Revision$
8 
9  Copyright (c) Brigham and Women's Hospital (BWH) All Rights Reserved.
10 
11  See License.txt or http://www.slicer.org/copyright/copyright.txt for details.
12 
13 ==========================================================================*/
14 #ifndef itkDiffusionTensor3DNearestCorrection_h
15 #define itkDiffusionTensor3DNearestCorrection_h
16 
17 #include "itkUnaryFunctorImageFilter.h"
18 #include "vnl/vnl_math.h"
19 #include <itkMatrix.h>
21 
22 namespace itk
23 {
24 
40 namespace Functor
41 {
42 
43 template <class TInput, class TOutput>
45 {
46 public:
48  {
49  }
51  {
52  }
53  bool operator!=( const DiffusionTensor3DNearest & other ) const
54  {
55  return *this != other;
56  }
57 
58  bool operator==( const DiffusionTensor3DNearest & other ) const
59  {
60  return !( *this != other );
61  }
62 
63  inline DiffusionTensor3D<TOutput> operator()
64  ( const DiffusionTensor3D<TInput> & tensorA )
65  {
66  DiffusionTensor3DExtended<double> tensorDouble( tensorA );
67  Matrix<double, 3, 3> B;
68  Matrix<double, 3, 3> A;
69  Matrix<double, 3, 3> transpose;
70  Matrix<double, 3, 3> H;
71  Matrix<double, 3, 3> mat;
72  A = tensorDouble.GetTensor2Matrix();
73  transpose = A.GetTranspose();
74  B = (A + transpose) / 2;
75  transpose = B.GetTranspose();
76  H = transpose * B;
77  tensorDouble.SetTensorFromMatrix(H);
78 
81  tensorDouble.ComputeEigenAnalysis( eigenValues, eigenVectors );
82  for( int i = 0; i < 3; i++ )
83  {
84  mat[i][i] = sqrt(eigenValues[i]);
85  }
86  eigenVectors = eigenVectors.GetTranspose();
87  H = eigenVectors * mat * eigenVectors.GetInverse();
88  mat = (B + H) / 2;
89  tensorDouble.SetTensorFromMatrix( mat );
90  tensorDouble.ComputeEigenAnalysis( eigenValues, eigenVectors ); // sometimes very small negative eigenvalues
91  // appear; we suppress them
92  mat.Fill(0);
93  for( int i = 0; i < 3; i++ )
94  {
95  mat[i][i] = ( eigenValues[i] <= 0 ? ITK_DIFFUSION_TENSOR_3D_ZERO : eigenValues[i] );
96  }
97  eigenVectors = eigenVectors.GetTranspose();
98  tensorDouble.SetTensorFromMatrix<double>( eigenVectors * mat * eigenVectors.GetInverse() );
99 
100  DiffusionTensor3D<TOutput> tensor;
101  for( int i = 0; i < 6; i++ )
102  {
103  tensor[i] = ( TOutput ) tensorDouble[i];
104  }
105  return tensor;
106  }
107 
108 };
109 } // end of Functor namespace
110 
111 template <class TInputImage, class TOutputImage>
113  public
114  UnaryFunctorImageFilter<TInputImage, TOutputImage,
115  Functor::DiffusionTensor3DNearest<
116  typename TInputImage::PixelType::ComponentType,
117  typename TOutputImage::PixelType::ComponentType> >
118 {
119 public:
122  typedef UnaryFunctorImageFilter<TInputImage, TOutputImage,
123  Functor::DiffusionTensor3DNearest<typename TInputImage::PixelType,
124  typename TOutputImage::PixelType> > Superclass;
125  typedef SmartPointer<Self> Pointer;
126  typedef SmartPointer<const Self> ConstPointer;
127 
129  itkTypeMacro(DiffusionTensor3DNearestCorrectionFilter, UnaryFunctorImageFilter);
130 
132  itkNewMacro( Self );
133 
134 #ifdef ITK_USE_CONCEPT_CHECKING
135 
136  itkConceptMacro( InputTensorTypeCheck,
137  ( Concept::SameType<DiffusionTensor3D<typename TInputImage::PixelType::ComponentType>,
138  typename TInputImage::PixelType> ) );
139  itkConceptMacro( OutputTensorTypeCheck,
140  ( Concept::SameType<DiffusionTensor3D<typename TOutputImage::PixelType::ComponentType>,
141  typename TOutputImage::PixelType> ) );
142 
144 #endif
145 protected:
147  {
148  }
150  {
151  }
152 private:
153  DiffusionTensor3DNearestCorrectionFilter( const Self & ); // purposely not implemented
154  void operator=( const Self & ); // purposely not implemented
155 
156 };
157 
158 } // end namespace itk
159 
160 #endif
void SetTensorFromMatrix(Matrix< C, 3, 3 > matrix)
bool operator!=(const DiffusionTensor3DNearest &other) const
UnaryFunctorImageFilter< TInputImage, TOutputImage, Functor::DiffusionTensor3DNearest< typename TInputImage::PixelType, typename TOutputImage::PixelType > > Superclass
Simplified inverse ITK transforms.
* itkConceptMacro(OutputEqualityComparableCheck, *(Concept::EqualityComparable< OutputImagePixelType >))
itkTypeMacro(DiffusionTensor3DNearestCorrectionFilter, UnaryFunctorImageFilter)
#define ITK_DIFFUSION_TENSOR_3D_ZERO
bool operator==(const DiffusionTensor3DNearest &other) const