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>
106 this->Superclass::Initialize();
107 this->Superclass::MultiThreadingInitialize();
109 if(m_ThreaderMSE != NULL) {
110 delete [] m_ThreaderMSE;
112 #if ITK_VERSION_MAJOR <= 4
113 m_ThreaderMSE = new double[this->m_NumberOfThreads];
115 m_ThreaderMSE = new double[this->m_NumberOfWorkUnits];
118 if(m_ThreaderMSEDerivatives != NULL) {
119 delete [] m_ThreaderMSEDerivatives;
121 #if ITK_VERSION_MAJOR <= 4
122 m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfThreads];
123 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
125 m_ThreaderMSEDerivatives = new DerivativeType[this->m_NumberOfWorkUnits];
126 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
128 m_ThreaderMSEDerivatives[threadID].SetSize( this->m_NumberOfParameters );
132 template < class TFixedImage, class TMovingImage >
134 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
135 ::GetValueThreadProcessSample( unsigned int threadID,
136 unsigned long fixedImageSample,
137 const MovingImagePointType & itkNotUsed(mappedPoint),
138 double movingImageValue) const
140 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
142 m_ThreaderMSE[threadID] += diff*diff;
147 template < class TFixedImage, class TMovingImage >
148 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
150 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
151 ::GetValue( const ParametersType & parameters ) const
153 itkDebugMacro("GetValue( " << parameters << " ) ");
155 if( !this->m_FixedImage ) {
156 itkExceptionMacro( << "Fixed image has not been assigned" );
159 #if ITK_VERSION_MAJOR <= 4
160 memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
162 memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
165 // Set up the parameters in the transform
166 this->m_Transform->SetParameters( parameters );
168 // MUST BE CALLED TO INITIATE PROCESSING
169 this->GetValueMultiThreadedInitiate();
171 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
172 << this->m_NumberOfPixelsCounted << " / "
173 << this->m_NumberOfFixedImageSamples
176 if( this->m_NumberOfPixelsCounted <
177 this->m_NumberOfFixedImageSamples / 4 ) {
178 itkExceptionMacro( "Too many samples map outside moving image buffer: "
179 << this->m_NumberOfPixelsCounted << " / "
180 << this->m_NumberOfFixedImageSamples
184 double mse = m_ThreaderMSE[0];
185 #if ITK_VERSION_MAJOR <= 4
186 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
188 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
190 mse += m_ThreaderMSE[t];
192 mse /= this->m_NumberOfPixelsCounted;
198 template < class TFixedImage, class TMovingImage >
200 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
201 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
202 unsigned long fixedImageSample,
203 const MovingImagePointType & itkNotUsed(mappedPoint),
204 double movingImageValue,
205 const ImageDerivativesType &
206 movingImageGradientValue ) const
208 double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
210 m_ThreaderMSE[threadID] += diff*diff;
213 //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
215 // Need to use one of the threader transforms if we're
218 // Use a raw pointer here to avoid the overhead of smart pointers.
219 // For instance, Register and UnRegister have mutex locks around
220 // the reference counts.
221 TransformType* transform;
224 transform = this->m_ThreaderTransform[threadID - 1];
226 transform = this->m_Transform;
229 // Jacobian should be evaluated at the unmapped (fixed image) point.
230 TransformJacobianType jacobian;
231 transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian);
233 unsigned int par, dim;
234 for( par=0; par<this->m_NumberOfParameters; par+=3) {
236 // for(unsigned int dim=0; dim<MovingImageDimension; dim++)
238 // sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
240 // m_ThreaderMSEDerivatives[threadID][par] += sum;
242 // JV only for 3D BLUT FFD
243 if (jacobian( 0, par ) )
244 for( dim=0; dim<3; dim++)
245 m_ThreaderMSEDerivatives[threadID][par+dim ] += 2.0 * diff* jacobian( dim, par+dim ) * movingImageGradientValue[dim];
252 * Get the both Value and Derivative Measure
254 template < class TFixedImage, class TMovingImage >
256 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
257 ::GetValueAndDerivative( const ParametersType & parameters,
259 DerivativeType & derivative) const
262 if( !this->m_FixedImage ) {
263 itkExceptionMacro( << "Fixed image has not been assigned" );
266 // Set up the parameters in the transform
267 this->m_Transform->SetParameters( parameters );
269 // Reset the joint pdfs to zero
270 #if ITK_VERSION_MAJOR <= 4
271 memset( m_ThreaderMSE, 0, this->m_NumberOfThreads * sizeof(MeasureType) );
273 memset( m_ThreaderMSE, 0, this->m_NumberOfWorkUnits * sizeof(MeasureType) );
276 // Set output values to zero
277 if(derivative.GetSize() != this->m_NumberOfParameters) {
278 derivative = DerivativeType( this->m_NumberOfParameters );
280 memset( derivative.data_block(),
282 this->m_NumberOfParameters * sizeof(double) );
284 #if ITK_VERSION_MAJOR <= 4
285 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
287 for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
289 memset( m_ThreaderMSEDerivatives[threadID].data_block(),
291 this->m_NumberOfParameters * sizeof(double) );
294 // MUST BE CALLED TO INITIATE PROCESSING
295 this->GetValueAndDerivativeMultiThreadedInitiate();
297 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
298 << this->m_NumberOfPixelsCounted << " / "
299 << this->m_NumberOfFixedImageSamples
302 if( this->m_NumberOfPixelsCounted <
303 this->m_NumberOfFixedImageSamples / 4 ) {
304 itkExceptionMacro( "Too many samples map outside moving image buffer: "
305 << this->m_NumberOfPixelsCounted << " / "
306 << this->m_NumberOfFixedImageSamples
311 #if ITK_VERSION_MAJOR <= 4
312 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
314 for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
316 value += m_ThreaderMSE[t];
317 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
319 derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
323 value /= this->m_NumberOfPixelsCounted;
324 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
326 derivative[parameter] /= this->m_NumberOfPixelsCounted;
328 //DD(parameter<<"\t"<<derivative[parameter]);
335 * Get the match measure derivative
337 template < class TFixedImage, class TMovingImage >
339 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
340 ::GetDerivative( const ParametersType & parameters,
341 DerivativeType & derivative ) const
343 if( !this->m_FixedImage ) {
344 itkExceptionMacro( << "Fixed image has not been assigned" );
348 // call the combined version
349 this->GetValueAndDerivative( parameters, value, derivative );
352 } // end namespace itk