ANIMA  4.0
animaMCMMeanSquaresImageToImageMetric.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <vnl/vnl_matrix_fixed.h>
5 #include <animaBaseTensorTools.h>
7 #include <itkImageRegionConstIteratorWithIndex.h>
8 
9 namespace anima
10 {
11 
12 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
15 {
16  m_FixedImagePoints.clear();
17  m_FixedImageValues.clear();
18 
20  mcmCreator.SetNumberOfCompartments(0);
21  mcmCreator.SetModelWithStationaryWaterComponent(true);
22  mcmCreator.SetModelWithFreeWaterComponent(false);
23 
24  m_ZeroDiffusionModel = mcmCreator.GetNewMultiCompartmentModel();
25  m_ZeroDiffusionVector = m_ZeroDiffusionModel->GetModelVector();
26 
27  m_L2DistanceComputer = anima::MCML2DistanceComputer::New();
28 }
29 
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  itkExceptionMacro("Fixed image has not been assigned");
39 
40  if (this->m_NumberOfPixelsCounted == 0)
41  return 0;
42 
43  this->SetTransformParameters(parameters);
44 
45  PixelType movingValue;
46 
47  OutputPointType transformedPoint;
48  ContinuousIndexType transformedIndex;
49 
50  MovingImageType *movingImage = const_cast <MovingImageType *> (this->GetMovingImage());
51  MCModelPointer currentMovingValue = movingImage->GetDescriptionModel()->Clone();
52 
53  double measure = 0;
54 
55  for (unsigned int i = 0;i < this->m_NumberOfPixelsCounted;++i)
56  {
57  transformedPoint = this->m_Transform->TransformPoint(m_FixedImagePoints[i]);
58  this->m_Interpolator->GetInputImage()->TransformPhysicalPointToContinuousIndex(transformedPoint,transformedIndex);
59 
60  if( this->m_Interpolator->IsInsideBuffer(transformedIndex))
61  {
62  movingValue = this->m_Interpolator->EvaluateAtContinuousIndex(transformedIndex);
63 
64  if (!isZero(movingValue))
65  {
66  currentMovingValue->SetModelVector(movingValue);
67 
68  if (this->GetModelRotation() != Superclass::NONE)
69  currentMovingValue->Reorient(this->m_OrientationMatrix, (this->GetModelRotation() == Superclass::PPD));
70 
71  // Now compute actual measure, depends on model compartment types
72  measure += m_L2DistanceComputer->ComputeDistance(m_FixedImageValues[i],currentMovingValue);
73  }
74  else
75  measure += m_L2DistanceComputer->ComputeDistance(m_FixedImageValues[i],m_ZeroDiffusionModel);
76  }
77  else
78  measure += m_L2DistanceComputer->ComputeDistance(m_FixedImageValues[i],m_ZeroDiffusionModel);
79  }
80 
81  measure /= this->m_NumberOfPixelsCounted;
82  return measure;
83 }
84 
85 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
86 bool
88 ::isZero(PixelType &vector) const
89 {
90  unsigned int ndim = vector.GetSize();
91 
92  for (unsigned int i = 0;i < ndim;++i)
93  {
94  if (vector[i] != 0)
95  return false;
96  }
97 
98  return true;
99 }
100 
101 template < class TFixedImagePixelType, class TMovingImagePixelType, unsigned int ImageDimension >
102 void
105 {
106  if(!this->m_FixedImage)
107  itkExceptionMacro( << "Fixed image has not been assigned" );
108 
109  FixedImageType *fixedImage = const_cast <FixedImageType *> (this->GetFixedImage());
110  this->m_NumberOfPixelsCounted = this->GetFixedImageRegion().GetNumberOfPixels();
111  typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
112  FixedIteratorType ti(fixedImage, this->GetFixedImageRegion());
113  typename FixedImageType::IndexType index;
114 
115  m_FixedImagePoints.resize(this->m_NumberOfPixelsCounted);
116  m_FixedImageValues.resize(this->m_NumberOfPixelsCounted);
117 
118  InputPointType inputPoint;
119 
120  unsigned int pos = 0;
121  PixelType fixedValue;
122  while(!ti.IsAtEnd())
123  {
124  index = ti.GetIndex();
125  fixedImage->TransformIndexToPhysicalPoint(index, inputPoint);
126 
127  m_FixedImagePoints[pos] = inputPoint;
128  fixedValue = ti.Get();
129 
130  if (!isZero(fixedValue))
131  {
132  m_FixedImageValues[pos] = fixedImage->GetDescriptionModel()->Clone();
133  m_FixedImageValues[pos]->SetModelVector(fixedValue);
134  }
135  else
136  {
137  m_FixedImageValues[pos] = m_ZeroDiffusionModel->Clone();
138  m_FixedImageValues[pos]->SetModelVector(m_ZeroDiffusionVector);
139  }
140 
141  ++ti;
142  ++pos;
143  }
144 }
145 
146 } // end namespace anima
MeasureType GetValue(const TransformParametersType &parameters) const ITK_OVERRIDE
Superclass::TransformParametersType TransformParametersType
itk::ContinuousIndex< double, ImageDimension > ContinuousIndexType
Really this class is some simplified factory that creates the MCM that it knows.