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 __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_txx
20 #define __clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD_txx
22 #include "clitkOptNormalizedCorrelationImageToImageMetricFor3DBLUTFFD.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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
38 ::NormalizedCorrelationImageToImageMetricFor3DBLUTFFD()
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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
62 ::~NormalizedCorrelationImageToImageMetricFor3DBLUTFFD()
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 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
106 ::PrintSelf(std::ostream& os, itk::Indent indent) const
109 Superclass::PrintSelf(os, indent);
117 template <class TFixedImage, class TMovingImage>
119 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
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 #if ITK_VERSION_MAJOR <= 4
134 m_ThreaderSFF = new double[this->m_NumberOfThreads];
136 m_ThreaderSFF = new double[this->m_NumberOfWorkUnits];
140 if(m_ThreaderSMM != NULL) {
141 delete [] m_ThreaderSMM;
143 #if ITK_VERSION_MAJOR <= 4
144 m_ThreaderSMM = new double[this->m_NumberOfThreads];
146 m_ThreaderSMM = new double[this->m_NumberOfWorkUnits];
149 if(m_ThreaderSFM != NULL) {
150 delete [] m_ThreaderSFM;
152 #if ITK_VERSION_MAJOR <= 4
153 m_ThreaderSFM = new double[this->m_NumberOfThreads];
155 m_ThreaderSFM = new double[this->m_NumberOfWorkUnits];
158 if(this->m_SubtractMean) {
159 if(m_ThreaderSF != NULL) {
160 delete [] m_ThreaderSF;
162 #if ITK_VERSION_MAJOR <= 4
163 m_ThreaderSF = new double[this->m_NumberOfThreads];
165 m_ThreaderSF = new double[this->m_NumberOfWorkUnits];
168 if(m_ThreaderSM != NULL) {
169 delete [] m_ThreaderSM;
171 #if ITK_VERSION_MAJOR <= 4
172 m_ThreaderSM = new double[this->m_NumberOfThreads];
174 m_ThreaderSM = new double[this->m_NumberOfWorkUnits];
178 if(m_ThreaderDerivativeF != NULL) {
179 delete [] m_ThreaderDerivativeF;
181 #if ITK_VERSION_MAJOR <= 4
182 m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfThreads];
183 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
185 m_ThreaderDerivativeF = new DerivativeType[this->m_NumberOfWorkUnits];
186 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
188 m_ThreaderDerivativeF[threadID].SetSize( this->m_NumberOfParameters );
191 if(m_ThreaderDerivativeM != NULL) {
192 delete [] m_ThreaderDerivativeM;
194 #if ITK_VERSION_MAJOR <= 4
195 m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfThreads];
196 for(unsigned int threadID=0; threadID<this->m_NumberOfThreads; threadID++) {
198 m_ThreaderDerivativeM = new DerivativeType[this->m_NumberOfWorkUnits];
199 for(unsigned int threadID=0; threadID<this->m_NumberOfWorkUnits; threadID++) {
201 m_ThreaderDerivativeM[threadID].SetSize( this->m_NumberOfParameters );
206 template < class TFixedImage, class TMovingImage >
208 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
209 ::GetValueThreadProcessSample(
210 unsigned int threadID,
211 unsigned long fixedImageSample,
212 const MovingImagePointType & itkNotUsed(mappedPoint),
213 double movingImageValue) const
215 const RealType fixedImageValue= this->m_FixedImageSamples[fixedImageSample].value;
216 m_ThreaderSFF[threadID] += fixedImageValue * fixedImageValue;
217 m_ThreaderSMM[threadID] += movingImageValue * movingImageValue;
218 m_ThreaderSFM[threadID] += fixedImageValue * movingImageValue;
219 if ( this->m_SubtractMean ) {
220 m_ThreaderSF[threadID] += fixedImageValue;
221 m_ThreaderSM[threadID] += movingImageValue;
227 template < class TFixedImage, class TMovingImage >
228 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
230 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
231 ::GetValue( const ParametersType & parameters ) const
233 itkDebugMacro("GetValue( " << parameters << " ) ");
235 if( !this->m_FixedImage ) {
236 itkExceptionMacro( << "Fixed image has not been assigned" );
240 //Reset the accumulators
241 #if ITK_VERSION_MAJOR <= 4
242 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
243 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
244 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
246 memset( m_ThreaderSFF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
247 memset( m_ThreaderSMM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
248 memset( m_ThreaderSFM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
250 if(this->m_SubtractMean) {
251 #if ITK_VERSION_MAJOR <= 4
252 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
253 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
255 memset( m_ThreaderSF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
256 memset( m_ThreaderSM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
261 // Set up the parameters in the transform
262 this->m_Transform->SetParameters( parameters );
264 // MUST BE CALLED TO INITIATE PROCESSING
265 this->GetValueMultiThreadedInitiate();
267 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
268 << this->m_NumberOfPixelsCounted << " / "
269 << this->m_NumberOfFixedImageSamples
272 if( this->m_NumberOfPixelsCounted <
273 this->m_NumberOfFixedImageSamples / 4 ) {
274 itkExceptionMacro( "Too many samples map outside moving image buffer: "
275 << this->m_NumberOfPixelsCounted << " / "
276 << this->m_NumberOfFixedImageSamples
280 // Accumulate the threads
281 AccumulateType sff, smm, sfm, sf, sm;
282 sff = m_ThreaderSFF[0];
283 smm = m_ThreaderSMM[0];
284 sfm = m_ThreaderSFM[0];
285 sf = m_ThreaderSF[0];
286 sm = m_ThreaderSM[0];
288 #if ITK_VERSION_MAJOR <= 4
289 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
291 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
293 sff += m_ThreaderSFF[t];
294 smm += m_ThreaderSMM[t];
295 sfm += m_ThreaderSFM[t];
296 if ( this->m_SubtractMean ) {
297 sf += m_ThreaderSF[t];
298 sm += m_ThreaderSM[t];
302 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
303 sff -= ( sf * sf / this->m_NumberOfPixelsCounted );
304 smm -= ( sm * sm / this->m_NumberOfPixelsCounted );
305 sfm -= ( sf * sm / this->m_NumberOfPixelsCounted );
309 const RealType denom = -1.0 * vcl_sqrt(sff * smm );
311 if( this->m_NumberOfPixelsCounted > 0 && denom != 0.0) {
312 measure = sfm / denom;
314 measure = itk::NumericTraits< MeasureType >::Zero;
322 template < class TFixedImage, class TMovingImage >
323 typename NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
325 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
326 ::ComputeSums( const ParametersType & parameters ) const
328 //No checking for the fixed image, done in the caller
329 //Reset the accumulators
330 #if ITK_VERSION_MAJOR <= 4
331 memset( m_ThreaderSFF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
332 memset( m_ThreaderSMM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
333 memset( m_ThreaderSFM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
335 memset( m_ThreaderSFF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
336 memset( m_ThreaderSMM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
337 memset( m_ThreaderSFM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
339 if(this->m_SubtractMean) {
340 #if ITK_VERSION_MAJOR <= 4
341 memset( m_ThreaderSF, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
342 memset( m_ThreaderSM, 0, this->m_NumberOfThreads * sizeof(AccumulateType) );
344 memset( m_ThreaderSF, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
345 memset( m_ThreaderSM, 0, this->m_NumberOfWorkUnits * sizeof(AccumulateType) );
350 // Set up the parameters in the transform
351 this->m_Transform->SetParameters( parameters );
353 // MUST BE CALLED TO INITIATE PROCESSING
354 this->GetValueMultiThreadedInitiate();
356 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
357 << this->m_NumberOfPixelsCounted << " / "
358 << this->m_NumberOfFixedImageSamples
361 if( this->m_NumberOfPixelsCounted <
362 this->m_NumberOfFixedImageSamples / 4 ) {
363 itkExceptionMacro( "Too many samples map outside moving image buffer: "
364 << this->m_NumberOfPixelsCounted << " / "
365 << this->m_NumberOfFixedImageSamples
369 // Accumulate the threads
370 m_SFF = m_ThreaderSFF[0];
371 m_SMM = m_ThreaderSMM[0];
372 m_SFM = m_ThreaderSFM[0];
373 m_SF = m_ThreaderSF[0];
374 m_SM = m_ThreaderSM[0];
376 #if ITK_VERSION_MAJOR <= 4
377 for(unsigned int t=1; t<this->m_NumberOfThreads; t++) {
379 for(unsigned int t=1; t<this->m_NumberOfWorkUnits; t++) {
381 m_SFF += m_ThreaderSFF[t];
382 m_SMM += m_ThreaderSMM[t];
383 m_SFM += m_ThreaderSFM[t];
384 if ( this->m_SubtractMean ) {
385 m_SF += m_ThreaderSF[t];
386 m_SM += m_ThreaderSM[t];
390 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
391 m_SFF -= ( m_SF * m_SF / this->m_NumberOfPixelsCounted );
392 m_SMM -= ( m_SM * m_SM / this->m_NumberOfPixelsCounted );
393 m_SFM -= ( m_SF * m_SM / this->m_NumberOfPixelsCounted );
394 m_FixedMean=m_SF / this->m_NumberOfPixelsCounted;
395 m_MovingMean=m_SM / this->m_NumberOfPixelsCounted;
399 m_Denom = -1.0 * vcl_sqrt(m_SFF * m_SMM );
401 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
402 measure = m_SFM / m_Denom;
404 measure = itk::NumericTraits< MeasureType >::Zero;
411 template < class TFixedImage, class TMovingImage >
413 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
414 ::GetValueAndDerivativeThreadProcessSample(
415 unsigned int threadID,
416 unsigned long fixedImageSample,
417 const MovingImagePointType & itkNotUsed(mappedPoint),
418 double movingImageValue,
419 const ImageDerivativesType &
420 movingImageGradientValue
424 const RealType fixedImageValue=this->m_FixedImageSamples[fixedImageSample].value;
425 const FixedImagePointType fixedImagePoint = this->m_FixedImageSamples[fixedImageSample].point;
427 // Need to use one of the threader transforms if we're
430 // Use a raw pointer here to avoid the overhead of smart pointers.
431 // For instance, Register and UnRegister have mutex locks around
432 // the reference counts.
433 TransformType* transform;
436 transform = this->m_ThreaderTransform[threadID - 1];
438 transform = this->m_Transform;
441 // Jacobian should be evaluated at the unmapped (fixed image) point.
442 TransformJacobianType jacobian;
443 transform->ComputeJacobianWithRespectToParameters( fixedImagePoint, jacobian );
445 // for(unsigned int par=0; par<this->m_NumberOfParameters; par++)
447 // RealType sumF = itk::NumericTraits< RealType >::Zero;
448 // RealType sumM = itk::NumericTraits< RealType >::Zero;
449 // for(unsigned int dim=0; dim<MovingImageDimension; dim++)
451 // const RealType differential = jacobian( dim, par ) * movingImageGradientValue[dim];
452 // if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 )
454 // sumF += (fixedImageValue-m_FixedMean) * differential;
455 // sumM += (movingImageValue-m_MovingMean) * differential;
459 // sumF += differential * fixedImageValue;
460 // sumM += differential * movingImageValue;
463 // m_ThreaderDerivativeF[threadID][par] += sumF;
464 // m_ThreaderDerivativeM[threadID][par] += sumM;
468 unsigned int par, dim;
469 RealType differential;
470 for( par=0; par<this->m_NumberOfParameters; par+=3) {
471 // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
472 if (jacobian( 0, par ) ) {
473 for(dim=0; dim<3; dim++) {
474 differential = jacobian( dim, par+dim ) * movingImageGradientValue[dim];
476 if ( this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0 ) {
477 m_ThreaderDerivativeF[threadID][par+dim]+= (fixedImageValue-m_FixedMean) * differential;
478 m_ThreaderDerivativeM[threadID][par+dim]+= (movingImageValue-m_MovingMean) * differential;
480 m_ThreaderDerivativeF[threadID][par+dim]+= differential * fixedImageValue;
481 m_ThreaderDerivativeM[threadID][par+dim]+= differential * movingImageValue;
492 * Get the both Value and Derivative Measure
494 template < class TFixedImage, class TMovingImage >
496 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
497 ::GetValueAndDerivative( const ParametersType & parameters,
499 DerivativeType & derivative) const
502 if( !this->m_FixedImage ) {
503 itkExceptionMacro( << "Fixed image has not been assigned" );
506 // Set up the parameters in the transform
507 this->m_Transform->SetParameters( parameters );
509 //We need the sums and the value to be calculated first
510 value=this->ComputeSums(parameters);
512 //Set output values to zero
513 if(derivative.GetSize() != this->m_NumberOfParameters) {
514 derivative = DerivativeType( this->m_NumberOfParameters );
516 memset( derivative.data_block(),
518 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
520 #if ITK_VERSION_MAJOR <= 4
521 for( unsigned int threadID = 0; threadID<this->m_NumberOfThreads; threadID++ ) {
523 for( unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits; threadID++ ) {
525 memset( m_ThreaderDerivativeF[threadID].data_block(),
527 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
529 memset( m_ThreaderDerivativeM[threadID].data_block(),
531 this->m_NumberOfParameters * sizeof(typename DerivativeType::ValueType) );
534 // MUST BE CALLED TO INITIATE PROCESSING
535 this->GetValueAndDerivativeMultiThreadedInitiate();
537 // Accumulate over the threads
538 DerivativeType derivativeF(this->m_NumberOfParameters), derivativeM(this->m_NumberOfParameters);
539 #if ITK_VERSION_MAJOR <= 4
540 for(unsigned int t=0; t<this->m_NumberOfThreads; t++) {
542 for(unsigned int t=0; t<this->m_NumberOfWorkUnits; t++) {
544 for(unsigned int parameter = 0; parameter < this->m_NumberOfParameters; parameter++) {
545 derivativeF[parameter] += m_ThreaderDerivativeF[t][parameter];
546 derivativeM[parameter] += m_ThreaderDerivativeM[t][parameter];
550 //Compute derivatives
551 if( this->m_NumberOfPixelsCounted > 0 && m_Denom != 0.0) {
552 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
553 derivative[i] = ( derivativeF[i] - (m_SFM/m_SMM)* derivativeM[i] ) / m_Denom;
556 for(unsigned int i=0; i<this->m_NumberOfParameters; i++) {
557 derivative[i] = itk::NumericTraits< MeasureType >::Zero;
565 * Get the match measure derivative
567 template < class TFixedImage, class TMovingImage >
569 NormalizedCorrelationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
570 ::GetDerivative( const ParametersType & parameters,
571 DerivativeType & derivative ) const
573 if( !this->m_FixedImage ) {
574 itkExceptionMacro( << "Fixed image has not been assigned" );
578 // call the combined version
579 this->GetValueAndDerivative( parameters, value, derivative );
582 } // end namespace clitk