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 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
 121 ::Initialize(void)
 122 #else
 123 ::Initialize(void) throw ( itk::ExceptionObject )
 124 #endif
 125 {
 126
 127   this->Superclass::Initialize();
 128   this->Superclass::MultiThreadingInitialize();
 129
 130
 131   /**
 132    * Allocate memory for the accumulators (set to zero in GetValue)
 133    */
 134   if(m_ThreaderSFF != NULL) {
 135     delete [] m_ThreaderSFF;
 136   }
 137 #if ITK_VERSION_MAJOR <= 4
 138   m_ThreaderSFF = new double[this->m_NumberOfThreads];
 139 #else
 140   m_ThreaderSFF = new double[this->m_NumberOfWorkUnits];
 141 #endif
 142
 143
 144   if(m_ThreaderSMM != NULL) {
 145     delete [] m_ThreaderSMM;
 146   }
 147 #if ITK_VERSION_MAJOR <= 4
 148   m_ThreaderSMM = new double[this->m_NumberOfThreads];
 149 #else
 150   m_ThreaderSMM = new double[this->m_NumberOfWorkUnits];
 151 #endif
 152
 153   if(m_ThreaderSFM != NULL) {
 154     delete [] m_ThreaderSFM;
 155   }
 156 #if ITK_VERSION_MAJOR <= 4
 157   m_ThreaderSFM = new double[this->m_NumberOfThreads];
 158 #else
 159   m_ThreaderSFM = new double[this->m_NumberOfWorkUnits];
 160 #endif
 161
 162   if(this->m_SubtractMean) {
 163     if(m_ThreaderSF != NULL) {
 164       delete [] m_ThreaderSF;
 165     }
 166 #if ITK_VERSION_MAJOR <= 4
 167     m_ThreaderSF = new double[this->m_NumberOfThreads];
 168 #else
 169     m_ThreaderSF = new double[this->m_NumberOfWorkUnits];
 170 #endif
 171
 172     if(m_ThreaderSM != NULL) {
 173       delete [] m_ThreaderSM;
 174     }
 175 #if ITK_VERSION_MAJOR <= 4
 176     m_ThreaderSM = new double[this->m_NumberOfThreads];
 177 #else
 178     m_ThreaderSM = new double[this->m_NumberOfWorkUnits];
 179 #endif
 180   }
 181
 182   if(m_ThreaderDerivativeF != NULL) {
 183     delete [] m_ThreaderDerivativeF;
 184   }
 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++) {
 188 #else
 189   m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfWorkUnits];
 190   for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
 191 #endif
 192     m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
 193   }
 194
 195   if(m_ThreaderDerivativeM != NULL) {
 196     delete [] m_ThreaderDerivativeM;
 197   }
 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++) {
 201 #else
 202   m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfWorkUnits];
 203   for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
 204 #endif
 205     m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
 206   }
 207 }
 208
 209
 210 template < class TFixedImage, class TMovingImage  >
 211 inline bool
 212 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 213 ::GetValueThreadProcessSample(
 214   unsigned int threadID,
 215   unsigned long fixedImageSample,
 216   const MovingImagePointType & itkNotUsed(mappedPoint),
 217   double movingImageValue) const
 218 {
 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;
 226   }
 227
 228   return true;
 229 }
 230
 231 template < class TFixedImage, class TMovingImage  >
 232 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 233 ::MeasureType
 234 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 235 ::GetValue( const ParametersType & parameters ) const
 236 {
 237   itkDebugMacro("GetValue( " << parameters << " ) ");
 238
 239   if( !this->m_FixedImage ) {
 240     itkExceptionMacro( << "Fixed image has not been assigned" );
 241   }
 242
 243
 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) );
 249 #else
 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) );
 253 #endif
 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) );
 258 #else
 259     memset( m_ThreaderSF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 260     memset( m_ThreaderSM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 261 #endif
 262   }
 263
 264
 265   // Set up the parameters in the transform
 266   this->m_Transform->SetParameters( parameters );
 267
 268   // MUST BE CALLED TO INITIATE PROCESSING
 269   this->GetValueMultiThreadedInitiate();
 270
 271   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 272                  << this->m_NumberOfPixelsCounted << " / "
 273                  << this->m_NumberOfFixedImageSamples
 274                  << std::endl );
 275
 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
 281                        << std::endl );
 282   }
 283
 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];
 291
 292 #if ITK_VERSION_MAJOR <= 4
 293   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
 294 #else
 295   for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
 296 #endif
 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];
 303     }
 304   }
 305
 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 );
 310   }
 311
 312
 313   const RealType denom = -1.0 * std::sqrt(sff * smm );
 314   MeasureType measure;
 315   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 316     measure = sfm / denom;
 317   } else {
 318     measure = itk::NumericTraits< MeasureType >::Zero;
 319   }
 320
 321
 322   return measure;
 323 }
 324
 325
 326 template < class TFixedImage, class TMovingImage  >
 327 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 328 ::MeasureType
 329 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 330 ::ComputeSums( const ParametersType & parameters ) const
 331 {
 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) );
 338 #else
 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) );
 342 #endif
 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) );
 347 #else
 348     memset( m_ThreaderSF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 349     memset( m_ThreaderSM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 350 #endif
 351   }
 352
 353
 354   // Set up the parameters in the transform
 355   this->m_Transform->SetParameters( parameters );
 356
 357   // MUST BE CALLED TO INITIATE PROCESSING
 358   this->GetValueMultiThreadedInitiate();
 359
 360   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 361                  << this->m_NumberOfPixelsCounted << " / "
 362                  << this->m_NumberOfFixedImageSamples
 363                  << std::endl );
 364
 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
 370                        << std::endl );
 371   }
 372
 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];
 379
 380 #if ITK_VERSION_MAJOR <= 4
 381   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
 382 #else
 383   for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
 384 #endif
 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];
 391     }
 392   }
 393
 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;
 400   }
 401
 402
 403   m_Denom = -1.0 * std::sqrt(m_SFF * m_SMM );
 404   MeasureType measure;
 405   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
 406     measure = m_SFM / m_Denom;
 407   } else {
 408     measure = itk::NumericTraits< MeasureType >::Zero;
 409   }
 410
 411   return measure;
 412 }
 413
 414
 415 template < class TFixedImage, class TMovingImage  >
 416 inline bool
 417 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<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
 425 ) const
 426 {
 427
 428   const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
 429   const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
 430
 431   // Need to use one of the threader transforms if we're
 432   // not in thread 0.
 433   //
 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;
 438
 439   if (threadID > 0) {
 440     transform = this->m_ThreaderTransform[threadID - 1];
 441   } else {
 442     transform = this->m_Transform;
 443   }
 444
 445   // Jacobian should be evaluated at the unmapped (fixed image) point.
 446   TransformJacobianType jacobian;
 447   transform->ComputeJacobianWithRespectToParameters( fixedImagePoint, jacobian );
 448
 449 //          for(unsigned int par=0; par<this->m_NumberOfParameters; par++)
 450 //            {
 451 //       RealType sumF = itk::NumericTraits< RealType >::Zero;
 452 //       RealType sumM = itk::NumericTraits< RealType >::Zero;
 453 //       for(unsigned int dim=0; dim<MovingImageDimension; dim++)
 454 //         {
 455 //           const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
 456 //           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 )
 457 //             {
 458 //               sumF += (fixedImageValue-m_FixedMean)  * differential;
 459 //               sumM += (movingImageValue-m_MovingMean) * differential;
 460 //             }
 461 //           else
 462 //             {
 463 //               sumF += differential * fixedImageValue;
 464 //               sumM += differential * movingImageValue;
 465 //             }
 466 //         }
 467 //       m_ThreaderDerivativeF[threadID][par] += sumF;
 468 //       m_ThreaderDerivativeM[threadID][par] += sumM;
 469 //            }
 470
 471   // JV
 472   unsigned int par, dim;
 473   RealType differential;
 474   for( par=0; par<this->m_NumberOfParameters; par+=3) {
 475     // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
 476     if (jacobian( 0, par ) ) {
 477       for(dim=0; dim<3; dim++) {
 478         differential = jacobian( dim, par+dim ) * movingImageGradientValue[dim];
 479
 480         if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 481           m_ThreaderDerivativeF[threadID][par+dim]+= (fixedImageValue-m_FixedMean)  * differential;
 482           m_ThreaderDerivativeM[threadID][par+dim]+= (movingImageValue-m_MovingMean) * differential;
 483         } else {
 484           m_ThreaderDerivativeF[threadID][par+dim]+= differential * fixedImageValue;
 485           m_ThreaderDerivativeM[threadID][par+dim]+= differential * movingImageValue;
 486         }
 487       }
 488     }
 489   }
 490
 491   return true;
 492 }
 493
 494
 495 /**
 496  * Get the both Value and Derivative Measure
 497  */
 498 template < class TFixedImage, class TMovingImage  >
 499 void
 500 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 501 ::GetValueAndDerivative( const ParametersType & parameters,
 502                          MeasureType & value,
 503                          DerivativeType & derivative) const
 504 {
 505
 506   if( !this->m_FixedImage ) {
 507     itkExceptionMacro( << "Fixed image has not been assigned" );
 508   }
 509
 510   // Set up the parameters in the transform
 511   this->m_Transform->SetParameters( parameters );
 512
 513   //We need the sums and the value to be calculated first
 514   value=this->ComputeSums(parameters);
 515
 516   //Set output values to zero
 517   if(derivative.GetSize() != this->m_NumberOfParameters) {
 518     derivative = DerivativeType( this->m_NumberOfParameters );
 519   }
 520   memset( derivative.data_block(),
 521           0,
 522           this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 523
 524 #if ITK_VERSION_MAJOR <= 4
 525   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
 526 #else
 527   for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
 528 #endif
 529     memset( m_ThreaderDerivativeF[threadID].data_block(),
 530             0,
 531             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 532
 533     memset( m_ThreaderDerivativeM[threadID].data_block(),
 534             0,
 535             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 536   }
 537
 538   // MUST BE CALLED TO INITIATE PROCESSING
 539   this->GetValueAndDerivativeMultiThreadedInitiate();
 540
 541   // Accumulate over the threads
 542   DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
 543 #if ITK_VERSION_MAJOR <= 4
 544   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
 545 #else
 546   for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
 547 #endif
 548     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
 549       derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
 550       derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
 551     }
 552   }
 553
 554   //Compute derivatives
 555   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
 556     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
 557       derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
 558     }
 559   } else {
 560     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
 561       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 562     }
 563   }
 564
 565 }
 566
 567
 568 /**
 569  * Get the match measure derivative
 570  */
 571 template < class TFixedImage, class TMovingImage  >
 572 void
 573 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 574 ::GetDerivative( const ParametersType & parameters,
 575                  DerivativeType & derivative ) const
 576 {
 577   if( !this->m_FixedImage ) {
 578     itkExceptionMacro( << "Fixed image has not been assigned" );
 579   }
 580
 581   MeasureType value;
 582   // call the combined version
 583   this->GetValueAndDerivative( parameters, value, derivative );
 584 }
 585
 586 } // end namespace clitk
 587
 588
 589 #endif