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 /*=========================================================================
21 Program: Insight Segmentation & Registration Toolkit
22 Module: $RCSfile: itkMeanSquaresImageToImageMetricFor3DBLUTFFD.txx,v $
24 Date: $Date: 2010/06/14 17:32:07 $
25 Version: $Revision: 1.1 $
27 Copyright (c) Insight Software Consortium. All rights reserved.
28 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
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.
34 =========================================================================*/
35 #ifndef __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
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"
43 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx"
45 #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
46 #include "itkImageRegionConstIteratorWithIndex.h"
54 template <class TFixedImage, class TMovingImage>
55 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
56 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
58 itkDebugMacro("Constructor");
62 * Get the match Measure
64 template <class TFixedImage, class TMovingImage>
65 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
66 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
67 ::GetValue( const TransformParametersType & parameters ) const
70 itkDebugMacro("GetValue( " << parameters << " ) ");
72 FixedImageConstPointer fixedImage = this->m_FixedImage;
75 itkExceptionMacro( << "Fixed image has not been assigned" );
78 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
81 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
83 typename FixedImageType::IndexType index;
85 MeasureType measure = NumericTraits< MeasureType >::Zero;
87 this->m_NumberOfPixelsCounted = 0;
89 this->SetTransformParameters( parameters );
91 while(!ti.IsAtEnd()) {
93 index = ti.GetIndex();
95 InputPointType inputPoint;
96 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
98 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
103 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
105 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
110 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
111 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
112 const RealType fixedValue = ti.Get();
113 this->m_NumberOfPixelsCounted++;
114 const RealType diff = movingValue - fixedValue;
115 measure += diff * diff;
121 if( !this->m_NumberOfPixelsCounted ) {
122 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
124 measure /= this->m_NumberOfPixelsCounted;
132 * Get the Derivative Measure
134 template < class TFixedImage, class TMovingImage>
136 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
137 ::GetDerivative( const TransformParametersType & parameters,
138 DerivativeType & derivative ) const
141 itkDebugMacro("GetDerivative( " << parameters << " ) ");
143 if( !this->GetGradientImage() ) {
144 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
147 FixedImageConstPointer fixedImage = this->m_FixedImage;
150 itkExceptionMacro( << "Fixed image has not been assigned" );
153 const unsigned int ImageDimension = FixedImageType::ImageDimension;
156 typedef itk::ImageRegionConstIteratorWithIndex<
157 FixedImageType> FixedIteratorType;
159 typedef itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
162 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
164 typename FixedImageType::IndexType index;
166 this->m_NumberOfPixelsCounted = 0;
168 this->SetTransformParameters( parameters );
170 const unsigned int ParametersDimension = this->GetNumberOfParameters();
171 derivative = DerivativeType( ParametersDimension );
172 derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
176 while(!ti.IsAtEnd()) {
178 index = ti.GetIndex();
180 InputPointType inputPoint;
181 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
183 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
188 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
190 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
195 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
196 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
198 TransformJacobianType jacobian;
199 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
201 const RealType fixedValue = ti.Value();
202 this->m_NumberOfPixelsCounted++;
203 const RealType diff = movingValue - fixedValue;
205 // Get the gradient by NearestNeighboorInterpolation:
206 // which is equivalent to round up the point components.
207 typedef typename OutputPointType::CoordRepType CoordRepType;
208 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
209 MovingImageContinuousIndexType;
211 MovingImageContinuousIndexType tempIndex;
212 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
214 typename MovingImageType::IndexType mappedIndex;
215 mappedIndex.CopyWithRound( tempIndex );
217 const GradientPixelType gradient =
218 this->GetGradientImage()->GetPixel( mappedIndex );
220 for(unsigned int par=0; par<ParametersDimension; par++) {
221 RealType sum = NumericTraits< RealType >::Zero;
222 for(unsigned int dim=0; dim<ImageDimension; dim++) {
223 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
225 derivative[par] += sum;
232 if( !this->m_NumberOfPixelsCounted ) {
233 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
235 for(unsigned int i=0; i<ParametersDimension; i++) {
236 derivative[i] /= this->m_NumberOfPixelsCounted;
244 * Get both the match Measure and theDerivative Measure
246 template <class TFixedImage, class TMovingImage>
248 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
249 ::GetValueAndDerivative(const TransformParametersType & parameters,
250 MeasureType & value, DerivativeType & derivative) const
253 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
255 if( !this->GetGradientImage() ) {
256 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
259 FixedImageConstPointer fixedImage = this->m_FixedImage;
262 itkExceptionMacro( << "Fixed image has not been assigned" );
265 const unsigned int ImageDimension = FixedImageType::ImageDimension;
267 typedef itk::ImageRegionConstIteratorWithIndex<
268 FixedImageType> FixedIteratorType;
270 typedef itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
273 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
275 typename FixedImageType::IndexType index;
277 MeasureType measure = NumericTraits< MeasureType >::Zero;
279 this->m_NumberOfPixelsCounted = 0;
281 this->SetTransformParameters( parameters );
283 const unsigned int ParametersDimension = this->GetNumberOfParameters();
284 derivative = DerivativeType( ParametersDimension );
285 derivative.Fill( NumericTraits<typename DerivativeType::ValueType>::Zero );
289 while(!ti.IsAtEnd()) {
291 index = ti.GetIndex();
293 InputPointType inputPoint;
294 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
296 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
301 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
303 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
308 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
309 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
311 TransformJacobianType jacobian;
312 this->m_Transform->ComputeJacobianWithRespectToParameters( inputPoint, jacobian );
314 const RealType fixedValue = ti.Value();
315 this->m_NumberOfPixelsCounted++;
317 const RealType diff = movingValue - fixedValue;
319 measure += diff * diff;
321 // Get the gradient by NearestNeighboorInterpolation:
322 // which is equivalent to round up the point components.
323 typedef typename OutputPointType::CoordRepType CoordRepType;
324 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
325 MovingImageContinuousIndexType;
327 MovingImageContinuousIndexType tempIndex;
328 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
330 typename MovingImageType::IndexType mappedIndex;
331 mappedIndex.CopyWithRound( tempIndex );
333 const GradientPixelType gradient =
334 this->GetGradientImage()->GetPixel( mappedIndex );
336 for(unsigned int par=0; par<ParametersDimension; par++) {
337 RealType sum = NumericTraits< RealType >::Zero;
338 for(unsigned int dim=0; dim<ImageDimension; dim++) {
339 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
341 derivative[par] += sum;
348 if( !this->m_NumberOfPixelsCounted ) {
349 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
351 for(unsigned int i=0; i<ParametersDimension; i++) {
352 derivative[i] /= this->m_NumberOfPixelsCounted;
354 measure /= this->m_NumberOfPixelsCounted;
361 } // end namespace itk