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://oncora1.lyon.fnclcc.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 ======================================================================-====*/
18 #ifndef _clitkCorrelationRatioImageToImageMetric_txx
19 #define _clitkCorrelationRatioImageToImageMetric_txx
22 * @file clitkCorrelationRatioImageToImageMetric.txx
23 * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
24 * @date July 30 18:14:53 2007
26 * @brief Compute the correlation ratio between 2 images
31 #include "clitkCorrelationRatioImageToImageMetric.h"
32 #include "itkImageRegionConstIteratorWithIndex.h"
33 #include "itkImageRegionConstIterator.h"
34 #include "itkImageRegionIterator.h"
42 template <class TFixedImage, class TMovingImage>
43 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
44 ::CorrelationRatioImageToImageMetric()
50 template <class TFixedImage, class TMovingImage>
52 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
53 ::Initialize(void) throw ( ExceptionObject )
56 this->Superclass::Initialize();
58 // Compute the minimum and maximum for the FixedImage over the FixedImageRegion.
59 // We can't use StatisticsImageFilter to do this because the filter computes the min/max for the largest possible region
60 double fixedImageMin = NumericTraits<double>::max();
61 double fixedImageMax = NumericTraits<double>::NonpositiveMin();
63 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
64 FixedIteratorType fixedImageIterator(
65 this->m_FixedImage, this->GetFixedImageRegion() );
67 for ( fixedImageIterator.GoToBegin();
68 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator )
71 double sample = static_cast<double>( fixedImageIterator.Get() );
73 if ( sample < fixedImageMin )
75 fixedImageMin = sample;
78 if ( sample > fixedImageMax )
80 fixedImageMax = sample;
84 // Compute binsize for the fixedImage
85 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
86 m_FixedImageMin=fixedImageMin;
87 //Allocate mempry and initialise the fixed image bin
88 m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
89 m_mMSVPB.resize( m_NumberOfBins, 0.0 );
90 m_mSMVPB.resize( m_NumberOfBins, 0.0 );
95 * Get the match Measure
97 template <class TFixedImage, class TMovingImage>
98 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
99 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
100 ::GetValue( const TransformParametersType & parameters ) const
103 itkDebugMacro("GetValue( " << parameters << " ) ");
105 FixedImageConstPointer fixedImage = this->m_FixedImage;
109 itkExceptionMacro( << "Fixed image has not been assigned" );
112 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
115 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
117 typename FixedImageType::IndexType index;
119 MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
121 this->m_NumberOfPixelsCounted = 0;
122 this->SetTransformParameters( parameters );
125 //temporary measures for the calculation
132 index = ti.GetIndex();
134 typename Superclass::InputPointType inputPoint;
135 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
137 // Verify that the point is in the fixed Image Mask
138 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
144 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
146 //Verify that the point is in the moving Image Mask
147 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
153 // Verify is the interpolated value is in the buffer
154 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
156 //Accumulate calculations for the correlation ratio
157 //For each pixel the is in both masks and the buffer we adapt the following measures:
158 //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
159 //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
160 //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
162 //get the value of the moving image
163 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
164 // for the variance of the overlapping moving image we accumulate the following measures
165 const RealType movingSquaredValue=movingValue*movingValue;
166 mMSV+=movingSquaredValue;
169 //get the fixed value
170 const RealType fixedValue = ti.Get();
172 //check in which bin the fixed value belongs, get the index
173 const double fixedImageBinTerm = (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
174 const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
175 //adapt the measures per bin
176 this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
177 this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
178 //increase the fixed image bin and the total pixel count
179 this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
180 this->m_NumberOfPixelsCounted++;
186 if( !this->m_NumberOfPixelsCounted )
188 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
193 //apdapt the measures per bin
194 for (unsigned int i=0; i< m_NumberOfBins; i++ ){
195 if (this->m_NumberOfPixelsCountedPerBin[i]>0){
196 measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
200 //Normalize with the global measures
201 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
212 * Get the Derivative Measure
214 template < class TFixedImage, class TMovingImage>
216 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
217 ::GetDerivative( const TransformParametersType & parameters,
218 DerivativeType & derivative ) const
221 itkDebugMacro("GetDerivative( " << parameters << " ) ");
223 if( !this->GetGradientImage() )
225 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
228 FixedImageConstPointer fixedImage = this->m_FixedImage;
232 itkExceptionMacro( << "Fixed image has not been assigned" );
235 const unsigned int ImageDimension = FixedImageType::ImageDimension;
238 typedef itk::ImageRegionConstIteratorWithIndex<
239 FixedImageType> FixedIteratorType;
241 typedef itk::ImageRegionConstIteratorWithIndex<
242 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
245 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
247 typename FixedImageType::IndexType index;
249 this->m_NumberOfPixelsCounted = 0;
251 this->SetTransformParameters( parameters );
253 const unsigned int ParametersDimension = this->GetNumberOfParameters();
254 derivative = DerivativeType( ParametersDimension );
255 derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
262 index = ti.GetIndex();
264 typename Superclass::InputPointType inputPoint;
265 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
267 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
273 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
275 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
281 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
283 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
285 const TransformJacobianType & jacobian =
286 this->m_Transform->GetJacobian( inputPoint );
289 const RealType fixedValue = ti.Value();
290 this->m_NumberOfPixelsCounted++;
291 const RealType diff = movingValue - fixedValue;
293 // Get the gradient by NearestNeighboorInterpolation:
294 // which is equivalent to round up the point components.
295 typedef typename Superclass::OutputPointType OutputPointType;
296 typedef typename OutputPointType::CoordRepType CoordRepType;
297 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
298 MovingImageContinuousIndexType;
300 MovingImageContinuousIndexType tempIndex;
301 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
303 typename MovingImageType::IndexType mappedIndex;
304 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
306 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
309 const GradientPixelType gradient =
310 this->GetGradientImage()->GetPixel( mappedIndex );
312 for(unsigned int par=0; par<ParametersDimension; par++)
314 RealType sum = NumericTraits< RealType >::Zero;
315 for(unsigned int dim=0; dim<ImageDimension; dim++)
317 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
319 derivative[par] += sum;
326 if( !this->m_NumberOfPixelsCounted )
328 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
332 for(unsigned int i=0; i<ParametersDimension; i++)
334 derivative[i] /= this->m_NumberOfPixelsCounted;
342 * Get both the match Measure and the Derivative Measure
344 template <class TFixedImage, class TMovingImage>
346 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
347 ::GetValueAndDerivative(const TransformParametersType & parameters,
348 MeasureType & value, DerivativeType & derivative) const
351 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
353 if( !this->GetGradientImage() )
355 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
358 FixedImageConstPointer fixedImage = this->m_FixedImage;
362 itkExceptionMacro( << "Fixed image has not been assigned" );
365 const unsigned int ImageDimension = FixedImageType::ImageDimension;
367 typedef itk::ImageRegionConstIteratorWithIndex<
368 FixedImageType> FixedIteratorType;
370 typedef itk::ImageRegionConstIteratorWithIndex<
371 ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
374 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
376 typename FixedImageType::IndexType index;
378 MeasureType measure = NumericTraits< MeasureType >::Zero;
380 this->m_NumberOfPixelsCounted = 0;
382 this->SetTransformParameters( parameters );
384 const unsigned int ParametersDimension = this->GetNumberOfParameters();
385 derivative = DerivativeType( ParametersDimension );
386 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
393 index = ti.GetIndex();
395 typename Superclass::InputPointType inputPoint;
396 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
398 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) )
404 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
406 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) )
412 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
414 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
416 const TransformJacobianType & jacobian =
417 this->m_Transform->GetJacobian( inputPoint );
420 const RealType fixedValue = ti.Value();
421 this->m_NumberOfPixelsCounted++;
423 const RealType diff = movingValue - fixedValue;
425 measure += diff * diff;
427 // Get the gradient by NearestNeighboorInterpolation:
428 // which is equivalent to round up the point components.
429 typedef typename Superclass::OutputPointType OutputPointType;
430 typedef typename OutputPointType::CoordRepType CoordRepType;
431 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
432 MovingImageContinuousIndexType;
434 MovingImageContinuousIndexType tempIndex;
435 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
437 typename MovingImageType::IndexType mappedIndex;
438 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
440 mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
443 const GradientPixelType gradient =
444 this->GetGradientImage()->GetPixel( mappedIndex );
446 for(unsigned int par=0; par<ParametersDimension; par++)
448 RealType sum = NumericTraits< RealType >::Zero;
449 for(unsigned int dim=0; dim<ImageDimension; dim++)
451 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
453 derivative[par] += sum;
460 if( !this->m_NumberOfPixelsCounted )
462 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
466 for(unsigned int i=0; i<ParametersDimension; i++)
468 derivative[i] /= this->m_NumberOfPixelsCounted;
470 measure /= this->m_NumberOfPixelsCounted;
477 } // end namespace itk