]> Creatis software - clitk.git/blob - registration/clitkNormalizedCorrelationImageToImageMetric.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkNormalizedCorrelationImageToImageMetric.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 #ifndef _clitkNormalizedCorrelationImageToImageMetric_txx
20 #define _clitkNormalizedCorrelationImageToImageMetric_txx
21
22 // First make sure that the configuration is available.
23 // This line can be removed once the optimized versions
24 // gets integrated into the main directories.
25 #include "itkConfigure.h"
26
27 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
28 #include "clitkOptNormalizedCorrelationImageToImageMetric.txx"
29 #else
30
31
32 #include "clitkNormalizedCorrelationImageToImageMetric.h"
33
34 #include "itkImageRegionConstIteratorWithIndex.h"
35
36
37 namespace clitk
38 {
39
40 /*
41  * Constructor
42  */
43 template <class TFixedImage, class TMovingImage>
44 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
45 ::NormalizedCorrelationImageToImageMetric()
46 {
47   m_SubtractMean = false;
48 }
49
50 /*
51  * Get the match Measure
52  */
53 template <class TFixedImage, class TMovingImage>
54 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
55 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
56 ::GetValue( const TransformParametersType & parameters ) const
57 {
58
59   FixedImageConstPointer fixedImage = this->m_FixedImage;
60
61   if( !fixedImage ) {
62     itkExceptionMacro( << "Fixed image has not been assigned" );
63   }
64
65   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
66
67   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
68
69   typename FixedImageType::IndexType index;
70
71   MeasureType measure;
72
73   this->m_NumberOfPixelsCounted = 0;
74
75   this->SetTransformParameters( parameters );
76
77   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
78
79   AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero;
80   AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero;
81   AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero;
82   AccumulateType sf  = itk::NumericTraits< AccumulateType >::Zero;
83   AccumulateType sm  = itk::NumericTraits< AccumulateType >::Zero;
84
85   while(!ti.IsAtEnd()) {
86
87     index = ti.GetIndex();
88
89     InputPointType inputPoint;
90     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
91
92     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
93       ++ti;
94       continue;
95     }
96
97     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
98
99     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
100       ++ti;
101       continue;
102     }
103
104     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
105       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
106       const RealType fixedValue   = ti.Get();
107       sff += fixedValue  * fixedValue;
108       smm += movingValue * movingValue;
109       sfm += fixedValue  * movingValue;
110       if ( this->m_SubtractMean ) {
111         sf += fixedValue;
112         sm += movingValue;
113       }
114       this->m_NumberOfPixelsCounted++;
115     }
116
117     ++ti;
118   }
119
120   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
121     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
122     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
123     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
124   }
125
126   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
127
128   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
129     measure = sfm / denom;
130   } else {
131     measure = itk::NumericTraits< MeasureType >::Zero;
132   }
133
134   return measure;
135
136 }
137
138
139
140
141
142 /*
143  * Get the Derivative Measure
144  */
145 template < class TFixedImage, class TMovingImage>
146 void
147 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
148 ::GetDerivative( const TransformParametersType & parameters,
149                  DerivativeType & derivative ) const
150 {
151
152   if( !this->GetGradientImage() ) {
153     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
154   }
155
156   FixedImageConstPointer fixedImage = this->m_FixedImage;
157
158   if( !fixedImage ) {
159     itkExceptionMacro( << "Fixed image has not been assigned" );
160   }
161
162   const unsigned int dimension = FixedImageType::ImageDimension;
163
164   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
165
166   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
167
168   typename FixedImageType::IndexType index;
169
170   this->m_NumberOfPixelsCounted = 0;
171
172   this->SetTransformParameters( parameters );
173
174   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
175
176   AccumulateType sff  = itk::NumericTraits< AccumulateType >::Zero;
177   AccumulateType smm  = itk::NumericTraits< AccumulateType >::Zero;
178   AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero;
179   AccumulateType sf  = itk::NumericTraits< AccumulateType >::Zero;
180   AccumulateType sm  = itk::NumericTraits< AccumulateType >::Zero;
181
182   const unsigned int ParametersDimension = this->GetNumberOfParameters();
183   derivative = DerivativeType( ParametersDimension );
184   derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
185
186   DerivativeType derivativeF = DerivativeType( ParametersDimension );
187   derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
188
189   DerivativeType derivativeM = DerivativeType( ParametersDimension );
190   derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
191
192   ti.GoToBegin();
193   // First compute the sums
194   while(!ti.IsAtEnd()) {
195
196     index = ti.GetIndex();
197
198     InputPointType inputPoint;
199     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
200
201     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
202       ++ti;
203       continue;
204     }
205
206     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
207
208     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
209       ++ti;
210       continue;
211     }
212
213     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
214       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
215       const RealType fixedValue   = ti.Get();
216       sff += fixedValue  * fixedValue;
217       smm += movingValue * movingValue;
218       sfm += fixedValue  * movingValue;
219       if ( this->m_SubtractMean ) {
220         sf += fixedValue;
221         sm += movingValue;
222       }
223       this->m_NumberOfPixelsCounted++;
224     }
225
226     ++ti;
227   }
228
229   // Compute contributions to derivatives
230   ti.GoToBegin();
231   while(!ti.IsAtEnd()) {
232
233     index = ti.GetIndex();
234
235     InputPointType inputPoint;
236     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
237
238     if ( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
239       ++ti;
240       continue;
241     }
242
243     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
244
245     if ( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
246       ++ti;
247       continue;
248     }
249
250     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
251       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
252       const RealType fixedValue   = ti.Get();
253
254 #if ITK_VERSION_MAJOR >= 4
255       TransformJacobianType jacobian;
256       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
257 #else
258       const TransformJacobianType & jacobian =
259         this->m_Transform->GetJacobian( inputPoint );
260 #endif
261
262       // Get the gradient by NearestNeighboorInterpolation:
263       // which is equivalent to round up the point components.
264       typedef typename OutputPointType::CoordRepType CoordRepType;
265       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
266       MovingImageContinuousIndexType;
267
268       MovingImageContinuousIndexType tempIndex;
269       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
270
271       typename MovingImageType::IndexType mappedIndex;
272       mappedIndex.CopyWithRound( tempIndex );
273
274       const GradientPixelType gradient =
275         this->GetGradientImage()->GetPixel( mappedIndex );
276
277       for(unsigned int par=0; par<ParametersDimension; par++) {
278         RealType sumF = itk::NumericTraits< RealType >::Zero;
279         RealType sumM = itk::NumericTraits< RealType >::Zero;
280         for(unsigned int dim=0; dim<dimension; dim++) {
281           const RealType differential = jacobian( dim, par ) * gradient[dim];
282           sumF += fixedValue  * differential;
283           sumM += movingValue * differential;
284           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
285             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
286             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
287           }
288         }
289         derivativeF[par] += sumF;
290         derivativeM[par] += sumM;
291       }
292     }
293
294     ++ti;
295   }
296
297   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
298     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
299     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
300     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
301   }
302
303   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
304
305   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
306     for(unsigned int i=0; i<ParametersDimension; i++) {
307       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
308     }
309   } else {
310     for(unsigned int i=0; i<ParametersDimension; i++) {
311       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
312     }
313   }
314
315 }
316
317
318 /*
319  * Get both the match Measure and theDerivative Measure
320  */
321 template <class TFixedImage, class TMovingImage>
322 void
323 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
324 ::GetValueAndDerivative(const TransformParametersType & parameters,
325                         MeasureType & value, DerivativeType  & derivative) const
326 {
327
328
329   if( !this->GetGradientImage() ) {
330     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
331   }
332
333   FixedImageConstPointer fixedImage = this->m_FixedImage;
334
335   if( !fixedImage ) {
336     itkExceptionMacro( << "Fixed image has not been assigned" );
337   }
338
339   const unsigned int dimension = FixedImageType::ImageDimension;
340
341   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
342
343   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
344
345   typename FixedImageType::IndexType index;
346
347   this->m_NumberOfPixelsCounted = 0;
348
349   this->SetTransformParameters( parameters );
350
351   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
352
353   AccumulateType sff  = itk::NumericTraits< AccumulateType >::Zero;
354   AccumulateType smm  = itk::NumericTraits< AccumulateType >::Zero;
355   AccumulateType sfm  = itk::NumericTraits< AccumulateType >::Zero;
356   AccumulateType sf   = itk::NumericTraits< AccumulateType >::Zero;
357   AccumulateType sm   = itk::NumericTraits< AccumulateType >::Zero;
358
359   const unsigned int ParametersDimension = this->GetNumberOfParameters();
360   derivative = DerivativeType( ParametersDimension );
361   derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
362
363   DerivativeType derivativeF = DerivativeType( ParametersDimension );
364   derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
365
366   DerivativeType derivativeM = DerivativeType( ParametersDimension );
367   derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
368
369   DerivativeType derivativeM1 = DerivativeType( ParametersDimension );
370   derivativeM1.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
371
372   ti.GoToBegin();
373   // First compute the sums
374   while(!ti.IsAtEnd()) {
375
376     index = ti.GetIndex();
377
378     InputPointType inputPoint;
379     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
380
381     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
382       ++ti;
383       continue;
384     }
385
386     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
387
388     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
389       ++ti;
390       continue;
391     }
392
393     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
394       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
395       const RealType fixedValue   = ti.Get();
396       sff += fixedValue  * fixedValue;
397       smm += movingValue * movingValue;
398       sfm += fixedValue  * movingValue;
399       if ( this->m_SubtractMean ) {
400         sf += fixedValue;
401         sm += movingValue;
402       }
403       this->m_NumberOfPixelsCounted++;
404     }
405
406     ++ti;
407   }
408
409
410   // Compute contributions to derivatives
411   ti.GoToBegin();
412   while(!ti.IsAtEnd()) {
413
414     index = ti.GetIndex();
415
416     InputPointType inputPoint;
417     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
418
419     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
420       ++ti;
421       continue;
422     }
423
424     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
425
426     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
427       ++ti;
428       continue;
429     }
430
431     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
432       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
433       const RealType fixedValue     = ti.Get();
434
435 #if ITK_VERSION_MAJOR >= 4
436       TransformJacobianType jacobian;
437       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
438 #else
439       const TransformJacobianType & jacobian =
440         this->m_Transform->GetJacobian( inputPoint );
441 #endif
442
443       // Get the gradient by NearestNeighboorInterpolation:
444       // which is equivalent to round up the point components.
445       typedef typename OutputPointType::CoordRepType CoordRepType;
446       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
447       MovingImageContinuousIndexType;
448
449       MovingImageContinuousIndexType tempIndex;
450       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
451
452       typename MovingImageType::IndexType mappedIndex;
453       mappedIndex.CopyWithRound( tempIndex );
454
455       const GradientPixelType gradient =
456         this->GetGradientImage()->GetPixel( mappedIndex );
457
458       for(unsigned int par=0; par<ParametersDimension; par++) {
459         RealType sumF = itk::NumericTraits< RealType >::Zero;
460         RealType sumM = itk::NumericTraits< RealType >::Zero;
461         for(unsigned int dim=0; dim<dimension; dim++) {
462           const RealType differential = jacobian( dim, par ) * gradient[dim];
463           sumF += fixedValue  * differential;
464           sumM += movingValue * differential;
465           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
466             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
467             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
468           }
469         }
470         derivativeF[par] += sumF;
471         derivativeM[par] += sumM;
472       }
473     }
474     ++ti;
475   }
476
477   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
478     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
479     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
480     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
481   }
482
483   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
484
485   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
486     for(unsigned int i=0; i<ParametersDimension; i++) {
487       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
488     }
489     value = sfm / denom;
490   } else {
491     for(unsigned int i=0; i<ParametersDimension; i++) {
492       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
493     }
494     value = itk::NumericTraits< MeasureType >::Zero;
495   }
496
497
498
499 }
500
501 template < class TFixedImage, class TMovingImage>
502 void
503 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
504 ::PrintSelf(std::ostream& os, itk::Indent indent) const
505 {
506   Superclass::PrintSelf(os, indent);
507   os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
508 }
509
510 } // end namespace itk
511
512
513 #endif // opt
514
515 #endif // _clitkNormalizedCorrelationImageToImageMetric.txx