]> Creatis software - clitk.git/blob - registration/clitkOptNormalizedCorrelationImageToImageMetric.txx
small bugs corrections
[clitk.git] / registration / clitkOptNormalizedCorrelationImageToImageMetric.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://oncora1.lyon.fnclcc.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 __clitkOptNormalizedCorrelationImageToImageMetric_txx
20 #define __clitkOptNormalizedCorrelationImageToImageMetric_txx
21
22 #include "clitkOptNormalizedCorrelationImageToImageMetric.h"
23 #include "itkCovariantVector.h"
24 #include "itkImageRandomConstIteratorWithIndex.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionIterator.h"
27 #include "itkImageIterator.h"
28 #include "vnl/vnl_math.h"
29
30 namespace clitk
31 {
32
33 /**
34  * Constructor
35  */
36 template < class TFixedImage, class TMovingImage >
37 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
38 ::NormalizedCorrelationImageToImageMetric()
39 {
40   this->SetComputeGradient(true);
41
42   m_ThreaderSFF = NULL;
43   m_ThreaderSMM = NULL;
44   m_ThreaderSFM = NULL;
45   m_ThreaderSF = NULL;
46   m_ThreaderSM = NULL;
47   m_ThreaderDerivativeF = NULL;
48   m_ThreaderDerivativeM = NULL;
49   this->m_WithinThreadPreProcess = false;
50   this->m_WithinThreadPostProcess = false;
51
52   //  For backward compatibility, the default behavior is to use all the pixels
53   //  in the fixed image.
54   this->UseAllPixelsOn();
55
56   m_SubtractMean=false;
57
58 }
59
60 template < class TFixedImage, class TMovingImage >
61 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
62 ::~NormalizedCorrelationImageToImageMetric()
63 {
64   if(m_ThreaderSFF != NULL) {
65     delete [] m_ThreaderSFF;
66   }
67   m_ThreaderSFF = NULL;
68
69   if(m_ThreaderSMM != NULL) {
70     delete [] m_ThreaderSMM;
71   }
72   m_ThreaderSMM = NULL;
73
74   if(m_ThreaderSFM != NULL) {
75     delete [] m_ThreaderSFM;
76   }
77   m_ThreaderSFM = NULL;
78
79   if(m_ThreaderSF != NULL) {
80     delete [] m_ThreaderSF;
81   }
82   m_ThreaderSF = NULL;
83
84   if(m_ThreaderSM != NULL) {
85     delete [] m_ThreaderSM;
86   }
87   m_ThreaderSM = NULL;
88
89   if(m_ThreaderDerivativeF != NULL) {
90     delete [] m_ThreaderDerivativeF;
91   }
92   m_ThreaderDerivativeF = NULL;
93
94   if(m_ThreaderDerivativeM != NULL) {
95     delete [] m_ThreaderDerivativeM;
96   }
97   m_ThreaderDerivativeM = NULL;
98 }
99
100 /**
101  * Print out internal information about this class
102  */
103 template < class TFixedImage, class TMovingImage  >
104 void
105 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
106 ::PrintSelf(std::ostream& os, itk::Indent indent) const
107 {
108
109   Superclass::PrintSelf(os, indent);
110
111 }
112
113
114 /**
115  * Initialize
116  */
117 template <class TFixedImage, class TMovingImage>
118 void
119 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
120 ::Initialize(void) throw ( itk::ExceptionObject )
121 {
122
123   this->Superclass::Initialize();
124   this->Superclass::MultiThreadingInitialize();
125
126
127   /**
128    * Allocate memory for the accumulators (set to zero in GetValue)
129    */
130   if(m_ThreaderSFF != NULL) {
131     delete [] m_ThreaderSFF;
132   }
133   m_ThreaderSFF = new double[this->m_NumberOfThreads];
134
135
136   if(m_ThreaderSMM != NULL) {
137     delete [] m_ThreaderSMM;
138   }
139   m_ThreaderSMM = new double[this->m_NumberOfThreads];
140
141   if(m_ThreaderSFM != NULL) {
142     delete [] m_ThreaderSFM;
143   }
144   m_ThreaderSFM = new double[this->m_NumberOfThreads];
145
146   if(this->m_SubtractMean) {
147     if(m_ThreaderSF != NULL) {
148       delete [] m_ThreaderSF;
149     }
150     m_ThreaderSF = new double[this->m_NumberOfThreads];
151
152     if(m_ThreaderSM != NULL) {
153       delete [] m_ThreaderSM;
154     }
155     m_ThreaderSM = new double[this->m_NumberOfThreads];
156   }
157
158   if(m_ThreaderDerivativeF != NULL) {
159     delete [] m_ThreaderDerivativeF;
160   }
161   m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfThreads];
162   for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
163     m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
164   }
165
166   if(m_ThreaderDerivativeM != NULL) {
167     delete [] m_ThreaderDerivativeM;
168   }
169   m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfThreads];
170   for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
171     m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
172   }
173 }
174
175
176 template < class TFixedImage, class TMovingImage  >
177 inline bool
178 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
179 ::GetValueThreadProcessSample(
180   unsigned int threadID,
181   unsigned long fixedImageSample,
182   const MovingImagePointType & itkNotUsed(mappedPoint),
183   double movingImageValue) const
184 {
185   const RealType fixedImageValue= this->m_FixedImageSamples[fixedImageSample].value;
186   m_ThreaderSFF[threadID] += fixedImageValue  * fixedImageValue;
187   m_ThreaderSMM[threadID] += movingImageValue * movingImageValue;
188   m_ThreaderSFM[threadID] += fixedImageValue  * movingImageValue;
189   if ( this->m_SubtractMean ) {
190     m_ThreaderSF[threadID] += fixedImageValue;
191     m_ThreaderSM[threadID] += movingImageValue;
192   }
193
194   return true;
195 }
196
197 template < class TFixedImage, class TMovingImage  >
198 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
199 ::MeasureType
200 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
201 ::GetValue( const ParametersType & parameters ) const
202 {
203   itkDebugMacro("GetValue( " << parameters << " ) ");
204
205   if( !this->m_FixedImage ) {
206     itkExceptionMacro( << "Fixed image has not been assigned" );
207   }
208
209
210   //Reset the accumulators
211   memset( m_ThreaderSFF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
212   memset( m_ThreaderSMM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
213   memset( m_ThreaderSFM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
214   if(this->m_SubtractMean) {
215     memset( m_ThreaderSF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
216     memset( m_ThreaderSM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
217   }
218
219
220   // Set up the parameters in the transform
221   this->m_Transform->SetParameters( parameters );
222   this->m_Parameters = parameters;
223
224   // MUST BE CALLED TO INITIATE PROCESSING
225   this->GetValueMultiThreadedInitiate();
226
227   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
228                  << this->m_NumberOfPixelsCounted << " / "
229                  << this->m_NumberOfFixedImageSamples
230                  << std::endl );
231
232   if( this->m_NumberOfPixelsCounted <
233       this->m_NumberOfFixedImageSamples / 4 ) {
234     itkExceptionMacro( "Too many samples map outside moving image buffer: "
235                        << this->m_NumberOfPixelsCounted << " / "
236                        << this->m_NumberOfFixedImageSamples
237                        << std::endl );
238   }
239
240   // Accumulate the threads
241   AccumulateType sff, smm, sfm, sf, sm;
242   sff = m_ThreaderSFF[0];
243   smm = m_ThreaderSMM[0];
244   sfm = m_ThreaderSFM[0];
245   sf  = m_ThreaderSF[0];
246   sm  = m_ThreaderSM[0];
247
248   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
249     sff +=  m_ThreaderSFF[t];
250     smm +=  m_ThreaderSMM[t];
251     sfm +=  m_ThreaderSFM[t];
252     if ( this->m_SubtractMean ) {
253       sf +=  m_ThreaderSF[t];
254       sm +=  m_ThreaderSM[t];
255     }
256   }
257
258   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
259     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
260     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
261     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
262   }
263
264
265   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
266   MeasureType measure;
267   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
268     measure = sfm / denom;
269   } else {
270     measure = itk::NumericTraits< MeasureType >::Zero;
271   }
272
273
274   return measure;
275 }
276
277
278 template < class TFixedImage, class TMovingImage  >
279 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
280 ::MeasureType
281 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
282 ::ComputeSums( const ParametersType & parameters ) const
283 {
284   //No checking for the fixed image,  done in the caller
285   //Reset the accumulators
286   memset( m_ThreaderSFF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
287   memset( m_ThreaderSMM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
288   memset( m_ThreaderSFM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
289   if(this->m_SubtractMean) {
290     memset( m_ThreaderSF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
291     memset( m_ThreaderSM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
292   }
293
294
295   // Set up the parameters in the transform
296   this->m_Transform->SetParameters( parameters );
297   this->m_Parameters = parameters;
298
299   // MUST BE CALLED TO INITIATE PROCESSING
300   this->GetValueMultiThreadedInitiate();
301
302   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
303                  << this->m_NumberOfPixelsCounted << " / "
304                  << this->m_NumberOfFixedImageSamples
305                  << std::endl );
306
307   if( this->m_NumberOfPixelsCounted <
308       this->m_NumberOfFixedImageSamples / 4 ) {
309     itkExceptionMacro( "Too many samples map outside moving image buffer: "
310                        << this->m_NumberOfPixelsCounted << " / "
311                        << this->m_NumberOfFixedImageSamples
312                        << std::endl );
313   }
314
315   // Accumulate the threads
316   m_SFF = m_ThreaderSFF[0];
317   m_SMM = m_ThreaderSMM[0];
318   m_SFM = m_ThreaderSFM[0];
319   m_SF  = m_ThreaderSF[0];
320   m_SM  = m_ThreaderSM[0];
321
322   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
323     m_SFF +=  m_ThreaderSFF[t];
324     m_SMM +=  m_ThreaderSMM[t];
325     m_SFM +=  m_ThreaderSFM[t];
326     if ( this->m_SubtractMean ) {
327       m_SF +=  m_ThreaderSF[t];
328       m_SM +=  m_ThreaderSM[t];
329     }
330   }
331
332   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
333     m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
334     m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
335     m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
336     m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
337     m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
338   }
339
340
341   m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
342   MeasureType measure;
343   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
344     measure = m_SFM / m_Denom;
345   } else {
346     measure = itk::NumericTraits< MeasureType >::Zero;
347   }
348
349   return measure;
350 }
351
352
353 template < class TFixedImage, class TMovingImage  >
354 inline bool
355 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
356 ::GetValueAndDerivativeThreadProcessSample(
357   unsigned int threadID,
358   unsigned long fixedImageSample,
359   const MovingImagePointType & itkNotUsed(mappedPoint),
360   double movingImageValue,
361   const ImageDerivativesType &
362   movingImageGradientValue
363 ) const
364 {
365
366   const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
367   const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
368
369   // Need to use one of the threader transforms if we're
370   // not in thread 0.
371   //
372   // Use a raw pointer here to avoid the overhead of smart pointers.
373   // For instance, Register and UnRegister have mutex locks around
374   // the reference counts.
375   TransformType* transform;
376
377   if (threadID > 0) {
378     transform = this->m_ThreaderTransform[threadID - 1];
379   } else {
380     transform = this->m_Transform;
381   }
382
383   // Jacobian should be evaluated at the unmapped (fixed image) point.
384   const TransformJacobianType & jacobian = transform
385       ->GetJacobian( fixedImagePoint );
386
387   for(unsigned int par=0; par<this->m_NumberOfParameters; par++) {
388     RealType sumF = itk::NumericTraits< RealType >::Zero;
389     RealType sumM = itk::NumericTraits< RealType >::Zero;
390     for(unsigned int dim=0; dim<MovingImageDimension; dim++) {
391       const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
392       if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
393         sumF += (fixedImageValue-m_FixedMean)  * differential;
394         sumM += (movingImageValue-m_MovingMean) * differential;
395       } else {
396         sumF += differential * fixedImageValue;
397         sumM += differential * movingImageValue;
398       }
399     }
400     m_ThreaderDerivativeF[threadID][par] += sumF;
401     m_ThreaderDerivativeM[threadID][par] += sumM;
402   }
403
404   return true;
405 }
406
407
408 /**
409  * Get the both Value and Derivative Measure
410  */
411 template < class TFixedImage, class TMovingImage  >
412 void
413 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
414 ::GetValueAndDerivative( const ParametersType & parameters,
415                          MeasureType & value,
416                          DerivativeType & derivative) const
417 {
418
419   if( !this->m_FixedImage ) {
420     itkExceptionMacro( << "Fixed image has not been assigned" );
421   }
422
423   // Set up the parameters in the transform
424   this->m_Transform->SetParameters( parameters );
425   this->m_Parameters = parameters;
426
427   //We need the sums and the value to be calculated first
428   value=this->ComputeSums(parameters);
429
430   //Set output values to zero
431   if(derivative.GetSize() != this->m_NumberOfParameters) {
432     derivative = DerivativeType( this->m_NumberOfParameters );
433   }
434   memset( derivative.data_block(),
435           0,
436           this->m_NumberOfParameters * sizeof(ITK_TYPENAME DerivativeType::ValueType) );
437
438   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
439     memset( m_ThreaderDerivativeF[threadID].data_block(),
440             0,
441             this->m_NumberOfParameters * sizeof(ITK_TYPENAME DerivativeType::ValueType) );
442
443     memset( m_ThreaderDerivativeM[threadID].data_block(),
444             0,
445             this->m_NumberOfParameters * sizeof(ITK_TYPENAME DerivativeType::ValueType) );
446   }
447
448   // MUST BE CALLED TO INITIATE PROCESSING
449   this->GetValueAndDerivativeMultiThreadedInitiate();
450
451   // Accumulate over the threads
452   DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
453   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
454     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
455       derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
456       derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
457     }
458   }
459
460   //Compute derivatives
461   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
462     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
463       derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
464     }
465   } else {
466     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
467       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
468     }
469   }
470
471 }
472
473
474 /**
475  * Get the match measure derivative
476  */
477 template < class TFixedImage, class TMovingImage  >
478 void
479 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
480 ::GetDerivative( const ParametersType & parameters,
481                  DerivativeType & derivative ) const
482 {
483   if( !this->m_FixedImage ) {
484     itkExceptionMacro( << "Fixed image has not been assigned" );
485   }
486
487   MeasureType value;
488   // call the combined version
489   this->GetValueAndDerivative( parameters, value, derivative );
490 }
491
492 } // end namespace clitk
493
494
495 #endif