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