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 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
44 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx"
47 #include "itkMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
48 #include "itkImageRegionConstIteratorWithIndex.h"
56 template <class TFixedImage, class TMovingImage>
57 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
58 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
60 itkDebugMacro("Constructor");
64 * Get the match Measure
66 template <class TFixedImage, class TMovingImage>
67 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>::MeasureType
68 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
69 ::GetValue( const TransformParametersType & parameters ) const
72 itkDebugMacro("GetValue( " << parameters << " ) ");
74 FixedImageConstPointer fixedImage = this->m_FixedImage;
77 itkExceptionMacro( << "Fixed image has not been assigned" );
80 typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
83 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
85 typename FixedImageType::IndexType index;
87 MeasureType measure = NumericTraits< MeasureType >::Zero;
89 this->m_NumberOfPixelsCounted = 0;
91 this->SetTransformParameters( parameters );
93 while(!ti.IsAtEnd()) {
95 index = ti.GetIndex();
97 InputPointType inputPoint;
98 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
100 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
105 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
107 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
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;
123 if( !this->m_NumberOfPixelsCounted ) {
124 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
126 measure /= this->m_NumberOfPixelsCounted;
134 * Get the Derivative Measure
136 template < class TFixedImage, class TMovingImage>
138 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
139 ::GetDerivative( const TransformParametersType & parameters,
140 DerivativeType & derivative ) const
143 itkDebugMacro("GetDerivative( " << parameters << " ) ");
145 if( !this->GetGradientImage() ) {
146 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
149 FixedImageConstPointer fixedImage = this->m_FixedImage;
152 itkExceptionMacro( << "Fixed image has not been assigned" );
155 const unsigned int ImageDimension = FixedImageType::ImageDimension;
158 typedef itk::ImageRegionConstIteratorWithIndex<
159 FixedImageType> FixedIteratorType;
161 typedef itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
164 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
166 typename FixedImageType::IndexType index;
168 this->m_NumberOfPixelsCounted = 0;
170 this->SetTransformParameters( parameters );
172 const unsigned int ParametersDimension = this->GetNumberOfParameters();
173 derivative = DerivativeType( ParametersDimension );
174 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
178 while(!ti.IsAtEnd()) {
180 index = ti.GetIndex();
182 InputPointType inputPoint;
183 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
185 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
190 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
192 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
197 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
198 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
200 const TransformJacobianType & jacobian =
201 this->m_Transform->GetJacobian( inputPoint );
204 const RealType fixedValue = ti.Value();
205 this->m_NumberOfPixelsCounted++;
206 const RealType diff = movingValue - fixedValue;
208 // Get the gradient by NearestNeighboorInterpolation:
209 // which is equivalent to round up the point components.
210 typedef typename OutputPointType::CoordRepType CoordRepType;
211 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
212 MovingImageContinuousIndexType;
214 MovingImageContinuousIndexType tempIndex;
215 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
217 typename MovingImageType::IndexType mappedIndex;
218 mappedIndex.CopyWithRound( tempIndex );
220 const GradientPixelType gradient =
221 this->GetGradientImage()->GetPixel( mappedIndex );
223 for(unsigned int par=0; par<ParametersDimension; par++) {
224 RealType sum = NumericTraits< RealType >::Zero;
225 for(unsigned int dim=0; dim<ImageDimension; dim++) {
226 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
228 derivative[par] += sum;
235 if( !this->m_NumberOfPixelsCounted ) {
236 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
238 for(unsigned int i=0; i<ParametersDimension; i++) {
239 derivative[i] /= this->m_NumberOfPixelsCounted;
247 * Get both the match Measure and theDerivative Measure
249 template <class TFixedImage, class TMovingImage>
251 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
252 ::GetValueAndDerivative(const TransformParametersType & parameters,
253 MeasureType & value, DerivativeType & derivative) const
256 itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
258 if( !this->GetGradientImage() ) {
259 itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
262 FixedImageConstPointer fixedImage = this->m_FixedImage;
265 itkExceptionMacro( << "Fixed image has not been assigned" );
268 const unsigned int ImageDimension = FixedImageType::ImageDimension;
270 typedef itk::ImageRegionConstIteratorWithIndex<
271 FixedImageType> FixedIteratorType;
273 typedef itk::ImageRegionConstIteratorWithIndex<GradientImageType> GradientIteratorType;
276 FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
278 typename FixedImageType::IndexType index;
280 MeasureType measure = NumericTraits< MeasureType >::Zero;
282 this->m_NumberOfPixelsCounted = 0;
284 this->SetTransformParameters( parameters );
286 const unsigned int ParametersDimension = this->GetNumberOfParameters();
287 derivative = DerivativeType( ParametersDimension );
288 derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
292 while(!ti.IsAtEnd()) {
294 index = ti.GetIndex();
296 InputPointType inputPoint;
297 fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
299 if( this->m_FixedImageMask && !this->m_FixedImageMask->IsInside( inputPoint ) ) {
304 OutputPointType transformedPoint = this->m_Transform->TransformPoint( inputPoint );
306 if( this->m_MovingImageMask && !this->m_MovingImageMask->IsInside( transformedPoint ) ) {
311 if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) ) {
312 const RealType movingValue = this->m_Interpolator->Evaluate( transformedPoint );
314 const TransformJacobianType & jacobian =
315 this->m_Transform->GetJacobian( inputPoint );
318 const RealType fixedValue = ti.Value();
319 this->m_NumberOfPixelsCounted++;
321 const RealType diff = movingValue - fixedValue;
323 measure += diff * diff;
325 // Get the gradient by NearestNeighboorInterpolation:
326 // which is equivalent to round up the point components.
327 typedef typename OutputPointType::CoordRepType CoordRepType;
328 typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
329 MovingImageContinuousIndexType;
331 MovingImageContinuousIndexType tempIndex;
332 this->m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
334 typename MovingImageType::IndexType mappedIndex;
335 mappedIndex.CopyWithRound( tempIndex );
337 const GradientPixelType gradient =
338 this->GetGradientImage()->GetPixel( mappedIndex );
340 for(unsigned int par=0; par<ParametersDimension; par++) {
341 RealType sum = NumericTraits< RealType >::Zero;
342 for(unsigned int dim=0; dim<ImageDimension; dim++) {
343 sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
345 derivative[par] += sum;
352 if( !this->m_NumberOfPixelsCounted ) {
353 itkExceptionMacro(<<"All the points mapped to outside of the moving image");
355 for(unsigned int i=0; i<ParametersDimension; i++) {
356 derivative[i] /= this->m_NumberOfPixelsCounted;
358 measure /= this->m_NumberOfPixelsCounted;
365 } // end namespace itk