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 #ifndef __clitkOptNormalizedCorrelationImageToImageMetric_txx
20 #define __clitkOptNormalizedCorrelationImageToImageMetric_txx
22 #include "clitkOptNormalizedCorrelationImageToImageMetric.h"
23 #include "itkCovariantVector.h"
24 #include "itkImageRandomConstIteratorWithIndex.h"
25 #include "itkImageRegionConstIterator.h"
26 #include "itkImageRegionIterator.h"
27 #include "itkImageIterator.h"
28 #include "vnl/vnl_math.h"
36 template < class TFixedImage, class TMovingImage >
37 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
38 ::NormalizedCorrelationImageToImageMetric()
40 this->SetComputeGradient(true);
47 m_ThreaderDerivativeF = NULL;
48 m_ThreaderDerivativeM = NULL;
49 this->m_WithinThreadPreProcess = false;
50 this->m_WithinThreadPostProcess = false;
52 // For backward compatibility, the default behavior is to use all the pixels
53 // in the fixed image.
54 this->UseAllPixelsOn();
60 template < class TFixedImage, class TMovingImage >
61 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
62 ::~NormalizedCorrelationImageToImageMetric()
64 if(m_ThreaderSFF != NULL) {
65 delete [] m_ThreaderSFF;
69 if(m_ThreaderSMM != NULL) {
70 delete [] m_ThreaderSMM;
74 if(m_ThreaderSFM != NULL) {
75 delete [] m_ThreaderSFM;
79 if(m_ThreaderSF != NULL) {
80 delete [] m_ThreaderSF;
84 if(m_ThreaderSM != NULL) {
85 delete [] m_ThreaderSM;
89 if(m_ThreaderDerivativeF != NULL) {
90 delete [] m_ThreaderDerivativeF;
92 m_ThreaderDerivativeF = NULL;
94 if(m_ThreaderDerivativeM != NULL) {
95 delete [] m_ThreaderDerivativeM;
97 m_ThreaderDerivativeM = NULL;
101 * Print out internal information about this class
103 template < class TFixedImage, class TMovingImage >
105 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
106 ::PrintSelf(std::ostream& os, itk::Indent indent) const
109 Superclass::PrintSelf(os, indent);
117 template <class TFixedImage, class TMovingImage>
119 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
120 ::Initialize(void) throw ( itk::ExceptionObject )
123 this->Superclass::Initialize();
124 this->Superclass::MultiThreadingInitialize();
128 * Allocate memory for the accumulators (set to zero in GetValue)
130 if(m_ThreaderSFF != NULL) {
131 delete [] m_ThreaderSFF;
133 m_ThreaderSFF = new double[this->m_NumberOfThreads];
136 if(m_ThreaderSMM != NULL) {
137 delete [] m_ThreaderSMM;
139 m_ThreaderSMM = new double[this->m_NumberOfThreads];
141 if(m_ThreaderSFM != NULL) {
142 delete [] m_ThreaderSFM;
144 m_ThreaderSFM = new double[this->m_NumberOfThreads];
146 if(this->m_SubtractMean) {
147 if(m_ThreaderSF != NULL) {
148 delete [] m_ThreaderSF;
150 m_ThreaderSF = new double[this->m_NumberOfThreads];
152 if(m_ThreaderSM != NULL) {
153 delete [] m_ThreaderSM;
155 m_ThreaderSM = new double[this->m_NumberOfThreads];
158 if(m_ThreaderDerivativeF != NULL) {
159 delete [] m_ThreaderDerivativeF;
161 m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfThreads];
162 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
163 m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
166 if(m_ThreaderDerivativeM != NULL) {
167 delete [] m_ThreaderDerivativeM;
169 m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfThreads];
170 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
171 m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
176 template < class TFixedImage, class TMovingImage >
178 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
179 ::GetValueThreadProcessSample(
180 unsigned int threadID,
181 unsigned long fixedImageSample,
182 const MovingImagePointType & itkNotUsed(mappedPoint),
183 double movingImageValue) const
185 const RealType fixedImageValue= this->m_FixedImageSamples[fixedImageSample].value;
186 m_ThreaderSFF[threadID] += fixedImageValue * fixedImageValue;
187 m_ThreaderSMM[threadID] += movingImageValue * movingImageValue;
188 m_ThreaderSFM[threadID] += fixedImageValue * movingImageValue;
189 if ( this->m_SubtractMean ) {
190 m_ThreaderSF[threadID] += fixedImageValue;
191 m_ThreaderSM[threadID] += movingImageValue;
197 template < class TFixedImage, class TMovingImage >
198 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
200 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
201 ::GetValue( const ParametersType & parameters ) const
203 itkDebugMacro("GetValue( " << parameters << " ) ");
205 if( !this->m_FixedImage ) {
206 itkExceptionMacro( << "Fixed image has not been assigned" );
210 //Reset the accumulators
211 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
212 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
213 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
214 if(this->m_SubtractMean) {
215 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
216 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
220 // Set up the parameters in the transform
221 this->m_Transform->SetParameters( parameters );
223 // MUST BE CALLED TO INITIATE PROCESSING
224 this->GetValueMultiThreadedInitiate();
226 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
227 << this->m_NumberOfPixelsCounted << " / "
228 << this->m_NumberOfFixedImageSamples
231 if( this->m_NumberOfPixelsCounted <
232 this->m_NumberOfFixedImageSamples / 4 ) {
233 itkExceptionMacro( "Too many samples map outside moving image buffer: "
234 << this->m_NumberOfPixelsCounted << " / "
235 << this->m_NumberOfFixedImageSamples
239 // Accumulate the threads
240 AccumulateType sff, smm, sfm, sf, sm;
241 sff = m_ThreaderSFF[0];
242 smm = m_ThreaderSMM[0];
243 sfm = m_ThreaderSFM[0];
244 sf = m_ThreaderSF[0];
245 sm = m_ThreaderSM[0];
247 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
248 sff += m_ThreaderSFF[t];
249 smm += m_ThreaderSMM[t];
250 sfm += m_ThreaderSFM[t];
251 if ( this->m_SubtractMean ) {
252 sf += m_ThreaderSF[t];
253 sm += m_ThreaderSM[t];
257 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
258 sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
259 smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
260 sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
264 const RealType denom = -1.0 * vcl_sqrt(sff * smm );
266 if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
267 measure = sfm / denom;
269 measure = itk::NumericTraits< MeasureType >::Zero;
277 template < class TFixedImage, class TMovingImage >
278 typename NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
280 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
281 ::ComputeSums( const ParametersType & parameters ) const
283 //No checking for the fixed image, done in the caller
284 //Reset the accumulators
285 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
286 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
287 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
288 if(this->m_SubtractMean) {
289 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
290 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
294 // Set up the parameters in the transform
295 this->m_Transform->SetParameters( parameters );
297 // MUST BE CALLED TO INITIATE PROCESSING
298 this->GetValueMultiThreadedInitiate();
300 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
301 << this->m_NumberOfPixelsCounted << " / "
302 << this->m_NumberOfFixedImageSamples
305 if( this->m_NumberOfPixelsCounted <
306 this->m_NumberOfFixedImageSamples / 4 ) {
307 itkExceptionMacro( "Too many samples map outside moving image buffer: "
308 << this->m_NumberOfPixelsCounted << " / "
309 << this->m_NumberOfFixedImageSamples
313 // Accumulate the threads
314 m_SFF = m_ThreaderSFF[0];
315 m_SMM = m_ThreaderSMM[0];
316 m_SFM = m_ThreaderSFM[0];
317 m_SF = m_ThreaderSF[0];
318 m_SM = m_ThreaderSM[0];
320 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
321 m_SFF += m_ThreaderSFF[t];
322 m_SMM += m_ThreaderSMM[t];
323 m_SFM += m_ThreaderSFM[t];
324 if ( this->m_SubtractMean ) {
325 m_SF += m_ThreaderSF[t];
326 m_SM += m_ThreaderSM[t];
330 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
331 m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
332 m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
333 m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
334 m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
335 m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
339 m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
341 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
342 measure = m_SFM / m_Denom;
344 measure = itk::NumericTraits< MeasureType >::Zero;
351 template < class TFixedImage, class TMovingImage >
353 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
354 ::GetValueAndDerivativeThreadProcessSample(
355 unsigned int threadID,
356 unsigned long fixedImageSample,
357 const MovingImagePointType & itkNotUsed(mappedPoint),
358 double movingImageValue,
359 const ImageDerivativesType &
360 movingImageGradientValue
364 const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
365 const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
367 // Need to use one of the threader transforms if we're
370 // Use a raw pointer here to avoid the overhead of smart pointers.
371 // For instance, Register and UnRegister have mutex locks around
372 // the reference counts.
373 TransformType* transform;
376 transform = this->m_ThreaderTransform[threadID - 1];
378 transform = this->m_Transform;
381 // Jacobian should be evaluated at the unmapped (fixed image) point.
382 TransformJacobianType jacobian;
383 transform->ComputeJacobianWithRespectToParameters(fixedImagePoint, jacobian);
385 for(unsigned int par=0; par<this->m_NumberOfParameters; par++) {
386 RealType sumF = itk::NumericTraits< RealType >::Zero;
387 RealType sumM = itk::NumericTraits< RealType >::Zero;
388 for(unsigned int dim=0; dim<MovingImageDimension; dim++) {
389 const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
390 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
391 sumF += (fixedImageValue-m_FixedMean) * differential;
392 sumM += (movingImageValue-m_MovingMean) * differential;
394 sumF += differential * fixedImageValue;
395 sumM += differential * movingImageValue;
398 m_ThreaderDerivativeF[threadID][par] += sumF;
399 m_ThreaderDerivativeM[threadID][par] += sumM;
407 * Get the both Value and Derivative Measure
409 template < class TFixedImage, class TMovingImage >
411 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
412 ::GetValueAndDerivative( const ParametersType & parameters,
414 DerivativeType & derivative) const
417 if( !this->m_FixedImage ) {
418 itkExceptionMacro( << "Fixed image has not been assigned" );
421 // Set up the parameters in the transform
422 this->m_Transform->SetParameters( parameters );
424 //We need the sums and the value to be calculated first
425 value=this->ComputeSums(parameters);
427 //Set output values to zero
428 if(derivative.GetSize() != this->m_NumberOfParameters) {
429 derivative = DerivativeType( this->m_NumberOfParameters );
431 memset( derivative.data_block(),
433 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
435 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
436 memset( m_ThreaderDerivativeF[threadID].data_block(),
438 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
440 memset( m_ThreaderDerivativeM[threadID].data_block(),
442 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
445 // MUST BE CALLED TO INITIATE PROCESSING
446 this->GetValueAndDerivativeMultiThreadedInitiate();
448 // Accumulate over the threads
449 DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
450 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
451 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
452 derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
453 derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
457 //Compute derivatives
458 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
459 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
460 derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
463 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
464 derivative[i] = itk::NumericTraits< MeasureType >::Zero;
472 * Get the match measure derivative
474 template < class TFixedImage, class TMovingImage >
476 NormalizedCorrelationImageToImageMetric<TFixedImage,TMovingImage>
477 ::GetDerivative( const ParametersType & parameters,
478 DerivativeType & derivative ) const
480 if( !this->m_FixedImage ) {
481 itkExceptionMacro( << "Fixed image has not been assigned" );
485 // call the combined version
486 this->GetValueAndDerivative( parameters, value, derivative );
489 } // end namespace clitk