]> Creatis software - clitk.git/blob - registration/clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.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 __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_txx
20 #define __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_txx
21
22 #include "clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
38 ::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD()
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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
62 ::~NormalizedCorrelationImageToImageMetricFor3DBLUTFFD()
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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
199 ::MeasureType
200 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<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 #if ITK_VERSION_MAJOR < 4
223   this->m_Parameters = parameters;
224 #endif
225
226   // MUST BE CALLED TO INITIATE PROCESSING
227   this->GetValueMultiThreadedInitiate();
228
229   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
230                  << this->m_NumberOfPixelsCounted << " / "
231                  << this->m_NumberOfFixedImageSamples
232                  << std::endl );
233
234   if( this->m_NumberOfPixelsCounted <
235       this->m_NumberOfFixedImageSamples / 4 ) {
236     itkExceptionMacro( "Too many samples map outside moving image buffer: "
237                        << this->m_NumberOfPixelsCounted << " / "
238                        << this->m_NumberOfFixedImageSamples
239                        << std::endl );
240   }
241
242   // Accumulate the threads
243   AccumulateType sff, smm, sfm, sf, sm;
244   sff = m_ThreaderSFF[0];
245   smm = m_ThreaderSMM[0];
246   sfm = m_ThreaderSFM[0];
247   sf  = m_ThreaderSF[0];
248   sm  = m_ThreaderSM[0];
249
250   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
251     sff +=  m_ThreaderSFF[t];
252     smm +=  m_ThreaderSMM[t];
253     sfm +=  m_ThreaderSFM[t];
254     if ( this->m_SubtractMean ) {
255       sf +=  m_ThreaderSF[t];
256       sm +=  m_ThreaderSM[t];
257     }
258   }
259
260   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
261     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
262     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
263     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
264   }
265
266
267   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
268   MeasureType measure;
269   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
270     measure = sfm / denom;
271   } else {
272     measure = itk::NumericTraits< MeasureType >::Zero;
273   }
274
275
276   return measure;
277 }
278
279
280 template < class TFixedImage, class TMovingImage  >
281 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
282 ::MeasureType
283 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
284 ::ComputeSums( const ParametersType & parameters ) const
285 {
286   //No checking for the fixed image,  done in the caller
287   //Reset the accumulators
288   memset( m_ThreaderSFF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
289   memset( m_ThreaderSMM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
290   memset( m_ThreaderSFM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
291   if(this->m_SubtractMean) {
292     memset( m_ThreaderSF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
293     memset( m_ThreaderSM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
294   }
295
296
297   // Set up the parameters in the transform
298   this->m_Transform->SetParameters( parameters );
299 #if ITK_VERSION_MAJOR < 4
300   this->m_Parameters = parameters;
301 #endif
302
303   // MUST BE CALLED TO INITIATE PROCESSING
304   this->GetValueMultiThreadedInitiate();
305
306   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
307                  << this->m_NumberOfPixelsCounted << " / "
308                  << this->m_NumberOfFixedImageSamples
309                  << std::endl );
310
311   if( this->m_NumberOfPixelsCounted <
312       this->m_NumberOfFixedImageSamples / 4 ) {
313     itkExceptionMacro( "Too many samples map outside moving image buffer: "
314                        << this->m_NumberOfPixelsCounted << " / "
315                        << this->m_NumberOfFixedImageSamples
316                        << std::endl );
317   }
318
319   // Accumulate the threads
320   m_SFF = m_ThreaderSFF[0];
321   m_SMM = m_ThreaderSMM[0];
322   m_SFM = m_ThreaderSFM[0];
323   m_SF  = m_ThreaderSF[0];
324   m_SM  = m_ThreaderSM[0];
325
326   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
327     m_SFF +=  m_ThreaderSFF[t];
328     m_SMM +=  m_ThreaderSMM[t];
329     m_SFM +=  m_ThreaderSFM[t];
330     if ( this->m_SubtractMean ) {
331       m_SF +=  m_ThreaderSF[t];
332       m_SM +=  m_ThreaderSM[t];
333     }
334   }
335
336   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
337     m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
338     m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
339     m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
340     m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
341     m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
342   }
343
344
345   m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
346   MeasureType measure;
347   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
348     measure = m_SFM / m_Denom;
349   } else {
350     measure = itk::NumericTraits< MeasureType >::Zero;
351   }
352
353   return measure;
354 }
355
356
357 template < class TFixedImage, class TMovingImage  >
358 inline bool
359 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
360 ::GetValueAndDerivativeThreadProcessSample(
361   unsigned int threadID,
362   unsigned long fixedImageSample,
363   const MovingImagePointType & itkNotUsed(mappedPoint),
364   double movingImageValue,
365   const ImageDerivativesType &
366   movingImageGradientValue
367 ) const
368 {
369
370   const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
371   const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
372
373   // Need to use one of the threader transforms if we're
374   // not in thread 0.
375   //
376   // Use a raw pointer here to avoid the overhead of smart pointers.
377   // For instance, Register and UnRegister have mutex locks around
378   // the reference counts.
379   TransformType* transform;
380
381   if (threadID > 0) {
382     transform = this->m_ThreaderTransform[threadID - 1];
383   } else {
384     transform = this->m_Transform;
385   }
386
387   // Jacobian should be evaluated at the unmapped (fixed image) point.
388 #if ITK_VERSION_MAJOR >= 4
389   TransformJacobianType jacobian;
390   transform->ComputeJacobianWithRespectToParameters( fixedImagePoint, jacobian );
391 #else
392   const TransformJacobianType & jacobian = transform->GetJacobian( fixedImagePoint );
393 #endif
394
395 //          for(unsigned int par=0; par<this->m_NumberOfParameters; par++)
396 //            {
397 //       RealType sumF = itk::NumericTraits< RealType >::Zero;
398 //       RealType sumM = itk::NumericTraits< RealType >::Zero;
399 //       for(unsigned int dim=0; dim<MovingImageDimension; dim++)
400 //         {
401 //           const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
402 //           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 )
403 //             {
404 //               sumF += (fixedImageValue-m_FixedMean)  * differential;
405 //               sumM += (movingImageValue-m_MovingMean) * differential;
406 //             }
407 //           else
408 //             {
409 //               sumF += differential * fixedImageValue;
410 //               sumM += differential * movingImageValue;
411 //             }
412 //         }
413 //       m_ThreaderDerivativeF[threadID][par] += sumF;
414 //       m_ThreaderDerivativeM[threadID][par] += sumM;
415 //            }
416
417   // JV
418   unsigned int par, dim;
419   RealType differential;
420   for( par=0; par<this->m_NumberOfParameters; par+=3) {
421     // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
422     if (jacobian( 0, par ) ) {
423       for(dim=0; dim<3; dim++) {
424         differential = jacobian( dim, par+dim ) * movingImageGradientValue[dim];
425
426         if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
427           m_ThreaderDerivativeF[threadID][par+dim]+= (fixedImageValue-m_FixedMean)  * differential;
428           m_ThreaderDerivativeM[threadID][par+dim]+= (movingImageValue-m_MovingMean) * differential;
429         } else {
430           m_ThreaderDerivativeF[threadID][par+dim]+= differential * fixedImageValue;
431           m_ThreaderDerivativeM[threadID][par+dim]+= differential * movingImageValue;
432         }
433       }
434     }
435   }
436
437   return true;
438 }
439
440
441 /**
442  * Get the both Value and Derivative Measure
443  */
444 template < class TFixedImage, class TMovingImage  >
445 void
446 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
447 ::GetValueAndDerivative( const ParametersType & parameters,
448                          MeasureType & value,
449                          DerivativeType & derivative) const
450 {
451
452   if( !this->m_FixedImage ) {
453     itkExceptionMacro( << "Fixed image has not been assigned" );
454   }
455
456   // Set up the parameters in the transform
457   this->m_Transform->SetParameters( parameters );
458 #if ITK_VERSION_MAJOR < 4
459   this->m_Parameters = parameters;
460 #endif
461
462   //We need the sums and the value to be calculated first
463   value=this->ComputeSums(parameters);
464
465   //Set output values to zero
466   if(derivative.GetSize() != this->m_NumberOfParameters) {
467     derivative = DerivativeType( this->m_NumberOfParameters );
468   }
469   memset( derivative.data_block(),
470           0,
471           this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
472
473   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
474     memset( m_ThreaderDerivativeF[threadID].data_block(),
475             0,
476             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
477
478     memset( m_ThreaderDerivativeM[threadID].data_block(),
479             0,
480             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
481   }
482
483   // MUST BE CALLED TO INITIATE PROCESSING
484   this->GetValueAndDerivativeMultiThreadedInitiate();
485
486   // Accumulate over the threads
487   DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
488   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
489     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
490       derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
491       derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
492     }
493   }
494
495   //Compute derivatives
496   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
497     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
498       derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
499     }
500   } else {
501     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
502       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
503     }
504   }
505
506 }
507
508
509 /**
510  * Get the match measure derivative
511  */
512 template < class TFixedImage, class TMovingImage  >
513 void
514 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
515 ::GetDerivative( const ParametersType & parameters,
516                  DerivativeType & derivative ) const
517 {
518   if( !this->m_FixedImage ) {
519     itkExceptionMacro( << "Fixed image has not been assigned" );
520   }
521
522   MeasureType value;
523   // call the combined version
524   this->GetValueAndDerivative( parameters, value, derivative );
525 }
526
527 } // end namespace clitk
528
529
530 #endif