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)
 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 #if ITK_VERSION_MAJOR <= 4
 134   m_ThreaderSFF = new double[this->m_NumberOfThreads];
 135 #else
 136   m_ThreaderSFF = new double[this->m_NumberOfWorkUnits];
 137 #endif
 138
 139
 140   if(m_ThreaderSMM != NULL) {
 141     delete [] m_ThreaderSMM;
 142   }
 143 #if ITK_VERSION_MAJOR <= 4
 144   m_ThreaderSMM = new double[this->m_NumberOfThreads];
 145 #else
 146   m_ThreaderSMM = new double[this->m_NumberOfWorkUnits];
 147 #endif
 148
 149   if(m_ThreaderSFM != NULL) {
 150     delete [] m_ThreaderSFM;
 151   }
 152 #if ITK_VERSION_MAJOR <= 4
 153   m_ThreaderSFM = new double[this->m_NumberOfThreads];
 154 #else
 155   m_ThreaderSFM = new double[this->m_NumberOfWorkUnits];
 156 #endif
 157
 158   if(this->m_SubtractMean) {
 159     if(m_ThreaderSF != NULL) {
 160       delete [] m_ThreaderSF;
 161     }
 162 #if ITK_VERSION_MAJOR <= 4
 163     m_ThreaderSF = new double[this->m_NumberOfThreads];
 164 #else
 165     m_ThreaderSF = new double[this->m_NumberOfWorkUnits];
 166 #endif
 167
 168     if(m_ThreaderSM != NULL) {
 169       delete [] m_ThreaderSM;
 170     }
 171 #if ITK_VERSION_MAJOR <= 4
 172     m_ThreaderSM = new double[this->m_NumberOfThreads];
 173 #else
 174     m_ThreaderSM = new double[this->m_NumberOfWorkUnits];
 175 #endif
 176   }
 177
 178   if(m_ThreaderDerivativeF != NULL) {
 179     delete [] m_ThreaderDerivativeF;
 180   }
 181 #if ITK_VERSION_MAJOR <= 4
 182   m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfThreads];
 183   for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
 184 #else
 185   m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfWorkUnits];
 186   for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
 187 #endif
 188     m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
 189   }
 190
 191   if(m_ThreaderDerivativeM != NULL) {
 192     delete [] m_ThreaderDerivativeM;
 193   }
 194 #if ITK_VERSION_MAJOR <= 4
 195   m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfThreads];
 196   for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
 197 #else
 198   m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfWorkUnits];
 199   for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
 200 #endif
 201     m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
 202   }
 203 }
 204
 205
 206 template < class TFixedImage, class TMovingImage  >
 207 inline bool
 208 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 209 ::GetValueThreadProcessSample(
 210   unsigned int threadID,
 211   unsigned long fixedImageSample,
 212   const MovingImagePointType & itkNotUsed(mappedPoint),
 213   double movingImageValue) const
 214 {
 215   const RealType fixedImageValue= this->m_FixedImageSamples[fixedImageSample].value;
 216   m_ThreaderSFF[threadID] += fixedImageValue  * fixedImageValue;
 217   m_ThreaderSMM[threadID] += movingImageValue * movingImageValue;
 218   m_ThreaderSFM[threadID] += fixedImageValue  * movingImageValue;
 219   if ( this->m_SubtractMean ) {
 220     m_ThreaderSF[threadID] += fixedImageValue;
 221     m_ThreaderSM[threadID] += movingImageValue;
 222   }
 223
 224   return true;
 225 }
 226
 227 template < class TFixedImage, class TMovingImage  >
 228 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 229 ::MeasureType
 230 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 231 ::GetValue( const ParametersType & parameters ) const
 232 {
 233   itkDebugMacro("GetValue( " << parameters << " ) ");
 234
 235   if( !this->m_FixedImage ) {
 236     itkExceptionMacro( << "Fixed image has not been assigned" );
 237   }
 238
 239
 240   //Reset the accumulators
 241 #if ITK_VERSION_MAJOR <= 4
 242   memset( m_ThreaderSFF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 243   memset( m_ThreaderSMM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 244   memset( m_ThreaderSFM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 245 #else
 246   memset( m_ThreaderSFF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 247   memset( m_ThreaderSMM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 248   memset( m_ThreaderSFM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 249 #endif
 250   if(this->m_SubtractMean) {
 251 #if ITK_VERSION_MAJOR <= 4
 252     memset( m_ThreaderSF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 253     memset( m_ThreaderSM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 254 #else
 255     memset( m_ThreaderSF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 256     memset( m_ThreaderSM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 257 #endif
 258   }
 259
 260
 261   // Set up the parameters in the transform
 262   this->m_Transform->SetParameters( parameters );
 263
 264   // MUST BE CALLED TO INITIATE PROCESSING
 265   this->GetValueMultiThreadedInitiate();
 266
 267   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 268                  << this->m_NumberOfPixelsCounted << " / "
 269                  << this->m_NumberOfFixedImageSamples
 270                  << std::endl );
 271
 272   if( this->m_NumberOfPixelsCounted <
 273       this->m_NumberOfFixedImageSamples / 4 ) {
 274     itkExceptionMacro( "Too many samples map outside moving image buffer: "
 275                        << this->m_NumberOfPixelsCounted << " / "
 276                        << this->m_NumberOfFixedImageSamples
 277                        << std::endl );
 278   }
 279
 280   // Accumulate the threads
 281   AccumulateType sff, smm, sfm, sf, sm;
 282   sff = m_ThreaderSFF[0];
 283   smm = m_ThreaderSMM[0];
 284   sfm = m_ThreaderSFM[0];
 285   sf  = m_ThreaderSF[0];
 286   sm  = m_ThreaderSM[0];
 287
 288 #if ITK_VERSION_MAJOR <= 4
 289   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
 290 #else
 291   for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
 292 #endif
 293     sff +=  m_ThreaderSFF[t];
 294     smm +=  m_ThreaderSMM[t];
 295     sfm +=  m_ThreaderSFM[t];
 296     if ( this->m_SubtractMean ) {
 297       sf +=  m_ThreaderSF[t];
 298       sm +=  m_ThreaderSM[t];
 299     }
 300   }
 301
 302   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 303     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 304     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 305     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 306   }
 307
 308
 309   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 310   MeasureType measure;
 311   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 312     measure = sfm / denom;
 313   } else {
 314     measure = itk::NumericTraits< MeasureType >::Zero;
 315   }
 316
 317
 318   return measure;
 319 }
 320
 321
 322 template < class TFixedImage, class TMovingImage  >
 323 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 324 ::MeasureType
 325 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 326 ::ComputeSums( const ParametersType & parameters ) const
 327 {
 328   //No checking for the fixed image,  done in the caller
 329   //Reset the accumulators
 330 #if ITK_VERSION_MAJOR <= 4
 331   memset( m_ThreaderSFF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 332   memset( m_ThreaderSMM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 333   memset( m_ThreaderSFM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 334 #else
 335   memset( m_ThreaderSFF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 336   memset( m_ThreaderSMM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 337   memset( m_ThreaderSFM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 338 #endif
 339   if(this->m_SubtractMean) {
 340 #if ITK_VERSION_MAJOR <= 4
 341     memset( m_ThreaderSF,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 342     memset( m_ThreaderSM,  0,  this->m_NumberOfThreads * sizeof(AccumulateType) );
 343 #else
 344     memset( m_ThreaderSF,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 345     memset( m_ThreaderSM,  0,  this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
 346 #endif
 347   }
 348
 349
 350   // Set up the parameters in the transform
 351   this->m_Transform->SetParameters( parameters );
 352
 353   // MUST BE CALLED TO INITIATE PROCESSING
 354   this->GetValueMultiThreadedInitiate();
 355
 356   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
 357                  << this->m_NumberOfPixelsCounted << " / "
 358                  << this->m_NumberOfFixedImageSamples
 359                  << std::endl );
 360
 361   if( this->m_NumberOfPixelsCounted <
 362       this->m_NumberOfFixedImageSamples / 4 ) {
 363     itkExceptionMacro( "Too many samples map outside moving image buffer: "
 364                        << this->m_NumberOfPixelsCounted << " / "
 365                        << this->m_NumberOfFixedImageSamples
 366                        << std::endl );
 367   }
 368
 369   // Accumulate the threads
 370   m_SFF = m_ThreaderSFF[0];
 371   m_SMM = m_ThreaderSMM[0];
 372   m_SFM = m_ThreaderSFM[0];
 373   m_SF  = m_ThreaderSF[0];
 374   m_SM  = m_ThreaderSM[0];
 375
 376 #if ITK_VERSION_MAJOR <= 4
 377   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
 378 #else
 379   for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
 380 #endif
 381     m_SFF +=  m_ThreaderSFF[t];
 382     m_SMM +=  m_ThreaderSMM[t];
 383     m_SFM +=  m_ThreaderSFM[t];
 384     if ( this->m_SubtractMean ) {
 385       m_SF +=  m_ThreaderSF[t];
 386       m_SM +=  m_ThreaderSM[t];
 387     }
 388   }
 389
 390   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 391     m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
 392     m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
 393     m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
 394     m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
 395     m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
 396   }
 397
 398
 399   m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
 400   MeasureType measure;
 401   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
 402     measure = m_SFM / m_Denom;
 403   } else {
 404     measure = itk::NumericTraits< MeasureType >::Zero;
 405   }
 406
 407   return measure;
 408 }
 409
 410
 411 template < class TFixedImage, class TMovingImage  >
 412 inline bool
 413 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 414 ::GetValueAndDerivativeThreadProcessSample(
 415   unsigned int threadID,
 416   unsigned long fixedImageSample,
 417   const MovingImagePointType & itkNotUsed(mappedPoint),
 418   double movingImageValue,
 419   const ImageDerivativesType &
 420   movingImageGradientValue
 421 ) const
 422 {
 423
 424   const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
 425   const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
 426
 427   // Need to use one of the threader transforms if we're
 428   // not in thread 0.
 429   //
 430   // Use a raw pointer here to avoid the overhead of smart pointers.
 431   // For instance, Register and UnRegister have mutex locks around
 432   // the reference counts.
 433   TransformType* transform;
 434
 435   if (threadID > 0) {
 436     transform = this->m_ThreaderTransform[threadID - 1];
 437   } else {
 438     transform = this->m_Transform;
 439   }
 440
 441   // Jacobian should be evaluated at the unmapped (fixed image) point.
 442   TransformJacobianType jacobian;
 443   transform->ComputeJacobianWithRespectToParameters( fixedImagePoint, jacobian );
 444
 445 //          for(unsigned int par=0; par<this->m_NumberOfParameters; par++)
 446 //            {
 447 //       RealType sumF = itk::NumericTraits< RealType >::Zero;
 448 //       RealType sumM = itk::NumericTraits< RealType >::Zero;
 449 //       for(unsigned int dim=0; dim<MovingImageDimension; dim++)
 450 //         {
 451 //           const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
 452 //           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 )
 453 //             {
 454 //               sumF += (fixedImageValue-m_FixedMean)  * differential;
 455 //               sumM += (movingImageValue-m_MovingMean) * differential;
 456 //             }
 457 //           else
 458 //             {
 459 //               sumF += differential * fixedImageValue;
 460 //               sumM += differential * movingImageValue;
 461 //             }
 462 //         }
 463 //       m_ThreaderDerivativeF[threadID][par] += sumF;
 464 //       m_ThreaderDerivativeM[threadID][par] += sumM;
 465 //            }
 466
 467   // JV
 468   unsigned int par, dim;
 469   RealType differential;
 470   for( par=0; par<this->m_NumberOfParameters; par+=3) {
 471     // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
 472     if (jacobian( 0, par ) ) {
 473       for(dim=0; dim<3; dim++) {
 474         differential = jacobian( dim, par+dim ) * movingImageGradientValue[dim];
 475
 476         if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 477           m_ThreaderDerivativeF[threadID][par+dim]+= (fixedImageValue-m_FixedMean)  * differential;
 478           m_ThreaderDerivativeM[threadID][par+dim]+= (movingImageValue-m_MovingMean) * differential;
 479         } else {
 480           m_ThreaderDerivativeF[threadID][par+dim]+= differential * fixedImageValue;
 481           m_ThreaderDerivativeM[threadID][par+dim]+= differential * movingImageValue;
 482         }
 483       }
 484     }
 485   }
 486
 487   return true;
 488 }
 489
 490
 491 /**
 492  * Get the both Value and Derivative Measure
 493  */
 494 template < class TFixedImage, class TMovingImage  >
 495 void
 496 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 497 ::GetValueAndDerivative( const ParametersType & parameters,
 498                          MeasureType & value,
 499                          DerivativeType & derivative) const
 500 {
 501
 502   if( !this->m_FixedImage ) {
 503     itkExceptionMacro( << "Fixed image has not been assigned" );
 504   }
 505
 506   // Set up the parameters in the transform
 507   this->m_Transform->SetParameters( parameters );
 508
 509   //We need the sums and the value to be calculated first
 510   value=this->ComputeSums(parameters);
 511
 512   //Set output values to zero
 513   if(derivative.GetSize() != this->m_NumberOfParameters) {
 514     derivative = DerivativeType( this->m_NumberOfParameters );
 515   }
 516   memset( derivative.data_block(),
 517           0,
 518           this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 519
 520 #if ITK_VERSION_MAJOR <= 4
 521   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
 522 #else
 523   for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
 524 #endif
 525     memset( m_ThreaderDerivativeF[threadID].data_block(),
 526             0,
 527             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 528
 529     memset( m_ThreaderDerivativeM[threadID].data_block(),
 530             0,
 531             this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
 532   }
 533
 534   // MUST BE CALLED TO INITIATE PROCESSING
 535   this->GetValueAndDerivativeMultiThreadedInitiate();
 536
 537   // Accumulate over the threads
 538   DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
 539 #if ITK_VERSION_MAJOR <= 4
 540   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
 541 #else
 542   for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
 543 #endif
 544     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
 545       derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
 546       derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
 547     }
 548   }
 549
 550   //Compute derivatives
 551   if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
 552     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
 553       derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
 554     }
 555   } else {
 556     for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
 557       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 558     }
 559   }
 560
 561 }
 562
 563
 564 /**
 565  * Get the match measure derivative
 566  */
 567 template < class TFixedImage, class TMovingImage  >
 568 void
 569 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
 570 ::GetDerivative( const ParametersType & parameters,
 571                  DerivativeType & derivative ) const
 572 {
 573   if( !this->m_FixedImage ) {
 574     itkExceptionMacro( << "Fixed image has not been assigned" );
 575   }
 576
 577   MeasureType value;
 578   // call the combined version
 579   this->GetValueAndDerivative( parameters, value, derivative );
 580 }
 581
 582 } // end namespace clitk
 583
 584
 585 #endif