]> Creatis software - clitk.git/blob - registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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: itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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 __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
37
38 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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 #include "itkStatisticsImageFilter.h"
46
47 #include "vnl/vnl_vector.txx"
48 #include "vnl/vnl_c_vector.txx"
49
50 namespace itk
51 {
52
53 /**
54  * Constructor
55  */
56 template < class TFixedImage, class TMovingImage >
57 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
58 ::MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
59 {
60   m_NumberOfHistogramBins = 50;
61
62   this->SetComputeGradient(false); // don't use the default gradient for now
63
64   // Initialize PDFs to NULL
65   m_JointPDF = NULL;
66   m_JointPDFDerivatives = NULL;
67
68   // Initialize memory
69   m_FixedImageMarginalPDF = NULL;
70   m_MovingImageMarginalPDF = NULL;
71
72   m_MovingImageNormalizedMin = 0.0;
73   m_FixedImageNormalizedMin = 0.0;
74   m_MovingImageTrueMin = 0.0;
75   m_MovingImageTrueMax = 0.0;
76   m_FixedImageBinSize = 0.0;
77   m_MovingImageBinSize = 0.0;
78
79   m_CubicBSplineDerivativeKernel = NULL;
80
81   // For multi-threading the metric
82   m_ThreaderFixedImageMarginalPDF = NULL;
83   m_ThreaderJointPDF = NULL;
84   m_ThreaderJointPDFDerivatives = NULL;
85   m_ThreaderJointPDFStartBin = NULL;
86   m_ThreaderJointPDFEndBin = NULL;
87   m_ThreaderJointPDFSum = NULL;
88   this->m_WithinThreadPreProcess = true;
89   this->m_WithinThreadPostProcess = false;
90
91   this->m_ThreaderMetricDerivative = NULL;
92   this->m_UseExplicitPDFDerivatives = true;
93   this->m_ImplicitDerivativesSecondPass = false;
94 }
95
96 template < class TFixedImage, class TMovingImage >
97 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
98 ::~MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
99 {
100
101   if(m_FixedImageMarginalPDF != NULL) {
102     delete [] m_FixedImageMarginalPDF;
103   }
104   m_FixedImageMarginalPDF = NULL;
105   if(m_MovingImageMarginalPDF != NULL) {
106     delete [] m_MovingImageMarginalPDF;
107   }
108   m_MovingImageMarginalPDF = NULL;
109   if(m_ThreaderJointPDF != NULL) {
110     delete [] m_ThreaderJointPDF;
111   }
112   m_ThreaderJointPDF = NULL;
113
114   if(m_ThreaderJointPDFDerivatives != NULL) {
115     delete [] m_ThreaderJointPDFDerivatives;
116   }
117   m_ThreaderJointPDFDerivatives = NULL;
118
119   if(m_ThreaderFixedImageMarginalPDF != NULL) {
120     delete [] m_ThreaderFixedImageMarginalPDF;
121   }
122   m_ThreaderFixedImageMarginalPDF = NULL;
123
124   if(m_ThreaderJointPDFStartBin != NULL) {
125     delete [] m_ThreaderJointPDFStartBin;
126   }
127   m_ThreaderJointPDFStartBin = NULL;
128
129   if(m_ThreaderJointPDFEndBin != NULL) {
130     delete [] m_ThreaderJointPDFEndBin;
131   }
132   m_ThreaderJointPDFEndBin = NULL;
133
134   if(m_ThreaderJointPDFSum != NULL) {
135     delete [] m_ThreaderJointPDFSum;
136   }
137   m_ThreaderJointPDFSum = NULL;
138
139   if( this->m_ThreaderMetricDerivative != NULL ) {
140     delete [] this->m_ThreaderMetricDerivative;
141   }
142   this->m_ThreaderMetricDerivative = NULL;
143 }
144
145 /**
146  * Print out internal information about this class
147  */
148 template < class TFixedImage, class TMovingImage  >
149 void
150 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
151 ::PrintSelf(std::ostream& os, Indent indent) const
152 {
153
154   Superclass::PrintSelf(os, indent);
155
156   os << indent << "NumberOfHistogramBins: ";
157   os << this->m_NumberOfHistogramBins << std::endl;
158
159   // Debugging information
160   os << indent << "FixedImageNormalizedMin: ";
161   os << this->m_FixedImageNormalizedMin << std::endl;
162   os << indent << "MovingImageNormalizedMin: ";
163   os << this->m_MovingImageNormalizedMin << std::endl;
164   os << indent << "MovingImageTrueMin: ";
165   os << this->m_MovingImageTrueMin << std::endl;
166   os << indent << "MovingImageTrueMax: ";
167   os << this->m_MovingImageTrueMax << std::endl;
168   os << indent << "FixedImageBinSize: ";
169   os << this->m_FixedImageBinSize << std::endl;
170   os << indent << "MovingImageBinSize: ";
171   os << this->m_MovingImageBinSize << std::endl;
172   os << indent << "UseExplicitPDFDerivatives: ";
173   os << this->m_UseExplicitPDFDerivatives << std::endl;
174   os << indent << "ImplicitDerivativesSecondPass: ";
175   os << this->m_ImplicitDerivativesSecondPass << std::endl;
176
177 }
178
179
180 /**
181  * Initialize
182  */
183 template <class TFixedImage, class TMovingImage>
184 void
185 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
186 ::Initialize(void) throw ( ExceptionObject )
187 {
188   this->Superclass::Initialize();
189   this->Superclass::MultiThreadingInitialize();
190
191   typedef StatisticsImageFilter<FixedImageType> FixedImageStatisticsFilterType;
192   typename FixedImageStatisticsFilterType::Pointer fixedImageStats =
193     FixedImageStatisticsFilterType::New();
194   fixedImageStats->SetInput( this->m_FixedImage );
195   fixedImageStats->SetNumberOfThreads( this->m_NumberOfThreads );
196   fixedImageStats->Update();
197
198   m_FixedImageTrueMin = fixedImageStats->GetMinimum();
199   m_FixedImageTrueMax = fixedImageStats->GetMaximum();
200   double fixedImageMin = m_FixedImageTrueMin;
201   double fixedImageMax = m_FixedImageTrueMax;
202
203   typedef StatisticsImageFilter<MovingImageType>
204   MovingImageStatisticsFilterType;
205   typename MovingImageStatisticsFilterType::Pointer movingImageStats =
206     MovingImageStatisticsFilterType::New();
207   movingImageStats->SetInput( this->m_MovingImage );
208   movingImageStats->SetNumberOfThreads( this->m_NumberOfThreads );
209   movingImageStats->Update();
210
211   m_MovingImageTrueMin = movingImageStats->GetMinimum();
212   m_MovingImageTrueMax = movingImageStats->GetMaximum();
213   double movingImageMin = m_MovingImageTrueMin;
214   double movingImageMax = m_MovingImageTrueMax;
215
216   itkDebugMacro( " FixedImageMin: " << fixedImageMin <<
217                  " FixedImageMax: " << fixedImageMax << std::endl );
218   itkDebugMacro( " MovingImageMin: " << movingImageMin <<
219                  " MovingImageMax: " << movingImageMax << std::endl );
220
221
222   /**
223    * Compute binsize for the histograms.
224    *
225    * The binsize for the image intensities needs to be adjusted so that
226    * we can avoid dealing with boundary conditions using the cubic
227    * spline as the Parzen window.  We do this by increasing the size
228    * of the bins so that the joint histogram becomes "padded" at the
229    * borders. Because we are changing the binsize,
230    * we also need to shift the minimum by the padded amount in order to
231    * avoid minimum values filling in our padded region.
232    *
233    * Note that there can still be non-zero bin values in the padded region,
234    * it's just that these bins will never be a central bin for the Parzen
235    * window.
236    *
237    */
238   const int padding = 2;  // this will pad by 2 bins
239
240   m_FixedImageBinSize = ( fixedImageMax - fixedImageMin )
241                         / static_cast<double>( m_NumberOfHistogramBins
242                             - 2 * padding );
243   m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize
244                               - static_cast<double>( padding );
245
246   m_MovingImageBinSize = ( movingImageMax - movingImageMin )
247                          / static_cast<double>( m_NumberOfHistogramBins
248                              - 2 * padding );
249   m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize
250                                - static_cast<double>( padding );
251
252   itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin );
253   itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin );
254   itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize );
255   itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize );
256
257
258   /**
259    * Allocate memory for the marginal PDF and initialize values
260    * to zero. The marginal PDFs are stored as std::vector.
261    */
262   if(m_FixedImageMarginalPDF != NULL) {
263     delete [] m_FixedImageMarginalPDF;
264   }
265   m_FixedImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
266   if(m_MovingImageMarginalPDF != NULL) {
267     delete [] m_MovingImageMarginalPDF;
268   }
269   m_MovingImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
270
271   /**
272    * Allocate memory for the joint PDF and joint PDF derivatives.
273    * The joint PDF and joint PDF derivatives are store as itk::Image.
274    */
275   m_JointPDF = JointPDFType::New();
276   m_JointPDFDerivatives = JointPDFDerivativesType::New();
277
278   // Instantiate a region, index, size
279   JointPDFRegionType             jointPDFRegion;
280   JointPDFIndexType              jointPDFIndex;
281   JointPDFSizeType               jointPDFSize;
282
283   // Deallocate the memory that may have been allocated for
284   // previous runs of the metric.
285   this->m_JointPDFDerivatives = NULL;  // by destroying the dynamic array
286   this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones
287   this->m_MetricDerivative = DerivativeType( 1 );
288
289   JointPDFDerivativesRegionType  jointPDFDerivativesRegion;
290
291   //
292   // Now allocate memory according to the user-selected method.
293   //
294   if( this->m_UseExplicitPDFDerivatives ) {
295     this->m_JointPDFDerivatives = JointPDFDerivativesType::New();
296
297     JointPDFDerivativesIndexType  jointPDFDerivativesIndex;
298     JointPDFDerivativesSizeType    jointPDFDerivativesSize;
299
300     // For the derivatives of the joint PDF define a region starting from {0,0,0}
301     // with size {m_NumberOfParameters,m_NumberOfHistogramBins,
302     // m_NumberOfHistogramBins}. The dimension represents transform parameters,
303     // fixed image parzen window index and moving image parzen window index,
304     // respectively.
305     jointPDFDerivativesIndex.Fill( 0 );
306     jointPDFDerivativesSize[0] = this->m_NumberOfParameters;
307     jointPDFDerivativesSize[1] = this->m_NumberOfHistogramBins;
308     jointPDFDerivativesSize[2] = this->m_NumberOfHistogramBins;
309
310     jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex );
311     jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize );
312
313     // Set the regions and allocate
314     m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion );
315     m_JointPDFDerivatives->Allocate();
316     m_JointPDFDerivativesBufferSize = jointPDFDerivativesSize[0]
317                                       * jointPDFDerivativesSize[1]
318                                       * jointPDFDerivativesSize[2]
319                                       * sizeof(JointPDFDerivativesValueType);
320
321   } else {
322     /** Allocate memory for helper array that will contain the pRatios
323      *  for each bin of the joint histogram. This is part of the effort
324      *  for flattening the computation of the PDF Jacobians.
325      */
326     this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins );
327     this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() );
328   }
329
330   // For the joint PDF define a region starting from {0,0}
331   // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}.
332   // The dimension represents fixed image parzen window index
333   // and moving image parzen window index, respectively.
334   jointPDFIndex.Fill( 0 );
335   jointPDFSize.Fill( m_NumberOfHistogramBins );
336
337   jointPDFRegion.SetIndex( jointPDFIndex );
338   jointPDFRegion.SetSize( jointPDFSize );
339
340   // Set the regions and allocate
341   m_JointPDF->SetRegions( jointPDFRegion );
342   m_JointPDF->Allocate();
343
344   m_JointPDFBufferSize = jointPDFSize[0] * jointPDFSize[1] * sizeof(PDFValueType);
345
346
347   /**
348    * Setup the kernels used for the Parzen windows.
349    */
350   m_CubicBSplineKernel = CubicBSplineFunctionType::New();
351   m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New();
352
353   /**
354    * Pre-compute the fixed image parzen window index for
355    * each point of the fixed image sample points list.
356    */
357   this->ComputeFixedImageParzenWindowIndices( this->m_FixedImageSamples );
358
359
360   if(m_ThreaderFixedImageMarginalPDF != NULL) {
361     delete [] m_ThreaderFixedImageMarginalPDF;
362   }
363   // Assumes number of threads doesn't change between calls to Initialize
364   m_ThreaderFixedImageMarginalPDF = new
365   PDFValueType[(this->m_NumberOfThreads-1)
366                * m_NumberOfHistogramBins];
367
368   if(m_ThreaderJointPDF != NULL) {
369     delete [] m_ThreaderJointPDF;
370   }
371   m_ThreaderJointPDF = new typename
372   JointPDFType::Pointer[this->m_NumberOfThreads-1];
373
374   if(m_ThreaderJointPDFStartBin != NULL) {
375     delete [] m_ThreaderJointPDFStartBin;
376   }
377   m_ThreaderJointPDFStartBin = new int[this->m_NumberOfThreads];
378
379   if(m_ThreaderJointPDFEndBin != NULL) {
380     delete [] m_ThreaderJointPDFEndBin;
381   }
382   m_ThreaderJointPDFEndBin = new int[this->m_NumberOfThreads];
383
384   if(m_ThreaderJointPDFSum != NULL) {
385     delete [] m_ThreaderJointPDFSum;
386   }
387   m_ThreaderJointPDFSum = new double[this->m_NumberOfThreads];
388
389   unsigned int threadID;
390
391   int binRange = m_NumberOfHistogramBins / this->m_NumberOfThreads;
392
393   for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
394     m_ThreaderJointPDF[threadID] = JointPDFType::New();
395     m_ThreaderJointPDF[threadID]->SetRegions( jointPDFRegion );
396     m_ThreaderJointPDF[threadID]->Allocate();
397
398     m_ThreaderJointPDFStartBin[threadID] = threadID * binRange;
399     m_ThreaderJointPDFEndBin[threadID] = (threadID + 1) * binRange - 1;
400   }
401
402   m_ThreaderJointPDFStartBin[this->m_NumberOfThreads-1] =
403     (this->m_NumberOfThreads - 1 ) * binRange;
404
405   m_ThreaderJointPDFEndBin[this->m_NumberOfThreads-1] = m_NumberOfHistogramBins - 1;
406
407   // Release memory of arrays that may have been used for
408   // previous executions of this metric with different settings
409   // of the memory caching flags.
410   if(m_ThreaderJointPDFDerivatives != NULL) {
411     delete [] m_ThreaderJointPDFDerivatives;
412   }
413   m_ThreaderJointPDFDerivatives = NULL;
414
415   if(m_ThreaderMetricDerivative != NULL) {
416     delete [] m_ThreaderMetricDerivative;
417   }
418   m_ThreaderMetricDerivative = NULL;
419
420
421   if( this->m_UseExplicitPDFDerivatives ) {
422     m_ThreaderJointPDFDerivatives = new typename
423     JointPDFDerivativesType::Pointer[this->m_NumberOfThreads-1];
424
425     for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
426       m_ThreaderJointPDFDerivatives[threadID] = JointPDFDerivativesType::New();
427       m_ThreaderJointPDFDerivatives[threadID]->SetRegions(
428         jointPDFDerivativesRegion );
429       m_ThreaderJointPDFDerivatives[threadID]->Allocate();
430     }
431   } else {
432     m_ThreaderMetricDerivative = new DerivativeType[this->m_NumberOfThreads-1];
433
434     for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
435       this->m_ThreaderMetricDerivative[threadID] = DerivativeType( this->GetNumberOfParameters() );
436     }
437   }
438 }
439
440 /**
441  * Uniformly sample the fixed image domain using a random walk
442  */
443 template < class TFixedImage, class TMovingImage >
444 void
445 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
446 ::ComputeFixedImageParzenWindowIndices(
447   FixedImageSampleContainer& samples )
448 {
449
450   typename FixedImageSampleContainer::iterator iter;
451   typename FixedImageSampleContainer::const_iterator end=samples.end();
452
453   for( iter=samples.begin(); iter != end; ++iter ) {
454
455     // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
456     double windowTerm = static_cast<double>( (*iter).value )
457                         / m_FixedImageBinSize
458                         - m_FixedImageNormalizedMin;
459     unsigned int pindex = static_cast<unsigned int>( windowTerm );
460
461     // Make sure the extreme values are in valid bins
462     if ( pindex < 2 ) {
463       pindex = 2;
464     } else if ( pindex > (m_NumberOfHistogramBins - 3) ) {
465       pindex = m_NumberOfHistogramBins - 3;
466     }
467
468     (*iter).valueIndex = pindex;
469   }
470
471 }
472
473 template < class TFixedImage, class TMovingImage  >
474 inline void
475 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
476 ::GetValueThreadPreProcess( unsigned int threadID,
477                             bool withinSampleThread ) const
478 {
479
480   this->Superclass::GetValueThreadPreProcess( threadID, withinSampleThread );
481
482   if(threadID > 0) {
483     memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
484             0,
485             m_JointPDFBufferSize );
486     memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
487               *m_NumberOfHistogramBins]),
488             0,
489             m_NumberOfHistogramBins*sizeof(PDFValueType) );
490   } else {
491     // zero-th thread uses the variables directly
492     memset( m_JointPDF->GetBufferPointer(),
493             0,
494             m_JointPDFBufferSize );
495     memset( m_FixedImageMarginalPDF,
496             0,
497             m_NumberOfHistogramBins*sizeof(PDFValueType) );
498   }
499 }
500
501 template < class TFixedImage, class TMovingImage  >
502 inline bool
503 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
504 ::GetValueThreadProcessSample( unsigned int threadID,
505                                unsigned long fixedImageSample,
506                                const MovingImagePointType & itkNotUsed(mappedPoint),
507                                double movingImageValue) const
508 {
509   /**
510    * Compute this sample's contribution to the marginal and
511    * joint distributions.
512    *
513    */
514
515   if(movingImageValue < m_MovingImageTrueMin) {
516     return false;
517   } else if(movingImageValue > m_MovingImageTrueMax) {
518     return false;
519   }
520
521   // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
522   double movingImageParzenWindowTerm = movingImageValue
523                                        / m_MovingImageBinSize
524                                        - m_MovingImageNormalizedMin;
525   // Same as floor
526   unsigned int movingImageParzenWindowIndex =
527     static_cast<unsigned int>( movingImageParzenWindowTerm );
528   if( movingImageParzenWindowIndex < 2 ) {
529     movingImageParzenWindowIndex = 2;
530   } else if( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
531     movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
532   }
533
534   unsigned int fixedImageParzenWindowIndex =
535     this->m_FixedImageSamples[fixedImageSample].valueIndex;
536   if(threadID > 0) {
537     m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
538                                     + fixedImageParzenWindowIndex] += 1;
539   } else {
540     m_FixedImageMarginalPDF[fixedImageParzenWindowIndex] += 1;
541   }
542
543   // Pointer to affected bin to be updated
544   JointPDFValueType *pdfPtr;
545   if(threadID > 0) {
546     pdfPtr = m_ThreaderJointPDF[threadID-1]->GetBufferPointer() +
547              ( fixedImageParzenWindowIndex
548                * m_ThreaderJointPDF[threadID-1]
549                ->GetOffsetTable()[1] );
550   } else {
551     pdfPtr = m_JointPDF->GetBufferPointer() +
552              ( fixedImageParzenWindowIndex
553                * m_JointPDF->GetOffsetTable()[1] );
554   }
555
556   // Move the pointer to the first affected bin
557   int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
558   pdfPtr += pdfMovingIndex;
559   int pdfMovingIndexMax = static_cast<int>(movingImageParzenWindowIndex) + 2;
560
561   double movingImageParzenWindowArg =
562     static_cast<double>( pdfMovingIndex )
563     - movingImageParzenWindowTerm;
564
565   while( pdfMovingIndex <= pdfMovingIndexMax ) {
566     *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
567                    ->Evaluate(
568                      movingImageParzenWindowArg ) );
569     movingImageParzenWindowArg += 1;
570     ++pdfMovingIndex;
571   }
572
573   return true;
574 }
575
576 template < class TFixedImage, class TMovingImage  >
577 inline void
578 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
579 ::GetValueThreadPostProcess( unsigned int threadID,
580                              bool itkNotUsed(withinSampleThread) ) const
581 {
582   unsigned int t;
583   int i;
584   int maxI;
585   maxI = m_NumberOfHistogramBins
586          * ( m_ThreaderJointPDFEndBin[threadID]
587              - m_ThreaderJointPDFStartBin[threadID] + 1);
588   JointPDFValueType *pdfPtr;
589   JointPDFValueType *pdfPtrStart;
590   pdfPtrStart = m_JointPDF->GetBufferPointer()
591                 + ( m_ThreaderJointPDFStartBin[threadID]
592                     * m_JointPDF->GetOffsetTable()[1] );
593   JointPDFValueType *tPdfPtr;
594   JointPDFValueType *tPdfPtrEnd;
595   unsigned int       tPdfPtrOffset;
596   tPdfPtrOffset = ( m_ThreaderJointPDFStartBin[threadID]
597                     * m_JointPDF->GetOffsetTable()[1] );
598   for(t=0; t<this->m_NumberOfThreads-1; t++) {
599     pdfPtr = pdfPtrStart;
600     tPdfPtr = m_ThreaderJointPDF[t]->GetBufferPointer() + tPdfPtrOffset;
601     tPdfPtrEnd = tPdfPtr + maxI;
602     //for(i=0; i < maxI; i++)
603     while(tPdfPtr < tPdfPtrEnd) {
604       *(pdfPtr++) += *(tPdfPtr++);
605     }
606     for(i = m_ThreaderJointPDFStartBin[threadID];
607         i <= m_ThreaderJointPDFEndBin[threadID];
608         i++) {
609       m_FixedImageMarginalPDF[i] += m_ThreaderFixedImageMarginalPDF[
610                                       (t*m_NumberOfHistogramBins) + i];
611     }
612   }
613   double jointPDFSum = 0.0;
614   pdfPtr = pdfPtrStart;
615   for(i = 0; i < maxI; i++) {
616     jointPDFSum += *(pdfPtr++);
617   }
618   if(threadID > 0) {
619     m_ThreaderJointPDFSum[threadID-1] = jointPDFSum;
620   } else {
621     m_JointPDFSum = jointPDFSum;
622   }
623 }
624
625 template < class TFixedImage, class TMovingImage  >
626 typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
627 ::MeasureType
628 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
629 ::GetValue( const ParametersType & parameters ) const
630 {
631   // Set up the parameters in the transform
632   this->m_Transform->SetParameters( parameters );
633 #if ITK_VERSION_MAJOR < 4
634   this->m_Parameters = parameters;
635 #endif
636
637   // MUST BE CALLED TO INITIATE PROCESSING
638   this->GetValueMultiThreadedInitiate();
639
640   // MUST BE CALLED TO INITIATE PROCESSING
641   this->GetValueMultiThreadedPostProcessInitiate();
642
643   for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
644     m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
645   }
646   if ( m_JointPDFSum == 0.0 ) {
647     itkExceptionMacro( "Joint PDF summed to zero" );
648   }
649
650   memset( m_MovingImageMarginalPDF,
651           0,
652           m_NumberOfHistogramBins*sizeof(PDFValueType) );
653
654   JointPDFValueType * pdfPtr;
655   PDFValueType * movingMarginalPtr;
656   unsigned int i, j;
657   double fixedPDFSum = 0.0;
658   double nFactor = 1.0 / m_JointPDFSum;
659   pdfPtr = m_JointPDF->GetBufferPointer();
660   for(i=0; i<m_NumberOfHistogramBins; i++) {
661     fixedPDFSum += m_FixedImageMarginalPDF[i];
662     movingMarginalPtr = m_MovingImageMarginalPDF;
663     for(j=0; j<m_NumberOfHistogramBins; j++) {
664       *(pdfPtr) *= nFactor;
665       *(movingMarginalPtr++) += *(pdfPtr++);
666     }
667   }
668
669   if( this->m_NumberOfPixelsCounted <
670       this->m_NumberOfFixedImageSamples / 16 ) {
671     itkExceptionMacro( "Too many samples map outside moving image buffer: "
672                        << this->m_NumberOfPixelsCounted << " / "
673                        << this->m_NumberOfFixedImageSamples
674                        << std::endl );
675   }
676
677   // Normalize the fixed image marginal PDF
678   if ( fixedPDFSum == 0.0 ) {
679     itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
680   }
681   for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
682     m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
683   }
684
685   /**
686    * Compute the metric by double summation over histogram.
687    */
688
689   // Setup pointer to point to the first bin
690   JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
691
692   double sum = 0.0;
693
694   for( unsigned int fixedIndex = 0;
695        fixedIndex < m_NumberOfHistogramBins;
696        ++fixedIndex ) {
697     double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
698
699     for( unsigned int movingIndex = 0;
700          movingIndex < m_NumberOfHistogramBins;
701          ++movingIndex, jointPDFPtr++ ) {
702       double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
703       double jointPDFValue = *(jointPDFPtr);
704
705       // check for non-zero bin contribution
706       if( jointPDFValue > 1e-16 &&  movingImagePDFValue > 1e-16 ) {
707
708         double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
709         if( fixedImagePDFValue > 1e-16) {
710           sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
711         }
712
713       }  // end if-block to check non-zero bin contribution
714     }  // end for-loop over moving index
715   }  // end for-loop over fixed index
716
717   return static_cast<MeasureType>( -1.0 * sum );
718
719 }
720
721
722 template < class TFixedImage, class TMovingImage  >
723 inline void
724 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
725 ::GetValueAndDerivativeThreadPreProcess( unsigned int threadID,
726     bool itkNotUsed(withinSampleThread) ) const
727 {
728   if(threadID > 0) {
729     memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
730             0,
731             m_JointPDFBufferSize );
732
733     memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
734               * m_NumberOfHistogramBins]),
735             0,
736             m_NumberOfHistogramBins*sizeof(PDFValueType) );
737
738     if( this->m_UseExplicitPDFDerivatives ) {
739       memset( m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer(),
740               0,
741               m_JointPDFDerivativesBufferSize );
742     }
743   } else {
744     memset( m_JointPDF->GetBufferPointer(),
745             0,
746             m_JointPDFBufferSize );
747     memset( m_FixedImageMarginalPDF,
748             0,
749             m_NumberOfHistogramBins*sizeof(PDFValueType) );
750
751     if( this->m_UseExplicitPDFDerivatives ) {
752       memset( m_JointPDFDerivatives->GetBufferPointer(),
753               0,
754               m_JointPDFDerivativesBufferSize );
755     }
756   }
757 }
758
759 template < class TFixedImage, class TMovingImage  >
760 inline bool
761 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
762 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
763     unsigned long fixedImageSample,
764     const MovingImagePointType & itkNotUsed(mappedPoint),
765     double movingImageValue,
766     const ImageDerivativesType &
767     movingImageGradientValue) const
768 {
769   /**
770    * Compute this sample's contribution to the marginal
771    *   and joint distributions.
772    *
773    */
774   if(movingImageValue < m_MovingImageTrueMin) {
775     return false;
776   } else if(movingImageValue > m_MovingImageTrueMax) {
777     return false;
778   }
779
780   unsigned int fixedImageParzenWindowIndex =
781     this->m_FixedImageSamples[fixedImageSample].valueIndex;
782
783   // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
784   double movingImageParzenWindowTerm = movingImageValue
785                                        / m_MovingImageBinSize
786                                        - m_MovingImageNormalizedMin;
787   unsigned int movingImageParzenWindowIndex =
788     static_cast<unsigned int>( movingImageParzenWindowTerm );
789
790   // Make sure the extreme values are in valid bins
791   if ( movingImageParzenWindowIndex < 2 ) {
792     movingImageParzenWindowIndex = 2;
793   } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
794     movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
795   }
796
797   // Since a zero-order BSpline (box car) kernel is used for
798   // the fixed image marginal pdf, we need only increment the
799   // fixedImageParzenWindowIndex by value of 1.0.
800   if(threadID > 0) {
801     ++m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
802                                       + fixedImageParzenWindowIndex];
803   } else {
804     ++m_FixedImageMarginalPDF[fixedImageParzenWindowIndex];
805   }
806
807   /**
808    * The region of support of the parzen window determines which bins
809    * of the joint PDF are effected by the pair of image values.
810    * Since we are using a cubic spline for the moving image parzen
811    * window, four bins are effected.  The fixed image parzen window is
812    * a zero-order spline (box car) and thus effects only one bin.
813    *
814    *  The PDF is arranged so that moving image bins corresponds to the
815    * zero-th (column) dimension and the fixed image bins corresponds
816    * to the first (row) dimension.
817    *
818    */
819
820   // Pointer to affected bin to be updated
821   JointPDFValueType *pdfPtr;
822   if(threadID > 0) {
823     pdfPtr = m_ThreaderJointPDF[threadID-1]
824              ->GetBufferPointer() +
825              ( fixedImageParzenWindowIndex
826                * m_NumberOfHistogramBins );
827   } else {
828     pdfPtr = m_JointPDF->GetBufferPointer() +
829              ( fixedImageParzenWindowIndex
830                * m_NumberOfHistogramBins );
831   }
832
833   // Move the pointer to the fist affected bin
834   int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
835   pdfPtr += pdfMovingIndex;
836   int pdfMovingIndexMax = static_cast<int>(movingImageParzenWindowIndex) + 2;
837
838   double movingImageParzenWindowArg = static_cast<double>( pdfMovingIndex )
839                                       - static_cast<double>( movingImageParzenWindowTerm );
840
841   while( pdfMovingIndex <= pdfMovingIndexMax ) {
842     *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
843                    ->Evaluate(
844                      movingImageParzenWindowArg ) );
845
846     if( this->m_UseExplicitPDFDerivatives || this->m_ImplicitDerivativesSecondPass ) {
847       // Compute the cubicBSplineDerivative for later repeated use.
848       double cubicBSplineDerivativeValue =
849         m_CubicBSplineDerivativeKernel->Evaluate( movingImageParzenWindowArg );
850
851       // Compute PDF derivative contribution.
852       this->ComputePDFDerivatives( threadID,
853                                    fixedImageSample,
854                                    pdfMovingIndex,
855                                    movingImageGradientValue,
856                                    cubicBSplineDerivativeValue );
857     }
858
859     movingImageParzenWindowArg += 1;
860     ++pdfMovingIndex;
861   }
862
863   return true;
864 }
865
866 template < class TFixedImage, class TMovingImage  >
867 inline void
868 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
869 ::GetValueAndDerivativeThreadPostProcess( unsigned int threadID,
870     bool withinSampleThread ) const
871 {
872   this->GetValueThreadPostProcess( threadID, withinSampleThread );
873
874   if( this->m_UseExplicitPDFDerivatives ) {
875     const unsigned int rowSize = this->m_NumberOfParameters * m_NumberOfHistogramBins;
876
877     const unsigned int maxI =
878       rowSize * ( m_ThreaderJointPDFEndBin[threadID]
879                   - m_ThreaderJointPDFStartBin[threadID] + 1 );
880
881     JointPDFDerivativesValueType *pdfDPtr;
882     JointPDFDerivativesValueType *pdfDPtrStart;
883     pdfDPtrStart = m_JointPDFDerivatives->GetBufferPointer()
884                    + ( m_ThreaderJointPDFStartBin[threadID] * rowSize );
885     JointPDFDerivativesValueType *tPdfDPtr;
886     JointPDFDerivativesValueType *tPdfDPtrEnd;
887     unsigned int tPdfDPtrOffset;
888     tPdfDPtrOffset = m_ThreaderJointPDFStartBin[threadID] *  rowSize;
889     for(unsigned int t=0; t<this->m_NumberOfThreads-1; t++) {
890       pdfDPtr = pdfDPtrStart;
891       tPdfDPtr = m_ThreaderJointPDFDerivatives[t]->GetBufferPointer()
892                  + tPdfDPtrOffset;
893       tPdfDPtrEnd = tPdfDPtr + maxI;
894       // for(i = 0; i < maxI; i++)
895       while(tPdfDPtr < tPdfDPtrEnd) {
896         *(pdfDPtr++) += *(tPdfDPtr++);
897       }
898     }
899
900     double nFactor = 1.0 / (m_MovingImageBinSize
901                             * this->m_NumberOfPixelsCounted);
902
903     pdfDPtr = pdfDPtrStart;
904     tPdfDPtrEnd = pdfDPtrStart + maxI;
905     //for(int i = 0; i < maxI; i++)
906     while(pdfDPtr < tPdfDPtrEnd) {
907       *(pdfDPtr++) *= nFactor;
908     }
909   }
910 }
911
912 /**
913  * Get the both Value and Derivative Measure
914  */
915 template < class TFixedImage, class TMovingImage  >
916 void
917 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
918 ::GetValueAndDerivative( const ParametersType & parameters,
919                          MeasureType & value,
920                          DerivativeType & derivative) const
921 {
922   // Set output values to zero
923   value = NumericTraits< MeasureType >::Zero;
924
925   if( this->m_UseExplicitPDFDerivatives ) {
926     // Set output values to zero
927     if(derivative.GetSize() != this->m_NumberOfParameters) {
928       derivative = DerivativeType( this->m_NumberOfParameters );
929     }
930     memset( derivative.data_block(),
931             0,
932             this->m_NumberOfParameters * sizeof(double) );
933   } else {
934     this->m_PRatioArray.Fill( 0.0 );
935     this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero );
936     for(unsigned int threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++ ) {
937       this->m_ThreaderMetricDerivative[threadID].Fill( NumericTraits< MeasureType >::Zero );
938     }
939     this->m_ImplicitDerivativesSecondPass = false;
940   }
941
942   // Set up the parameters in the transform
943   this->m_Transform->SetParameters( parameters );
944 #if ITK_VERSION_MAJOR < 4
945   this->m_Parameters = parameters;
946 #endif
947
948   // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
949   this->GetValueAndDerivativeMultiThreadedInitiate();
950
951   // CALL IF DOING THREADED POST PROCESSING
952   this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
953
954   for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
955     m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
956   }
957   if ( m_JointPDFSum == 0.0 ) {
958     itkExceptionMacro( "Joint PDF summed to zero" );
959   }
960
961   memset( m_MovingImageMarginalPDF,
962           0,
963           m_NumberOfHistogramBins*sizeof(PDFValueType) );
964
965   JointPDFValueType * pdfPtr;
966   PDFValueType * movingMarginalPtr;
967   unsigned int i, j;
968   double fixedPDFSum = 0.0;
969   const double normalizationFactor = 1.0 / m_JointPDFSum;
970
971   pdfPtr = m_JointPDF->GetBufferPointer();
972   for(i=0; i<m_NumberOfHistogramBins; i++) {
973     fixedPDFSum += m_FixedImageMarginalPDF[i];
974     movingMarginalPtr = m_MovingImageMarginalPDF;
975     for(j=0; j<m_NumberOfHistogramBins; j++) {
976       *(pdfPtr) *= normalizationFactor;
977       *(movingMarginalPtr++) += *(pdfPtr++);
978     }
979   }
980
981   if( this->m_NumberOfPixelsCounted <
982       this->m_NumberOfFixedImageSamples / 16 ) {
983     itkExceptionMacro( "Too many samples map outside moving image buffer: "
984                        << this->m_NumberOfPixelsCounted << " / "
985                        << this->m_NumberOfFixedImageSamples
986                        << std::endl );
987   }
988
989   // Normalize the fixed image marginal PDF
990   if ( fixedPDFSum == 0.0 ) {
991     itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
992   }
993   for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
994     m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
995   }
996
997   /**
998    * Compute the metric by double summation over histogram.
999    */
1000
1001   // Setup pointer to point to the first bin
1002   JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
1003
1004   // Initialize sum to zero
1005   double sum = 0.0;
1006
1007   const double nFactor = 1.0 / (m_MovingImageBinSize
1008                                 * this->m_NumberOfPixelsCounted);
1009
1010   for( unsigned int fixedIndex = 0;
1011        fixedIndex < m_NumberOfHistogramBins;
1012        ++fixedIndex ) {
1013     double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
1014
1015     for( unsigned int movingIndex = 0;
1016          movingIndex < m_NumberOfHistogramBins;
1017          ++movingIndex, jointPDFPtr++ ) {
1018       double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
1019       double jointPDFValue = *(jointPDFPtr);
1020
1021       // check for non-zero bin contribution
1022       if( jointPDFValue > 1e-16 &&  movingImagePDFValue > 1e-16 ) {
1023
1024         double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
1025
1026         if( fixedImagePDFValue > 1e-16) {
1027           sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
1028         }
1029
1030         if( this->m_UseExplicitPDFDerivatives ) {
1031           // move joint pdf derivative pointer to the right position
1032           JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1033                                          + ( fixedIndex  * m_JointPDFDerivatives->GetOffsetTable()[2] )
1034                                          + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1035
1036           for( unsigned int parameter=0; parameter < this->m_NumberOfParameters; ++parameter, derivPtr++ ) {
1037
1038             // Ref: eqn 23 of Thevenaz & Unser paper [3]
1039             derivative[parameter] -= (*derivPtr) * pRatio;
1040
1041           }  // end for-loop over parameters
1042         } else {
1043           this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor;
1044         }
1045       }  // end if-block to check non-zero bin contribution
1046     }  // end for-loop over moving index
1047   }  // end for-loop over fixed index
1048
1049   if( !(this->m_UseExplicitPDFDerivatives ) ) {
1050     // Second pass: This one is done for accumulating the contributions
1051     //              to the derivative array.
1052     //
1053     this->m_ImplicitDerivativesSecondPass = true;
1054     //
1055     // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
1056     this->GetValueAndDerivativeMultiThreadedInitiate();
1057
1058     // CALL IF DOING THREADED POST PROCESSING
1059     this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
1060
1061     // Consolidate the contributions from each one of the threads to the total
1062     // derivative.
1063     for(unsigned int t = 0; t < this->m_NumberOfThreads-1; t++ ) {
1064       DerivativeType * source = &(this->m_ThreaderMetricDerivative[t]);
1065       for(unsigned int pp=0; pp < this->m_NumberOfParameters; pp++ ) {
1066         this->m_MetricDerivative[pp] += (*source)[pp];
1067       }
1068     }
1069
1070     derivative = this->m_MetricDerivative;
1071   }
1072
1073   value = static_cast<MeasureType>( -1.0 * sum );
1074
1075 }
1076
1077 /**
1078  * Get the match measure derivative
1079  */
1080 template < class TFixedImage, class TMovingImage  >
1081 void
1082 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1083 ::GetDerivative( const ParametersType & parameters,
1084                  DerivativeType & derivative ) const
1085 {
1086   MeasureType value;
1087   // call the combined version
1088   this->GetValueAndDerivative( parameters, value, derivative );
1089 }
1090
1091
1092 /**
1093  * Compute PDF derivatives contribution for each parameter
1094  */
1095 template < class TFixedImage, class TMovingImage >
1096 void
1097 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1098 ::ComputePDFDerivatives( unsigned int threadID,
1099                          unsigned int sampleNumber,
1100                          int pdfMovingIndex,
1101                          const ImageDerivativesType & movingImageGradientValue,
1102                          double cubicBSplineDerivativeValue ) const
1103 {
1104   // Update bins in the PDF derivatives for the current intensity pair
1105   // Could pre-compute
1106   JointPDFDerivativesValueType * derivPtr;
1107
1108   double precomputedWeight = 0.0;
1109
1110   const int pdfFixedIndex = this->m_FixedImageSamples[sampleNumber].valueIndex;
1111
1112   DerivativeType * derivativeHelperArray = NULL;
1113
1114   if( this->m_UseExplicitPDFDerivatives ) {
1115     if(threadID > 0) {
1116       derivPtr = m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer()
1117                  + ( pdfFixedIndex  * m_JointPDFDerivatives->GetOffsetTable()[2] )
1118                  + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1119     } else {
1120       derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1121                  + ( pdfFixedIndex  * m_JointPDFDerivatives->GetOffsetTable()[2] )
1122                  + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1123     }
1124   } else {
1125     derivPtr = 0;
1126     // Recover the precomputed weight for this specific PDF bin
1127     precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex];
1128     if(threadID > 0) {
1129       derivativeHelperArray = &(this->m_ThreaderMetricDerivative[threadID-1]);
1130     } else {
1131       derivativeHelperArray = &(this->m_MetricDerivative);
1132     }
1133
1134   }
1135
1136   if( !this->m_TransformIsBSpline ) {
1137     /**
1138      * Generic version which works for all transforms.
1139      */
1140
1141     // Compute the transform Jacobian.
1142     // Should pre-compute
1143     typedef typename TransformType::JacobianType JacobianType;
1144
1145     // Need to use one of the threader transforms if we're
1146     // not in thread 0.
1147     //
1148     // Use a raw pointer here to avoid the overhead of smart pointers.
1149     // For instance, Register and UnRegister have mutex locks around
1150     // the reference counts.
1151     TransformType* transform;
1152
1153     if (threadID > 0) {
1154       transform = this->m_ThreaderTransform[threadID - 1];
1155     } else {
1156       transform = this->m_Transform;
1157     }
1158
1159 #if ITK_VERSION_MAJOR >= 4
1160     JacobianType jacobian;
1161     transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[sampleNumber].point, jacobian);
1162 #else
1163     const JacobianType& jacobian =
1164       transform->GetJacobian( this->m_FixedImageSamples[sampleNumber].point );
1165 #endif
1166
1167     //     for ( unsigned int mu = 0; mu < this->m_NumberOfParameters; mu++ )
1168     //       {
1169     //       double innerProduct = 0.0;
1170     //       for ( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ )
1171     //         {
1172     //         innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1173     //         }
1174
1175     //       const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1176
1177     //       if( this->m_UseExplicitPDFDerivatives )
1178     //         {
1179     //         *(derivPtr) -= derivativeContribution;
1180     //         ++derivPtr;
1181     //         }
1182     //       else
1183     //         {
1184     //         (*derivativeHelperArray)[mu] += precomputedWeight * derivativeContribution;
1185     //         }
1186     //       }
1187     //    }
1188     // JV
1189     unsigned int mu, dim;
1190     double innerProduct,derivativeContribution;
1191     for ( mu = 0; mu < this->m_NumberOfParameters; mu+=3 ) {
1192       // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
1193       if (jacobian[0][mu])
1194         for ( dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1195           innerProduct = jacobian[dim][mu+dim] * movingImageGradientValue[dim];
1196           derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1197
1198           if( this->m_UseExplicitPDFDerivatives ) {
1199             *(derivPtr) -= derivativeContribution;
1200             ++derivPtr;
1201           } else {
1202             (*derivativeHelperArray)[mu+dim] += precomputedWeight * derivativeContribution;
1203           }
1204         }
1205       else  if( this->m_UseExplicitPDFDerivatives ) derivPtr+=3;
1206     }
1207   } else {
1208     const WeightsValueType * weights = NULL;
1209     const IndexValueType   * indices = NULL;
1210
1211
1212     BSplineTransformWeightsType    * weightsHelper = NULL;
1213     BSplineTransformIndexArrayType * indicesHelper = NULL;
1214
1215     if( this->m_UseCachingOfBSplineWeights ) {
1216       //
1217       // If the transform is of type BSplineDeformableTransform, we can obtain
1218       // a speed up by only processing the affected parameters.  Note that
1219       // these pointers are just pointing to pre-allocated rows of the caching
1220       // arrays. There is therefore, no need to free this memory.
1221       //
1222       weights = this->m_BSplineTransformWeightsArray[sampleNumber];
1223       indices = this->m_BSplineTransformIndicesArray[sampleNumber];
1224     } else {
1225       if( threadID > 0 ) {
1226         weightsHelper = &(this->m_ThreaderBSplineTransformWeights[threadID-1]);
1227         indicesHelper = &(this->m_ThreaderBSplineTransformIndices[threadID-1]);
1228       } else {
1229         weightsHelper = &(this->m_BSplineTransformWeights);
1230         indicesHelper = &(this->m_BSplineTransformIndices);
1231       }
1232
1233 #if ITK_VERSION_MAJOR >= 4
1234       this->m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition(
1235         this->m_FixedImageSamples[sampleNumber].point,
1236         *weightsHelper, *indicesHelper );
1237 #else
1238       this->m_BSplineTransform->GetJacobian(
1239         this->m_FixedImageSamples[sampleNumber].point,
1240         *weightsHelper, *indicesHelper );
1241 #endif
1242     }
1243
1244     for( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1245
1246       double innerProduct;
1247       int parameterIndex;
1248
1249       for( unsigned int mu = 0; mu < this->m_NumBSplineWeights; mu++ ) {
1250
1251         /* The array weights contains the Jacobian values in a 1-D array
1252          * (because for each parameter the Jacobian is non-zero in only 1 of the
1253          * possible dimensions) which is multiplied by the moving image
1254          * gradient. */
1255         if( this->m_UseCachingOfBSplineWeights ) {
1256           innerProduct = movingImageGradientValue[dim] * weights[mu];
1257           parameterIndex = indices[mu] + this->m_BSplineParametersOffset[dim];
1258         } else {
1259           innerProduct = movingImageGradientValue[dim] * (*weightsHelper)[mu];
1260           parameterIndex = (*indicesHelper)[mu] + this->m_BSplineParametersOffset[dim];
1261         }
1262
1263         const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1264
1265         if( this->m_UseExplicitPDFDerivatives ) {
1266           JointPDFValueType * ptr = derivPtr + parameterIndex;
1267           *(ptr) -= derivativeContribution;
1268         } else {
1269           (*derivativeHelperArray)[parameterIndex] += precomputedWeight * derivativeContribution;
1270         }
1271
1272       } //end mu for loop
1273     } //end dim for loop
1274
1275   } // end if-block transform is BSpline
1276
1277 }
1278
1279
1280 } // end namespace itk
1281
1282
1283 #endif