tools/clitkCorrelationRatioImageToImageMetric.txx

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