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 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
106 ::Initialize(void) throw ( ExceptionObject )
110 this->Superclass::Initialize();
111 this->Superclass::MultiThreadingInitialize();
113 if(m_ThreaderMSE != NULL) {
114 delete [] m_ThreaderMSE;
116 #if ITK_VERSION_MAJOR <= 4
117 m_ThreaderMSE = new double[this->m_NumberOfThreads];
119 m_ThreaderMSE = new double[this->m_NumberOfWorkUnits];
122 if(m_ThreaderMSEDerivatives != NULL) {
123 delete [] m_ThreaderMSEDerivatives;
125 #if ITK_VERSION_MAJOR <= 4
126 m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfThreads];
127 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
129 m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfWorkUnits];
130 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
132 m_ThreaderMSEDerivatives[threadID].SetSize( this->m_NumberOfParameters );
136 template < class TFixedImage, class TMovingImage >
138 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
139 ::GetValueThreadProcessSample( unsigned int threadID,
140 unsigned long fixedImageSample,
141 const MovingImagePointType & itkNotUsed(mappedPoint),
142 double movingImageValue) const
144 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
146 m_ThreaderMSE[threadID] += diff*diff;
151 template < class TFixedImage, class TMovingImage >
152 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
154 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
155 ::GetValue( const ParametersType & parameters ) const
157 itkDebugMacro("GetValue( " << parameters << " ) ");
159 if( !this->m_FixedImage ) {
160 itkExceptionMacro( << "Fixed image has not been assigned" );
163 #if ITK_VERSION_MAJOR <= 4
164 memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
166 memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
169 // Set up the parameters in the transform
170 this->m_Transform->SetParameters( parameters );
172 // MUST BE CALLED TO INITIATE PROCESSING
173 this->GetValueMultiThreadedInitiate();
175 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
176 << this->m_NumberOfPixelsCounted << " / "
177 << this->m_NumberOfFixedImageSamples
180 if( this->m_NumberOfPixelsCounted <
181 this->m_NumberOfFixedImageSamples / 4 ) {
182 itkExceptionMacro( "Too many samples map outside moving image buffer: "
183 << this->m_NumberOfPixelsCounted << " / "
184 << this->m_NumberOfFixedImageSamples
188 double mse = m_ThreaderMSE[0];
189 #if ITK_VERSION_MAJOR <= 4
190 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
192 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
194 mse += m_ThreaderMSE[t];
196 mse /= this->m_NumberOfPixelsCounted;
202 template < class TFixedImage, class TMovingImage >
204 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
205 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
206 unsigned long fixedImageSample,
207 const MovingImagePointType & itkNotUsed(mappedPoint),
208 double movingImageValue,
209 const ImageDerivativesType &
210 movingImageGradientValue ) const
212 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
214 m_ThreaderMSE[threadID] += diff*diff;
217 //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
219 // Need to use one of the threader transforms if we're
222 // Use a raw pointer here to avoid the overhead of smart pointers.
223 // For instance, Register and UnRegister have mutex locks around
224 // the reference counts.
225 TransformType* transform;
228 transform = this->m_ThreaderTransform[threadID - 1];
230 transform = this->m_Transform;
233 // Jacobian should be evaluated at the unmapped (fixed image) point.
234 TransformJacobianType jacobian;
235 transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian);
237 unsigned int par, dim;
238 for( par=0; par<this->m_NumberOfParameters; par+=3) {
240 // for(unsigned int dim=0; dim<MovingImageDimension; dim++)
242 // sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
244 // m_ThreaderMSEDerivatives[threadID][par] += sum;
246 // JV only for 3D BLUT FFD
247 if (jacobian( 0, par ) )
248 for( dim=0; dim<3; dim++)
249 m_ThreaderMSEDerivatives[threadID][par+dim ] += 2.0 * diff* jacobian( dim, par+dim ) * movingImageGradientValue[dim];
256 * Get the both Value and Derivative Measure
258 template < class TFixedImage, class TMovingImage >
260 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
261 ::GetValueAndDerivative( const ParametersType & parameters,
263 DerivativeType & derivative) const
266 if( !this->m_FixedImage ) {
267 itkExceptionMacro( << "Fixed image has not been assigned" );
270 // Set up the parameters in the transform
271 this->m_Transform->SetParameters( parameters );
273 // Reset the joint pdfs to zero
274 #if ITK_VERSION_MAJOR <= 4
275 memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
277 memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
280 // Set output values to zero
281 if(derivative.GetSize() != this->m_NumberOfParameters) {
282 derivative = DerivativeType( this->m_NumberOfParameters );
284 memset( derivative.data_block(),
286 this->m_NumberOfParameters * sizeof(double) );
288 #if ITK_VERSION_MAJOR <= 4
289 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
291 for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
293 memset( m_ThreaderMSEDerivatives[threadID].data_block(),
295 this->m_NumberOfParameters * sizeof(double) );
298 // MUST BE CALLED TO INITIATE PROCESSING
299 this->GetValueAndDerivativeMultiThreadedInitiate();
301 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
302 << this->m_NumberOfPixelsCounted << " / "
303 << this->m_NumberOfFixedImageSamples
306 if( this->m_NumberOfPixelsCounted <
307 this->m_NumberOfFixedImageSamples / 4 ) {
308 itkExceptionMacro( "Too many samples map outside moving image buffer: "
309 << this->m_NumberOfPixelsCounted << " / "
310 << this->m_NumberOfFixedImageSamples
315 #if ITK_VERSION_MAJOR <= 4
316 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
318 for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
320 value += m_ThreaderMSE[t];
321 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
323 derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
327 value /= this->m_NumberOfPixelsCounted;
328 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
330 derivative[parameter] /= this->m_NumberOfPixelsCounted;
332 //DD(parameter<<"\t"<<derivative[parameter]);
339 * Get the match measure derivative
341 template < class TFixedImage, class TMovingImage >
343 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
344 ::GetDerivative( const ParametersType & parameters,
345 DerivativeType & derivative ) const
347 if( !this->m_FixedImage ) {
348 itkExceptionMacro( << "Fixed image has not been assigned" );
352 // call the combined version
353 this->GetValueAndDerivative( parameters, value, derivative );
356 } // end namespace itk