ANIMA  4.0
animaMTPairingCorrelationImageToImageMetric.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <vnl/vnl_matrix_fixed.h>
6 #include <itkImageRegionConstIteratorWithIndex.h>
7 
8 #include <boost/math/special_functions/factorials.hpp>
9 
10 namespace anima
11 {
12 
13 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
16 {
17  m_FixedImagePoints.clear();
18  m_FixedImageCompartmentWeights.clear();
19  m_FixedImageLogTensors.clear();
20  m_NumberOfFixedCompartments = 1;
21 
23  mcmCreator.SetNumberOfCompartments(0);
24  mcmCreator.SetModelWithStationaryWaterComponent(true);
25  mcmCreator.SetModelWithFreeWaterComponent(false);
26 
27  m_ZeroDiffusionModel = mcmCreator.GetNewMultiCompartmentModel();
28 
29  m_leCalculator = LECalculatorType::New();
30 }
31 
32 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
33 bool
36 {
37  FixedImageType *fixedImage = const_cast <FixedImageType *> (this->GetFixedImage());
38  MCModelPointer val = fixedImage->GetDescriptionModel();
39  for (unsigned int i = 0;i < val->GetNumberOfCompartments();++i)
40  {
41  if (!val->GetCompartment(i)->GetTensorCompatible())
42  return false;
43  }
44 
45  MovingImageType *movingImage = const_cast <MovingImageType *> (this->GetFixedImage());
46  val = movingImage->GetDescriptionModel();
47  for (unsigned int i = 0;i < val->GetNumberOfCompartments();++i)
48  {
49  if (!val->GetCompartment(i)->GetTensorCompatible())
50  return false;
51  }
52 
53  return true;
54 }
55 
56 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
59 ::GetValue(const TransformParametersType & parameters) const
60 {
61  FixedImageConstPointer fixedImage = this->m_FixedImage;
62 
63  if( !fixedImage )
64  itkExceptionMacro("Fixed image has not been assigned");
65 
66  bool tensorCompatibilityCondition = this->CheckTensorCompatibility();
67  if (!tensorCompatibilityCondition)
68  itkExceptionMacro("Only tensor compatible models handled")
69 
70  if (this->m_NumberOfPixelsCounted == 0)
71  return 0;
72 
73  this->SetTransformParameters( parameters );
74 
75  PixelType movingValue, workValue;
76 
77  OutputPointType transformedPoint;
78  ContinuousIndexType transformedIndex;
79 
80  MovingImageType *movingImage = const_cast <MovingImageType *> (this->GetMovingImage());
81  MCModelPointer currentMovingValue = movingImage->GetDescriptionModel()->Clone();
82 
83  std::vector < std::vector <double> > movingImageCompartmentWeights(this->m_NumberOfPixelsCounted);
84  std::vector < std::vector <PixelType> > movingImageLogTensors(this->m_NumberOfPixelsCounted);
85  std::vector <double> tmpWeights;
86  vnl_matrix <double> workLogMatrix(3,3);
87 
88  // Getting moving values
89  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
90  {
91  transformedPoint = this->m_Transform->TransformPoint( m_FixedImagePoints[i] );
92  this->m_Interpolator->GetInputImage()->TransformPhysicalPointToContinuousIndex(transformedPoint,transformedIndex);
93 
94  if( this->m_Interpolator->IsInsideBuffer( transformedIndex ) )
95  {
96  movingValue = this->m_Interpolator->EvaluateAtContinuousIndex( transformedIndex );
97 
98  if (!isZero(movingValue))
99  {
100  currentMovingValue->SetModelVector(movingValue);
101  if (this->GetModelRotation() != Superclass::NONE)
102  currentMovingValue->Reorient(this->m_OrientationMatrix, (this->GetModelRotation() == Superclass::PPD));
103 
104  tmpWeights = currentMovingValue->GetCompartmentWeights();
105  unsigned int internalCounter = 0;
106  for (unsigned int j = 0;j < currentMovingValue->GetNumberOfCompartments();++j)
107  {
108  if (tmpWeights[internalCounter] <= 0)
109  tmpWeights.erase(tmpWeights.begin() + internalCounter);
110  else
111  {
112  ++internalCounter;
113  m_leCalculator->GetTensorLogarithm(currentMovingValue->GetCompartment(j)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),workLogMatrix);
114  anima::GetVectorRepresentation(workLogMatrix,workValue,6,true);
115  movingImageLogTensors[i].push_back(workValue);
116  }
117  }
118 
119  movingImageCompartmentWeights[i] = tmpWeights;
120  }
121  else
122  {
123  movingImageLogTensors[i].resize(1);
124  anima::GetVectorRepresentation(m_ZeroDiffusionModel->GetCompartment(0)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),movingImageLogTensors[i][0],6,true);
125  movingImageCompartmentWeights[i].resize(1);
126  movingImageCompartmentWeights[i][0] = 1.0;
127  }
128  }
129  else
130  {
131  movingImageLogTensors[i].resize(1);
132  anima::GetVectorRepresentation(m_ZeroDiffusionModel->GetCompartment(0)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),movingImageLogTensors[i][0],6,true);
133  movingImageCompartmentWeights[i].resize(1);
134  movingImageCompartmentWeights[i][0] = 1.0;
135  }
136  }
137 
138  double mRS = this->ComputeMapping(m_FixedImageCompartmentWeights,m_FixedImageLogTensors,movingImageCompartmentWeights,movingImageLogTensors);
139  double mRR = this->ComputeMapping(m_FixedImageCompartmentWeights,m_FixedImageLogTensors,m_FixedImageCompartmentWeights,m_FixedImageLogTensors);
140  double mSS = this->ComputeMapping(movingImageCompartmentWeights,movingImageLogTensors,movingImageCompartmentWeights,movingImageLogTensors);
141  double mRT = 0;
142  double mST = 0;
143  double numMaxCompartments = std::max(movingImage->GetDescriptionModel()->GetNumberOfCompartments(),m_NumberOfFixedCompartments);
144  double epsilon = std::sqrt(numMaxCompartments / (3.0 * this->m_NumberOfPixelsCounted)) / numMaxCompartments;
145 
146  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
147  {
148  for (unsigned int k = 0;k < m_FixedImageCompartmentWeights[i].size();++k)
149  {
150  double dotProd = 0;
151  for (unsigned int l = 0;l < 3;++l)
152  {
153  unsigned int index = l * (l + 1) / 2 - 1;
154  dotProd += m_FixedImageLogTensors[i][k][index];
155  }
156 
157  mRT += m_FixedImageCompartmentWeights[i][k] * dotProd;
158  }
159 
160  for (unsigned int k = 0;k < movingImageCompartmentWeights[i].size();++k)
161  {
162  double dotProd = 0;
163  for (unsigned int l = 0;l < 3;++l)
164  {
165  unsigned int index = l * (l + 1) / 2 - 1;
166  dotProd += movingImageLogTensors[i][k][index];
167  }
168 
169  mST += movingImageCompartmentWeights[i][k] * dotProd;
170  }
171  }
172 
173  mRT *= epsilon;
174  mST *= epsilon;
175 
176  // Now computing the measure itself, going for some one to one pairing
177  double measure = (mRS - mRT * mST) * (mRS - mRT * mST);
178  double denom = (mRR - mRT * mRT) * (mSS - mST * mST);
179 
180  if (denom > 0)
181  measure /= denom;
182  else
183  measure = 0;
184 
185  return measure;
186 }
187 
188 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
189 double
191 ::ComputeMapping(const std::vector < std::vector <double> > &refImageCompartmentWeights, const std::vector < std::vector <PixelType> > &refImageLogTensors,
192  const std::vector < std::vector <double> > &movingImageCompartmentWeights, const std::vector < std::vector <PixelType> > &movingImageLogTensors) const
193 {
194  std::vector <unsigned int> currentPermutation;
195  double mappingDistanceValue = 0;
196  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
197  {
198  double bestValue = 0.0;
199  unsigned int fixedNumCompartments = refImageCompartmentWeights[i].size();
200  unsigned int movingNumCompartments = movingImageCompartmentWeights[i].size();
201 
202  unsigned int minCompartmentsNumber = fixedNumCompartments;
203  int rest = movingNumCompartments - minCompartmentsNumber;
204  bool fixedMin = true;
205  if (rest < 0)
206  {
207  fixedMin = false;
208  minCompartmentsNumber = movingNumCompartments;
209  rest *= -1;
210  }
211 
212  unsigned int maxCompartmentsNumber = minCompartmentsNumber + rest;
213 
214  currentPermutation.resize(maxCompartmentsNumber);
215  for (unsigned int j = 0;j < maxCompartmentsNumber;++j)
216  currentPermutation[j] = j;
217 
218  // Test all permutations
219  do
220  {
221  double distValue = 0;
222 
223  for (unsigned int j = 0;j < maxCompartmentsNumber;++j)
224  {
225  unsigned int firstIndex;
226  unsigned int secondIndex;
227 
228  if (fixedMin)
229  {
230  firstIndex = currentPermutation[j];
231  secondIndex = j;
232  }
233  else
234  {
235  firstIndex = j;
236  secondIndex = currentPermutation[j];
237  }
238 
239  if ((firstIndex >= refImageCompartmentWeights[i].size())||(secondIndex >= movingImageCompartmentWeights[i].size()))
240  continue;
241 
242  double dist = 0;
243  for (unsigned int k = 0;k < refImageLogTensors[i][firstIndex].GetSize();++k)
244  dist += refImageLogTensors[i][firstIndex][k] * movingImageLogTensors[i][secondIndex][k];
245 
246  distValue += refImageCompartmentWeights[i][firstIndex] * movingImageCompartmentWeights[i][secondIndex] * dist;
247  }
248 
249  if (std::abs(distValue) > std::abs(bestValue))
250  bestValue = distValue;
251 
252  } while(std::next_permutation(currentPermutation.begin(),currentPermutation.end()));
253 
254  mappingDistanceValue += bestValue;
255  }
256 
257  return mappingDistanceValue;
258 }
259 
260 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
261 bool
263 ::isZero(PixelType &vector) const
264 {
265  unsigned int ndim = vector.GetSize();
266 
267  for (unsigned int i = 0;i < ndim;++i)
268  {
269  if (vector[i] != 0)
270  return false;
271  }
272 
273  return true;
274 }
275 
276 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
277 void
280 {
281  if(!this->m_FixedImage)
282  itkExceptionMacro( << "Fixed image has not been assigned" );
283 
284  FixedImageType *fixedImage = const_cast <FixedImageType *> (this->GetFixedImage());
285 
286  this->m_NumberOfPixelsCounted = this->GetFixedImageRegion().GetNumberOfPixels();
287  typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
288  FixedIteratorType ti(fixedImage, this->GetFixedImageRegion());
289  typename FixedImageType::IndexType index;
290 
291  m_FixedImagePoints.resize(this->m_NumberOfPixelsCounted);
292  m_FixedImageCompartmentWeights.resize(this->m_NumberOfPixelsCounted);
293  m_FixedImageLogTensors.resize(this->m_NumberOfPixelsCounted);
294 
295  InputPointType inputPoint;
296  MCModelPointer fixedMCM = fixedImage->GetDescriptionModel()->Clone();
297  m_NumberOfFixedCompartments = fixedMCM->GetNumberOfCompartments();
298 
299  unsigned int pos = 0;
300  PixelType fixedValue, workValue;
301  std::vector <double> tmpWeights;
302  vnl_matrix <double> workLogMatrix(3,3);
303 
304  while(!ti.IsAtEnd())
305  {
306  index = ti.GetIndex();
307  fixedImage->TransformIndexToPhysicalPoint(index, inputPoint);
308 
309  m_FixedImagePoints[pos] = inputPoint;
310  fixedValue = ti.Get();
311 
312  if (!isZero(fixedValue))
313  {
314  fixedMCM->SetModelVector(fixedValue);
315  tmpWeights = fixedMCM->GetCompartmentWeights();
316  unsigned int internalCounter = 0;
317  for (unsigned int i = 0;i < fixedMCM->GetNumberOfCompartments();++i)
318  {
319  if (tmpWeights[internalCounter] <= 0)
320  tmpWeights.erase(tmpWeights.begin() + internalCounter);
321  else
322  {
323  ++internalCounter;
324  m_leCalculator->GetTensorLogarithm(fixedMCM->GetCompartment(i)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),workLogMatrix);
325  anima::GetVectorRepresentation(workLogMatrix,workValue,6,true);
326  m_FixedImageLogTensors[pos].push_back(workValue);
327  }
328  }
329 
330  m_FixedImageCompartmentWeights[pos] = tmpWeights;
331  }
332  else
333  {
334  m_FixedImageLogTensors[pos].resize(1);
335  anima::GetVectorRepresentation(m_ZeroDiffusionModel->GetCompartment(0)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),m_FixedImageLogTensors[pos][0],6,true);
336  m_FixedImageCompartmentWeights[pos].resize(1);
337  m_FixedImageCompartmentWeights[pos][0] = 1.0;
338  }
339 
340  ++ti;
341  ++pos;
342  }
343 }
344 
345 } // end 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)
MeasureType GetValue(const TransformParametersType &parameters) const ITK_OVERRIDE
double ComputeMapping(const std::vector< std::vector< double > > &refImageCompartmentWeights, const std::vector< std::vector< PixelType > > &refImageLogTensors, const std::vector< std::vector< double > > &movingImageCompartmentWeights, const std::vector< std::vector< PixelType > > &movingImageLogTensors) const
Really this class is some simplified factory that creates the MCM that it knows.