1 /*=========================================================================
4 Module: $RCSfile: clitkCorrelationRatioImageToImageMetric.h
6 Date: $Date: 2010/01/06 13:31:57 $
7 Version: $Revision: 1.1 $
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.
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.
17 =========================================================================*/
19 #ifndef _clitkCorrelationRatioImageToImageMetric_txx
20 #define _clitkCorrelationRatioImageToImageMetric_txx
23 * @file clitkCorrelationRatioImageToImageMetric.txx
24 * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
25 * @date July 30 18:14:53 2007
27 * @brief Compute the correlation ratio between 2 images
32 #include "clitkCorrelationRatioImageToImageMetric.h"
33 #include "itkImageRegionConstIteratorWithIndex.h"
34 #include "itkImageRegionConstIterator.h"
35 #include "itkImageRegionIterator.h"
43 template <class TFixedImage, class TMovingImage>
44 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
45 ::CorrelationRatioImageToImageMetric()
51 template <class TFixedImage, class TMovingImage>
53 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
54 ::Initialize(void) throw ( ExceptionObject )
57 this->Superclass::Initialize();
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();
64 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
65 FixedIteratorType fixedImageIterator(
66 this->m_FixedImage, this->GetFixedImageRegion() );
68 for ( fixedImageIterator.GoToBegin();
69 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator )
72 double sample = static_cast<double>( fixedImageIterator.Get() );
74 if ( sample < fixedImageMin )
76 fixedImageMin = sample;
79 if ( sample > fixedImageMax )
81 fixedImageMax = sample;
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 );
96 * Get the match Measure
98 template <class TFixedImage, class TMovingImage>
99 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
100 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
101 ::GetValue( const TransformParametersType & parameters ) const
104 itkDebugMacro("GetValue( " << parameters << " ) ");
106 FixedImageConstPointer fixedImage = this->m_FixedImage;
110 itkExceptionMacro( << "Fixed image has not been assigned" );
113 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
116 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
118 typename FixedImageType::IndexType index;
120 MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
122 this->m_NumberOfPixelsCounted = 0;
123 this->SetTransformParameters( parameters );
126 //temporary measures for the calculation
133 index = ti.GetIndex();
135 typename Superclass::InputPointType inputPoint;
136 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
138 // Verify that the point is in the fixed Image Mask
139 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
145 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
147 //Verify that the point is in the moving Image Mask
148 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
154 // Verify is the interpolated value is in the buffer
155 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
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]
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;
170 //get the fixed value
171 const RealType fixedValue = ti.Get();
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++;
187 if( !this->m_NumberOfPixelsCounted )
189 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
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]));
201 //Normalize with the global measures
202 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
213 * Get the Derivative Measure
215 template < class TFixedImage, class TMovingImage>
217 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
218 ::GetDerivative( const TransformParametersType & parameters,
219 DerivativeType & derivative ) const
222 itkDebugMacro("GetDerivative( " << parameters << " ) ");
224 if( !this->GetGradientImage() )
226 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
229 FixedImageConstPointer fixedImage = this->m_FixedImage;
233 itkExceptionMacro( << "Fixed image has not been assigned" );
236 const unsigned int ImageDimension = FixedImageType::ImageDimension;
239 typedef itk::ImageRegionConstIteratorWithIndex<
240 FixedImageType> FixedIteratorType;
242 typedef itk::ImageRegionConstIteratorWithIndex<
243 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
246 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
248 typename FixedImageType::IndexType index;
250 this->m_NumberOfPixelsCounted = 0;
252 this->SetTransformParameters( parameters );
254 const unsigned int ParametersDimension = this->GetNumberOfParameters();
255 derivative = DerivativeType( ParametersDimension );
256 derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
263 index = ti.GetIndex();
265 typename Superclass::InputPointType inputPoint;
266 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
268 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
274 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
276 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
282 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
284 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
286 const TransformJacobianType & jacobian =
287 this->m_Transform->GetJacobian( inputPoint );
290 const RealType fixedValue = ti.Value();
291 this->m_NumberOfPixelsCounted++;
292 const RealType diff = movingValue - fixedValue;
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;
301 MovingImageContinuousIndexType tempIndex;
302 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
304 typename MovingImageType::IndexType mappedIndex;
305 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
307 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
310 const GradientPixelType gradient =
311 this->GetGradientImage()->GetPixel( mappedIndex );
313 for(unsigned int par=0; par<ParametersDimension; par++)
315 RealType sum = NumericTraits< RealType >::Zero;
316 for(unsigned int dim=0; dim<ImageDimension; dim++)
318 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
320 derivative[par] += sum;
327 if( !this->m_NumberOfPixelsCounted )
329 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
333 for(unsigned int i=0; i<ParametersDimension; i++)
335 derivative[i] /= this->m_NumberOfPixelsCounted;
343 * Get both the match Measure and the Derivative Measure
345 template <class TFixedImage, class TMovingImage>
347 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
348 ::GetValueAndDerivative(const TransformParametersType & parameters,
349 MeasureType & value, DerivativeType & derivative) const
352 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
354 if( !this->GetGradientImage() )
356 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
359 FixedImageConstPointer fixedImage = this->m_FixedImage;
363 itkExceptionMacro( << "Fixed image has not been assigned" );
366 const unsigned int ImageDimension = FixedImageType::ImageDimension;
368 typedef itk::ImageRegionConstIteratorWithIndex<
369 FixedImageType> FixedIteratorType;
371 typedef itk::ImageRegionConstIteratorWithIndex<
372 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
375 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
377 typename FixedImageType::IndexType index;
379 MeasureType measure = NumericTraits< MeasureType >::Zero;
381 this->m_NumberOfPixelsCounted = 0;
383 this->SetTransformParameters( parameters );
385 const unsigned int ParametersDimension = this->GetNumberOfParameters();
386 derivative = DerivativeType( ParametersDimension );
387 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
394 index = ti.GetIndex();
396 typename Superclass::InputPointType inputPoint;
397 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
399 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
405 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
407 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
413 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
415 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
417 const TransformJacobianType & jacobian =
418 this->m_Transform->GetJacobian( inputPoint );
421 const RealType fixedValue = ti.Value();
422 this->m_NumberOfPixelsCounted++;
424 const RealType diff = movingValue - fixedValue;
426 measure += diff * diff;
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;
435 MovingImageContinuousIndexType tempIndex;
436 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
438 typename MovingImageType::IndexType mappedIndex;
439 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
441 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
444 const GradientPixelType gradient =
445 this->GetGradientImage()->GetPixel( mappedIndex );
447 for(unsigned int par=0; par<ParametersDimension; par++)
449 RealType sum = NumericTraits< RealType >::Zero;
450 for(unsigned int dim=0; dim<ImageDimension; dim++)
452 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
454 derivative[par] += sum;
461 if( !this->m_NumberOfPixelsCounted )
463 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
467 for(unsigned int i=0; i<ParametersDimension; i++)
469 derivative[i] /= this->m_NumberOfPixelsCounted;
471 measure /= this->m_NumberOfPixelsCounted;
478 } // end namespace itk