1 #ifndef _clitkCorrelationRatioImageToImageMetric_txx
2 #define _clitkCorrelationRatioImageToImageMetric_txx
5 * @file clitkCorrelationRatioImageToImageMetric.txx
6 * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
7 * @date July 30 18:14:53 2007
9 * @brief Compute the correlation ratio between 2 images
14 #include "clitkCorrelationRatioImageToImageMetric.h"
15 #include "itkImageRegionConstIteratorWithIndex.h"
16 #include "itkImageRegionConstIterator.h"
17 #include "itkImageRegionIterator.h"
25 template <class TFixedImage, class TMovingImage>
26 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
27 ::CorrelationRatioImageToImageMetric()
33 template <class TFixedImage, class TMovingImage>
35 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
36 ::Initialize(void) throw ( ExceptionObject )
39 this->Superclass::Initialize();
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();
46 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
47 FixedIteratorType fixedImageIterator(
48 this->m_FixedImage, this->GetFixedImageRegion() );
50 for ( fixedImageIterator.GoToBegin();
51 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator )
54 double sample = static_cast<double>( fixedImageIterator.Get() );
56 if ( sample < fixedImageMin )
58 fixedImageMin = sample;
61 if ( sample > fixedImageMax )
63 fixedImageMax = sample;
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 );
78 * Get the match Measure
80 template <class TFixedImage, class TMovingImage>
81 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
82 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
83 ::GetValue( const TransformParametersType & parameters ) const
86 itkDebugMacro("GetValue( " << parameters << " ) ");
88 FixedImageConstPointer fixedImage = this->m_FixedImage;
92 itkExceptionMacro( << "Fixed image has not been assigned" );
95 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
98 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
100 typename FixedImageType::IndexType index;
102 MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
104 this->m_NumberOfPixelsCounted = 0;
105 this->SetTransformParameters( parameters );
108 //temporary measures for the calculation
115 index = ti.GetIndex();
117 typename Superclass::InputPointType inputPoint;
118 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
120 // Verify that the point is in the fixed Image Mask
121 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
127 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
129 //Verify that the point is in the moving Image Mask
130 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
136 // Verify is the interpolated value is in the buffer
137 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
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]
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;
152 //get the fixed value
153 const RealType fixedValue = ti.Get();
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++;
169 if( !this->m_NumberOfPixelsCounted )
171 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
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]));
183 //Normalize with the global measures
184 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
195 * Get the Derivative Measure
197 template < class TFixedImage, class TMovingImage>
199 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
200 ::GetDerivative( const TransformParametersType & parameters,
201 DerivativeType & derivative ) const
204 itkDebugMacro("GetDerivative( " << parameters << " ) ");
206 if( !this->GetGradientImage() )
208 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
211 FixedImageConstPointer fixedImage = this->m_FixedImage;
215 itkExceptionMacro( << "Fixed image has not been assigned" );
218 const unsigned int ImageDimension = FixedImageType::ImageDimension;
221 typedef itk::ImageRegionConstIteratorWithIndex<
222 FixedImageType> FixedIteratorType;
224 typedef itk::ImageRegionConstIteratorWithIndex<
225 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
228 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
230 typename FixedImageType::IndexType index;
232 this->m_NumberOfPixelsCounted = 0;
234 this->SetTransformParameters( parameters );
236 const unsigned int ParametersDimension = this->GetNumberOfParameters();
237 derivative = DerivativeType( ParametersDimension );
238 derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
245 index = ti.GetIndex();
247 typename Superclass::InputPointType inputPoint;
248 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
250 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
256 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
258 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
264 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
266 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
268 const TransformJacobianType & jacobian =
269 this->m_Transform->GetJacobian( inputPoint );
272 const RealType fixedValue = ti.Value();
273 this->m_NumberOfPixelsCounted++;
274 const RealType diff = movingValue - fixedValue;
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;
283 MovingImageContinuousIndexType tempIndex;
284 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
286 typename MovingImageType::IndexType mappedIndex;
287 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
289 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
292 const GradientPixelType gradient =
293 this->GetGradientImage()->GetPixel( mappedIndex );
295 for(unsigned int par=0; par<ParametersDimension; par++)
297 RealType sum = NumericTraits< RealType >::Zero;
298 for(unsigned int dim=0; dim<ImageDimension; dim++)
300 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
302 derivative[par] += sum;
309 if( !this->m_NumberOfPixelsCounted )
311 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
315 for(unsigned int i=0; i<ParametersDimension; i++)
317 derivative[i] /= this->m_NumberOfPixelsCounted;
325 * Get both the match Measure and the Derivative Measure
327 template <class TFixedImage, class TMovingImage>
329 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
330 ::GetValueAndDerivative(const TransformParametersType & parameters,
331 MeasureType & value, DerivativeType & derivative) const
334 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
336 if( !this->GetGradientImage() )
338 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
341 FixedImageConstPointer fixedImage = this->m_FixedImage;
345 itkExceptionMacro( << "Fixed image has not been assigned" );
348 const unsigned int ImageDimension = FixedImageType::ImageDimension;
350 typedef itk::ImageRegionConstIteratorWithIndex<
351 FixedImageType> FixedIteratorType;
353 typedef itk::ImageRegionConstIteratorWithIndex<
354 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
357 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
359 typename FixedImageType::IndexType index;
361 MeasureType measure = NumericTraits< MeasureType >::Zero;
363 this->m_NumberOfPixelsCounted = 0;
365 this->SetTransformParameters( parameters );
367 const unsigned int ParametersDimension = this->GetNumberOfParameters();
368 derivative = DerivativeType( ParametersDimension );
369 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
376 index = ti.GetIndex();
378 typename Superclass::InputPointType inputPoint;
379 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
381 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
387 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
389 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
395 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
397 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
399 const TransformJacobianType & jacobian =
400 this->m_Transform->GetJacobian( inputPoint );
403 const RealType fixedValue = ti.Value();
404 this->m_NumberOfPixelsCounted++;
406 const RealType diff = movingValue - fixedValue;
408 measure += diff * diff;
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;
417 MovingImageContinuousIndexType tempIndex;
418 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
420 typename MovingImageType::IndexType mappedIndex;
421 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
423 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
426 const GradientPixelType gradient =
427 this->GetGradientImage()->GetPixel( mappedIndex );
429 for(unsigned int par=0; par<ParametersDimension; par++)
431 RealType sum = NumericTraits< RealType >::Zero;
432 for(unsigned int dim=0; dim<ImageDimension; dim++)
434 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
436 derivative[par] += sum;
443 if( !this->m_NumberOfPixelsCounted )
445 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
449 for(unsigned int i=0; i<ParametersDimension; i++)
451 derivative[i] /= this->m_NumberOfPixelsCounted;
453 measure /= this->m_NumberOfPixelsCounted;
460 } // end namespace itk