ANIMA  4.0
animaMCMPairingMeanSquaresImageToImageMetric.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_FixedImageValues.clear();
19 
21  mcmCreator.SetNumberOfCompartments(0);
22  mcmCreator.SetModelWithStationaryWaterComponent(true);
23  mcmCreator.SetModelWithFreeWaterComponent(false);
24 
25  m_ZeroDiffusionModel = mcmCreator.GetNewMultiCompartmentModel();
26  m_ZeroDiffusionVector = m_ZeroDiffusionModel->GetModelVector();
27 
28  m_OneToOneMapping = false;
29 
30  m_leCalculator = LECalculatorType::New();
31 }
32 
33 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
34 bool
37 {
38  FixedImageType *fixedImage = const_cast <FixedImageType *> (this->GetFixedImage());
39  MCModelPointer val = fixedImage->GetDescriptionModel();
40  for (unsigned int i = 0;i < val->GetNumberOfCompartments();++i)
41  {
42  if (!val->GetCompartment(i)->GetTensorCompatible())
43  return false;
44  }
45 
46  MovingImageType *movingImage = const_cast <MovingImageType *> (this->GetFixedImage());
47  val = movingImage->GetDescriptionModel();
48  for (unsigned int i = 0;i < val->GetNumberOfCompartments();++i)
49  {
50  if (!val->GetCompartment(i)->GetTensorCompatible())
51  return false;
52  }
53 
54  return true;
55 }
56 
57 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
60 ::GetValue(const TransformParametersType & parameters) const
61 {
62  FixedImageConstPointer fixedImage = this->m_FixedImage;
63 
64  if( !fixedImage )
65  itkExceptionMacro("Fixed image has not been assigned");
66 
67  if (this->m_NumberOfPixelsCounted == 0)
68  return 0;
69 
70  this->SetTransformParameters( parameters );
71 
72  PixelType movingValue;
73 
74  OutputPointType transformedPoint;
75  ContinuousIndexType transformedIndex;
76 
77  MovingImageType *movingImage = const_cast <MovingImageType *> (this->GetMovingImage());
78  MCModelPointer currentMovingValue = movingImage->GetDescriptionModel()->Clone();
79 
80  double measure = 0;
81  bool tensorCompatibilityCondition = this->CheckTensorCompatibility();
82 
83  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
84  {
85  transformedPoint = this->m_Transform->TransformPoint( m_FixedImagePoints[i] );
86  this->m_Interpolator->GetInputImage()->TransformPhysicalPointToContinuousIndex(transformedPoint,transformedIndex);
87 
88  if( this->m_Interpolator->IsInsideBuffer( transformedIndex ) )
89  {
90  movingValue = this->m_Interpolator->EvaluateAtContinuousIndex( transformedIndex );
91 
92  if (!isZero(movingValue))
93  {
94  currentMovingValue->SetModelVector(movingValue);
95 
96  if (this->GetModelRotation() != Superclass::NONE)
97  currentMovingValue->Reorient(this->m_OrientationMatrix, (this->GetModelRotation() == Superclass::PPD));
98 
99  // Now compute actual measure, depends on model compartment types
100  if (tensorCompatibilityCondition)
101  measure += this->ComputeTensorBasedMetricPart(i,currentMovingValue);
102  else
103  measure += this->ComputeNonTensorBasedMetricPart(i,currentMovingValue);
104  }
105  else
106  {
107  if (tensorCompatibilityCondition)
108  measure += this->ComputeTensorBasedMetricPart(i,m_ZeroDiffusionModel);
109  else
110  measure += this->ComputeNonTensorBasedMetricPart(i,m_ZeroDiffusionModel);
111  }
112  }
113  else
114  {
115  if (tensorCompatibilityCondition)
116  measure += this->ComputeTensorBasedMetricPart(i,m_ZeroDiffusionModel);
117  else
118  measure += this->ComputeNonTensorBasedMetricPart(i,m_ZeroDiffusionModel);
119  }
120 
121  }
122 
123  if (measure <= 0)
124  measure = 0;
125 
126  measure /= this->m_NumberOfPixelsCounted;
127  return measure;
128 }
129 
130 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
131 double
133 ::ComputeTensorBasedMetricPart(unsigned int index, const MCModelPointer &movingValue) const
134 {
135  unsigned int fixedNumCompartments = m_FixedImageValues[index]->GetNumberOfCompartments();
136  unsigned int movingNumCompartments = movingValue->GetNumberOfCompartments();
137 
138  typedef itk::VariableLengthVector <double> LogVectorType;
139  std::vector <LogVectorType> fixedLogVectors(fixedNumCompartments);
140  std::vector <LogVectorType> movingLogVectors(movingNumCompartments);
141  std::vector <double> fixedWeights(fixedNumCompartments);
142  std::vector <double> movingWeights(movingNumCompartments);
143  vnl_matrix <double> workLogMatrix(3,3);
144 
145  unsigned int pos = 0;
146  for (unsigned int i = 0;i < fixedNumCompartments;++i)
147  {
148  if (m_FixedImageValues[index]->GetCompartmentWeight(i) == 0)
149  continue;
150 
151  m_leCalculator->GetTensorLogarithm(m_FixedImageValues[index]->GetCompartment(i)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),
152  workLogMatrix);
153  anima::GetVectorRepresentation(workLogMatrix,fixedLogVectors[pos],6,true);
154 
155  fixedWeights[pos] = m_FixedImageValues[index]->GetCompartmentWeight(i);
156  ++pos;
157  }
158 
159  fixedNumCompartments = pos;
160 
161  if (fixedNumCompartments == 0)
162  return 0;
163  fixedLogVectors.resize(fixedNumCompartments);
164  fixedWeights.resize(fixedNumCompartments);
165 
166  pos = 0;
167  for (unsigned int i = 0;i < movingNumCompartments;++i)
168  {
169  if (movingValue->GetCompartmentWeight(i) == 0)
170  continue;
171 
172  m_leCalculator->GetTensorLogarithm(movingValue->GetCompartment(i)->GetDiffusionTensor().GetVnlMatrix().as_matrix(),
173  workLogMatrix);
174  anima::GetVectorRepresentation(workLogMatrix,movingLogVectors[pos],6,true);
175 
176  movingWeights[pos] = movingValue->GetCompartmentWeight(i);
177  ++pos;
178  }
179 
180  movingNumCompartments = pos;
181 
182  if (movingNumCompartments == 0)
183  return 0;
184  movingLogVectors.resize(movingNumCompartments);
185  movingWeights.resize(movingNumCompartments);
186 
187  double bestMetricValue = -1;
188 
189  unsigned int minCompartmentsNumber = fixedNumCompartments;
190  int rest = movingNumCompartments - minCompartmentsNumber;
191  bool fixedMin = true;
192  if (rest < 0)
193  {
194  fixedMin = false;
195  minCompartmentsNumber = movingNumCompartments;
196  rest *= -1;
197  }
198 
199  unsigned int maxCompartmentsNumber = minCompartmentsNumber + rest;
200  unsigned int totalNumPairingVectors = 1;
201  if (!m_OneToOneMapping)
202  totalNumPairingVectors = boost::math::factorial<double>(maxCompartmentsNumber-1)
203  / (boost::math::factorial<double>(minCompartmentsNumber-1)
204  * boost::math::factorial<double>(maxCompartmentsNumber - minCompartmentsNumber));
205 
206  std::vector < std::vector <unsigned int> > numPairingsVectors(totalNumPairingVectors);
207 
208  if (!m_OneToOneMapping)
209  {
210  std::vector <unsigned int> initialPairingsNumber(maxCompartmentsNumber-1,0);
211  std::vector <unsigned int> pairing(minCompartmentsNumber,1);
212  for (unsigned int i = minCompartmentsNumber-1;i < maxCompartmentsNumber-1;++i)
213  initialPairingsNumber[i] = 1;
214 
215  unsigned int countVector = 0;
216  do
217  {
218  unsigned int pos = 0;
219  std::fill(pairing.begin(),pairing.end(),1);
220  for (unsigned int i = 0;i < maxCompartmentsNumber-1;++i)
221  {
222  if (initialPairingsNumber[i] == 0)
223  {
224  ++pos;
225  continue;
226  }
227 
228  ++pairing[pos];
229  }
230 
231  numPairingsVectors[countVector] = pairing;
232  ++countVector;
233  } while (std::next_permutation(initialPairingsNumber.begin(),initialPairingsNumber.end()));
234  }
235  else
236  {
237  unsigned int numCompartmentUsed = minCompartmentsNumber;
238  if (rest > 0)
239  ++numCompartmentUsed;
240  std::vector <unsigned int> initialPairingsNumber(numCompartmentUsed,1);
241  if (rest > 0)
242  initialPairingsNumber[minCompartmentsNumber] = rest;
243  numPairingsVectors[0] = initialPairingsNumber;
244  }
245 
246  // Loop on all possible numbers of pairings
247  std::vector <unsigned int> currentPermutation(maxCompartmentsNumber);
248 
249  for (unsigned int l = 0;l < totalNumPairingVectors;++l)
250  {
251  currentPermutation.resize(maxCompartmentsNumber);
252 
253  pos = 0;
254  for (unsigned int i = 0;i < numPairingsVectors[l].size();++i)
255  for (unsigned int j = 0;j < numPairingsVectors[l][i];++j)
256  {
257  currentPermutation[pos] = i;
258  ++pos;
259  }
260 
261  do
262  {
263  double metricValue = 0;
264 
265  for (unsigned int j = 0;j < maxCompartmentsNumber;++j)
266  {
267  unsigned int firstIndex;
268  unsigned int secondIndex;
269 
270  if (fixedMin)
271  {
272  firstIndex = currentPermutation[j];
273  secondIndex = j;
274  }
275  else
276  {
277  firstIndex = j;
278  secondIndex = currentPermutation[j];
279  }
280 
281  if ((firstIndex >= fixedLogVectors.size())||(secondIndex >= movingLogVectors.size()))
282  continue;
283 
284  double dist = 0;
285  for (unsigned int k = 0;k < fixedLogVectors[firstIndex].GetSize();++k)
286  dist += (fixedLogVectors[firstIndex][k] - movingLogVectors[secondIndex][k]) * (fixedLogVectors[firstIndex][k] - movingLogVectors[secondIndex][k]);
287 
288  if (!m_OneToOneMapping)
289  dist /= numPairingsVectors[l][currentPermutation[j]];
290 
291  metricValue += fixedWeights[firstIndex] * movingWeights[secondIndex] * dist;
292  }
293 
294  if ((metricValue < bestMetricValue)||(bestMetricValue < 0))
295  bestMetricValue = metricValue;
296  } while(std::next_permutation(currentPermutation.begin(),currentPermutation.end()));
297  }
298 
299  return bestMetricValue;
300 }
301 
302 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
303 double
305 ::ComputeNonTensorBasedMetricPart(unsigned int index, const MCModelPointer &movingValue) const
306 {
307  itkExceptionMacro("DDI basic metric not implemented yet");
308  return 0;
309 }
310 
311 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
312 bool
314 ::isZero(PixelType &vector) const
315 {
316  unsigned int ndim = vector.GetSize();
317 
318  for (unsigned int i = 0;i < ndim;++i)
319  {
320  if (vector[i] != 0)
321  return false;
322  }
323 
324  return true;
325 }
326 
327 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
328 void
331 {
332  if(!this->m_FixedImage)
333  itkExceptionMacro( << "Fixed image has not been assigned" );
334 
335  FixedImageType *fixedImage = const_cast <FixedImageType *> (this->GetFixedImage());
336 
337  this->m_NumberOfPixelsCounted = this->GetFixedImageRegion().GetNumberOfPixels();
338  typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
339  FixedIteratorType ti(fixedImage, this->GetFixedImageRegion());
340  typename FixedImageType::IndexType index;
341 
342  m_FixedImagePoints.resize(this->m_NumberOfPixelsCounted);
343  m_FixedImageValues.resize(this->m_NumberOfPixelsCounted);
344 
345  InputPointType inputPoint;
346 
347  unsigned int pos = 0;
348  PixelType fixedValue;
349  while(!ti.IsAtEnd())
350  {
351  index = ti.GetIndex();
352  fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
353 
354  m_FixedImagePoints[pos] = inputPoint;
355  fixedValue = ti.Get();
356 
357  if (!isZero(fixedValue))
358  {
359  m_FixedImageValues[pos] = fixedImage->GetDescriptionModel()->Clone();
360  m_FixedImageValues[pos]->SetModelVector(fixedValue);
361  }
362  else
363  {
364  m_FixedImageValues[pos] = m_ZeroDiffusionModel->Clone();
365  m_FixedImageValues[pos]->SetModelVector(m_ZeroDiffusionVector);
366  }
367 
368  ++ti;
369  ++pos;
370  }
371 }
372 
373 } // end namespace anima
void GetVectorRepresentation(const vnl_matrix< T1 > &tensor, itk::VariableLengthVector< T2 > &vector, unsigned int vecDim=0, bool scale=false)
itk::ContinuousIndex< double, ImageDimension > ContinuousIndexType
double ComputeNonTensorBasedMetricPart(unsigned int index, const MCModelPointer &movingValue) const
double ComputeTensorBasedMetricPart(unsigned int index, const MCModelPointer &movingValue) const
Really this class is some simplified factory that creates the MCM that it knows.
MeasureType GetValue(const TransformParametersType &parameters) const ITK_OVERRIDE