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: itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.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 __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
38 #include "itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.h"
39 #include "itkCovariantVector.h"
40 #include "itkImageRandomConstIteratorWithIndex.h"
41 #include "itkImageRegionConstIterator.h"
42 #include "itkImageRegionIterator.h"
43 #include "itkImageIterator.h"
44 #include "vnl/vnl_math.h"
52 template < class TFixedImage, class TMovingImage >
53 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
54 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
56 this->SetComputeGradient(true);
59 m_ThreaderMSEDerivatives = NULL;
60 this->m_WithinThreadPreProcess = false;
61 this->m_WithinThreadPostProcess = false;
63 // For backward compatibility, the default behavior is to use all the pixels
64 // in the fixed image.
65 this->UseAllPixelsOn();
68 template < class TFixedImage, class TMovingImage >
69 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
70 ::~MeanSquaresImageToImageMetricFor3DBLUTFFD()
72 if(m_ThreaderMSE != NULL) {
73 delete [] m_ThreaderMSE;
77 if(m_ThreaderMSEDerivatives != NULL) {
78 delete [] m_ThreaderMSEDerivatives;
80 m_ThreaderMSEDerivatives = NULL;
84 * Print out internal information about this class
86 template < class TFixedImage, class TMovingImage >
88 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
89 ::PrintSelf(std::ostream& os, Indent indent) const
92 Superclass::PrintSelf(os, indent);
100 template <class TFixedImage, class TMovingImage>
102 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
103 ::Initialize(void) throw ( ExceptionObject )
106 this->Superclass::Initialize();
107 this->Superclass::MultiThreadingInitialize();
109 if(m_ThreaderMSE != NULL) {
110 delete [] m_ThreaderMSE;
112 m_ThreaderMSE = new double[this->m_NumberOfThreads];
114 if(m_ThreaderMSEDerivatives != NULL) {
115 delete [] m_ThreaderMSEDerivatives;
117 m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfThreads];
118 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
119 m_ThreaderMSEDerivatives[threadID].SetSize( this->m_NumberOfParameters );
123 template < class TFixedImage, class TMovingImage >
125 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
126 ::GetValueThreadProcessSample( unsigned int threadID,
127 unsigned long fixedImageSample,
128 const MovingImagePointType & itkNotUsed(mappedPoint),
129 double movingImageValue) const
131 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
133 m_ThreaderMSE[threadID] += diff*diff;
138 template < class TFixedImage, class TMovingImage >
139 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
141 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
142 ::GetValue( const ParametersType & parameters ) const
144 itkDebugMacro("GetValue( " << parameters << " ) ");
146 if( !this->m_FixedImage ) {
147 itkExceptionMacro( << "Fixed image has not been assigned" );
150 memset( m_ThreaderMSE,
152 this->m_NumberOfThreads * sizeof(MeasureType) );
154 // Set up the parameters in the transform
155 this->m_Transform->SetParameters( parameters );
157 // MUST BE CALLED TO INITIATE PROCESSING
158 this->GetValueMultiThreadedInitiate();
160 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
161 << this->m_NumberOfPixelsCounted << " / "
162 << this->m_NumberOfFixedImageSamples
165 if( this->m_NumberOfPixelsCounted <
166 this->m_NumberOfFixedImageSamples / 4 ) {
167 itkExceptionMacro( "Too many samples map outside moving image buffer: "
168 << this->m_NumberOfPixelsCounted << " / "
169 << this->m_NumberOfFixedImageSamples
173 double mse = m_ThreaderMSE[0];
174 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
175 mse += m_ThreaderMSE[t];
177 mse /= this->m_NumberOfPixelsCounted;
183 template < class TFixedImage, class TMovingImage >
185 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
186 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
187 unsigned long fixedImageSample,
188 const MovingImagePointType & itkNotUsed(mappedPoint),
189 double movingImageValue,
190 const ImageDerivativesType &
191 movingImageGradientValue ) const
193 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
195 m_ThreaderMSE[threadID] += diff*diff;
198 //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
200 // Need to use one of the threader transforms if we're
203 // Use a raw pointer here to avoid the overhead of smart pointers.
204 // For instance, Register and UnRegister have mutex locks around
205 // the reference counts.
206 TransformType* transform;
209 transform = this->m_ThreaderTransform[threadID - 1];
211 transform = this->m_Transform;
214 // Jacobian should be evaluated at the unmapped (fixed image) point.
215 TransformJacobianType jacobian;
216 transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian);
218 unsigned int par, dim;
219 for( par=0; par<this->m_NumberOfParameters; par+=3) {
221 // for(unsigned int dim=0; dim<MovingImageDimension; dim++)
223 // sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
225 // m_ThreaderMSEDerivatives[threadID][par] += sum;
227 // JV only for 3D BLUT FFD
228 if (jacobian( 0, par ) )
229 for( dim=0; dim<3; dim++)
230 m_ThreaderMSEDerivatives[threadID][par+dim ] += 2.0 * diff* jacobian( dim, par+dim ) * movingImageGradientValue[dim];
237 * Get the both Value and Derivative Measure
239 template < class TFixedImage, class TMovingImage >
241 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
242 ::GetValueAndDerivative( const ParametersType & parameters,
244 DerivativeType & derivative) const
247 if( !this->m_FixedImage ) {
248 itkExceptionMacro( << "Fixed image has not been assigned" );
251 // Set up the parameters in the transform
252 this->m_Transform->SetParameters( parameters );
254 // Reset the joint pdfs to zero
255 memset( m_ThreaderMSE,
257 this->m_NumberOfThreads * sizeof(MeasureType) );
259 // Set output values to zero
260 if(derivative.GetSize() != this->m_NumberOfParameters) {
261 derivative = DerivativeType( this->m_NumberOfParameters );
263 memset( derivative.data_block(),
265 this->m_NumberOfParameters * sizeof(double) );
267 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
268 memset( m_ThreaderMSEDerivatives[threadID].data_block(),
270 this->m_NumberOfParameters * sizeof(double) );
273 // MUST BE CALLED TO INITIATE PROCESSING
274 this->GetValueAndDerivativeMultiThreadedInitiate();
276 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
277 << this->m_NumberOfPixelsCounted << " / "
278 << this->m_NumberOfFixedImageSamples
281 if( this->m_NumberOfPixelsCounted <
282 this->m_NumberOfFixedImageSamples / 4 ) {
283 itkExceptionMacro( "Too many samples map outside moving image buffer: "
284 << this->m_NumberOfPixelsCounted << " / "
285 << this->m_NumberOfFixedImageSamples
290 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
291 value += m_ThreaderMSE[t];
292 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
294 derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
298 value /= this->m_NumberOfPixelsCounted;
299 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
301 derivative[parameter] /= this->m_NumberOfPixelsCounted;
303 //DD(parameter<<"\t"<<derivative[parameter]);
310 * Get the match measure derivative
312 template < class TFixedImage, class TMovingImage >
314 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
315 ::GetDerivative( const ParametersType & parameters,
316 DerivativeType & derivative ) const
318 if( !this->m_FixedImage ) {
319 itkExceptionMacro( << "Fixed image has not been assigned" );
323 // call the combined version
324 this->GetValueAndDerivative( parameters, value, derivative );
327 } // end namespace itk