]> Creatis software - clitk.git/blob - registration/itkOptMeanSquaresImageToImageMetricFor3DBLUTFFD.txx
min max (to remove)
[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://oncora1.lyon.fnclcc.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   this->m_Parameters = parameters;
157
158   // MUST BE CALLED TO INITIATE PROCESSING
159   this->GetValueMultiThreadedInitiate();
160
161   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
162                  << this->m_NumberOfPixelsCounted << " / "
163                  << this->m_NumberOfFixedImageSamples
164                  << std::endl );
165
166   if( this->m_NumberOfPixelsCounted <
167       this->m_NumberOfFixedImageSamples / 4 ) {
168     itkExceptionMacro( "Too many samples map outside moving image buffer: "
169                        << this->m_NumberOfPixelsCounted << " / "
170                        << this->m_NumberOfFixedImageSamples
171                        << std::endl );
172   }
173
174   double mse = m_ThreaderMSE[0];
175   for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
176     mse += m_ThreaderMSE[t];
177   }
178   mse /= this->m_NumberOfPixelsCounted;
179
180   return mse;
181 }
182
183
184 template < class TFixedImage, class TMovingImage  >
185 inline bool
186 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
187 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
188     unsigned long fixedImageSample,
189     const MovingImagePointType & itkNotUsed(mappedPoint),
190     double movingImageValue,
191     const ImageDerivativesType &
192     movingImageGradientValue ) const
193 {
194   double diff = movingImageValue - this->m_FixedImageSamples[fixedImageSample].value;
195
196   m_ThreaderMSE[threadID] += diff*diff;
197
198   //JV
199   //FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
200
201   // Need to use one of the threader transforms if we're
202   // not in thread 0.
203   //
204   // Use a raw pointer here to avoid the overhead of smart pointers.
205   // For instance, Register and UnRegister have mutex locks around
206   // the reference counts.
207   TransformType* transform;
208
209   if (threadID > 0) {
210     transform = this->m_ThreaderTransform[threadID - 1];
211   } else {
212     transform = this->m_Transform;
213   }
214
215   // Jacobian should be evaluated at the unmapped (fixed image) point.
216   const TransformJacobianType & jacobian = transform ->GetJacobian( this->m_FixedImageSamples[fixedImageSample].point );
217   //double sum;
218   unsigned int par, dim;
219   for( par=0; par<this->m_NumberOfParameters; par+=3) {
220     //     double sum = 0.0;
221     //     for(unsigned int dim=0; dim<MovingImageDimension; dim++)
222     //       {
223     //       sum += 2.0 * diff * jacobian( dim, par ) * movingImageGradientValue[dim];
224     //       }
225     //     m_ThreaderMSEDerivatives[threadID][par] += sum;
226
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];
231
232   }
233   return true;
234 }
235
236 /**
237  * Get the both Value and Derivative Measure
238  */
239 template < class TFixedImage, class TMovingImage  >
240 void
241 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
242 ::GetValueAndDerivative( const ParametersType & parameters,
243                          MeasureType & value,
244                          DerivativeType & derivative) const
245 {
246
247   if( !this->m_FixedImage ) {
248     itkExceptionMacro( << "Fixed image has not been assigned" );
249   }
250
251   // Set up the parameters in the transform
252   this->m_Transform->SetParameters( parameters );
253   this->m_Parameters = parameters;
254
255   // Reset the joint pdfs to zero
256   memset( m_ThreaderMSE,
257           0,
258           this->m_NumberOfThreads * sizeof(MeasureType) );
259
260   // Set output values to zero
261   if(derivative.GetSize() != this->m_NumberOfParameters) {
262     derivative = DerivativeType( this->m_NumberOfParameters );
263   }
264   memset( derivative.data_block(),
265           0,
266           this->m_NumberOfParameters * sizeof(double) );
267
268   for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
269     memset( m_ThreaderMSEDerivatives[threadID].data_block(),
270             0,
271             this->m_NumberOfParameters * sizeof(double) );
272   }
273
274   // MUST BE CALLED TO INITIATE PROCESSING
275   this->GetValueAndDerivativeMultiThreadedInitiate();
276
277   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
278                  << this->m_NumberOfPixelsCounted << " / "
279                  << this->m_NumberOfFixedImageSamples
280                  << std::endl );
281
282   if( this->m_NumberOfPixelsCounted <
283       this->m_NumberOfFixedImageSamples / 4 ) {
284     itkExceptionMacro( "Too many samples map outside moving image buffer: "
285                        << this->m_NumberOfPixelsCounted << " / "
286                        << this->m_NumberOfFixedImageSamples
287                        << std::endl );
288   }
289
290   value = 0;
291   for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
292     value += m_ThreaderMSE[t];
293     for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
294         parameter++) {
295       derivative[parameter] += m_ThreaderMSEDerivatives[t][parameter];
296     }
297   }
298
299   value /= this->m_NumberOfPixelsCounted;
300   for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters;
301       parameter++) {
302     derivative[parameter] /= this->m_NumberOfPixelsCounted;
303     // JV
304     //DD(parameter<<"\t"<<derivative[parameter]);
305   }
306
307 }
308
309
310 /**
311  * Get the match measure derivative
312  */
313 template < class TFixedImage, class TMovingImage  >
314 void
315 MeanSquaresImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
316 ::GetDerivative( const ParametersType & parameters,
317                  DerivativeType & derivative ) const
318 {
319   if( !this->m_FixedImage ) {
320     itkExceptionMacro( << "Fixed image has not been assigned" );
321   }
322
323   MeasureType value;
324   // call the combined version
325   this->GetValueAndDerivative( parameters, value, derivative );
326 }
327
328 } // end namespace itk
329
330
331 #endif