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://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 _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)
  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   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<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       TransformJacobianType jacobian;
 270       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint , jacobian);
 271
 272       const RealType fixedValue     = ti.Value();
 273       this->m_NumberOfPixelsCounted++;
 274       const RealType diff = movingValue - fixedValue;
 275
 276       // Get the gradient by NearestNeighboorInterpolation:
 277       // which is equivalent to round up the point components.
 278       typedef typename Superclass::OutputPointType OutputPointType;
 279       typedef typename OutputPointType::CoordRepType CoordRepType;
 280       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 281       MovingImageContinuousIndexType;
 282
 283       MovingImageContinuousIndexType tempIndex;
 284       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 285
 286       typename MovingImageType::IndexType mappedIndex;
 287       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
 288         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
 289       }
 290
 291       const GradientPixelType gradient =
 292         this->GetGradientImage()->GetPixel( mappedIndex );
 293
 294       for(unsigned int par=0; par<ParametersDimension; par++) {
 295         RealType sum = NumericTraits< RealType >::Zero;
 296         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 297           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 298         }
 299         derivative[par] += sum;
 300       }
 301     }
 302
 303     ++ti;
 304   }
 305
 306   if( !this->m_NumberOfPixelsCounted ) {
 307     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 308   } else {
 309     for(unsigned int i=0; i<ParametersDimension; i++) {
 310       derivative[i] /= this->m_NumberOfPixelsCounted;
 311     }
 312   }
 313
 314 }
 315
 316
 317 /*
 318  * Get both the match Measure and the Derivative Measure
 319  */
 320 template <class TFixedImage, class TMovingImage>
 321 void
 322 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
 323 ::GetValueAndDerivative(const TransformParametersType & parameters,
 324                         MeasureType & value, DerivativeType  & derivative) const
 325 {
 326
 327   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
 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 ImageDimension = FixedImageType::ImageDimension;
 340
 341   typedef  itk::ImageRegionConstIteratorWithIndex<
 342   FixedImageType> FixedIteratorType;
 343
 344   typedef  itk::ImageRegionConstIteratorWithIndex<
 345   typename Superclass::GradientImageType> GradientIteratorType;
 346
 347
 348   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
 349
 350   typename FixedImageType::IndexType index;
 351
 352   MeasureType measure = NumericTraits< MeasureType >::Zero;
 353
 354   this->m_NumberOfPixelsCounted = 0;
 355
 356   this->SetTransformParameters( parameters );
 357
 358   const unsigned int ParametersDimension = this->GetNumberOfParameters();
 359   derivative = DerivativeType( ParametersDimension );
 360   derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
 361
 362   ti.GoToBegin();
 363
 364   while(!ti.IsAtEnd()) {
 365
 366     index = ti.GetIndex();
 367
 368     typename Superclass::InputPointType inputPoint;
 369     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 370
 371     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
 372       ++ti;
 373       continue;
 374     }
 375
 376     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
 377
 378     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
 379       ++ti;
 380       continue;
 381     }
 382
 383     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
 384       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
 385
 386       TransformJacobianType jacobian;
 387         this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
 388
 389       const RealType fixedValue     = ti.Value();
 390       this->m_NumberOfPixelsCounted++;
 391
 392       const RealType diff = movingValue - fixedValue;
 393
 394       measure += diff * diff;
 395
 396       // Get the gradient by NearestNeighboorInterpolation:
 397       // which is equivalent to round up the point components.
 398       typedef typename Superclass::OutputPointType OutputPointType;
 399       typedef typename OutputPointType::CoordRepType CoordRepType;
 400       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
 401       MovingImageContinuousIndexType;
 402
 403       MovingImageContinuousIndexType tempIndex;
 404       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
 405
 406       typename MovingImageType::IndexType mappedIndex;
 407       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
 408         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
 409       }
 410
 411       const GradientPixelType gradient =
 412         this->GetGradientImage()->GetPixel( mappedIndex );
 413
 414       for(unsigned int par=0; par<ParametersDimension; par++) {
 415         RealType sum = NumericTraits< RealType >::Zero;
 416         for(unsigned int dim=0; dim<ImageDimension; dim++) {
 417           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
 418         }
 419         derivative[par] += sum;
 420       }
 421     }
 422
 423     ++ti;
 424   }
 425
 426   if( !this->m_NumberOfPixelsCounted ) {
 427     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
 428   } else {
 429     for(unsigned int i=0; i<ParametersDimension; i++) {
 430       derivative[i] /= this->m_NumberOfPixelsCounted;
 431     }
 432     measure /= this->m_NumberOfPixelsCounted;
 433   }
 434
 435   value = measure;
 436
 437 }
 438
 439 } // end namespace itk
 440
 441
 442 #endif
 443