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(ThreadIdType 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(ThreadIdType 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   ThreadIdType 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   ThreadIdType 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( ThreadIdType 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