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 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
57 ::Initialize(void) throw ( ExceptionObject )
61 this->Superclass::Initialize();
63 // Compute the minimum and maximum for the FixedImage over the FixedImageRegion.
64 // We can't use StatisticsImageFilter to do this because the filter computes the min/max for the largest possible region
65 double fixedImageMin = NumericTraits<double>::max();
66 double fixedImageMax = NumericTraits<double>::NonpositiveMin();
68 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
69 FixedIteratorType fixedImageIterator(
70 this->m_FixedImage, this->GetFixedImageRegion() );
72 for ( fixedImageIterator.GoToBegin();
73 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
75 double sample = static_cast<double>( fixedImageIterator.Get() );
77 if ( sample < fixedImageMin ) {
78 fixedImageMin = sample;
81 if ( sample > fixedImageMax ) {
82 fixedImageMax = sample;
86 // Compute binsize for the fixedImage
87 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
88 m_FixedImageMin=fixedImageMin;
89 //Allocate mempry and initialise the fixed image bin
90 m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
91 m_mMSVPB.resize( m_NumberOfBins, 0.0 );
92 m_mSMVPB.resize( m_NumberOfBins, 0.0 );
97 * Get the match Measure
99 template <class TFixedImage, class TMovingImage>
100 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
101 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
102 ::GetValue( const TransformParametersType & parameters ) const
105 itkDebugMacro("GetValue( " << parameters << " ) ");
107 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
130 while(!ti.IsAtEnd()) {
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 ) ) {
143 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
145 //Verify that the point is in the moving Image Mask
146 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 ) ) {
153 //Accumulate calculations for the correlation ratio
154 //For each pixel the is in both masks and the buffer we adapt the following measures:
155 //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
156 //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
157 //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
159 //get the value of the moving image
160 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
161 // for the variance of the overlapping moving image we accumulate the following measures
162 const RealType movingSquaredValue=movingValue*movingValue;
163 mMSV+=movingSquaredValue;
166 //get the fixed value
167 const RealType fixedValue = ti.Get();
169 //check in which bin the fixed value belongs, get the index
170 const double fixedImageBinTerm = (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
171 const unsigned int fixedImageBinIndex = static_cast<unsigned int>( std::floor(fixedImageBinTerm ) );
172 //adapt the measures per bin
173 this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
174 this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
175 //increase the fixed image bin and the total pixel count
176 this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
177 this->m_NumberOfPixelsCounted++;
183 if( !this->m_NumberOfPixelsCounted ) {
184 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
187 //apdapt the measures per bin
188 for (unsigned int i=0; i< m_NumberOfBins; i++ ) {
189 if (this->m_NumberOfPixelsCountedPerBin[i]>0) {
190 measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
194 //Normalize with the global measures
195 measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
206 * Get the Derivative Measure
208 template < class TFixedImage, class TMovingImage>
210 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
211 ::GetDerivative( const TransformParametersType & parameters,
212 DerivativeType & derivative ) const
215 itkDebugMacro("GetDerivative( " << parameters << " ) ");
217 if( !this->GetGradientImage() ) {
218 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
221 FixedImageConstPointer fixedImage = this->m_FixedImage;
224 itkExceptionMacro( << "Fixed image has not been assigned" );
227 const unsigned int ImageDimension = FixedImageType::ImageDimension;
230 typedef itk::ImageRegionConstIteratorWithIndex<
231 FixedImageType> FixedIteratorType;
233 typedef itk::ImageRegionConstIteratorWithIndex<
234 typename Superclass::GradientImageType> GradientIteratorType;
237 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
239 typename FixedImageType::IndexType index;
241 this->m_NumberOfPixelsCounted = 0;
243 this->SetTransformParameters( parameters );
245 const unsigned int ParametersDimension = this->GetNumberOfParameters();
246 derivative = DerivativeType( ParametersDimension );
247 derivative.Fill( itk::NumericTraits<typename DerivativeType::ValueType>::Zero );
251 while(!ti.IsAtEnd()) {
253 index = ti.GetIndex();
255 typename Superclass::InputPointType inputPoint;
256 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
258 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
263 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
265 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
270 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
271 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
273 TransformJacobianType jacobian;
274 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint , jacobian);
276 const RealType fixedValue = ti.Value();
277 this->m_NumberOfPixelsCounted++;
278 const RealType diff = movingValue - fixedValue;
280 // Get the gradient by NearestNeighboorInterpolation:
281 // which is equivalent to round up the point components.
282 typedef typename Superclass::OutputPointType OutputPointType;
283 typedef typename OutputPointType::CoordRepType CoordRepType;
284 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
285 MovingImageContinuousIndexType;
287 MovingImageContinuousIndexType tempIndex;
288 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
290 typename MovingImageType::IndexType mappedIndex;
291 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
292 mappedIndex[j] = static_cast<long>( std::round( tempIndex[j] ) );
295 const GradientPixelType gradient =
296 this->GetGradientImage()->GetPixel( mappedIndex );
298 for(unsigned int par=0; par<ParametersDimension; par++) {
299 RealType sum = NumericTraits< RealType >::Zero;
300 for(unsigned int dim=0; dim<ImageDimension; dim++) {
301 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
303 derivative[par] += sum;
310 if( !this->m_NumberOfPixelsCounted ) {
311 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
313 for(unsigned int i=0; i<ParametersDimension; i++) {
314 derivative[i] /= this->m_NumberOfPixelsCounted;
322 * Get both the match Measure and the Derivative Measure
324 template <class TFixedImage, class TMovingImage>
326 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
327 ::GetValueAndDerivative(const TransformParametersType & parameters,
328 MeasureType & value, DerivativeType & derivative) const
331 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
333 if( !this->GetGradientImage() ) {
334 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
337 FixedImageConstPointer fixedImage = this->m_FixedImage;
340 itkExceptionMacro( << "Fixed image has not been assigned" );
343 const unsigned int ImageDimension = FixedImageType::ImageDimension;
345 typedef itk::ImageRegionConstIteratorWithIndex<
346 FixedImageType> FixedIteratorType;
348 typedef itk::ImageRegionConstIteratorWithIndex<
349 typename Superclass::GradientImageType> GradientIteratorType;
352 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
354 typename FixedImageType::IndexType index;
356 MeasureType measure = NumericTraits< MeasureType >::Zero;
358 this->m_NumberOfPixelsCounted = 0;
360 this->SetTransformParameters( parameters );
362 const unsigned int ParametersDimension = this->GetNumberOfParameters();
363 derivative = DerivativeType( ParametersDimension );
364 derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
368 while(!ti.IsAtEnd()) {
370 index = ti.GetIndex();
372 typename Superclass::InputPointType inputPoint;
373 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
375 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
380 typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
382 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
387 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
388 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
390 TransformJacobianType jacobian;
391 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
393 const RealType fixedValue = ti.Value();
394 this->m_NumberOfPixelsCounted++;
396 const RealType diff = movingValue - fixedValue;
398 measure += diff * diff;
400 // Get the gradient by NearestNeighboorInterpolation:
401 // which is equivalent to round up the point components.
402 typedef typename Superclass::OutputPointType OutputPointType;
403 typedef typename OutputPointType::CoordRepType CoordRepType;
404 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
405 MovingImageContinuousIndexType;
407 MovingImageContinuousIndexType tempIndex;
408 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
410 typename MovingImageType::IndexType mappedIndex;
411 for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
412 mappedIndex[j] = static_cast<long>( std::round( tempIndex[j] ) );
415 const GradientPixelType gradient =
416 this->GetGradientImage()->GetPixel( mappedIndex );
418 for(unsigned int par=0; par<ParametersDimension; par++) {
419 RealType sum = NumericTraits< RealType >::Zero;
420 for(unsigned int dim=0; dim<ImageDimension; dim++) {
421 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
423 derivative[par] += sum;
430 if( !this->m_NumberOfPixelsCounted ) {
431 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
433 for(unsigned int i=0; i<ParametersDimension; i++) {
434 derivative[i] /= this->m_NumberOfPixelsCounted;
436 measure /= this->m_NumberOfPixelsCounted;
443 } // end namespace itk