]> Creatis software - clitk.git/blob - registration/itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx
Added more options to clitkImageUncertainty, default should be identical
[clitk.git] / registration / itkMeanSquaresImageToImageMetricFor3DBLUTFFD.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 /*=========================================================================
20
21   Program:   Insight Segmentation & Registration Toolkit
22   Module:    $RCSfile: itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx,v $
23   Language:  C++
24   Date:      $Date: 2010/06/14 17:32:07 $
25   Version:   $Revision: 1.1 $
26
27   Copyright (c) Insight Software Consortium. All rights reserved.
28   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
29
30      This software is distributed WITHOUT ANY WARRANTY; without even
31      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
32      PURPOSE.  See the above copyright notices for more information.
33
34 =========================================================================*/
35 #ifndef __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
37
38 // First make sure that the configuration is available.
39 // This line can be removed once the optimized versions
40 // gets integrated into the main directories.
41 #include "itkConfigure.h"
42
43 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
44 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx"
45 #else
46
47 #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
48 #include "itkImageRegionConstIteratorWithIndex.h"
49
50 namespace itk
51 {
52
53 /**
54  * Constructor
55  */
56 template <class TFixedImage, class TMovingImage>
57 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
58 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
59 {
60   itkDebugMacro("Constructor");
61 }
62
63 /**
64  * Get the match Measure
65  */
66 template <class TFixedImage, class TMovingImage>
67 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
68 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
69 ::GetValue( const TransformParametersType & parameters ) const
70 {
71
72   itkDebugMacro("GetValue( " << parameters << " ) ");
73
74   FixedImageConstPointer fixedImage = this->m_FixedImage;
75
76   if( !fixedImage ) {
77     itkExceptionMacro( << "Fixed image has not been assigned" );
78   }
79
80   typedef  itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
81
82
83   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
84
85   typename FixedImageType::IndexType index;
86
87   MeasureType measure = NumericTraits< MeasureType >::Zero;
88
89   this->m_NumberOfPixelsCounted = 0;
90
91   this->SetTransformParameters( parameters );
92
93   while(!ti.IsAtEnd()) {
94
95     index = ti.GetIndex();
96
97     InputPointType inputPoint;
98     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
99
100     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
101       ++ti;
102       continue;
103     }
104
105     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
106
107     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
108       ++ti;
109       continue;
110     }
111
112     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
113       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
114       const RealType fixedValue   = ti.Get();
115       this->m_NumberOfPixelsCounted++;
116       const RealType diff = movingValue - fixedValue;
117       measure += diff * diff;
118     }
119
120     ++ti;
121   }
122
123   if( !this->m_NumberOfPixelsCounted ) {
124     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
125   } else {
126     measure /= this->m_NumberOfPixelsCounted;
127   }
128
129   return measure;
130
131 }
132
133 /**
134  * Get the Derivative Measure
135  */
136 template < class TFixedImage, class TMovingImage>
137 void
138 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
139 ::GetDerivative( const TransformParametersType & parameters,
140                  DerivativeType & derivative  ) const
141 {
142
143   itkDebugMacro("GetDerivative( " << parameters << " ) ");
144
145   if( !this->GetGradientImage() ) {
146     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
147   }
148
149   FixedImageConstPointer fixedImage = this->m_FixedImage;
150
151   if( !fixedImage ) {
152     itkExceptionMacro( << "Fixed image has not been assigned" );
153   }
154
155   const unsigned int ImageDimension = FixedImageType::ImageDimension;
156
157
158   typedef  itk::ImageRegionConstIteratorWithIndex<
159   FixedImageType> FixedIteratorType;
160
161   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
162
163
164   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
165
166   typename FixedImageType::IndexType index;
167
168   this->m_NumberOfPixelsCounted = 0;
169
170   this->SetTransformParameters( parameters );
171
172   const unsigned int ParametersDimension = this->GetNumberOfParameters();
173   derivative = DerivativeType( ParametersDimension );
174   derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
175
176   ti.GoToBegin();
177
178   while(!ti.IsAtEnd()) {
179
180     index = ti.GetIndex();
181
182     InputPointType inputPoint;
183     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
184
185     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
186       ++ti;
187       continue;
188     }
189
190     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
191
192     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
193       ++ti;
194       continue;
195     }
196
197     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
198       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
199
200 #if ITK_VERSION_MAJOR >= 4
201       TransformJacobianType jacobian;
202       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
203 #else
204       const TransformJacobianType & jacobian =
205         this->m_Transform->GetJacobian( inputPoint );
206 #endif
207
208       const RealType fixedValue     = ti.Value();
209       this->m_NumberOfPixelsCounted++;
210       const RealType diff = movingValue - fixedValue;
211
212       // Get the gradient by NearestNeighboorInterpolation:
213       // which is equivalent to round up the point components.
214       typedef typename OutputPointType::CoordRepType CoordRepType;
215       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
216       MovingImageContinuousIndexType;
217
218       MovingImageContinuousIndexType tempIndex;
219       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
220
221       typename MovingImageType::IndexType mappedIndex;
222       mappedIndex.CopyWithRound( tempIndex );
223
224       const GradientPixelType gradient =
225         this->GetGradientImage()->GetPixel( mappedIndex );
226
227       for(unsigned int par=0; par<ParametersDimension; par++) {
228         RealType sum = NumericTraits< RealType >::Zero;
229         for(unsigned int dim=0; dim<ImageDimension; dim++) {
230           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
231         }
232         derivative[par] += sum;
233       }
234     }
235
236     ++ti;
237   }
238
239   if( !this->m_NumberOfPixelsCounted ) {
240     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
241   } else {
242     for(unsigned int i=0; i<ParametersDimension; i++) {
243       derivative[i] /= this->m_NumberOfPixelsCounted;
244     }
245   }
246
247 }
248
249
250 /*
251  * Get both the match Measure and theDerivative Measure
252  */
253 template <class TFixedImage, class TMovingImage>
254 void
255 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
256 ::GetValueAndDerivative(const TransformParametersType & parameters,
257                         MeasureType & value, DerivativeType  & derivative) const
258 {
259
260   itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
261
262   if( !this->GetGradientImage() ) {
263     itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
264   }
265
266   FixedImageConstPointer fixedImage = this->m_FixedImage;
267
268   if( !fixedImage ) {
269     itkExceptionMacro( << "Fixed image has not been assigned" );
270   }
271
272   const unsigned int ImageDimension = FixedImageType::ImageDimension;
273
274   typedef  itk::ImageRegionConstIteratorWithIndex<
275   FixedImageType> FixedIteratorType;
276
277   typedef  itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
278
279
280   FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
281
282   typename FixedImageType::IndexType index;
283
284   MeasureType measure = NumericTraits< MeasureType >::Zero;
285
286   this->m_NumberOfPixelsCounted = 0;
287
288   this->SetTransformParameters( parameters );
289
290   const unsigned int ParametersDimension = this->GetNumberOfParameters();
291   derivative = DerivativeType( ParametersDimension );
292   derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
293
294   ti.GoToBegin();
295
296   while(!ti.IsAtEnd()) {
297
298     index = ti.GetIndex();
299
300     InputPointType inputPoint;
301     fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
302
303     if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
304       ++ti;
305       continue;
306     }
307
308     OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
309
310     if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
311       ++ti;
312       continue;
313     }
314
315     if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
316       const RealType movingValue  = this->m_Interpolator->Evaluate( transformedPoint );
317
318 #if ITK_VERSION_MAJOR >= 4
319       TransformJacobianType jacobian;
320       this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
321 #else
322       const TransformJacobianType & jacobian =
323         this->m_Transform->GetJacobian( inputPoint );
324 #endif
325
326       const RealType fixedValue     = ti.Value();
327       this->m_NumberOfPixelsCounted++;
328
329       const RealType diff = movingValue - fixedValue;
330
331       measure += diff * diff;
332
333       // Get the gradient by NearestNeighboorInterpolation:
334       // which is equivalent to round up the point components.
335       typedef typename OutputPointType::CoordRepType CoordRepType;
336       typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
337       MovingImageContinuousIndexType;
338
339       MovingImageContinuousIndexType tempIndex;
340       this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
341
342       typename MovingImageType::IndexType mappedIndex;
343       mappedIndex.CopyWithRound( tempIndex );
344
345       const GradientPixelType gradient =
346         this->GetGradientImage()->GetPixel( mappedIndex );
347
348       for(unsigned int par=0; par<ParametersDimension; par++) {
349         RealType sum = NumericTraits< RealType >::Zero;
350         for(unsigned int dim=0; dim<ImageDimension; dim++) {
351           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
352         }
353         derivative[par] += sum;
354       }
355     }
356
357     ++ti;
358   }
359
360   if( !this->m_NumberOfPixelsCounted ) {
361     itkExceptionMacro(<<"All the points mapped to outside of the moving image");
362   } else {
363     for(unsigned int i=0; i<ParametersDimension; i++) {
364       derivative[i] /= this->m_NumberOfPixelsCounted;
365     }
366     measure /= this->m_NumberOfPixelsCounted;
367   }
368
369   value = measure;
370
371 }
372
373 } // end namespace itk
374
375
376 #endif
377
378 #endif