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