ANIMA  4.0
animaTensorCorrelationImageToImageMetric.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionConstIteratorWithIndex.h>
5 
6 #include <animaBaseTensorTools.h>
7 
8 namespace anima
9 {
10 
14 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
17 {
18  m_FixedImagePoints.clear();
19  m_FixedImageValues.clear();
20 
21  m_FixedTProduct = 0;
22  m_FixedDenominator = 0;
23 
24  m_LogEpsilon = 0;
25 }
26 
30 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
33 ::GetValue( const TransformParametersType & parameters ) const
34 {
35  FixedImageConstPointer fixedImage = this->m_FixedImage;
36 
37  if( !fixedImage )
38  {
39  itkExceptionMacro( << "Fixed image has not been assigned" );
40  }
41 
42  if ((this->m_NumberOfPixelsCounted <= 1)||(m_FixedDenominator == 0))
43  return 0;
44 
45  this->SetTransformParameters(parameters);
46 
47  unsigned int vectorSize = m_FixedImageValues[0].GetSize();
48  PixelType movingValue;
49 
50  OutputPointType transformedPoint;
51  ContinuousIndexType transformedIndex;
52 
53  unsigned int tensorDimension = floor((std::sqrt((double)(8 * vectorSize + 1)) - 1) / 2.0);
54 
55  vnl_matrix <double> tmpMat(tensorDimension, tensorDimension);
56  vnl_matrix <double> ppdOrientationMatrix(tensorDimension, tensorDimension);
57  typedef itk::Matrix <double, 3, 3> EigVecMatrixType;
58  typedef vnl_vector_fixed <double,3> EigValVectorType;
59  itk::SymmetricEigenAnalysis < EigVecMatrixType, EigValVectorType, EigVecMatrixType> eigenComputer(3);
60  EigVecMatrixType eigVecs;
61  EigValVectorType eigVals;
62 
63  vnl_matrix <double> currentTensor(tensorDimension, tensorDimension);
64  double mST = 0, mRS = 0, mSS = 0;
65 
66  double movingDenominator = 0;
67  double measureNumerator = 0;
68  unsigned int numOutside = 0;
69 
70  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
71  {
72  transformedPoint = this->m_Transform->TransformPoint(m_FixedImagePoints[i]);
73  this->m_Interpolator->GetInputImage()->TransformPhysicalPointToContinuousIndex(transformedPoint,transformedIndex);
74 
75  if( this->m_Interpolator->IsInsideBuffer(transformedIndex))
76  {
77  movingValue = this->m_Interpolator->EvaluateAtContinuousIndex(transformedIndex);
78 
79  if (this->GetModelRotation() != Superclass::NONE)
80  {
81  // Rotating tensor
82  anima::GetTensorFromVectorRepresentation(movingValue,tmpMat,tensorDimension,true);
83 
84  if (this->GetModelRotation() == Superclass::FINITE_STRAIN)
85  anima::RotateSymmetricMatrix(tmpMat,this->m_OrientationMatrix,currentTensor);
86  else
87  {
88  eigenComputer.ComputeEigenValuesAndVectors(tmpMat,eigVals,eigVecs);
89  anima::ExtractPPDRotationFromJacobianMatrix(this->m_OrientationMatrix,ppdOrientationMatrix,eigVecs);
90  anima::RotateSymmetricMatrix(tmpMat,ppdOrientationMatrix,currentTensor);
91  }
92 
93  anima::GetVectorRepresentation(currentTensor,movingValue,vectorSize,true);
94  }
95 
96  unsigned int pos_internal = 0;
97  for (unsigned int j = 0;j < tensorDimension;++j)
98  for (unsigned int k = 0;k <= j;++k)
99  {
100  if (i == j)
101  mST += movingValue[pos_internal] * m_LogEpsilon;
102 
103  mRS += m_FixedImageValues[i][pos_internal] * movingValue[pos_internal];
104  mSS += movingValue[pos_internal] * movingValue[pos_internal];
105 
106  ++pos_internal;
107  }
108  }
109  else
110  ++numOutside;
111  }
112 
113  if (this->m_NumberOfPixelsCounted / 2.0 < numOutside)
114  return 0.0;
115 
116  movingDenominator = mSS + mST * mST * (3.0 * this->m_NumberOfPixelsCounted * m_LogEpsilon * m_LogEpsilon - 2.0);
117 
118  if (movingDenominator == 0)
119  return 0.0;
120 
121  measureNumerator = mRS - mST * m_FixedTProduct;
122 
123  double measure = measureNumerator * measureNumerator / (m_FixedDenominator * movingDenominator);
124 
125  if (measure > 1)
126  measure = 1;
127  else if (measure < 0)
128  measure = 0;
129 
130  return measure;
131 }
132 
133 template <class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension>
134 void
137 {
138  FixedImageConstPointer fixedImage = this->m_FixedImage;
139 
140  if( !fixedImage )
141  {
142  itkExceptionMacro( << "Fixed image has not been assigned" );
143  }
144 
145  unsigned int vectorSize = fixedImage->GetNumberOfComponentsPerPixel();
146  unsigned int tensorDimension = floor((sqrt((double)(8 * vectorSize + 1)) - 1) / 2.0);
147 
148  this->m_NumberOfPixelsCounted = this->GetFixedImageRegion().GetNumberOfPixels();
149 
150  m_FixedDenominator = 0;
151  m_FixedTProduct = 0;
152 
153  // Choose log-epsilon wisely so that m(T,T) = 1
154  m_LogEpsilon = 1.0 / std::sqrt(3.0 * this->m_NumberOfPixelsCounted);
155 
156  typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
157 
158  FixedIteratorType ti(fixedImage, this->GetFixedImageRegion());
159  typename FixedImageType::IndexType index;
160 
161  m_FixedImagePoints.resize(this->m_NumberOfPixelsCounted);
162  m_FixedImageValues.resize(this->m_NumberOfPixelsCounted);
163 
164  InputPointType inputPoint;
165 
166  unsigned int pos = 0;
167  PixelType fixedValue;
168  while(!ti.IsAtEnd())
169  {
170  index = ti.GetIndex();
171  fixedImage->TransformIndexToPhysicalPoint(index, inputPoint);
172 
173  m_FixedImagePoints[pos] = inputPoint;
174  fixedValue = ti.Get();
175  m_FixedImageValues[pos] = fixedValue;
176 
177  unsigned int pos_internal = 0;
178  for (unsigned int i = 0;i < tensorDimension;++i)
179  for (unsigned int j = 0;j <= i;++j)
180  {
181  if (i == j)
182  m_FixedTProduct += fixedValue[pos_internal] * m_LogEpsilon;
183 
184  m_FixedDenominator += fixedValue[pos_internal] * fixedValue[pos_internal];
185 
186  ++pos_internal;
187  }
188 
189  ++ti;
190  ++pos;
191  }
192 
193  m_FixedDenominator += m_FixedTProduct * m_FixedTProduct * (3.0 * this->m_NumberOfPixelsCounted * m_LogEpsilon * m_LogEpsilon - 2.0);
194 }
195 
196 } // end of namespace anima
itk::ContinuousIndex< double, ImageDimension > ContinuousIndexType
void GetVectorRepresentation(const vnl_matrix< T1 > &tensor, itk::VariableLengthVector< T2 > &vector, unsigned int vecDim=0, bool scale=false)
void RotateSymmetricMatrix(T1 &tensor, T2 &rotationMatrix, T3 &rotated_tensor, unsigned int tensorDim)
Rotates a symmetric matrix by performing R^T D R where R is a rotation matrix and D the symmetric mat...
void GetTensorFromVectorRepresentation(const itk::VariableLengthVector< T1 > &vector, vnl_matrix< T2 > &tensor, unsigned int tensDim=0, bool scale=false)
MeasureType GetValue(const TransformParametersType &parameters) const ITK_OVERRIDE
void ExtractPPDRotationFromJacobianMatrix(vnl_matrix< RealType > &jacobianMatrix, vnl_matrix< RealType > &rotationMatrix, MatrixType &eigenVectors)