ANIMA  4.0
animaBaseAffineBlockMatcher.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 /* Transforms */
5 #include <itkTranslationTransform.h>
9 
10 namespace anima
11 {
12 
13 template <typename TInputImageType>
16 {
17  m_BlockTransformType = Translation;
18 
19  m_AngleMax = M_PI;
20  m_TranslateMax = 10;
21  m_ScaleMax = 3;
22 
23  m_SearchAngleRadius = 5;
24  m_SearchScaleRadius = 0.1;
25 
26  m_AffineDirection = 1;
27 }
28 
29 template <typename TInputImageType>
33 {
34  switch (m_BlockTransformType)
35  {
36  case Translation:
37  return AgregatorType::TRANSLATION;
38 
39  case Rigid:
40  return AgregatorType::RIGID;
41 
42  case Affine:
43  case Directional_Affine:
44  default:
45  return AgregatorType::AFFINE;
46  }
47 }
48 
49 template <typename TInputImageType>
53 {
54  BaseInputTransformPointer outputValue;
55 
56  switch(m_BlockTransformType)
57  {
58  case Translation:
59  {
60  typedef itk::TranslationTransform <double, InputImageType::ImageDimension> itkTransformType;
61  typename itkTransformType::Pointer tr = itkTransformType::New();
62  tr->SetIdentity();
63  outputValue = tr;
64 
65  break;
66  }
67 
68  case Rigid:
69  {
70  if (InputImageType::ImageDimension == 3)
71  {
72  typedef anima::LogRigid3DTransform<double> itkTransformType;
73  typename itkTransformType::Pointer tr = itkTransformType::New();
74  tr->SetIdentity();
75  tr->SetCenter(blockCenter);
76 
77  outputValue = tr;
78  }
79  else
80  std::cerr << "Only Rigid 3D transforms handled right now." << std::endl;
81  break;
82  }
83 
84  case Affine:
85  {
86  typedef anima::SplitAffine3DTransform <double> itkTransformType;
87  typename itkTransformType::Pointer tr = itkTransformType::New();
88  tr->SetIdentity();
89  tr->SetCenter(blockCenter);
90 
91  outputValue = tr;
92  break;
93  }
94 
95  case Directional_Affine:
96  default:
97  {
98  typedef anima::DirectionTransform <double> itkTransformType;
99  typename itkTransformType::Pointer tr = itkTransformType::New();
100  tr->SetIdentity();
101 
102  typename itkTransformType::HomogeneousMatrixType geometry;
103 
104  geometry.SetIdentity();
105  for (unsigned int i = 0;i < 3;++i)
106  for (unsigned int j = 0;j < 3;++j)
107  geometry(i,j) = this->GetReferenceImage()->GetDirection()(i,j) * this->GetReferenceImage()->GetSpacing()[j];
108 
109  for (unsigned int j = 0;j < 3;++j)
110  geometry(j,TInputImageType::ImageDimension) = blockCenter[j];
111 
112  tr->SetUniqueDirection(m_AffineDirection);
113  tr->SetGeometry(geometry, false);
114 
115  outputValue = tr;
116  break;
117  }
118  }
119 
120  return outputValue;
121 }
122 
123 template <typename TInputImageType>
124 void
126 ::BlockMatchingSetup(MetricPointer &metric, unsigned int block)
127 {
128  switch (m_BlockTransformType)
129  {
130  case Translation:
131  {
132  typedef itk::TranslationTransform <double, InputImageType::ImageDimension> itkTransformType;
133  itkTransformType *tr = dynamic_cast <itkTransformType *> (this->GetBlockTransformPointer(block).GetPointer());
134  tr->SetIdentity();
135  break;
136  }
137 
138  default:
139  {
140  typedef itk::MatrixOffsetTransformBase <double, InputImageType::ImageDimension> itkTransformType;
141  itkTransformType *tr = dynamic_cast <itkTransformType *> (this->GetBlockTransformPointer(block).GetPointer());
142  tr->SetIdentity();
143  break;
144  }
145  }
146 }
147 
148 template <typename TInputImageType>
149 void
152 {
153  if (this->GetOptimizerType() == Superclass::Exhaustive)
154  return;
155 
156  typedef anima::BobyqaOptimizer LocalOptimizerType;
157  LocalOptimizerType::ScalesType tmpScales(this->GetBlockTransformPointer(0)->GetNumberOfParameters());
158  LocalOptimizerType::ScalesType lowerBounds(this->GetBlockTransformPointer(0)->GetNumberOfParameters());
159  LocalOptimizerType::ScalesType upperBounds(this->GetBlockTransformPointer(0)->GetNumberOfParameters());
160  typename InputImageType::SpacingType fixedSpacing = this->GetReferenceImage()->GetSpacing();
161 
162  switch (m_BlockTransformType)
163  {
164  case Translation:
165  {
166  for (unsigned int i = 0;i < InputImageType::ImageDimension;++i)
167  {
168  tmpScales[i] = 1.0 / fixedSpacing[i];
169 
170  lowerBounds[i] = - m_TranslateMax * fixedSpacing[i];
171  upperBounds[i] = m_TranslateMax * fixedSpacing[i];
172  }
173 
174  break;
175  }
176 
177  case Rigid:
178  {
179  for (unsigned int i = 0;i < InputImageType::ImageDimension;++i)
180  {
181  tmpScales[i] = this->GetSearchRadius() * 180.0 / (m_SearchAngleRadius * M_PI);
182  lowerBounds[i] = - m_AngleMax * M_PI / 180.0;
183  upperBounds[i] = m_AngleMax * M_PI / 180.0;
184 
185  tmpScales[i+InputImageType::ImageDimension] = 1.0 / fixedSpacing[i];
186 
187  lowerBounds[i+InputImageType::ImageDimension] = - m_TranslateMax * fixedSpacing[i];
188  upperBounds[i+InputImageType::ImageDimension] = m_TranslateMax * fixedSpacing[i];
189  }
190 
191  break;
192  }
193 
194  case Affine:
195  {
196  // There are 12 parameters: 3 angles, 3 translations, 3 scales, 3 angles again
197  for (unsigned int i = 0;i < InputImageType::ImageDimension;++i)
198  {
199  // Angles
200  tmpScales[i] = this->GetSearchRadius() * 180.0 / (m_SearchAngleRadius * M_PI);
201  lowerBounds[i] = - m_AngleMax * M_PI / 180.0;
202  upperBounds[i] = m_AngleMax * M_PI / 180.0;
203 
204  // Translations
205  tmpScales[InputImageType::ImageDimension + i] = 1.0 / fixedSpacing[i];
206 
207  lowerBounds[InputImageType::ImageDimension + i] = - m_TranslateMax * fixedSpacing[i];
208  upperBounds[InputImageType::ImageDimension + i] = m_TranslateMax * fixedSpacing[i];
209 
210  // Scales
211  tmpScales[2 * InputImageType::ImageDimension + i] = this->GetSearchRadius() / m_SearchScaleRadius;
212 
213  lowerBounds[2 * InputImageType::ImageDimension + i] = - std::log (m_ScaleMax);
214  upperBounds[2 * InputImageType::ImageDimension + i] = std::log (m_ScaleMax);
215 
216  // Second angles
217  tmpScales[3 * InputImageType::ImageDimension + i] = this->GetSearchRadius() * 180.0 / (m_SearchAngleRadius * M_PI);
218 
219  lowerBounds[3 * InputImageType::ImageDimension + i] = - m_AngleMax * M_PI / 180.0;
220  upperBounds[3 * InputImageType::ImageDimension + i] = m_AngleMax * M_PI / 180.0;
221  }
222 
223  break;
224  }
225 
226  case Directional_Affine:
227  default:
228  {
229  // There is 1 parameter: 1 translation in voxel coordinates
230  tmpScales[0] = 1.0;
231 
232  lowerBounds[0] = - m_TranslateMax;
233  upperBounds[0] = m_TranslateMax;
234 
235  break;
236  }
237  }
238 
239  LocalOptimizerType * tmpOpt = dynamic_cast <LocalOptimizerType *> (optimizer.GetPointer());
240  tmpOpt->SetScales(tmpScales);
241  tmpOpt->SetLowerBounds(lowerBounds);
242  tmpOpt->SetUpperBounds(upperBounds);
243 }
244 
245 } // end namespace anima
AgregatorType::TRANSFORM_TYPE GetAgregatorInputTransformType()
Superclass::BaseInputTransformPointer BaseInputTransformPointer
Superclass::MetricPointer MetricPointer
Superclass::OptimizerPointer OptimizerPointer
virtual void BlockMatchingSetup(MetricPointer &metric, unsigned int block)
virtual BaseInputTransformPointer GetNewBlockTransform(PointType &blockCenter)
virtual void TransformDependantOptimizerSetup(OptimizerPointer &optimizer)