registration/clitkCorrelationRatioImageToImageMetric.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 _clitkCorrelationRatioImageToImageMetric_txx
  20 #define _clitkCorrelationRatioImageToImageMetric_txx
  21
  22 /**
  23  * @file   clitkCorrelationRatioImageToImageMetric.txx
  24  * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
  25  * @date   July 30  18:14:53 2007
  26  *
  27  * @brief  Compute the correlation ratio between 2 images
  28  *
  29  *
  30  */
  31
  32 #include "clitkCorrelationRatioImageToImageMetric.h"
  33 #include "itkImageRegionConstIteratorWithIndex.h"
  34 #include "itkImageRegionConstIterator.h"
  35 #include "itkImageRegionIterator.h"
  36
  37 namespace clitk
  38 {
  39
  40 /*
  41  * Constructor
  42  */
  43 template <class TFixedImage, class TMovingImage>
  44 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
  45 ::CorrelationRatioImageToImageMetric()
  46 {
  47   m_NumberOfBins = 50;
  48
  49 }
  50
  51 template <class TFixedImage, class TMovingImage>
  52 void
  53 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
  54 ::Initialize(void) throw ( ExceptionObject )
  55 {
  56
  57   this->Superclass::Initialize();
  58
  59   // Compute the minimum and maximum for the FixedImage over the FixedImageRegion.
  60   // We can't use StatisticsImageFilter to do this because the filter computes the min/max for the largest possible region
  61   double fixedImageMin = NumericTraits<double>::max();
  62   double fixedImageMax = NumericTraits<double>::NonpositiveMin();
  63
  64   typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
  65   FixedIteratorType fixedImageIterator(
  66     this->m_FixedImage, this->GetFixedImageRegion() );
  67
  68   for ( fixedImageIterator.GoToBegin();
  69         !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
  70
  71     double sample = static_cast<double>( fixedImageIterator.Get() );
  72
  73     if ( sample < fixedImageMin ) {
  74       fixedImageMin = sample;
  75     }
  76
  77     if ( sample > fixedImageMax ) {
  78       fixedImageMax = sample;
  79     }
  80   }
  81
  82   // Compute binsize for the fixedImage
  83   m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
  84   m_FixedImageMin=fixedImageMin;
  85   //Allocate mempry and initialise the fixed image bin
  86   m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
  87   m_mMSVPB.resize( m_NumberOfBins, 0.0 );
  88   m_mSMVPB.resize( m_NumberOfBins, 0.0 );
  89 }
  90
  91
  92 /*
  93  * Get the match Measure
  94  */
  95 template <class TFixedImage, class TMovingImage>
  96 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
  97 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
  98 ::GetValue( const TransformParametersType & parameters ) const
  99 {
 100
 101   itkDebugMacro("GetValue( " << parameters << " ) ");
 102
 103   FixedImageConstPointer fixedImage = this->m_FixedImage;
 104
 105   if( !fixedImage ) {
 106     itkExceptionMacro( << "Fixed image has not been assigned" );
 107   }
 108
 109   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
 110
 111
 112   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 113
 114   typename FixedImageType::IndexType index;
 115
 116   MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
 117
 118   this->m_NumberOfPixelsCounted = 0;
 119   this->SetTransformParameters( parameters );
 120
 121
 122   //temporary measures for the calculation
 123   RealType mSMV=0;
 124   RealType mMSV=0;
 125
 126   while(!ti.IsAtEnd()) {
 127
 128     index = ti.GetIndex();
 129
 130     typename Superclass::InputPointType inputPoint;
 131     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 132
 133     // Verify that the point is in the fixed Image Mask
 134     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 135       ++ti;
 136       continue;
 137     }
 138
 139     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 140
 141     //Verify that the point is in the moving Image Mask
 142     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 143       ++ti;
 144       continue;
 145     }
 146
 147     // Verify is the interpolated value is in the buffer
 148     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 149       //Accumulate calculations for the correlation ratio
 150       //For each pixel the is in both masks and the buffer we adapt the following measures:
 151       //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
 152       //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
 153       //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
 154
 155       //get the value of the moving image
 156       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 157       // for the variance of the overlapping moving image we accumulate the following measures
 158       const RealType movingSquaredValue=movingValue*movingValue;
 159       mMSV+=movingSquaredValue;
 160       mSMV+=movingValue;
 161
 162       //get the fixed value
 163       const RealType fixedValue   = ti.Get();
 164
 165       //check in which bin the fixed value belongs, get the index
 166       const double fixedImageBinTerm =        (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
 167       const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
 168       //adapt the measures per bin
 169       this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
 170       this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
 171       //increase the fixed image bin and the total pixel count
 172       this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
 173       this->m_NumberOfPixelsCounted++;
 174     }
 175
 176     ++ti;
 177   }
 178
 179   if( !this->m_NumberOfPixelsCounted ) {
 180     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 181   } else {
 182
 183     //apdapt the measures per bin
 184     for (unsigned int i=0; i< m_NumberOfBins; i++ ) {
 185       if (this->m_NumberOfPixelsCountedPerBin[i]>0) {
 186         measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
 187       }
 188     }
 189
 190     //Normalize with the global measures
 191     measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
 192     return measure;
 193
 194   }
 195 }
 196
 197
 198
 199
 200
 201 /*
 202  * Get the Derivative Measure
 203  */
 204 template < class TFixedImage, class TMovingImage>
 205 void
 206 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
 207 ::GetDerivative( const TransformParametersType & parameters,
 208                  DerivativeType & derivative  ) const
 209 {
 210
 211   itkDebugMacro("GetDerivative( " << parameters << " ) ");
 212
 213   if( !this->GetGradientImage() ) {
 214     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 215   }
 216
 217   FixedImageConstPointer fixedImage = this->m_FixedImage;
 218
 219   if( !fixedImage ) {
 220     itkExceptionMacro( << "Fixed image has not been assigned" );
 221   }
 222
 223   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 224
 225
 226   typedef  itk::ImageRegionConstIteratorWithIndex<
 227   FixedImageType> FixedIteratorType;
 228
 229   typedef  itk::ImageRegionConstIteratorWithIndex<
 230   ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
 231
 232
 233   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 234
 235   typename FixedImageType::IndexType index;
 236
 237   this->m_NumberOfPixelsCounted = 0;
 238
 239   this->SetTransformParameters( parameters );
 240
 241   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 242   derivative = DerivativeType( ParametersDimension );
 243   derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 244
 245   ti.GoToBegin();
 246
 247   while(!ti.IsAtEnd()) {
 248
 249     index = ti.GetIndex();
 250
 251     typename Superclass::InputPointType inputPoint;
 252     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 253
 254     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 255       ++ti;
 256       continue;
 257     }
 258
 259     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 260
 261     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 262       ++ti;
 263       continue;
 264     }
 265
 266     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 267       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 268
 269       const TransformJacobianType & jacobian =
 270         this->m_Transform->GetJacobian( inputPoint );
 271
 272
 273       const RealType fixedValue     = ti.Value();
 274       this->m_NumberOfPixelsCounted++;
 275       const RealType diff = movingValue - fixedValue;
 276
 277       // Get the gradient by NearestNeighboorInterpolation:
 278       // which is equivalent to round up the point components.
 279       typedef typename Superclass::OutputPointType OutputPointType;
 280       typedef typename OutputPointType::CoordRepType CoordRepType;
 281       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 282       MovingImageContinuousIndexType;
 283
 284       MovingImageContinuousIndexType tempIndex;
 285       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 286
 287       typename MovingImageType::IndexType mappedIndex;
 288       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
 289         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
 290       }
 291
 292       const GradientPixelType gradient =
 293         this->GetGradientImage()->GetPixel( mappedIndex );
 294
 295       for(unsigned int par=0; par<ParametersDimension; par++) {
 296         RealType sum = NumericTraits< RealType >::Zero;
 297         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 298           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 299         }
 300         derivative[par] += sum;
 301       }
 302     }
 303
 304     ++ti;
 305   }
 306
 307   if( !this->m_NumberOfPixelsCounted ) {
 308     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 309   } else {
 310     for(unsigned int i=0; i<ParametersDimension; i++) {
 311       derivative[i] /= this->m_NumberOfPixelsCounted;
 312     }
 313   }
 314
 315 }
 316
 317
 318 /*
 319  * Get both the match Measure and the Derivative Measure
 320  */
 321 template <class TFixedImage, class TMovingImage>
 322 void
 323 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
 324 ::GetValueAndDerivative(const TransformParametersType & parameters,
 325                         MeasureType & value, DerivativeType  & derivative) const
 326 {
 327
 328   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
 329
 330   if( !this->GetGradientImage() ) {
 331     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
 332   }
 333
 334   FixedImageConstPointer fixedImage = this->m_FixedImage;
 335
 336   if( !fixedImage ) {
 337     itkExceptionMacro( << "Fixed image has not been assigned" );
 338   }
 339
 340   const unsigned int ImageDimension = FixedImageType::ImageDimension;
 341
 342   typedef  itk::ImageRegionConstIteratorWithIndex<
 343   FixedImageType> FixedIteratorType;
 344
 345   typedef  itk::ImageRegionConstIteratorWithIndex<
 346   ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
 347
 348
 349   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 350
 351   typename FixedImageType::IndexType index;
 352
 353   MeasureType measure = NumericTraits< MeasureType >::Zero;
 354
 355   this->m_NumberOfPixelsCounted = 0;
 356
 357   this->SetTransformParameters( parameters );
 358
 359   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 360   derivative = DerivativeType( ParametersDimension );
 361   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
 362
 363   ti.GoToBegin();
 364
 365   while(!ti.IsAtEnd()) {
 366
 367     index = ti.GetIndex();
 368
 369     typename Superclass::InputPointType inputPoint;
 370     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 371
 372     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 373       ++ti;
 374       continue;
 375     }
 376
 377     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 378
 379     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 380       ++ti;
 381       continue;
 382     }
 383
 384     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 385       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 386
 387       const TransformJacobianType & jacobian =
 388         this->m_Transform->GetJacobian( inputPoint );
 389
 390
 391       const RealType fixedValue     = ti.Value();
 392       this->m_NumberOfPixelsCounted++;
 393
 394       const RealType diff = movingValue - fixedValue;
 395
 396       measure += diff * diff;
 397
 398       // Get the gradient by NearestNeighboorInterpolation:
 399       // which is equivalent to round up the point components.
 400       typedef typename Superclass::OutputPointType OutputPointType;
 401       typedef typename OutputPointType::CoordRepType CoordRepType;
 402       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 403       MovingImageContinuousIndexType;
 404
 405       MovingImageContinuousIndexType tempIndex;
 406       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 407
 408       typename MovingImageType::IndexType mappedIndex;
 409       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
 410         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
 411       }
 412
 413       const GradientPixelType gradient =
 414         this->GetGradientImage()->GetPixel( mappedIndex );
 415
 416       for(unsigned int par=0; par<ParametersDimension; par++) {
 417         RealType sum = NumericTraits< RealType >::Zero;
 418         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 419           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 420         }
 421         derivative[par] += sum;
 422       }
 423     }
 424
 425     ++ti;
 426   }
 427
 428   if( !this->m_NumberOfPixelsCounted ) {
 429     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 430   } else {
 431     for(unsigned int i=0; i<ParametersDimension; i++) {
 432       derivative[i] /= this->m_NumberOfPixelsCounted;
 433     }
 434     measure /= this->m_NumberOfPixelsCounted;
 435   }
 436
 437   value = measure;
 438
 439 }
 440
 441 } // end namespace itk
 442
 443
 444 #endif
 445