tools/clitkCorrelationRatioImageToImageMetric.txx

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