1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
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 ) {
71 double sample = static_cast<double>( fixedImageIterator.Get() );
73 if ( sample < fixedImageMin ) {
74 fixedImageMin = sample;
77 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;
106 itkExceptionMacro( << "Fixed image has not been assigned" );
109 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
112 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
114 typename FixedImageType::IndexType index;
116 MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
118 this->m_NumberOfPixelsCounted = 0;
119 this->SetTransformParameters( parameters );
122 //temporary measures for the calculation
126 while(!ti.IsAtEnd()) {
128 index = ti.GetIndex();
130 typename Superclass::InputPointType inputPoint;
131 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
133 // Verify that the point is in the fixed Image Mask
134 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
139 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
141 //Verify that the point is in the moving Image Mask
142 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
147 // Verify is the interpolated value is in the buffer
148 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
149 //Accumulate calculations for the correlation ratio
150 //For each pixel the is in both masks and the buffer we adapt the following measures:
151 //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
152 //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
153 //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
155 //get the value of the moving image
156 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
157 // for the variance of the overlapping moving image we accumulate the following measures
158 const RealType movingSquaredValue=movingValue*movingValue;
159 mMSV+=movingSquaredValue;
162 //get the fixed value
163 const RealType fixedValue = ti.Get();
165 //check in which bin the fixed value belongs, get the index
166 const double fixedImageBinTerm = (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
167 const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
168 //adapt the measures per bin
169 this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
170 this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
171 //increase the fixed image bin and the total pixel count
172 this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
173 this->m_NumberOfPixelsCounted++;
179 if( !this->m_NumberOfPixelsCounted ) {
180 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
183 //apdapt the measures per bin
184 for (unsigned int i=0; i< m_NumberOfBins; i++ ) {
185 if (this->m_NumberOfPixelsCountedPerBin[i]>0) {
186 measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
190 //Normalize with the global measures
191 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
202 * Get the Derivative Measure
204 template < class TFixedImage, class TMovingImage>
206 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
207 ::GetDerivative( const TransformParametersType & parameters,
208 DerivativeType & derivative ) const
211 itkDebugMacro("GetDerivative( " << parameters << " ) ");
213 if( !this->GetGradientImage() ) {
214 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
217 FixedImageConstPointer fixedImage = this->m_FixedImage;
220 itkExceptionMacro( << "Fixed image has not been assigned" );
223 const unsigned int ImageDimension = FixedImageType::ImageDimension;
226 typedef itk::ImageRegionConstIteratorWithIndex<
227 FixedImageType> FixedIteratorType;
229 typedef itk::ImageRegionConstIteratorWithIndex<
230 typename Superclass::GradientImageType> GradientIteratorType;
233 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
235 typename FixedImageType::IndexType index;
237 this->m_NumberOfPixelsCounted = 0;
239 this->SetTransformParameters( parameters );
241 const unsigned int ParametersDimension = this->GetNumberOfParameters();
242 derivative = DerivativeType( ParametersDimension );
243 derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
247 while(!ti.IsAtEnd()) {
249 index = ti.GetIndex();
251 typename Superclass::InputPointType inputPoint;
252 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
254 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
259 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
261 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
266 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
267 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
269 TransformJacobianType jacobian;
270 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint , jacobian);
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++ ) {
288 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
291 const GradientPixelType gradient =
292 this->GetGradientImage()->GetPixel( mappedIndex );
294 for(unsigned int par=0; par<ParametersDimension; par++) {
295 RealType sum = NumericTraits< RealType >::Zero;
296 for(unsigned int dim=0; dim<ImageDimension; dim++) {
297 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
299 derivative[par] += sum;
306 if( !this->m_NumberOfPixelsCounted ) {
307 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
309 for(unsigned int i=0; i<ParametersDimension; i++) {
310 derivative[i] /= this->m_NumberOfPixelsCounted;
318 * Get both the match Measure and the Derivative Measure
320 template <class TFixedImage, class TMovingImage>
322 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
323 ::GetValueAndDerivative(const TransformParametersType & parameters,
324 MeasureType & value, DerivativeType & derivative) const
327 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
329 if( !this->GetGradientImage() ) {
330 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
333 FixedImageConstPointer fixedImage = this->m_FixedImage;
336 itkExceptionMacro( << "Fixed image has not been assigned" );
339 const unsigned int ImageDimension = FixedImageType::ImageDimension;
341 typedef itk::ImageRegionConstIteratorWithIndex<
342 FixedImageType> FixedIteratorType;
344 typedef itk::ImageRegionConstIteratorWithIndex<
345 typename Superclass::GradientImageType> GradientIteratorType;
348 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
350 typename FixedImageType::IndexType index;
352 MeasureType measure = NumericTraits< MeasureType >::Zero;
354 this->m_NumberOfPixelsCounted = 0;
356 this->SetTransformParameters( parameters );
358 const unsigned int ParametersDimension = this->GetNumberOfParameters();
359 derivative = DerivativeType( ParametersDimension );
360 derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
364 while(!ti.IsAtEnd()) {
366 index = ti.GetIndex();
368 typename Superclass::InputPointType inputPoint;
369 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
371 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
376 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
378 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
383 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
384 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
386 TransformJacobianType jacobian;
387 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
389 const RealType fixedValue = ti.Value();
390 this->m_NumberOfPixelsCounted++;
392 const RealType diff = movingValue - fixedValue;
394 measure += diff * diff;
396 // Get the gradient by NearestNeighboorInterpolation:
397 // which is equivalent to round up the point components.
398 typedef typename Superclass::OutputPointType OutputPointType;
399 typedef typename OutputPointType::CoordRepType CoordRepType;
400 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
401 MovingImageContinuousIndexType;
403 MovingImageContinuousIndexType tempIndex;
404 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
406 typename MovingImageType::IndexType mappedIndex;
407 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
408 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
411 const GradientPixelType gradient =
412 this->GetGradientImage()->GetPixel( mappedIndex );
414 for(unsigned int par=0; par<ParametersDimension; par++) {
415 RealType sum = NumericTraits< RealType >::Zero;
416 for(unsigned int dim=0; dim<ImageDimension; dim++) {
417 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
419 derivative[par] += sum;
426 if( !this->m_NumberOfPixelsCounted ) {
427 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
429 for(unsigned int i=0; i<ParametersDimension; i++) {
430 derivative[i] /= this->m_NumberOfPixelsCounted;
432 measure /= this->m_NumberOfPixelsCounted;
439 } // end namespace itk