Slicer  5.0
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 "itkMath.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:
47  DiffusionTensor3DNearest() = default;
48  ~DiffusionTensor3DNearest() = default;
49  bool operator!=( const DiffusionTensor3DNearest & other ) const
50  {
51  return *this != other;
52  }
53 
54  bool operator==( const DiffusionTensor3DNearest & other ) const
55  {
56  return !( *this != other );
57  }
58 
59  inline DiffusionTensor3D<TOutput> operator()
60  ( const DiffusionTensor3D<TInput> & tensorA )
61  {
62  DiffusionTensor3DExtended<double> tensorDouble( tensorA );
63  Matrix<double, 3, 3> B;
64  Matrix<double, 3, 3> A;
65  Matrix<double, 3, 3> transpose;
66  Matrix<double, 3, 3> H;
67  Matrix<double, 3, 3> mat;
68  A = tensorDouble.GetTensor2Matrix();
69  transpose = A.GetTranspose();
70  B = (A + transpose) / 2;
71  transpose = B.GetTranspose();
72  H = transpose * B;
73  tensorDouble.SetTensorFromMatrix(H);
74 
77  tensorDouble.ComputeEigenAnalysis( eigenValues, eigenVectors );
78  for( int i = 0; i < 3; i++ )
79  {
80  mat[i][i] = sqrt(eigenValues[i]);
81  }
82  eigenVectors = eigenVectors.GetTranspose();
83  H = eigenVectors * mat * eigenVectors.GetInverse();
84  mat = (B + H) / 2;
85  tensorDouble.SetTensorFromMatrix( mat );
86  tensorDouble.ComputeEigenAnalysis( eigenValues, eigenVectors ); // sometimes very small negative eigenvalues
87  // appear; we suppress them
88  mat.Fill(0);
89  for( int i = 0; i < 3; i++ )
90  {
91  mat[i][i] = ( eigenValues[i] <= 0 ? ITK_DIFFUSION_TENSOR_3D_ZERO : eigenValues[i] );
92  }
93  eigenVectors = eigenVectors.GetTranspose();
94  tensorDouble.SetTensorFromMatrix<double>( eigenVectors * mat * eigenVectors.GetInverse() );
95 
96  DiffusionTensor3D<TOutput> tensor;
97  for( int i = 0; i < 6; i++ )
98  {
99  tensor[i] = ( TOutput ) tensorDouble[i];
100  }
101  return tensor;
102  }
103 
104 };
105 } // end of Functor namespace
106 
107 template <class TInputImage, class TOutputImage>
109  public
110  UnaryFunctorImageFilter<TInputImage, TOutputImage,
111  Functor::DiffusionTensor3DNearest<
112  typename TInputImage::PixelType::ComponentType,
113  typename TOutputImage::PixelType::ComponentType> >
114 {
115 public:
118  typedef UnaryFunctorImageFilter<TInputImage, TOutputImage,
119  Functor::DiffusionTensor3DNearest<typename TInputImage::PixelType,
120  typename TOutputImage::PixelType> > Superclass;
121  typedef SmartPointer<Self> Pointer;
122  typedef SmartPointer<const Self> ConstPointer;
123 
126  void operator=( const Self & ) = delete;
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  ~DiffusionTensor3DNearestCorrectionFilter() override = default;
148 };
149 
150 } // end namespace itk
151 
152 #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
~DiffusionTensor3DNearestCorrectionFilter() override=default
bool operator==(const DiffusionTensor3DNearest &other) const