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