registration/clitkNormalizedCorrelationImageToImageMetric.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 _clitkNormalizedCorrelationImageToImageMetric_txx
  20 #define _clitkNormalizedCorrelationImageToImageMetric_txx
  21
  22 // First make sure that the configuration is available.
  23 // This line can be removed once the optimized versions
  24 // gets integrated into the main directories.
  25 #include "itkConfigure.h"
  26
  27 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
  28 #include "clitkOptNormalizedCorrelationImageToImageMetric.txx"
  29 #else
  30
  31
  32 #include "clitkNormalizedCorrelationImageToImageMetric.h"
  33
  34 #include "itkImageRegionConstIteratorWithIndex.h"
  35
  36
  37 namespace clitk
  38 {
  39
  40 /*
  41  * Constructor
  42  */
  43 template <class TFixedImage, class TMovingImage>
  44 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
  45 ::NormalizedCorrelationImageToImageMetric()
  46 {
  47   m_SubtractMean = false;
  48 }
  49
  50 /*
  51  * Get the match Measure
  52  */
  53 template <class TFixedImage, class TMovingImage>
  54 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
  55 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
  56 ::GetValue( const TransformParametersType & parameters ) const
  57 {
  58
  59   FixedImageConstPointer fixedImage = this->m_FixedImage;
  60
  61   if( !fixedImage ) {
  62     itkExceptionMacro( << "Fixed image has not been assigned" );
  63   }
  64
  65   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
  66
  67   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
  68
  69   typename FixedImageType::IndexType index;
  70
  71   MeasureType measure;
  72
  73   this->m_NumberOfPixelsCounted = 0;
  74
  75   this->SetTransformParameters( parameters );
  76
  77   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
  78
  79   AccumulateType sff = itk::NumericTraits< AccumulateType >::Zero;
  80   AccumulateType smm = itk::NumericTraits< AccumulateType >::Zero;
  81   AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero;
  82   AccumulateType sf  = itk::NumericTraits< AccumulateType >::Zero;
  83   AccumulateType sm  = itk::NumericTraits< AccumulateType >::Zero;
  84
  85   while(!ti.IsAtEnd()) {
  86
  87     index = ti.GetIndex();
  88
  89     InputPointType inputPoint;
  90     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
  91
  92     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
  93       ++ti;
  94       continue;
  95     }
  96
  97     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
  98
  99     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 100       ++ti;
 101       continue;
 102     }
 103
 104     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 105       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 106       const RealType fixedValue   = ti.Get();
 107       sff += fixedValue  * fixedValue;
 108       smm += movingValue * movingValue;
 109       sfm += fixedValue  * movingValue;
 110       if ( this->m_SubtractMean ) {
 111         sf += fixedValue;
 112         sm += movingValue;
 113       }
 114       this->m_NumberOfPixelsCounted++;
 115     }
 116
 117     ++ti;
 118   }
 119
 120   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 121     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 122     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 123     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 124   }
 125
 126   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 127
 128   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 129     measure = sfm / denom;
 130   } else {
 131     measure = itk::NumericTraits< MeasureType >::Zero;
 132   }
 133
 134   return measure;
 135
 136 }
 137
 138
 139
 140
 141
 142 /*
 143  * Get the Derivative Measure
 144  */
 145 template < class TFixedImage, class TMovingImage>
 146 void
 147 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
 148 ::GetDerivative( const TransformParametersType & parameters,
 149                  DerivativeType & derivative ) const
 150 {
 151
 152   if( !this->GetGradientImage() ) {
 153     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 154   }
 155
 156   FixedImageConstPointer fixedImage = this->m_FixedImage;
 157
 158   if( !fixedImage ) {
 159     itkExceptionMacro( << "Fixed image has not been assigned" );
 160   }
 161
 162   const unsigned int dimension = FixedImageType::ImageDimension;
 163
 164   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
 165
 166   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 167
 168   typename FixedImageType::IndexType index;
 169
 170   this->m_NumberOfPixelsCounted = 0;
 171
 172   this->SetTransformParameters( parameters );
 173
 174   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
 175
 176   AccumulateType sff  = itk::NumericTraits< AccumulateType >::Zero;
 177   AccumulateType smm  = itk::NumericTraits< AccumulateType >::Zero;
 178   AccumulateType sfm = itk::NumericTraits< AccumulateType >::Zero;
 179   AccumulateType sf  = itk::NumericTraits< AccumulateType >::Zero;
 180   AccumulateType sm  = itk::NumericTraits< AccumulateType >::Zero;
 181
 182   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 183   derivative = DerivativeType( ParametersDimension );
 184   derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 185
 186   DerivativeType derivativeF = DerivativeType( ParametersDimension );
 187   derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 188
 189   DerivativeType derivativeM = DerivativeType( ParametersDimension );
 190   derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 191
 192   ti.GoToBegin();
 193   // First compute the sums
 194   while(!ti.IsAtEnd()) {
 195
 196     index = ti.GetIndex();
 197
 198     InputPointType inputPoint;
 199     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 200
 201     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 202       ++ti;
 203       continue;
 204     }
 205
 206     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 207
 208     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 209       ++ti;
 210       continue;
 211     }
 212
 213     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 214       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 215       const RealType fixedValue   = ti.Get();
 216       sff += fixedValue  * fixedValue;
 217       smm += movingValue * movingValue;
 218       sfm += fixedValue  * movingValue;
 219       if ( this->m_SubtractMean ) {
 220         sf += fixedValue;
 221         sm += movingValue;
 222       }
 223       this->m_NumberOfPixelsCounted++;
 224     }
 225
 226     ++ti;
 227   }
 228
 229   // Compute contributions to derivatives
 230   ti.GoToBegin();
 231   while(!ti.IsAtEnd()) {
 232
 233     index = ti.GetIndex();
 234
 235     InputPointType inputPoint;
 236     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 237
 238     if ( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 239       ++ti;
 240       continue;
 241     }
 242
 243     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 244
 245     if ( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 246       ++ti;
 247       continue;
 248     }
 249
 250     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 251       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 252       const RealType fixedValue   = ti.Get();
 253
 254 #if ITK_VERSION_MAJOR >= 4
 255       TransformJacobianType jacobian;
 256       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 257 #else
 258       const TransformJacobianType & jacobian =
 259         this->m_Transform->GetJacobian( inputPoint );
 260 #endif
 261
 262       // Get the gradient by NearestNeighboorInterpolation:
 263       // which is equivalent to round up the point components.
 264       typedef typename OutputPointType::CoordRepType CoordRepType;
 265       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 266       MovingImageContinuousIndexType;
 267
 268       MovingImageContinuousIndexType tempIndex;
 269       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 270
 271       typename MovingImageType::IndexType mappedIndex;
 272       mappedIndex.CopyWithRound( tempIndex );
 273
 274       const GradientPixelType gradient =
 275         this->GetGradientImage()->GetPixel( mappedIndex );
 276
 277       for(unsigned int par=0; par<ParametersDimension; par++) {
 278         RealType sumF = itk::NumericTraits< RealType >::Zero;
 279         RealType sumM = itk::NumericTraits< RealType >::Zero;
 280         for(unsigned int dim=0; dim<dimension; dim++) {
 281           const RealType differential = jacobian( dim, par ) * gradient[dim];
 282           sumF += fixedValue  * differential;
 283           sumM += movingValue * differential;
 284           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 285             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
 286             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
 287           }
 288         }
 289         derivativeF[par] += sumF;
 290         derivativeM[par] += sumM;
 291       }
 292     }
 293
 294     ++ti;
 295   }
 296
 297   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 298     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 299     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 300     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 301   }
 302
 303   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 304
 305   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 306     for(unsigned int i=0; i<ParametersDimension; i++) {
 307       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
 308     }
 309   } else {
 310     for(unsigned int i=0; i<ParametersDimension; i++) {
 311       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 312     }
 313   }
 314
 315 }
 316
 317
 318 /*
 319  * Get both the match Measure and theDerivative Measure
 320  */
 321 template <class TFixedImage, class TMovingImage>
 322 void
 323 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
 324 ::GetValueAndDerivative(const TransformParametersType & parameters,
 325                         MeasureType & value, DerivativeType  & derivative) const
 326 {
 327
 328
 329   if( !this->GetGradientImage() ) {
 330     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 331   }
 332
 333   FixedImageConstPointer fixedImage = this->m_FixedImage;
 334
 335   if( !fixedImage ) {
 336     itkExceptionMacro( << "Fixed image has not been assigned" );
 337   }
 338
 339   const unsigned int dimension = FixedImageType::ImageDimension;
 340
 341   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
 342
 343   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 344
 345   typename FixedImageType::IndexType index;
 346
 347   this->m_NumberOfPixelsCounted = 0;
 348
 349   this->SetTransformParameters( parameters );
 350
 351   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
 352
 353   AccumulateType sff  = itk::NumericTraits< AccumulateType >::Zero;
 354   AccumulateType smm  = itk::NumericTraits< AccumulateType >::Zero;
 355   AccumulateType sfm  = itk::NumericTraits< AccumulateType >::Zero;
 356   AccumulateType sf   = itk::NumericTraits< AccumulateType >::Zero;
 357   AccumulateType sm   = itk::NumericTraits< AccumulateType >::Zero;
 358
 359   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 360   derivative = DerivativeType( ParametersDimension );
 361   derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 362
 363   DerivativeType derivativeF = DerivativeType( ParametersDimension );
 364   derivativeF.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 365
 366   DerivativeType derivativeM = DerivativeType( ParametersDimension );
 367   derivativeM.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 368
 369   DerivativeType derivativeM1 = DerivativeType( ParametersDimension );
 370   derivativeM1.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
 371
 372   ti.GoToBegin();
 373   // First compute the sums
 374   while(!ti.IsAtEnd()) {
 375
 376     index = ti.GetIndex();
 377
 378     InputPointType inputPoint;
 379     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 380
 381     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 382       ++ti;
 383       continue;
 384     }
 385
 386     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 387
 388     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 389       ++ti;
 390       continue;
 391     }
 392
 393     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 394       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 395       const RealType fixedValue   = ti.Get();
 396       sff += fixedValue  * fixedValue;
 397       smm += movingValue * movingValue;
 398       sfm += fixedValue  * movingValue;
 399       if ( this->m_SubtractMean ) {
 400         sf += fixedValue;
 401         sm += movingValue;
 402       }
 403       this->m_NumberOfPixelsCounted++;
 404     }
 405
 406     ++ti;
 407   }
 408
 409
 410   // Compute contributions to derivatives
 411   ti.GoToBegin();
 412   while(!ti.IsAtEnd()) {
 413
 414     index = ti.GetIndex();
 415
 416     InputPointType inputPoint;
 417     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 418
 419     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 420       ++ti;
 421       continue;
 422     }
 423
 424     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 425
 426     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 427       ++ti;
 428       continue;
 429     }
 430
 431     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 432       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 433       const RealType fixedValue     = ti.Get();
 434
 435 #if ITK_VERSION_MAJOR >= 4
 436       TransformJacobianType jacobian;
 437       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 438 #else
 439       const TransformJacobianType & jacobian =
 440         this->m_Transform->GetJacobian( inputPoint );
 441 #endif
 442
 443       // Get the gradient by NearestNeighboorInterpolation:
 444       // which is equivalent to round up the point components.
 445       typedef typename OutputPointType::CoordRepType CoordRepType;
 446       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 447       MovingImageContinuousIndexType;
 448
 449       MovingImageContinuousIndexType tempIndex;
 450       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 451
 452       typename MovingImageType::IndexType mappedIndex;
 453       mappedIndex.CopyWithRound( tempIndex );
 454
 455       const GradientPixelType gradient =
 456         this->GetGradientImage()->GetPixel( mappedIndex );
 457
 458       for(unsigned int par=0; par<ParametersDimension; par++) {
 459         RealType sumF = itk::NumericTraits< RealType >::Zero;
 460         RealType sumM = itk::NumericTraits< RealType >::Zero;
 461         for(unsigned int dim=0; dim<dimension; dim++) {
 462           const RealType differential = jacobian( dim, par ) * gradient[dim];
 463           sumF += fixedValue  * differential;
 464           sumM += movingValue * differential;
 465           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 466             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
 467             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
 468           }
 469         }
 470         derivativeF[par] += sumF;
 471         derivativeM[par] += sumM;
 472       }
 473     }
 474     ++ti;
 475   }
 476
 477   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 478     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 479     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 480     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 481   }
 482
 483   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 484
 485   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 486     for(unsigned int i=0; i<ParametersDimension; i++) {
 487       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
 488     }
 489     value = sfm / denom;
 490   } else {
 491     for(unsigned int i=0; i<ParametersDimension; i++) {
 492       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 493     }
 494     value = itk::NumericTraits< MeasureType >::Zero;
 495   }
 496
 497
 498
 499 }
 500
 501 template < class TFixedImage, class TMovingImage>
 502 void
 503 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
 504 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 505 {
 506   Superclass::PrintSelf(os, indent);
 507   os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
 508 }
 509
 510 } // end namespace itk
 511
 512
 513 #endif // opt
 514
 515 #endif // _clitkNormalizedCorrelationImageToImageMetric.txx