]> Creatis software - clitk.git/blob - registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx
Merge branch 'VTK6_Qt5_Overlay4D' into VTK6_Qt5_Binarize
[clitk.git] / registration / itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
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
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 /*=========================================================================
20
21   Program:   Insight Segmentation & Registration Toolkit
22   Module:    $RCSfile: itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx,v $
23   Language:  C++
24   Date:      $Date: 2010/06/14 17:32:07 $
25   Version:   $Revision: 1.1 $
26
27   Copyright (c) Insight Software Consortium. All rights reserved.
28   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
29
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.
33
34 =========================================================================*/
35 #ifndef __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD_txx
37
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"
45
46 namespace itk
47 {
48
49 /**
50  * Constructor
51  */
52 template < class TFixedImage, class TMovingImage >
53 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
54 ::MeanSquaresImageToImageMetricFor3DBLUTFFD()
55 {
56   this->SetComputeGradient(true);
57
58   m_ThreaderMSE = NULL;
59   m_ThreaderMSEDerivatives = NULL;
60   this->m_WithinThreadPreProcess = false;
61   this->m_WithinThreadPostProcess = false;
62
63   //  For backward compatibility, the default behavior is to use all the pixels
64   //  in the fixed image.
65   this->UseAllPixelsOn();
66 }
67
68 template < class TFixedImage, class TMovingImage >
69 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
70 ::~MeanSquaresImageToImageMetricFor3DBLUTFFD()
71 {
72   if(m_ThreaderMSE != NULL) {
73     delete [] m_ThreaderMSE;
74   }
75   m_ThreaderMSE = NULL;
76
77   if(m_ThreaderMSEDerivatives != NULL) {
78     delete [] m_ThreaderMSEDerivatives;
79   }
80   m_ThreaderMSEDerivatives = NULL;
81 }
82
83 /**
84  * Print out internal information about this class
85  */
86 template < class TFixedImage, class TMovingImage  >
87 void
88 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
89 ::PrintSelf(std::ostream& os, Indent indent) const
90 {
91
92   Superclass::PrintSelf(os, indent);
93
94 }
95
96
97 /**
98  * Initialize
99  */
100 template <class TFixedImage, class TMovingImage>
101 void
102 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
103 ::Initialize(void) throw ( ExceptionObject )
104 {
105
106   this->Superclass::Initialize();
107   this->Superclass::MultiThreadingInitialize();
108
109   if(m_ThreaderMSE != NULL) {
110     delete [] m_ThreaderMSE;
111   }
112   m_ThreaderMSE = new double[this->m_NumberOfThreads];
113
114   if(m_ThreaderMSEDerivatives != NULL) {
115     delete [] m_ThreaderMSEDerivatives;
116   }
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 );
120   }
121 }
122
123 template < class TFixedImage, class TMovingImage  >
124 inline bool
125 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
126 ::GetValueThreadProcessSample( unsigned int threadID,
127                                unsigned long fixedImageSample,
128                                const MovingImagePointType & itkNotUsed(mappedPoint),
129                                double movingImageValue) const
130 {
131   double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
132
133   m_ThreaderMSE[threadID] += diff*diff;
134
135   return true;
136 }
137
138 template < class TFixedImage, class TMovingImage  >
139 typename MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
140 ::MeasureType
141 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
142 ::GetValue( const ParametersType & parameters ) const
143 {
144   itkDebugMacro("GetValue( " << parameters << " ) ");
145
146   if( !this->m_FixedImage ) {
147     itkExceptionMacro( << "Fixed image has not been assigned" );
148   }
149
150   memset( m_ThreaderMSE,
151           0,
152           this->m_NumberOfThreads * sizeof(MeasureType) );
153
154   // Set up the parameters in the transform
155   this->m_Transform->SetParameters( parameters );
156 #if ITK_VERSION_MAJOR < 4
157   this->m_Parameters = parameters;
158 #endif
159
160   // MUST BE CALLED TO INITIATE PROCESSING
161   this->GetValueMultiThreadedInitiate();
162
163   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
164                  << this->m_NumberOfPixelsCounted << " / "
165                  << this->m_NumberOfFixedImageSamples
166                  << std::endl );
167
168   if( this->m_NumberOfPixelsCounted <
169       this->m_NumberOfFixedImageSamples / 4 ) {
170     itkExceptionMacro( "Too many samples map outside moving image buffer: "
171                        << this->m_NumberOfPixelsCounted << " / "
172                        << this->m_NumberOfFixedImageSamples
173                        << std::endl );
174   }
175
176   double mse = m_ThreaderMSE[0];
177   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
178     mse += m_ThreaderMSE[t];
179   }
180   mse /= this->m_NumberOfPixelsCounted;
181
182   return mse;
183 }
184
185
186 template < class TFixedImage, class TMovingImage  >
187 inline bool
188 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
189 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
190     unsigned long fixedImageSample,
191     const MovingImagePointType & itkNotUsed(mappedPoint),
192     double movingImageValue,
193     const ImageDerivativesType &
194     movingImageGradientValue ) const
195 {
196   double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
197
198   m_ThreaderMSE[threadID] += diff*diff;
199
200   //JV
201   //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
202
203   // Need to use one of the threader transforms if we're
204   // not in thread 0.
205   //
206   // Use a raw pointer here to avoid the overhead of smart pointers.
207   // For instance, Register and UnRegister have mutex locks around
208   // the reference counts.
209   TransformType* transform;
210
211   if (threadID > 0) {
212     transform = this->m_ThreaderTransform[threadID - 1];
213   } else {
214     transform = this->m_Transform;
215   }
216
217   // Jacobian should be evaluated at the unmapped (fixed image) point.
218 #if ITK_VERSION_MAJOR >= 4
219   TransformJacobianType jacobian;
220   transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[fixedImageSample].point, jacobian);
221 #else
222   const TransformJacobianType & jacobian = transform ->GetJacobian( this->m_FixedImageSamples[fixedImageSample].point );
223 #endif
224   //double sum;
225   unsigned int par, dim;
226   for( par=0; par<this->m_NumberOfParameters; par+=3) {
227     //     double sum = 0.0;
228     //     for(unsigned int dim=0; dim<MovingImageDimension; dim++)
229     //       {
230     //       sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
231     //       }
232     //     m_ThreaderMSEDerivatives[threadID][par] += sum;
233
234     // JV only for 3D BLUT FFD
235     if (jacobian( 0, par ) )
236       for( dim=0; dim<3; dim++)
237         m_ThreaderMSEDerivatives[threadID][par+dim  ]  += 2.0 * diff* jacobian( dim, par+dim   ) * movingImageGradientValue[dim];
238
239   }
240   return true;
241 }
242
243 /**
244  * Get the both Value and Derivative Measure
245  */
246 template < class TFixedImage, class TMovingImage  >
247 void
248 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
249 ::GetValueAndDerivative( const ParametersType & parameters,
250                          MeasureType & value,
251                          DerivativeType & derivative) const
252 {
253
254   if( !this->m_FixedImage ) {
255     itkExceptionMacro( << "Fixed image has not been assigned" );
256   }
257
258   // Set up the parameters in the transform
259   this->m_Transform->SetParameters( parameters );
260 #if ITK_VERSION_MAJOR < 4
261   this->m_Parameters = parameters;
262 #endif
263
264   // Reset the joint pdfs to zero
265   memset( m_ThreaderMSE,
266           0,
267           this->m_NumberOfThreads * sizeof(MeasureType) );
268
269   // Set output values to zero
270   if(derivative.GetSize() != this->m_NumberOfParameters) {
271     derivative = DerivativeType( this->m_NumberOfParameters );
272   }
273   memset( derivative.data_block(),
274           0,
275           this->m_NumberOfParameters * sizeof(double) );
276
277   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
278     memset( m_ThreaderMSEDerivatives[threadID].data_block(),
279             0,
280             this->m_NumberOfParameters * sizeof(double) );
281   }
282
283   // MUST BE CALLED TO INITIATE PROCESSING
284   this->GetValueAndDerivativeMultiThreadedInitiate();
285
286   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
287                  << this->m_NumberOfPixelsCounted << " / "
288                  << this->m_NumberOfFixedImageSamples
289                  << std::endl );
290
291   if( this->m_NumberOfPixelsCounted <
292       this->m_NumberOfFixedImageSamples / 4 ) {
293     itkExceptionMacro( "Too many samples map outside moving image buffer: "
294                        << this->m_NumberOfPixelsCounted << " / "
295                        << this->m_NumberOfFixedImageSamples
296                        << std::endl );
297   }
298
299   value = 0;
300   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
301     value += m_ThreaderMSE[t];
302     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
303         parameter++) {
304       derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
305     }
306   }
307
308   value /= this->m_NumberOfPixelsCounted;
309   for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
310       parameter++) {
311     derivative[parameter] /= this->m_NumberOfPixelsCounted;
312     // JV
313     //DD(parameter<<"\t"<<derivative[parameter]);
314   }
315
316 }
317
318
319 /**
320  * Get the match measure derivative
321  */
322 template < class TFixedImage, class TMovingImage  >
323 void
324 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
325 ::GetDerivative( const ParametersType & parameters,
326                  DerivativeType & derivative ) const
327 {
328   if( !this->m_FixedImage ) {
329     itkExceptionMacro( << "Fixed image has not been assigned" );
330   }
331
332   MeasureType value;
333   // call the combined version
334   this->GetValueAndDerivative( parameters, value, derivative );
335 }
336
337 } // end namespace itk
338
339
340 #endif