]> Creatis software - clitk.git/blob - registration/clitkCorrelationRatioImageToImageMetric.txx
min max (to remove)
[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://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
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   ITK_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<ITK_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       const TransformJacobianType & jacobian =
270         this->m_Transform->GetJacobian( inputPoint );
271
272
273       const RealType fixedValue     = ti.Value();
274       this->m_NumberOfPixelsCounted++;
275       const RealType diff = movingValue - fixedValue;
276
277       // Get the gradient by NearestNeighboorInterpolation:
278       // which is equivalent to round up the point components.
279       typedef typename Superclass::OutputPointType OutputPointType;
280       typedef typename OutputPointType::CoordRepType CoordRepType;
281       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
282       MovingImageContinuousIndexType;
283
284       MovingImageContinuousIndexType tempIndex;
285       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
286
287       typename MovingImageType::IndexType mappedIndex;
288       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
289         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
290       }
291
292       const GradientPixelType gradient =
293         this->GetGradientImage()->GetPixel( mappedIndex );
294
295       for(unsigned int par=0; par<ParametersDimension; par++) {
296         RealType sum = NumericTraits< RealType >::Zero;
297         for(unsigned int dim=0; dim<ImageDimension; dim++) {
298           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
299         }
300         derivative[par] += sum;
301       }
302     }
303
304     ++ti;
305   }
306
307   if( !this->m_NumberOfPixelsCounted ) {
308     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
309   } else {
310     for(unsigned int i=0; i<ParametersDimension; i++) {
311       derivative[i] /= this->m_NumberOfPixelsCounted;
312     }
313   }
314
315 }
316
317
318 /*
319  * Get both the match Measure and the Derivative Measure
320  */
321 template <class TFixedImage, class TMovingImage>
322 void
323 CorrelationRatioImageToImageMetric<TFixedImage,TMovingImage>
324 ::GetValueAndDerivative(const TransformParametersType & parameters,
325                         MeasureType & value, DerivativeType  & derivative) const
326 {
327
328   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
329
330   if( !this->GetGradientImage() ) {
331     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
332   }
333
334   FixedImageConstPointer fixedImage = this->m_FixedImage;
335
336   if( !fixedImage ) {
337     itkExceptionMacro( << "Fixed image has not been assigned" );
338   }
339
340   const unsigned int ImageDimension = FixedImageType::ImageDimension;
341
342   typedef  itk::ImageRegionConstIteratorWithIndex<
343   FixedImageType> FixedIteratorType;
344
345   typedef  itk::ImageRegionConstIteratorWithIndex<
346   ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
347
348
349   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
350
351   typename FixedImageType::IndexType index;
352
353   MeasureType measure = NumericTraits< MeasureType >::Zero;
354
355   this->m_NumberOfPixelsCounted = 0;
356
357   this->SetTransformParameters( parameters );
358
359   const unsigned int ParametersDimension = this->GetNumberOfParameters();
360   derivative = DerivativeType( ParametersDimension );
361   derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
362
363   ti.GoToBegin();
364
365   while(!ti.IsAtEnd()) {
366
367     index = ti.GetIndex();
368
369     typename Superclass::InputPointType inputPoint;
370     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
371
372     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
373       ++ti;
374       continue;
375     }
376
377     typename Superclass::OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
378
379     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
380       ++ti;
381       continue;
382     }
383
384     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
385       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
386
387       const TransformJacobianType & jacobian =
388         this->m_Transform->GetJacobian( inputPoint );
389
390
391       const RealType fixedValue     = ti.Value();
392       this->m_NumberOfPixelsCounted++;
393
394       const RealType diff = movingValue - fixedValue;
395
396       measure += diff * diff;
397
398       // Get the gradient by NearestNeighboorInterpolation:
399       // which is equivalent to round up the point components.
400       typedef typename Superclass::OutputPointType OutputPointType;
401       typedef typename OutputPointType::CoordRepType CoordRepType;
402       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
403       MovingImageContinuousIndexType;
404
405       MovingImageContinuousIndexType tempIndex;
406       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
407
408       typename MovingImageType::IndexType mappedIndex;
409       for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ ) {
410         mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
411       }
412
413       const GradientPixelType gradient =
414         this->GetGradientImage()->GetPixel( mappedIndex );
415
416       for(unsigned int par=0; par<ParametersDimension; par++) {
417         RealType sum = NumericTraits< RealType >::Zero;
418         for(unsigned int dim=0; dim<ImageDimension; dim++) {
419           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
420         }
421         derivative[par] += sum;
422       }
423     }
424
425     ++ti;
426   }
427
428   if( !this->m_NumberOfPixelsCounted ) {
429     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
430   } else {
431     for(unsigned int i=0; i<ParametersDimension; i++) {
432       derivative[i] /= this->m_NumberOfPixelsCounted;
433     }
434     measure /= this->m_NumberOfPixelsCounted;
435   }
436
437   value = measure;
438
439 }
440
441 } // end namespace itk
442
443
444 #endif
445