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