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