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://oncora1.lyon.fnclcc.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 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
  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<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 185
 186   DerivativeType derivativeF = DerivativeType( ParametersDimension );
 187   derivativeF.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 188
 189   DerivativeType derivativeM = DerivativeType( ParametersDimension );
 190   derivativeM.Fill( itk::NumericTraits<ITK_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       const TransformJacobianType & jacobian =
 255         this->m_Transform->GetJacobian( inputPoint );
 256
 257       // Get the gradient by NearestNeighboorInterpolation:
 258       // which is equivalent to round up the point components.
 259       typedef typename OutputPointType::CoordRepType CoordRepType;
 260       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 261       MovingImageContinuousIndexType;
 262
 263       MovingImageContinuousIndexType tempIndex;
 264       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 265
 266       typename MovingImageType::IndexType mappedIndex;
 267       mappedIndex.CopyWithRound( tempIndex );
 268
 269       const GradientPixelType gradient =
 270         this->GetGradientImage()->GetPixel( mappedIndex );
 271
 272       for(unsigned int par=0; par<ParametersDimension; par++) {
 273         RealType sumF = itk::NumericTraits< RealType >::Zero;
 274         RealType sumM = itk::NumericTraits< RealType >::Zero;
 275         for(unsigned int dim=0; dim<dimension; dim++) {
 276           const RealType differential = jacobian( dim, par ) * gradient[dim];
 277           sumF += fixedValue  * differential;
 278           sumM += movingValue * differential;
 279           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 280             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
 281             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
 282           }
 283         }
 284         derivativeF[par] += sumF;
 285         derivativeM[par] += sumM;
 286       }
 287     }
 288
 289     ++ti;
 290   }
 291
 292   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 293     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 294     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 295     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 296   }
 297
 298   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 299
 300   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 301     for(unsigned int i=0; i<ParametersDimension; i++) {
 302       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
 303     }
 304   } else {
 305     for(unsigned int i=0; i<ParametersDimension; i++) {
 306       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 307     }
 308   }
 309
 310 }
 311
 312
 313 /*
 314  * Get both the match Measure and theDerivative Measure
 315  */
 316 template <class TFixedImage, class TMovingImage>
 317 void
 318 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
 319 ::GetValueAndDerivative(const TransformParametersType & parameters,
 320                         MeasureType & value, DerivativeType  & derivative) const
 321 {
 322
 323
 324   if( !this->GetGradientImage() ) {
 325     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 326   }
 327
 328   FixedImageConstPointer fixedImage = this->m_FixedImage;
 329
 330   if( !fixedImage ) {
 331     itkExceptionMacro( << "Fixed image has not been assigned" );
 332   }
 333
 334   const unsigned int dimension = FixedImageType::ImageDimension;
 335
 336   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
 337
 338   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 339
 340   typename FixedImageType::IndexType index;
 341
 342   this->m_NumberOfPixelsCounted = 0;
 343
 344   this->SetTransformParameters( parameters );
 345
 346   typedef  typename itk::NumericTraits< MeasureType >::AccumulateType AccumulateType;
 347
 348   AccumulateType sff  = itk::NumericTraits< AccumulateType >::Zero;
 349   AccumulateType smm  = itk::NumericTraits< AccumulateType >::Zero;
 350   AccumulateType sfm  = itk::NumericTraits< AccumulateType >::Zero;
 351   AccumulateType sf   = itk::NumericTraits< AccumulateType >::Zero;
 352   AccumulateType sm   = itk::NumericTraits< AccumulateType >::Zero;
 353
 354   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 355   derivative = DerivativeType( ParametersDimension );
 356   derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 357
 358   DerivativeType derivativeF = DerivativeType( ParametersDimension );
 359   derivativeF.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 360
 361   DerivativeType derivativeM = DerivativeType( ParametersDimension );
 362   derivativeM.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 363
 364   DerivativeType derivativeM1 = DerivativeType( ParametersDimension );
 365   derivativeM1.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 366
 367   ti.GoToBegin();
 368   // First compute the sums
 369   while(!ti.IsAtEnd()) {
 370
 371     index = ti.GetIndex();
 372
 373     InputPointType inputPoint;
 374     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 375
 376     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 377       ++ti;
 378       continue;
 379     }
 380
 381     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 382
 383     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 384       ++ti;
 385       continue;
 386     }
 387
 388     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 389       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 390       const RealType fixedValue   = ti.Get();
 391       sff += fixedValue  * fixedValue;
 392       smm += movingValue * movingValue;
 393       sfm += fixedValue  * movingValue;
 394       if ( this->m_SubtractMean ) {
 395         sf += fixedValue;
 396         sm += movingValue;
 397       }
 398       this->m_NumberOfPixelsCounted++;
 399     }
 400
 401     ++ti;
 402   }
 403
 404
 405   // Compute contributions to derivatives
 406   ti.GoToBegin();
 407   while(!ti.IsAtEnd()) {
 408
 409     index = ti.GetIndex();
 410
 411     InputPointType inputPoint;
 412     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 413
 414     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 415       ++ti;
 416       continue;
 417     }
 418
 419     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 420
 421     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 422       ++ti;
 423       continue;
 424     }
 425
 426     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 427       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 428       const RealType fixedValue     = ti.Get();
 429
 430       const TransformJacobianType & jacobian =
 431         this->m_Transform->GetJacobian( inputPoint );
 432
 433       // Get the gradient by NearestNeighboorInterpolation:
 434       // which is equivalent to round up the point components.
 435       typedef typename OutputPointType::CoordRepType CoordRepType;
 436       typedef itk::ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 437       MovingImageContinuousIndexType;
 438
 439       MovingImageContinuousIndexType tempIndex;
 440       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 441
 442       typename MovingImageType::IndexType mappedIndex;
 443       mappedIndex.CopyWithRound( tempIndex );
 444
 445       const GradientPixelType gradient =
 446         this->GetGradientImage()->GetPixel( mappedIndex );
 447
 448       for(unsigned int par=0; par<ParametersDimension; par++) {
 449         RealType sumF = itk::NumericTraits< RealType >::Zero;
 450         RealType sumM = itk::NumericTraits< RealType >::Zero;
 451         for(unsigned int dim=0; dim<dimension; dim++) {
 452           const RealType differential = jacobian( dim, par ) * gradient[dim];
 453           sumF += fixedValue  * differential;
 454           sumM += movingValue * differential;
 455           if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 456             sumF -= differential * sf / this->m_NumberOfPixelsCounted;
 457             sumM -= differential * sm / this->m_NumberOfPixelsCounted;
 458           }
 459         }
 460         derivativeF[par] += sumF;
 461         derivativeM[par] += sumM;
 462       }
 463     }
 464     ++ti;
 465   }
 466
 467   if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
 468     sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
 469     smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
 470     sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
 471   }
 472
 473   const RealType denom = -1.0 * vcl_sqrt(sff * smm );
 474
 475   if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
 476     for(unsigned int i=0; i<ParametersDimension; i++) {
 477       derivative[i] = ( derivativeF[i] - (sfm/smm)* derivativeM[i] ) / denom;
 478     }
 479     value = sfm / denom;
 480   } else {
 481     for(unsigned int i=0; i<ParametersDimension; i++) {
 482       derivative[i] = itk::NumericTraits< MeasureType >::Zero;
 483     }
 484     value = itk::NumericTraits< MeasureType >::Zero;
 485   }
 486
 487
 488
 489 }
 490
 491 template < class TFixedImage, class TMovingImage>
 492 void
 493 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
 494 ::PrintSelf(std::ostream& os, itk::Indent indent) const
 495 {
 496   Superclass::PrintSelf(os, indent);
 497   os << indent << "SubtractMean: " << m_SubtractMean << std::endl;
 498 }
 499
 500 } // end namespace itk
 501
 502
 503 #endif // opt
 504
 505 #endif // _clitkNormalizedCorrelationImageToImageMetric.txx