1 /*=========================================================================
6 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
7 l'Image). All rights reserved. See Doc/License.txt or
8 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notices for more information.
14 =========================================================================*/
16 #ifndef _clitkCorrelationRatioImageToImageMetric_txx
17 #define _clitkCorrelationRatioImageToImageMetric_txx
20 * @file clitkCorrelationRatioImageToImageMetric.txx
21 * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
22 * @date July 30 18:14:53 2007
24 * @brief Compute the correlation ratio between 2 images
29 #include "clitkCorrelationRatioImageToImageMetric.h"
30 #include "itkImageRegionConstIteratorWithIndex.h"
31 #include "itkImageRegionConstIterator.h"
32 #include "itkImageRegionIterator.h"
40 template <class TFixedImage, class TMovingImage>
41 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
42 ::CorrelationRatioImageToImageMetric()
48 template <class TFixedImage, class TMovingImage>
50 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
51 ::Initialize(void) throw ( ExceptionObject )
54 this->Superclass::Initialize();
56 // Compute the minimum and maximum for the FixedImage over the FixedImageRegion.
57 // We can't use StatisticsImageFilter to do this because the filter computes the min/max for the largest possible region
58 double fixedImageMin = NumericTraits<double>::max();
59 double fixedImageMax = NumericTraits<double>::NonpositiveMin();
61 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
62 FixedIteratorType fixedImageIterator(
63 this->m_FixedImage, this->GetFixedImageRegion() );
65 for ( fixedImageIterator.GoToBegin();
66 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator )
69 double sample = static_cast<double>( fixedImageIterator.Get() );
71 if ( sample < fixedImageMin )
73 fixedImageMin = sample;
76 if ( sample > fixedImageMax )
78 fixedImageMax = sample;
82 // Compute binsize for the fixedImage
83 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
84 m_FixedImageMin=fixedImageMin;
85 //Allocate mempry and initialise the fixed image bin
86 m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
87 m_mMSVPB.resize( m_NumberOfBins, 0.0 );
88 m_mSMVPB.resize( m_NumberOfBins, 0.0 );
93 * Get the match Measure
95 template <class TFixedImage, class TMovingImage>
96 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
97 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
98 ::GetValue( const TransformParametersType & parameters ) const
101 itkDebugMacro("GetValue( " << parameters << " ) ");
103 FixedImageConstPointer fixedImage = this->m_FixedImage;
107 itkExceptionMacro( << "Fixed image has not been assigned" );
110 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
113 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
115 typename FixedImageType::IndexType index;
117 MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
119 this->m_NumberOfPixelsCounted = 0;
120 this->SetTransformParameters( parameters );
123 //temporary measures for the calculation
130 index = ti.GetIndex();
132 typename Superclass::InputPointType inputPoint;
133 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
135 // Verify that the point is in the fixed Image Mask
136 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
142 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
144 //Verify that the point is in the moving Image Mask
145 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
151 // Verify is the interpolated value is in the buffer
152 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
154 //Accumulate calculations for the correlation ratio
155 //For each pixel the is in both masks and the buffer we adapt the following measures:
156 //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
157 //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
158 //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
160 //get the value of the moving image
161 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
162 // for the variance of the overlapping moving image we accumulate the following measures
163 const RealType movingSquaredValue=movingValue*movingValue;
164 mMSV+=movingSquaredValue;
167 //get the fixed value
168 const RealType fixedValue = ti.Get();
170 //check in which bin the fixed value belongs, get the index
171 const double fixedImageBinTerm = (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
172 const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
173 //adapt the measures per bin
174 this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
175 this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
176 //increase the fixed image bin and the total pixel count
177 this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
178 this->m_NumberOfPixelsCounted++;
184 if( !this->m_NumberOfPixelsCounted )
186 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
191 //apdapt the measures per bin
192 for (unsigned int i=0; i< m_NumberOfBins; i++ ){
193 if (this->m_NumberOfPixelsCountedPerBin[i]>0){
194 measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
198 //Normalize with the global measures
199 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
210 * Get the Derivative Measure
212 template < class TFixedImage, class TMovingImage>
214 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
215 ::GetDerivative( const TransformParametersType & parameters,
216 DerivativeType & derivative ) const
219 itkDebugMacro("GetDerivative( " << parameters << " ) ");
221 if( !this->GetGradientImage() )
223 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
226 FixedImageConstPointer fixedImage = this->m_FixedImage;
230 itkExceptionMacro( << "Fixed image has not been assigned" );
233 const unsigned int ImageDimension = FixedImageType::ImageDimension;
236 typedef itk::ImageRegionConstIteratorWithIndex<
237 FixedImageType> FixedIteratorType;
239 typedef itk::ImageRegionConstIteratorWithIndex<
240 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
243 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
245 typename FixedImageType::IndexType index;
247 this->m_NumberOfPixelsCounted = 0;
249 this->SetTransformParameters( parameters );
251 const unsigned int ParametersDimension = this->GetNumberOfParameters();
252 derivative = DerivativeType( ParametersDimension );
253 derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
260 index = ti.GetIndex();
262 typename Superclass::InputPointType inputPoint;
263 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
265 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
271 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
273 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
279 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
281 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
283 const TransformJacobianType & jacobian =
284 this->m_Transform->GetJacobian( inputPoint );
287 const RealType fixedValue = ti.Value();
288 this->m_NumberOfPixelsCounted++;
289 const RealType diff = movingValue - fixedValue;
291 // Get the gradient by NearestNeighboorInterpolation:
292 // which is equivalent to round up the point components.
293 typedef typename Superclass::OutputPointType OutputPointType;
294 typedef typename OutputPointType::CoordRepType CoordRepType;
295 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
296 MovingImageContinuousIndexType;
298 MovingImageContinuousIndexType tempIndex;
299 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
301 typename MovingImageType::IndexType mappedIndex;
302 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
304 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
307 const GradientPixelType gradient =
308 this->GetGradientImage()->GetPixel( mappedIndex );
310 for(unsigned int par=0; par<ParametersDimension; par++)
312 RealType sum = NumericTraits< RealType >::Zero;
313 for(unsigned int dim=0; dim<ImageDimension; dim++)
315 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
317 derivative[par] += sum;
324 if( !this->m_NumberOfPixelsCounted )
326 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
330 for(unsigned int i=0; i<ParametersDimension; i++)
332 derivative[i] /= this->m_NumberOfPixelsCounted;
340 * Get both the match Measure and the Derivative Measure
342 template <class TFixedImage, class TMovingImage>
344 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
345 ::GetValueAndDerivative(const TransformParametersType & parameters,
346 MeasureType & value, DerivativeType & derivative) const
349 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
351 if( !this->GetGradientImage() )
353 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
356 FixedImageConstPointer fixedImage = this->m_FixedImage;
360 itkExceptionMacro( << "Fixed image has not been assigned" );
363 const unsigned int ImageDimension = FixedImageType::ImageDimension;
365 typedef itk::ImageRegionConstIteratorWithIndex<
366 FixedImageType> FixedIteratorType;
368 typedef itk::ImageRegionConstIteratorWithIndex<
369 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
372 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
374 typename FixedImageType::IndexType index;
376 MeasureType measure = NumericTraits< MeasureType >::Zero;
378 this->m_NumberOfPixelsCounted = 0;
380 this->SetTransformParameters( parameters );
382 const unsigned int ParametersDimension = this->GetNumberOfParameters();
383 derivative = DerivativeType( ParametersDimension );
384 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
391 index = ti.GetIndex();
393 typename Superclass::InputPointType inputPoint;
394 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
396 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
402 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
404 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
410 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
412 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
414 const TransformJacobianType & jacobian =
415 this->m_Transform->GetJacobian( inputPoint );
418 const RealType fixedValue = ti.Value();
419 this->m_NumberOfPixelsCounted++;
421 const RealType diff = movingValue - fixedValue;
423 measure += diff * diff;
425 // Get the gradient by NearestNeighboorInterpolation:
426 // which is equivalent to round up the point components.
427 typedef typename Superclass::OutputPointType OutputPointType;
428 typedef typename OutputPointType::CoordRepType CoordRepType;
429 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
430 MovingImageContinuousIndexType;
432 MovingImageContinuousIndexType tempIndex;
433 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
435 typename MovingImageType::IndexType mappedIndex;
436 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
438 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
441 const GradientPixelType gradient =
442 this->GetGradientImage()->GetPixel( mappedIndex );
444 for(unsigned int par=0; par<ParametersDimension; par++)
446 RealType sum = NumericTraits< RealType >::Zero;
447 for(unsigned int dim=0; dim<ImageDimension; dim++)
449 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
451 derivative[par] += sum;
458 if( !this->m_NumberOfPixelsCounted )
460 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
464 for(unsigned int i=0; i<ParametersDimension; i++)
466 derivative[i] /= this->m_NumberOfPixelsCounted;
468 measure /= this->m_NumberOfPixelsCounted;
475 } // end namespace itk