]> Creatis software - clitk.git/blob - tools/clitkCorrelationRatioImageToImageMetric.txx
Replaced round function with portable rint function
[clitk.git] / tools / clitkCorrelationRatioImageToImageMetric.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
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
8
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.
12
13   It is distributed under dual licence
14
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
20
21 /**
22  * @file   clitkCorrelationRatioImageToImageMetric.txx
23  * @author Jef Vandemeulebroucke <jef@creatis.insa-lyon.fr>
24  * @date   July 30  18:14:53 2007
25  *
26  * @brief  Compute the correlation ratio between 2 images
27  *
28  *
29  */
30
31 #include "clitkCorrelationRatioImageToImageMetric.h"
32 #include "itkImageRegionConstIteratorWithIndex.h"
33 #include "itkImageRegionConstIterator.h"
34 #include "itkImageRegionIterator.h"
35
36 namespace clitk
37 {
38
39 /*
40  * Constructor
41  */
42 template <class TFixedImage, class TMovingImage>
43 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
44 ::CorrelationRatioImageToImageMetric()
45 {
46   m_NumberOfBins = 50;
47
48 }
49
50 template <class TFixedImage, class TMovingImage>
51 void
52 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
53 ::Initialize(void) throw ( ExceptionObject )
54 {
55
56   this->Superclass::Initialize();
57
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();
62
63   typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
64   FixedIteratorType fixedImageIterator(
65     this->m_FixedImage, this->GetFixedImageRegion() );
66
67   for ( fixedImageIterator.GoToBegin();
68         !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
69
70     double sample = static_cast<double>( fixedImageIterator.Get() );
71
72     if ( sample < fixedImageMin ) {
73       fixedImageMin = sample;
74     }
75
76     if ( sample > fixedImageMax ) {
77       fixedImageMax = sample;
78     }
79   }
80
81   // Compute binsize for the fixedImage
82   m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / m_NumberOfBins;
83   m_FixedImageMin=fixedImageMin;
84   //Allocate mempry and initialise the fixed image bin
85   m_NumberOfPixelsCountedPerBin.resize( m_NumberOfBins, 0 );
86   m_mMSVPB.resize( m_NumberOfBins, 0.0 );
87   m_mSMVPB.resize( m_NumberOfBins, 0.0 );
88 }
89
90
91 /*
92  * Get the match Measure
93  */
94 template <class TFixedImage, class TMovingImage>
95 typename CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>::MeasureType
96 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
97 ::GetValue( const TransformParametersType & parameters ) const
98 {
99
100   itkDebugMacro("GetValue( " << parameters << " ) ");
101
102   FixedImageConstPointer fixedImage = this->m_FixedImage;
103
104   if( !fixedImage ) {
105     itkExceptionMacro( << "Fixed image has not been assigned" );
106   }
107
108   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
109
110
111   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
112
113   typename FixedImageType::IndexType index;
114
115   MeasureType measure = itk::NumericTraits< MeasureType >::Zero;
116
117   this->m_NumberOfPixelsCounted = 0;
118   this->SetTransformParameters( parameters );
119
120
121   //temporary measures for the calculation
122   RealType mSMV=0;
123   RealType mMSV=0;
124
125   while(!ti.IsAtEnd()) {
126
127     index = ti.GetIndex();
128
129     typename Superclass::InputPointType inputPoint;
130     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
131
132     // Verify that the point is in the fixed Image Mask
133     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
134       ++ti;
135       continue;
136     }
137
138     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
139
140     //Verify that the point is in the moving Image Mask
141     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
142       ++ti;
143       continue;
144     }
145
146     // Verify is the interpolated value is in the buffer
147     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
148       //Accumulate calculations for the correlation ratio
149       //For each pixel the is in both masks and the buffer we adapt the following measures:
150       //movingMeanSquaredValue mMSV; movingSquaredMeanValue mSMV;
151       //movingMeanSquaredValuePerBin[i] mSMVPB; movingSquaredMeanValuePerBin[i] mSMVPB
152       //NumberOfPixelsCounted, NumberOfPixelsCountedPerBin[i]
153
154       //get the value of the moving image
155       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
156       // for the variance of the overlapping moving image we accumulate the following measures
157       const RealType movingSquaredValue=movingValue*movingValue;
158       mMSV+=movingSquaredValue;
159       mSMV+=movingValue;
160
161       //get the fixed value
162       const RealType fixedValue   = ti.Get();
163
164       //check in which bin the fixed value belongs, get the index
165       const double fixedImageBinTerm =        (fixedValue - m_FixedImageMin) / m_FixedImageBinSize;
166       const unsigned int fixedImageBinIndex = static_cast<unsigned int>( vcl_floor(fixedImageBinTerm ) );
167       //adapt the measures per bin
168       this->m_mMSVPB[fixedImageBinIndex]+=movingSquaredValue;
169       this->m_mSMVPB[fixedImageBinIndex]+=movingValue;
170       //increase the fixed image bin and the total pixel count
171       this->m_NumberOfPixelsCountedPerBin[fixedImageBinIndex]+=1;
172       this->m_NumberOfPixelsCounted++;
173     }
174
175     ++ti;
176   }
177
178   if( !this->m_NumberOfPixelsCounted ) {
179     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
180   } else {
181
182     //apdapt the measures per bin
183     for (unsigned int i=0; i< m_NumberOfBins; i++ ) {
184       if (this->m_NumberOfPixelsCountedPerBin[i]>0) {
185         measure+=(this->m_mMSVPB[i]-((this->m_mSMVPB[i]*this->m_mSMVPB[i])/this->m_NumberOfPixelsCountedPerBin[i]));
186       }
187     }
188
189     //Normalize with the global measures
190     measure /= (mMSV-((mSMV*mSMV)/ this->m_NumberOfPixelsCounted));
191     return measure;
192
193   }
194 }
195
196
197
198
199
200 /*
201  * Get the Derivative Measure
202  */
203 template < class TFixedImage, class TMovingImage>
204 void
205 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
206 ::GetDerivative( const TransformParametersType & parameters,
207                  DerivativeType & derivative  ) const
208 {
209
210   itkDebugMacro("GetDerivative( " << parameters << " ) ");
211
212   if( !this->GetGradientImage() ) {
213     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
214   }
215
216   FixedImageConstPointer fixedImage = this->m_FixedImage;
217
218   if( !fixedImage ) {
219     itkExceptionMacro( << "Fixed image has not been assigned" );
220   }
221
222   const unsigned int ImageDimension = FixedImageType::ImageDimension;
223
224
225   typedef  itk::ImageRegionConstIteratorWithIndex<
226   FixedImageType> FixedIteratorType;
227
228   typedef  itk::ImageRegionConstIteratorWithIndex<
229   ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
230
231
232   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
233
234   typename FixedImageType::IndexType index;
235
236   this->m_NumberOfPixelsCounted = 0;
237
238   this->SetTransformParameters( parameters );
239
240   const unsigned int ParametersDimension = this->GetNumberOfParameters();
241   derivative = DerivativeType( ParametersDimension );
242   derivative.Fill( itk::NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
243
244   ti.GoToBegin();
245
246   while(!ti.IsAtEnd()) {
247
248     index = ti.GetIndex();
249
250     typename Superclass::InputPointType inputPoint;
251     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
252
253     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
254       ++ti;
255       continue;
256     }
257
258     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
259
260     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
261       ++ti;
262       continue;
263     }
264
265     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
266       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
267
268       const TransformJacobianType & jacobian =
269         this->m_Transform->GetJacobian( inputPoint );
270
271
272       const RealType fixedValue     = ti.Value();
273       this->m_NumberOfPixelsCounted++;
274       const RealType diff = movingValue - fixedValue;
275
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;
282
283       MovingImageContinuousIndexType tempIndex;
284       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
285
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] ) );
289       }
290
291       const GradientPixelType gradient =
292         this->GetGradientImage()->GetPixel( mappedIndex );
293
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];
298         }
299         derivative[par] += sum;
300       }
301     }
302
303     ++ti;
304   }
305
306   if( !this->m_NumberOfPixelsCounted ) {
307     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
308   } else {
309     for(unsigned int i=0; i<ParametersDimension; i++) {
310       derivative[i] /= this->m_NumberOfPixelsCounted;
311     }
312   }
313
314 }
315
316
317 /*
318  * Get both the match Measure and the Derivative Measure
319  */
320 template <class TFixedImage, class TMovingImage>
321 void
322 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
323 ::GetValueAndDerivative(const TransformParametersType & parameters,
324                         MeasureType & value, DerivativeType  & derivative) const
325 {
326
327   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
328
329   if( !this->GetGradientImage() ) {
330     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
331   }
332
333   FixedImageConstPointer fixedImage = this->m_FixedImage;
334
335   if( !fixedImage ) {
336     itkExceptionMacro( << "Fixed image has not been assigned" );
337   }
338
339   const unsigned int ImageDimension = FixedImageType::ImageDimension;
340
341   typedef  itk::ImageRegionConstIteratorWithIndex<
342   FixedImageType> FixedIteratorType;
343
344   typedef  itk::ImageRegionConstIteratorWithIndex<
345   ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
346
347
348   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
349
350   typename FixedImageType::IndexType index;
351
352   MeasureType measure = NumericTraits< MeasureType >::Zero;
353
354   this->m_NumberOfPixelsCounted = 0;
355
356   this->SetTransformParameters( parameters );
357
358   const unsigned int ParametersDimension = this->GetNumberOfParameters();
359   derivative = DerivativeType( ParametersDimension );
360   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
361
362   ti.GoToBegin();
363
364   while(!ti.IsAtEnd()) {
365
366     index = ti.GetIndex();
367
368     typename Superclass::InputPointType inputPoint;
369     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
370
371     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
372       ++ti;
373       continue;
374     }
375
376     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
377
378     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
379       ++ti;
380       continue;
381     }
382
383     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
384       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
385
386       const TransformJacobianType & jacobian =
387         this->m_Transform->GetJacobian( inputPoint );
388
389
390       const RealType fixedValue     = ti.Value();
391       this->m_NumberOfPixelsCounted++;
392
393       const RealType diff = movingValue - fixedValue;
394
395       measure += diff * diff;
396
397       // Get the gradient by NearestNeighboorInterpolation:
398       // which is equivalent to round up the point components.
399       typedef typename Superclass::OutputPointType OutputPointType;
400       typedef typename OutputPointType::CoordRepType CoordRepType;
401       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
402       MovingImageContinuousIndexType;
403
404       MovingImageContinuousIndexType tempIndex;
405       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
406
407       typename MovingImageType::IndexType mappedIndex;
408       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
409         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
410       }
411
412       const GradientPixelType gradient =
413         this->GetGradientImage()->GetPixel( mappedIndex );
414
415       for(unsigned int par=0; par<ParametersDimension; par++) {
416         RealType sum = NumericTraits< RealType >::Zero;
417         for(unsigned int dim=0; dim<ImageDimension; dim++) {
418           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
419         }
420         derivative[par] += sum;
421       }
422     }
423
424     ++ti;
425   }
426
427   if( !this->m_NumberOfPixelsCounted ) {
428     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
429   } else {
430     for(unsigned int i=0; i<ParametersDimension; i++) {
431       derivative[i] /= this->m_NumberOfPixelsCounted;
432     }
433     measure /= this->m_NumberOfPixelsCounted;
434   }
435
436   value = measure;
437
438 }
439
440 } // end namespace itk
441
442
443 #endif
444