1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef __clitkOptNormalizedCorrelationImageToImageMetric_txx
20 #define __clitkOptNormalizedCorrelationImageToImageMetric_txx
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"
36 template < class TFixedImage, class TMovingImage >
37 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
38 ::NormalizedCorrelationImageToImageMetric()
40 this->SetComputeGradient(true);
47 m_ThreaderDerivativeF = NULL;
48 m_ThreaderDerivativeM = NULL;
49 this->m_WithinThreadPreProcess = false;
50 this->m_WithinThreadPostProcess = false;
52 // For backward compatibility, the default behavior is to use all the pixels
53 // in the fixed image.
54 this->UseAllPixelsOn();
60 template < class TFixedImage, class TMovingImage >
61 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
62 ::~NormalizedCorrelationImageToImageMetric()
64 if(m_ThreaderSFF != NULL) {
65 delete [] m_ThreaderSFF;
69 if(m_ThreaderSMM != NULL) {
70 delete [] m_ThreaderSMM;
74 if(m_ThreaderSFM != NULL) {
75 delete [] m_ThreaderSFM;
79 if(m_ThreaderSF != NULL) {
80 delete [] m_ThreaderSF;
84 if(m_ThreaderSM != NULL) {
85 delete [] m_ThreaderSM;
89 if(m_ThreaderDerivativeF != NULL) {
90 delete [] m_ThreaderDerivativeF;
92 m_ThreaderDerivativeF = NULL;
94 if(m_ThreaderDerivativeM != NULL) {
95 delete [] m_ThreaderDerivativeM;
97 m_ThreaderDerivativeM = NULL;
101 * Print out internal information about this class
103 template < class TFixedImage, class TMovingImage >
105 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
106 ::PrintSelf(std::ostream& os, itk::Indent indent) const
109 Superclass::PrintSelf(os, indent);
117 template <class TFixedImage, class TMovingImage>
119 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
120 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
123 ::Initialize(void) throw ( itk::ExceptionObject )
127 this->Superclass::Initialize();
128 this->Superclass::MultiThreadingInitialize();
132 * Allocate memory for the accumulators (set to zero in GetValue)
134 if(m_ThreaderSFF != NULL) {
135 delete [] m_ThreaderSFF;
137 #if ITK_VERSION_MAJOR <= 4
138 m_ThreaderSFF = new double[this->m_NumberOfThreads];
140 m_ThreaderSFF = new double[this->m_NumberOfWorkUnits];
144 if(m_ThreaderSMM != NULL) {
145 delete [] m_ThreaderSMM;
147 #if ITK_VERSION_MAJOR <= 4
148 m_ThreaderSMM = new double[this->m_NumberOfThreads];
150 m_ThreaderSMM = new double[this->m_NumberOfWorkUnits];
153 if(m_ThreaderSFM != NULL) {
154 delete [] m_ThreaderSFM;
156 #if ITK_VERSION_MAJOR <= 4
157 m_ThreaderSFM = new double[this->m_NumberOfThreads];
159 m_ThreaderSFM = new double[this->m_NumberOfWorkUnits];
162 if(this->m_SubtractMean) {
163 if(m_ThreaderSF != NULL) {
164 delete [] m_ThreaderSF;
166 #if ITK_VERSION_MAJOR <= 4
167 m_ThreaderSF = new double[this->m_NumberOfThreads];
169 m_ThreaderSF = new double[this->m_NumberOfWorkUnits];
172 if(m_ThreaderSM != NULL) {
173 delete [] m_ThreaderSM;
175 #if ITK_VERSION_MAJOR <= 4
176 m_ThreaderSM = new double[this->m_NumberOfThreads];
178 m_ThreaderSM = new double[this->m_NumberOfWorkUnits];
182 if(m_ThreaderDerivativeF != NULL) {
183 delete [] m_ThreaderDerivativeF;
185 #if ITK_VERSION_MAJOR <= 4
186 m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfThreads];
187 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
189 m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfWorkUnits];
190 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
192 m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
195 if(m_ThreaderDerivativeM != NULL) {
196 delete [] m_ThreaderDerivativeM;
198 #if ITK_VERSION_MAJOR <= 4
199 m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfThreads];
200 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
202 m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfWorkUnits];
203 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
205 m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
210 template < class TFixedImage, class TMovingImage >
212 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
213 ::GetValueThreadProcessSample(
214 unsigned int threadID,
215 unsigned long fixedImageSample,
216 const MovingImagePointType & itkNotUsed(mappedPoint),
217 double movingImageValue) const
219 const RealType fixedImageValue= this->m_FixedImageSamples[fixedImageSample].value;
220 m_ThreaderSFF[threadID] += fixedImageValue * fixedImageValue;
221 m_ThreaderSMM[threadID] += movingImageValue * movingImageValue;
222 m_ThreaderSFM[threadID] += fixedImageValue * movingImageValue;
223 if ( this->m_SubtractMean ) {
224 m_ThreaderSF[threadID] += fixedImageValue;
225 m_ThreaderSM[threadID] += movingImageValue;
231 template < class TFixedImage, class TMovingImage >
232 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
234 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
235 ::GetValue( const ParametersType & parameters ) const
237 itkDebugMacro("GetValue( " << parameters << " ) ");
239 if( !this->m_FixedImage ) {
240 itkExceptionMacro( << "Fixed image has not been assigned" );
244 //Reset the accumulators
245 #if ITK_VERSION_MAJOR <= 4
246 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
247 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
248 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
250 memset( m_ThreaderSFF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
251 memset( m_ThreaderSMM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
252 memset( m_ThreaderSFM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
254 if(this->m_SubtractMean) {
255 #if ITK_VERSION_MAJOR <= 4
256 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
257 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
259 memset( m_ThreaderSF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
260 memset( m_ThreaderSM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
265 // Set up the parameters in the transform
266 this->m_Transform->SetParameters( parameters );
268 // MUST BE CALLED TO INITIATE PROCESSING
269 this->GetValueMultiThreadedInitiate();
271 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
272 << this->m_NumberOfPixelsCounted << " / "
273 << this->m_NumberOfFixedImageSamples
276 if( this->m_NumberOfPixelsCounted <
277 this->m_NumberOfFixedImageSamples / 4 ) {
278 itkExceptionMacro( "Too many samples map outside moving image buffer: "
279 << this->m_NumberOfPixelsCounted << " / "
280 << this->m_NumberOfFixedImageSamples
284 // Accumulate the threads
285 AccumulateType sff, smm, sfm, sf, sm;
286 sff = m_ThreaderSFF[0];
287 smm = m_ThreaderSMM[0];
288 sfm = m_ThreaderSFM[0];
289 sf = m_ThreaderSF[0];
290 sm = m_ThreaderSM[0];
292 #if ITK_VERSION_MAJOR <= 4
293 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
295 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
297 sff += m_ThreaderSFF[t];
298 smm += m_ThreaderSMM[t];
299 sfm += m_ThreaderSFM[t];
300 if ( this->m_SubtractMean ) {
301 sf += m_ThreaderSF[t];
302 sm += m_ThreaderSM[t];
306 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
307 sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
308 smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
309 sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
313 const RealType denom = -1.0 * vcl_sqrt(sff * smm );
315 if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
316 measure = sfm / denom;
318 measure = itk::NumericTraits< MeasureType >::Zero;
326 template < class TFixedImage, class TMovingImage >
327 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
329 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
330 ::ComputeSums( const ParametersType & parameters ) const
332 //No checking for the fixed image, done in the caller
333 //Reset the accumulators
334 #if ITK_VERSION_MAJOR <= 4
335 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
336 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
337 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
339 memset( m_ThreaderSFF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
340 memset( m_ThreaderSMM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
341 memset( m_ThreaderSFM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
343 if(this->m_SubtractMean) {
344 #if ITK_VERSION_MAJOR <= 4
345 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
346 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
348 memset( m_ThreaderSF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
349 memset( m_ThreaderSM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
354 // Set up the parameters in the transform
355 this->m_Transform->SetParameters( parameters );
357 // MUST BE CALLED TO INITIATE PROCESSING
358 this->GetValueMultiThreadedInitiate();
360 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
361 << this->m_NumberOfPixelsCounted << " / "
362 << this->m_NumberOfFixedImageSamples
365 if( this->m_NumberOfPixelsCounted <
366 this->m_NumberOfFixedImageSamples / 4 ) {
367 itkExceptionMacro( "Too many samples map outside moving image buffer: "
368 << this->m_NumberOfPixelsCounted << " / "
369 << this->m_NumberOfFixedImageSamples
373 // Accumulate the threads
374 m_SFF = m_ThreaderSFF[0];
375 m_SMM = m_ThreaderSMM[0];
376 m_SFM = m_ThreaderSFM[0];
377 m_SF = m_ThreaderSF[0];
378 m_SM = m_ThreaderSM[0];
380 #if ITK_VERSION_MAJOR <= 4
381 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
383 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
385 m_SFF += m_ThreaderSFF[t];
386 m_SMM += m_ThreaderSMM[t];
387 m_SFM += m_ThreaderSFM[t];
388 if ( this->m_SubtractMean ) {
389 m_SF += m_ThreaderSF[t];
390 m_SM += m_ThreaderSM[t];
394 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
395 m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
396 m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
397 m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
398 m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
399 m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
403 m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
405 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
406 measure = m_SFM / m_Denom;
408 measure = itk::NumericTraits< MeasureType >::Zero;
415 template < class TFixedImage, class TMovingImage >
417 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
418 ::GetValueAndDerivativeThreadProcessSample(
419 unsigned int threadID,
420 unsigned long fixedImageSample,
421 const MovingImagePointType & itkNotUsed(mappedPoint),
422 double movingImageValue,
423 const ImageDerivativesType &
424 movingImageGradientValue
428 const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
429 const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
431 // Need to use one of the threader transforms if we're
434 // Use a raw pointer here to avoid the overhead of smart pointers.
435 // For instance, Register and UnRegister have mutex locks around
436 // the reference counts.
437 TransformType* transform;
440 transform = this->m_ThreaderTransform[threadID - 1];
442 transform = this->m_Transform;
445 // Jacobian should be evaluated at the unmapped (fixed image) point.
446 TransformJacobianType jacobian;
447 transform->ComputeJacobianWithRespectToParameters(fixedImagePoint, jacobian);
449 for(unsigned int par=0; par<this->m_NumberOfParameters; par++) {
450 RealType sumF = itk::NumericTraits< RealType >::Zero;
451 RealType sumM = itk::NumericTraits< RealType >::Zero;
452 for(unsigned int dim=0; dim<MovingImageDimension; dim++) {
453 const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
454 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
455 sumF += (fixedImageValue-m_FixedMean) * differential;
456 sumM += (movingImageValue-m_MovingMean) * differential;
458 sumF += differential * fixedImageValue;
459 sumM += differential * movingImageValue;
462 m_ThreaderDerivativeF[threadID][par] += sumF;
463 m_ThreaderDerivativeM[threadID][par] += sumM;
471 * Get the both Value and Derivative Measure
473 template < class TFixedImage, class TMovingImage >
475 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
476 ::GetValueAndDerivative( const ParametersType & parameters,
478 DerivativeType & derivative) const
481 if( !this->m_FixedImage ) {
482 itkExceptionMacro( << "Fixed image has not been assigned" );
485 // Set up the parameters in the transform
486 this->m_Transform->SetParameters( parameters );
488 //We need the sums and the value to be calculated first
489 value=this->ComputeSums(parameters);
491 //Set output values to zero
492 if(derivative.GetSize() != this->m_NumberOfParameters) {
493 derivative = DerivativeType( this->m_NumberOfParameters );
495 memset( derivative.data_block(),
497 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
498 #if ITK_VERSION_MAJOR <= 4
499 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
501 for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
503 memset( m_ThreaderDerivativeF[threadID].data_block(),
505 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
507 memset( m_ThreaderDerivativeM[threadID].data_block(),
509 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
512 // MUST BE CALLED TO INITIATE PROCESSING
513 this->GetValueAndDerivativeMultiThreadedInitiate();
515 // Accumulate over the threads
516 DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
517 #if ITK_VERSION_MAJOR <= 4
518 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
520 for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
522 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
523 derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
524 derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
528 //Compute derivatives
529 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
530 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
531 derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
534 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
535 derivative[i] = itk::NumericTraits< MeasureType >::Zero;
543 * Get the match measure derivative
545 template < class TFixedImage, class TMovingImage >
547 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
548 ::GetDerivative( const ParametersType & parameters,
549 DerivativeType & derivative ) const
551 if( !this->m_FixedImage ) {
552 itkExceptionMacro( << "Fixed image has not been assigned" );
556 // call the combined version
557 this->GetValueAndDerivative( parameters, value, derivative );
560 } // end namespace clitk