Orfeo Toolbox  4.0
otbNCCRegistrationFunction.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12  Some parts of this code are derived from ITK. See ITKCopyright.txt
13  for details.
14 
15 
16  This software is distributed WITHOUT ANY WARRANTY; without even
17  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18  PURPOSE. See the above copyright notices for more information.
19 
20 =========================================================================*/
21 #ifndef __otbNCCRegistrationFunction_txx
22 #define __otbNCCRegistrationFunction_txx
23 
25 
26 #include "otbMacro.h"
27 #include "vnl/vnl_math.h"
28 
29 namespace otb
30 {
31 
32 /*
33  * Default constructor
34  */
35 template <class TFixedImage, class TMovingImage, class TDisplacementField>
38 {
39 
40  RadiusType r;
41  unsigned int j;
42  for (j = 0; j < ImageDimension; ++j)
43  {
44  r[j] = 1;
45  }
46  this->SetRadius(r);
47  m_MetricTotal = 0.0;
48 
49  m_TimeStep = 1.0;
50  m_DenominatorThreshold = 1e-9;
51  m_IntensityDifferenceThreshold = 0.001;
52  this->SetMovingImage(NULL);
53  this->SetFixedImage(NULL);
54  m_FixedImageSpacing.Fill(1.0);
55  m_FixedImageOrigin.Fill(0.0);
56  m_FixedImageGradientCalculator = GradientCalculatorType::New();
57 
58  typename DefaultInterpolatorType::Pointer interp =
59  DefaultInterpolatorType::New();
60 
61  m_MovingImageInterpolator = static_cast<InterpolatorType*>(
62  interp.GetPointer());
63 
64 }
65 
66 /*
67  * Standard "PrintSelf" method.
68  */
69 template <class TFixedImage, class TMovingImage, class TDisplacementField>
70 void
72 ::PrintSelf(std::ostream& os, itk::Indent indent) const
73 {
74 
75  Superclass::PrintSelf(os, indent);
76  /*
77  os << indent << "MovingImageIterpolator: ";
78  os << m_MovingImageInterpolator.GetPointer() << std::endl;
79  os << indent << "FixedImageGradientCalculator: ";
80  os << m_FixedImageGradientCalculator.GetPointer() << std::endl;
81  os << indent << "DenominatorThreshold: ";
82  os << m_DenominatorThreshold << std::endl;
83  os << indent << "IntensityDifferenceThreshold: ";
84  os << m_IntensityDifferenceThreshold << std::endl;
85  */
86 }
87 
88 /*
89  * Set the function state values before each iteration
90  */
91 template <class TFixedImage, class TMovingImage, class TDisplacementField>
92 void
95 {
96  if (!this->m_MovingImage || !this->m_FixedImage || !m_MovingImageInterpolator)
97  {
98  itkExceptionMacro(<< "MovingImage, FixedImage and/or Interpolator not set");
99  }
100 
101  // cache fixed image information
102  m_FixedImageSpacing = this->m_FixedImage->GetSpacing();
103  m_FixedImageOrigin = this->m_FixedImage->GetOrigin();
104 
105  // setup gradient calculator
106  m_FixedImageGradientCalculator->SetInputImage(this->m_FixedImage);
107 
108  // setup moving image interpolator
109  m_MovingImageInterpolator->SetInputImage(this->m_MovingImage);
110 
111  otbMsgDevMacro( " total metric " << m_MetricTotal << " field size " <<
112  this->GetDisplacementField()->GetLargestPossibleRegion().GetSize() << " image size " <<
113  this->m_FixedImage->GetLargestPossibleRegion().GetSize() );
114  m_MetricTotal = 0.0;
115 
116 }
117 
118 /*
119  * Compute update at a non boundary neighbourhood
120  */
121 template <class TFixedImage, class TMovingImage, class TDisplacementField>
123 ::PixelType
125 ::ComputeUpdate(const NeighborhoodType& it, void * itkNotUsed(globalData),
126  const FloatOffsetType& itkNotUsed(offset))
127 {
128 
129  PixelType update;
130  update.Fill(0.0);
131  unsigned int j;
132 
133  IndexType oindex = it.GetIndex();
134 
135  typename FixedImageType::SizeType hradius = it.GetRadius();
136 
137  FixedImageType* img = const_cast<FixedImageType *>(this->m_FixedImage.GetPointer());
138 
139  typename FixedImageType::SizeType imagesize = img->GetLargestPossibleRegion().GetSize();
140 
142  hoodIt(hradius, img, img->GetRequestedRegion());
143  hoodIt.SetLocation(oindex);
144 
145  double sff = 0.0;
146  double smm = 0.0;
147  double sfm = 0.0;
148  double fixedValue;
149  double movingValue;
150 
151  double derivativeF[ImageDimension];
152  double derivativeM[ImageDimension];
153  for (j = 0; j < ImageDimension; ++j)
154  {
155  derivativeF[j] = 0;
156  derivativeM[j] = 0;
157  }
158 
159  unsigned int indct;
160  unsigned int hoodlen = hoodIt.Size();
161 
162  for (indct = 0; indct < hoodlen - 1; indct++)
163  {
164 
165  IndexType index = hoodIt.GetIndex(indct);
166 
167  bool inimage = true;
168  for (unsigned int dd = 0; dd < ImageDimension; dd++)
169  {
170  if (index[dd] < 0 || index[dd] >
171  static_cast<typename IndexType::IndexValueType>(imagesize[dd] - 1)) inimage = false;
172  }
173 
174  if (inimage)
175  {
176 
177  // Get fixed image related information
178 
179  CovariantVectorType fixedGradient;
180  double fixedGradientSquaredMagnitude = 0;
181 
182  // Note: no need to check the index is within
183  // fixed image buffer. This is done by the external filter.
184  fixedValue = (double) this->m_FixedImage->GetPixel(index);
185 
186  fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex(index);
187  for (j = 0; j < ImageDimension; ++j)
188  {
189  fixedGradientSquaredMagnitude += vnl_math_sqr(fixedGradient[j]) * m_FixedImageSpacing[j];
190  }
191 
192  // Get moving image related information
193 
194  PointType mappedPoint;
195 
196  typedef typename TDisplacementField::PixelType DisplacementPixelType;
197 
198  // Edited by OTB developers
199  DisplacementPixelType vec;
200  vec.Fill(0);
201  if (this->GetDisplacementField()->GetBufferedRegion().IsInside(index))
202  {
203  vec = this->GetDisplacementField()->GetPixel(index);
204  }
205  // End Edited by OTB developers
206 
207  for (j = 0; j < ImageDimension; ++j)
208  {
209  mappedPoint[j] = double(index[j]) * m_FixedImageSpacing[j] +
210  m_FixedImageOrigin[j];
211  mappedPoint[j] += vec[j];
212  }
213  if (m_MovingImageInterpolator->IsInsideBuffer(mappedPoint))
214  {
215  movingValue = m_MovingImageInterpolator->Evaluate(mappedPoint);
216  }
217  else
218  {
219  movingValue = 0.0;
220  }
221 
222  sff += fixedValue * fixedValue;
223  smm += movingValue * movingValue;
224  sfm += fixedValue * movingValue;
225 
226  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
227  {
228  double differential = fixedGradient[dim];
229  derivativeF[dim] += fixedValue * differential;
230  derivativeM[dim] += movingValue * differential;
231  }
232  }
233 
234  }
235 
236  double updatenorm = 0.0;
237  if ((sff * smm) != 0.0)
238  {
239  double factor = 1.0 / vcl_sqrt(sff * smm);
240  for (unsigned int i = 0; i < ImageDimension; ++i)
241  {
242  update[i] = factor * (derivativeF[i] - (sfm / smm) * derivativeM[i]);
243  updatenorm += (update[i] * update[i]);
244  }
245  updatenorm = vcl_sqrt(updatenorm);
246  m_MetricTotal += sfm * factor;
247  this->m_Energy += sfm * factor;
248  }
249  else
250  {
251  for (unsigned int i = 0; i < ImageDimension; ++i)
252  {
253  update[i] = 0.0;
254  }
255  updatenorm = 1.0;
256  }
257 
258 // if ( fixedValue > 0.40 && movingValue > 0.40) std::cout << " update norm " << updatenorm;
259 
260  if (this->GetNormalizeGradient() && updatenorm != 0.0)
261  {
262  update /= (updatenorm);
263  }
264 
265  return update * this->m_GradientStep;
266 }
267 
268 } // end namespace otbs
269 
270 #endif

Generated at Sat Mar 8 2014 16:10:57 for Orfeo Toolbox with doxygen 1.8.3.1