]> Creatis software - clitk.git/blob - registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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: itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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 __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
37
38
39 // First make sure that the configuration is available.
40 // This line can be removed once the optimized versions
41 // gets integrated into the main directories.
42 #include "itkConfigure.h"
43
44 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
45 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"
46 #else
47
48
49 #include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h"
50 #include "itkBSplineInterpolateImageFunction.h"
51 #include "itkCovariantVector.h"
52 #include "itkImageRandomConstIteratorWithIndex.h"
53 #include "itkImageRegionConstIterator.h"
54 #include "itkImageRegionIterator.h"
55 #include "itkImageIterator.h"
56 #include "vnl/vnl_math.h"
57 #include "itkBSplineDeformableTransform.h"
58
59 namespace itk
60 {
61
62
63 /**
64  * Constructor
65  */
66 template < class TFixedImage, class TMovingImage >
67 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
68 ::MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
69 {
70
71   m_NumberOfSpatialSamples = 500;
72   m_NumberOfHistogramBins = 50;
73
74   this->SetComputeGradient(false); // don't use the default gradient for now
75
76   m_InterpolatorIsBSpline = false;
77   m_TransformIsBSpline    = false;
78
79   // Initialize PDFs to NULL
80   m_JointPDF = NULL;
81   m_JointPDFDerivatives = NULL;
82
83   m_UseExplicitPDFDerivatives = true;
84
85   typename BSplineTransformType::Pointer transformer =
86     BSplineTransformType::New();
87   this->SetTransform (transformer);
88
89   typename BSplineInterpolatorType::Pointer interpolator =
90     BSplineInterpolatorType::New();
91   this->SetInterpolator (interpolator);
92
93   // Initialize memory
94   m_MovingImageNormalizedMin = 0.0;
95   m_FixedImageNormalizedMin = 0.0;
96   m_MovingImageTrueMin = 0.0;
97   m_MovingImageTrueMax = 0.0;
98   m_FixedImageBinSize = 0.0;
99   m_MovingImageBinSize = 0.0;
100   m_CubicBSplineDerivativeKernel = NULL;
101   m_BSplineInterpolator = NULL;
102   m_DerivativeCalculator = NULL;
103   m_NumParametersPerDim = 0;
104   m_NumBSplineWeights = 0;
105   m_BSplineTransform = NULL;
106   m_NumberOfParameters = 0;
107   m_UseAllPixels = false;
108   m_ReseedIterator = false;
109   m_RandomSeed = -1;
110   m_UseCachingOfBSplineWeights = true;
111 }
112
113
114 /**
115  * Print out internal information about this class
116  */
117 template < class TFixedImage, class TMovingImage  >
118 void
119 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
120 ::PrintSelf(std::ostream& os, Indent indent) const
121 {
122
123   Superclass::PrintSelf(os, indent);
124
125   os << indent << "NumberOfSpatialSamples: ";
126   os << m_NumberOfSpatialSamples << std::endl;
127   os << indent << "NumberOfHistogramBins: ";
128   os << m_NumberOfHistogramBins << std::endl;
129   os << indent << "UseAllPixels: ";
130   os << m_UseAllPixels << std::endl;
131
132   // Debugging information
133   os << indent << "NumberOfParameters: ";
134   os << m_NumberOfParameters << std::endl;
135   os << indent << "FixedImageNormalizedMin: ";
136   os << m_FixedImageNormalizedMin << std::endl;
137   os << indent << "MovingImageNormalizedMin: ";
138   os << m_MovingImageNormalizedMin << std::endl;
139   os << indent << "MovingImageTrueMin: ";
140   os << m_MovingImageTrueMin << std::endl;
141   os << indent << "MovingImageTrueMax: ";
142   os << m_MovingImageTrueMax << std::endl;
143   os << indent << "FixedImageBinSize: ";
144   os << m_FixedImageBinSize << std::endl;
145   os << indent << "MovingImageBinSize: ";
146   os << m_MovingImageBinSize << std::endl;
147   os << indent << "InterpolatorIsBSpline: ";
148   os << m_InterpolatorIsBSpline << std::endl;
149   os << indent << "TransformIsBSpline: ";
150   os << m_TransformIsBSpline << std::endl;
151   os << indent << "UseCachingOfBSplineWeights: ";
152   os << m_UseCachingOfBSplineWeights << std::endl;
153   os << indent << "UseExplicitPDFDerivatives: ";
154   os << m_UseExplicitPDFDerivatives << std::endl;
155
156 }
157
158
159 /**
160  * Initialize
161  */
162 template <class TFixedImage, class TMovingImage>
163 void
164 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
165 ::Initialize(void) throw ( ExceptionObject )
166 {
167   this->Superclass::Initialize();
168
169   // Cache the number of transformation parameters
170   m_NumberOfParameters = this->m_Transform->GetNumberOfParameters();
171
172   /**
173    * Compute the minimum and maximum for the FixedImage over
174    * the FixedImageRegion.
175    *
176    * NB: We can't use StatisticsImageFilter to do this because
177    * the filter computes the min/max for the largest possible region
178    */
179   double fixedImageMin = NumericTraits<double>::max();
180   double fixedImageMax = NumericTraits<double>::NonpositiveMin();
181
182   typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
183   FixedIteratorType fixedImageIterator(
184     this->m_FixedImage, this->GetFixedImageRegion() );
185
186   for ( fixedImageIterator.GoToBegin();
187         !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
188
189     double sample = static_cast<double>( fixedImageIterator.Get() );
190
191     if ( sample < fixedImageMin ) {
192       fixedImageMin = sample;
193     }
194
195     if ( sample > fixedImageMax ) {
196       fixedImageMax = sample;
197     }
198   }
199
200
201   /**
202    * Compute the minimum and maximum for the entire moving image
203    * in the buffer.
204    */
205   double movingImageMin = NumericTraits<double>::max();
206   double movingImageMax = NumericTraits<double>::NonpositiveMin();
207
208   typedef ImageRegionConstIterator<MovingImageType> MovingIteratorType;
209   MovingIteratorType movingImageIterator(
210     this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() );
211
212   for ( movingImageIterator.GoToBegin();
213         !movingImageIterator.IsAtEnd(); ++movingImageIterator) {
214     double sample = static_cast<double>( movingImageIterator.Get() );
215
216     if ( sample < movingImageMin ) {
217       movingImageMin = sample;
218     }
219
220     if ( sample > movingImageMax ) {
221       movingImageMax = sample;
222     }
223   }
224
225   m_MovingImageTrueMin = movingImageMin;
226   m_MovingImageTrueMax = movingImageMax;
227
228   itkDebugMacro( " FixedImageMin: " << fixedImageMin <<
229                  " FixedImageMax: " << fixedImageMax << std::endl );
230   itkDebugMacro( " MovingImageMin: " << movingImageMin <<
231                  " MovingImageMax: " << movingImageMax << std::endl );
232
233
234   /**
235    * Compute binsize for the histograms.
236    *
237    * The binsize for the image intensities needs to be adjusted so that
238    * we can avoid dealing with boundary conditions using the cubic
239    * spline as the Parzen window.  We do this by increasing the size
240    * of the bins so that the joint histogram becomes "padded" at the
241    * borders. Because we are changing the binsize,
242    * we also need to shift the minimum by the padded amount in order to
243    * avoid minimum values filling in our padded region.
244    *
245    * Note that there can still be non-zero bin values in the padded region,
246    * it's just that these bins will never be a central bin for the Parzen
247    * window.
248    *
249    */
250   const int padding = 2;  // this will pad by 2 bins
251
252   m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) /
253                         static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
254   m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize -
255                               static_cast<double>( padding );
256
257   m_MovingImageBinSize = ( movingImageMax - movingImageMin ) /
258                          static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
259   m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize -
260                                static_cast<double>( padding );
261
262
263   itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin );
264   itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin );
265   itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize );
266   itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize );
267
268   if( m_UseAllPixels ) {
269     m_NumberOfSpatialSamples =
270       this->GetFixedImageRegion().GetNumberOfPixels();
271   }
272
273   /**
274    * Allocate memory for the fixed image sample container.
275    */
276   m_FixedImageSamples.resize( m_NumberOfSpatialSamples );
277
278
279   /**
280    * Allocate memory for the marginal PDF and initialize values
281    * to zero. The marginal PDFs are stored as std::vector.
282    */
283   m_FixedImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 );
284   m_MovingImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 );
285
286   /**
287    * Allocate memory for the joint PDF and joint PDF derivatives.
288    * The joint PDF and joint PDF derivatives are store as itk::Image.
289    */
290   m_JointPDF = JointPDFType::New();
291
292   // Instantiate a region, index, size
293   JointPDFRegionType            jointPDFRegion;
294   JointPDFIndexType              jointPDFIndex;
295   JointPDFSizeType              jointPDFSize;
296
297   // Deallocate the memory that may have been allocated for
298   // previous runs of the metric.
299   this->m_JointPDFDerivatives = NULL;  // by destroying the dynamic array
300   this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones
301   this->m_MetricDerivative = DerivativeType( 1 );
302
303   //
304   // Now allocate memory according to the user-selected method.
305   //
306   if( this->m_UseExplicitPDFDerivatives ) {
307     this->m_JointPDFDerivatives = JointPDFDerivativesType::New();
308     JointPDFDerivativesRegionType  jointPDFDerivativesRegion;
309     JointPDFDerivativesIndexType  jointPDFDerivativesIndex;
310     JointPDFDerivativesSizeType    jointPDFDerivativesSize;
311
312     // For the derivatives of the joint PDF define a region starting from {0,0,0}
313     // with size {m_NumberOfParameters,m_NumberOfHistogramBins,
314     // m_NumberOfHistogramBins}. The dimension represents transform parameters,
315     // fixed image parzen window index and moving image parzen window index,
316     // respectively.
317     jointPDFDerivativesIndex.Fill( 0 );
318     jointPDFDerivativesSize[0] = m_NumberOfParameters;
319     jointPDFDerivativesSize[1] = m_NumberOfHistogramBins;
320     jointPDFDerivativesSize[2] = m_NumberOfHistogramBins;
321
322     jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex );
323     jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize );
324
325     // Set the regions and allocate
326     m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion );
327     m_JointPDFDerivatives->Allocate();
328   } else {
329     /** Allocate memory for helper array that will contain the pRatios
330      *  for each bin of the joint histogram. This is part of the effort
331      *  for flattening the computation of the PDF Jacobians.
332      */
333     this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins );
334     this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() );
335   }
336
337   // For the joint PDF define a region starting from {0,0}
338   // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}.
339   // The dimension represents fixed image parzen window index
340   // and moving image parzen window index, respectively.
341   jointPDFIndex.Fill( 0 );
342   jointPDFSize.Fill( m_NumberOfHistogramBins );
343
344   jointPDFRegion.SetIndex( jointPDFIndex );
345   jointPDFRegion.SetSize( jointPDFSize );
346
347   // Set the regions and allocate
348   m_JointPDF->SetRegions( jointPDFRegion );
349   m_JointPDF->Allocate();
350
351
352   /**
353    * Setup the kernels used for the Parzen windows.
354    */
355   m_CubicBSplineKernel = CubicBSplineFunctionType::New();
356   m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New();
357
358
359   if( m_UseAllPixels ) {
360     /**
361      * Take all the pixels within the fixed image region)
362      * to create the sample points list.
363      */
364     this->SampleFullFixedImageDomain( m_FixedImageSamples );
365   } else {
366     /**
367      * Uniformly sample the fixed image (within the fixed image region)
368      * to create the sample points list.
369      */
370     this->SampleFixedImageDomain( m_FixedImageSamples );
371   }
372
373   /**
374    * Pre-compute the fixed image parzen window index for
375    * each point of the fixed image sample points list.
376    */
377   this->ComputeFixedImageParzenWindowIndices( m_FixedImageSamples );
378
379   /**
380    * Check if the interpolator is of type BSplineInterpolateImageFunction.
381    * If so, we can make use of its EvaluateDerivatives method.
382    * Otherwise, we instantiate an external central difference
383    * derivative calculator.
384    *
385    * TODO: Also add it the possibility of using the default gradient
386    * provided by the superclass.
387    *
388    */
389   m_InterpolatorIsBSpline = true;
390
391   BSplineInterpolatorType * testPtr = dynamic_cast<BSplineInterpolatorType *>(
392                                         this->m_Interpolator.GetPointer() );
393   if ( !testPtr ) {
394     m_InterpolatorIsBSpline = false;
395
396     m_DerivativeCalculator = DerivativeFunctionType::New();
397
398 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
399     m_DerivativeCalculator->UseImageDirectionOn();
400 #endif
401
402     m_DerivativeCalculator->SetInputImage( this->m_MovingImage );
403
404     m_BSplineInterpolator = NULL;
405     itkDebugMacro( "Interpolator is not BSpline" );
406   } else {
407     m_BSplineInterpolator = testPtr;
408
409 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
410     m_BSplineInterpolator->UseImageDirectionOn();
411 #endif
412
413     m_DerivativeCalculator = NULL;
414     itkDebugMacro( "Interpolator is BSpline" );
415   }
416
417   /**
418    * Check if the transform is of type BSplineDeformableTransform.
419    *
420    * If so, several speed up features are implemented.
421    * [1] Precomputing the results of bulk transform for each sample point.
422    * [2] Precomputing the BSpline weights for each sample point,
423    *     to be used later to directly compute the deformation vector
424    * [3] Precomputing the indices of the parameters within the
425    *     the support region of each sample point.
426    */
427   m_TransformIsBSpline = true;
428
429   BSplineTransformType * testPtr2 = dynamic_cast<BSplineTransformType *>(
430                                       this->m_Transform.GetPointer() );
431   if( !testPtr2 ) {
432     m_TransformIsBSpline = false;
433     m_BSplineTransform = NULL;
434     itkDebugMacro( "Transform is not BSplineDeformable" );
435   } else {
436     m_BSplineTransform = testPtr2;
437     m_NumParametersPerDim =
438       m_BSplineTransform->GetNumberOfParametersPerDimension();
439     m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights();
440     itkDebugMacro( "Transform is BSplineDeformable" );
441   }
442
443   if( this->m_TransformIsBSpline ) {
444     // First, deallocate memory that may have been used from previous run of the Metric
445     this->m_BSplineTransformWeightsArray.SetSize( 1, 1 );
446     this->m_BSplineTransformIndicesArray.SetSize( 1, 1 );
447     this->m_PreTransformPointsArray.resize( 1 );
448     this->m_WithinSupportRegionArray.resize( 1 );
449     this->m_Weights.SetSize( 1 );
450     this->m_Indices.SetSize( 1 );
451
452
453     if( this->m_UseCachingOfBSplineWeights ) {
454       m_BSplineTransformWeightsArray.SetSize(
455         m_NumberOfSpatialSamples, m_NumBSplineWeights );
456       m_BSplineTransformIndicesArray.SetSize(
457         m_NumberOfSpatialSamples, m_NumBSplineWeights );
458       m_PreTransformPointsArray.resize( m_NumberOfSpatialSamples );
459       m_WithinSupportRegionArray.resize( m_NumberOfSpatialSamples );
460
461       this->PreComputeTransformValues();
462     } else {
463       this->m_Weights.SetSize( this->m_NumBSplineWeights );
464       this->m_Indices.SetSize( this->m_NumBSplineWeights );
465     }
466
467     for ( unsigned int j = 0; j < FixedImageDimension; j++ ) {
468       m_ParametersOffset[j] = j *
469                               m_BSplineTransform->GetNumberOfParametersPerDimension();
470     }
471   }
472
473 }
474
475
476 /**
477  * Uniformly sample the fixed image domain using a random walk
478  */
479 template < class TFixedImage, class TMovingImage >
480 void
481 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
482 ::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples )
483 {
484
485   // Set up a random interator within the user specified fixed image region.
486   typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
487   RandomIterator randIter( this->m_FixedImage, this->GetFixedImageRegion() );
488
489   randIter.SetNumberOfSamples( m_NumberOfSpatialSamples );
490   randIter.GoToBegin();
491
492   typename FixedImageSpatialSampleContainer::iterator iter;
493   typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
494
495   if( this->m_FixedImageMask ) {
496
497     InputPointType inputPoint;
498
499     iter=samples.begin();
500     int count = 0;
501     int samples_found = 0;
502     int maxcount = m_NumberOfSpatialSamples * 10;
503     while( iter != end ) {
504
505       if ( count > maxcount ) {
506 #if 0
507         itkExceptionMacro(
508           "Drew too many samples from the mask (is it too small?): "
509           << maxcount << std::endl );
510 #else
511 samples.resize(samples_found);
512 //        this->SetNumberOfSpatialSamples(sample_found);
513 break;
514 #endif
515       }
516       count++;
517
518       // Get sampled index
519       FixedImageIndexType index = randIter.GetIndex();
520       // Check if the Index is inside the mask, translate index to point
521       this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
522
523       // If not inside the mask, ignore the point
524       if( !this->m_FixedImageMask->IsInside( inputPoint ) ) {
525         ++randIter; // jump to another random position
526         continue;
527       }
528
529       // Get sampled fixed image value
530       (*iter).FixedImageValue = randIter.Get();
531       // Translate index to point
532       (*iter).FixedImagePointValue = inputPoint;
533       samples_found++;
534       // Jump to random position
535       ++randIter;
536       ++iter;
537     }
538   } else {
539     for( iter=samples.begin(); iter != end; ++iter ) {
540       // Get sampled index
541       FixedImageIndexType index = randIter.GetIndex();
542       // Get sampled fixed image value
543       (*iter).FixedImageValue = randIter.Get();
544       // Translate index to point
545       this->m_FixedImage->TransformIndexToPhysicalPoint( index,
546           (*iter).FixedImagePointValue );
547       // Jump to random position
548       ++randIter;
549
550     }
551   }
552 }
553
554 /**
555  * Sample the fixed image domain using all pixels in the Fixed image region
556  */
557 template < class TFixedImage, class TMovingImage >
558 void
559 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
560 ::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples )
561 {
562
563   // Set up a region interator within the user specified fixed image region.
564   typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
565   RegionIterator regionIter( this->m_FixedImage, this->GetFixedImageRegion() );
566
567   regionIter.GoToBegin();
568
569   typename FixedImageSpatialSampleContainer::iterator iter;
570   typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
571
572   if( this->m_FixedImageMask ) {
573     InputPointType inputPoint;
574
575     iter=samples.begin();
576     unsigned long nSamplesPicked = 0;
577
578     while( iter != end && !regionIter.IsAtEnd() ) {
579       // Get sampled index
580       FixedImageIndexType index = regionIter.GetIndex();
581       // Check if the Index is inside the mask, translate index to point
582       this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
583
584       // If not inside the mask, ignore the point
585       if( !this->m_FixedImageMask->IsInside( inputPoint ) ) {
586         ++regionIter; // jump to next pixel
587         continue;
588       }
589
590       // Get sampled fixed image value
591       (*iter).FixedImageValue = regionIter.Get();
592       // Translate index to point
593       (*iter).FixedImagePointValue = inputPoint;
594
595       ++regionIter;
596       ++iter;
597       ++nSamplesPicked;
598     }
599
600     // If we picked fewer samples than the desired number,
601     // resize the container
602     if (nSamplesPicked != this->m_NumberOfSpatialSamples) {
603       this->m_NumberOfSpatialSamples = nSamplesPicked;
604       samples.resize(this->m_NumberOfSpatialSamples);
605     }
606   } else { // not restricting sample throwing to a mask
607
608     // cannot sample more than the number of pixels in the image region
609     if (  this->m_NumberOfSpatialSamples
610           > this->GetFixedImageRegion().GetNumberOfPixels()) {
611       this->m_NumberOfSpatialSamples
612       = this->GetFixedImageRegion().GetNumberOfPixels();
613       samples.resize(this->m_NumberOfSpatialSamples);
614     }
615
616     for( iter=samples.begin(); iter != end; ++iter ) {
617       // Get sampled index
618       FixedImageIndexType index = regionIter.GetIndex();
619       // Get sampled fixed image value
620       (*iter).FixedImageValue = regionIter.Get();
621       // Translate index to point
622       this->m_FixedImage->TransformIndexToPhysicalPoint( index,
623           (*iter).FixedImagePointValue );
624       ++regionIter;
625     }
626   }
627 }
628
629 /**
630  * Uniformly sample the fixed image domain using a random walk
631  */
632 template < class TFixedImage, class TMovingImage >
633 void
634 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
635 ::ComputeFixedImageParzenWindowIndices(
636   FixedImageSpatialSampleContainer& samples )
637 {
638
639   typename FixedImageSpatialSampleContainer::iterator iter;
640   typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
641
642   for( iter=samples.begin(); iter != end; ++iter ) {
643
644     // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
645     double windowTerm =
646       static_cast<double>( (*iter).FixedImageValue ) / m_FixedImageBinSize -
647       m_FixedImageNormalizedMin;
648     unsigned int pindex = static_cast<unsigned int>( vcl_floor(windowTerm ) );
649
650     // Make sure the extreme values are in valid bins
651     if ( pindex < 2 ) {
652       pindex = 2;
653     } else if ( pindex > (m_NumberOfHistogramBins - 3) ) {
654       pindex = m_NumberOfHistogramBins - 3;
655     }
656
657     (*iter).FixedImageParzenWindowIndex = pindex;
658
659   }
660
661 }
662
663 /**
664  * Get the match Measure
665  */
666 template < class TFixedImage, class TMovingImage  >
667 typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
668 ::MeasureType
669 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
670 ::GetValue( const ParametersType& parameters ) const
671 {
672
673   // Reset marginal pdf to all zeros.
674   // Assumed the size has already been set to NumberOfHistogramBins
675   // in Initialize().
676   for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) {
677     m_FixedImageMarginalPDF[j]  = 0.0;
678     m_MovingImageMarginalPDF[j] = 0.0;
679   }
680
681   // Reset the joint pdfs to zero
682   m_JointPDF->FillBuffer( 0.0 );
683
684   // Set up the parameters in the transform
685   this->m_Transform->SetParameters( parameters );
686
687
688   // Declare iterators for iteration over the sample container
689   typename FixedImageSpatialSampleContainer::const_iterator fiter;
690   typename FixedImageSpatialSampleContainer::const_iterator fend =
691     m_FixedImageSamples.end();
692
693   unsigned long nSamples=0;
694   unsigned long nFixedImageSamples=0;
695
696
697   for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
698
699     // Get moving image value
700     MovingImagePointType mappedPoint;
701     bool sampleOk;
702     double movingImageValue;
703
704     this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
705                           sampleOk, movingImageValue );
706
707     ++nFixedImageSamples;
708
709     if( sampleOk ) {
710
711       ++nSamples;
712
713       /**
714        * Compute this sample's contribution to the marginal and
715        * joint distributions.
716        *
717        */
718
719       // Determine parzen window arguments (see eqn 6 of Mattes paper [2])
720       double movingImageParzenWindowTerm =
721         movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
722       unsigned int movingImageParzenWindowIndex =
723         static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
724
725       // Make sure the extreme values are in valid bins
726       if ( movingImageParzenWindowIndex < 2 ) {
727         movingImageParzenWindowIndex = 2;
728       } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
729         movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
730       }
731
732
733       // Since a zero-order BSpline (box car) kernel is used for
734       // the fixed image marginal pdf, we need only increment the
735       // fixedImageParzenWindowIndex by value of 1.0.
736       m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] +=
737         static_cast<PDFValueType>( 1 );
738
739       /**
740        * The region of support of the parzen window determines which bins
741        * of the joint PDF are effected by the pair of image values.
742        * Since we are using a cubic spline for the moving image parzen
743        * window, four bins are affected.  The fixed image parzen window is
744        * a zero-order spline (box car) and thus effects only one bin.
745        *
746        *  The PDF is arranged so that moving image bins corresponds to the
747        * zero-th (column) dimension and the fixed image bins corresponds
748        * to the first (row) dimension.
749        *
750        */
751
752       // Pointer to affected bin to be updated
753       JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
754                                   ( (*fiter).FixedImageParzenWindowIndex
755                                     * m_JointPDF->GetOffsetTable()[1] );
756
757       // Move the pointer to the first affected bin
758       int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
759       pdfPtr += pdfMovingIndex;
760
761       for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
762            + 2;
763            pdfMovingIndex++, pdfPtr++ ) {
764
765         // Update PDF for the current intensity pair
766         double movingImageParzenWindowArg =
767           static_cast<double>( pdfMovingIndex ) -
768           static_cast<double>( movingImageParzenWindowTerm );
769
770         *(pdfPtr) += static_cast<PDFValueType>(
771                        m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );
772
773       }  //end parzen windowing for loop
774
775     } //end if-block check sampleOk
776   } // end iterating over fixed image spatial sample container for loop
777
778   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
779                  << nSamples << " / " << m_NumberOfSpatialSamples
780                  << std::endl );
781
782   if( nSamples < m_NumberOfSpatialSamples / 16 ) {
783     itkExceptionMacro( "Too many samples map outside moving image buffer: "
784                        << nSamples << " / " << m_NumberOfSpatialSamples
785                        << std::endl );
786   }
787
788   this->m_NumberOfPixelsCounted = nSamples;
789
790
791   /**
792    * Normalize the PDFs, compute moving image marginal PDF
793    *
794    */
795   typedef ImageRegionIterator<JointPDFType> JointPDFIteratorType;
796   JointPDFIteratorType jointPDFIterator ( m_JointPDF,
797                                           m_JointPDF->GetBufferedRegion() );
798
799   jointPDFIterator.GoToBegin();
800
801   // Compute joint PDF normalization factor
802   // (to ensure joint PDF sum adds to 1.0)
803   double jointPDFSum = 0.0;
804
805   while( !jointPDFIterator.IsAtEnd() ) {
806     jointPDFSum += jointPDFIterator.Get();
807     ++jointPDFIterator;
808   }
809
810   if ( jointPDFSum == 0.0 ) {
811     itkExceptionMacro( "Joint PDF summed to zero" );
812   }
813
814
815   // Normalize the PDF bins
816   jointPDFIterator.GoToEnd();
817   while( !jointPDFIterator.IsAtBegin() ) {
818     --jointPDFIterator;
819     jointPDFIterator.Value() /= static_cast<PDFValueType>( jointPDFSum );
820   }
821
822
823   // Normalize the fixed image marginal PDF
824   double fixedPDFSum = 0.0;
825   for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) {
826     fixedPDFSum += m_FixedImageMarginalPDF[bin];
827   }
828
829   if ( fixedPDFSum == 0.0 ) {
830     itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
831   }
832
833   for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
834     m_FixedImageMarginalPDF[bin] /= static_cast<PDFValueType>( fixedPDFSum );
835   }
836
837
838   // Compute moving image marginal PDF by summing over fixed image bins.
839   typedef ImageLinearIteratorWithIndex<JointPDFType> JointPDFLinearIterator;
840   JointPDFLinearIterator linearIter( m_JointPDF,
841                                      m_JointPDF->GetBufferedRegion() );
842
843   linearIter.SetDirection( 1 );
844   linearIter.GoToBegin();
845   unsigned int movingIndex1 = 0;
846
847   while( !linearIter.IsAtEnd() ) {
848
849     double sum = 0.0;
850
851     while( !linearIter.IsAtEndOfLine() ) {
852       sum += linearIter.Get();
853       ++linearIter;
854     }
855
856     m_MovingImageMarginalPDF[movingIndex1] = static_cast<PDFValueType>(sum);
857
858     linearIter.NextLine();
859     ++movingIndex1;
860
861   }
862
863   /**
864    * Compute the metric by double summation over histogram.
865    */
866
867   // Setup pointer to point to the first bin
868   JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
869
870   // Initialze sum to zero
871   double sum = 0.0;
872
873   for( unsigned int fixedIndex = 0;
874        fixedIndex < m_NumberOfHistogramBins;
875        ++fixedIndex ) {
876
877     double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
878
879     for( unsigned int movingIndex = 0;
880          movingIndex < m_NumberOfHistogramBins;
881          ++movingIndex, jointPDFPtr++ ) {
882
883       double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
884       double jointPDFValue = *(jointPDFPtr);
885
886       // check for non-zero bin contribution
887       if( jointPDFValue > 1e-16 &&  movingImagePDFValue > 1e-16 ) {
888
889         double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
890         if( fixedImagePDFValue > 1e-16) {
891           sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
892         }
893
894       }  // end if-block to check non-zero bin contribution
895     }  // end for-loop over moving index
896   }  // end for-loop over fixed index
897
898   return static_cast<MeasureType>( -1.0 * sum );
899
900 }
901
902
903 /**
904  * Get the both Value and Derivative Measure
905  */
906 template < class TFixedImage, class TMovingImage  >
907 void
908 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
909 ::GetValueAndDerivative(
910   const ParametersType& parameters,
911   MeasureType& value,
912   DerivativeType& derivative) const
913 {
914
915   // Set output values to zero
916   value = NumericTraits< MeasureType >::Zero;
917
918   if( this->m_UseExplicitPDFDerivatives ) {
919     m_JointPDFDerivatives->FillBuffer( 0.0 );
920     derivative = DerivativeType( this->GetNumberOfParameters() );
921     derivative.Fill( NumericTraits< MeasureType >::Zero );
922   } else {
923     this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero );
924     this->m_PRatioArray.Fill( 0.0 );
925   }
926
927   // Reset marginal pdf to all zeros.
928   // Assumed the size has already been set to NumberOfHistogramBins
929   // in Initialize().
930   for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) {
931     m_FixedImageMarginalPDF[j]  = 0.0;
932     m_MovingImageMarginalPDF[j] = 0.0;
933   }
934
935   // Reset the joint pdfs to zero
936   m_JointPDF->FillBuffer( 0.0 );
937
938
939   // Set up the parameters in the transform
940   this->m_Transform->SetParameters( parameters );
941
942
943   // Declare iterators for iteration over the sample container
944   typename FixedImageSpatialSampleContainer::const_iterator fiter;
945   typename FixedImageSpatialSampleContainer::const_iterator fend =
946     m_FixedImageSamples.end();
947
948   unsigned long nSamples=0;
949   unsigned long nFixedImageSamples=0;
950
951   for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
952
953     // Get moving image value
954     MovingImagePointType mappedPoint;
955     bool sampleOk;
956     double movingImageValue;
957
958
959     this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
960                           sampleOk, movingImageValue );
961
962     if( sampleOk ) {
963       ++nSamples;
964
965       // Get moving image derivative at the mapped position
966       ImageDerivativesType movingImageGradientValue;
967       this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue );
968
969
970       /**
971        * Compute this sample's contribution to the marginal
972        *   and joint distributions.
973        *
974        */
975
976       // Determine parzen window arguments (see eqn 6 of Mattes paper [2])
977       double movingImageParzenWindowTerm =
978         movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
979       unsigned int movingImageParzenWindowIndex =
980         static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
981
982       // Make sure the extreme values are in valid bins
983       if ( movingImageParzenWindowIndex < 2 ) {
984         movingImageParzenWindowIndex = 2;
985       } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
986         movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
987       }
988
989
990       // Since a zero-order BSpline (box car) kernel is used for
991       // the fixed image marginal pdf, we need only increment the
992       // fixedImageParzenWindowIndex by value of 1.0.
993       m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] +=
994         static_cast<PDFValueType>( 1 );
995
996       /**
997        * The region of support of the parzen window determines which bins
998        * of the joint PDF are effected by the pair of image values.
999        * Since we are using a cubic spline for the moving image parzen
1000        * window, four bins are effected.  The fixed image parzen window is
1001        * a zero-order spline (box car) and thus effects only one bin.
1002        *
1003        *  The PDF is arranged so that moving image bins corresponds to the
1004        * zero-th (column) dimension and the fixed image bins corresponds
1005        * to the first (row) dimension.
1006        *
1007        */
1008
1009       // Pointer to affected bin to be updated
1010       JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
1011                                   ( (*fiter).FixedImageParzenWindowIndex * m_NumberOfHistogramBins );
1012
1013       // Move the pointer to the fist affected bin
1014       int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
1015       pdfPtr += pdfMovingIndex;
1016
1017       for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
1018            + 2;
1019            pdfMovingIndex++, pdfPtr++ ) {
1020         // Update PDF for the current intensity pair
1021         double movingImageParzenWindowArg =
1022           static_cast<double>( pdfMovingIndex ) -
1023           static_cast<double>(movingImageParzenWindowTerm);
1024
1025         *(pdfPtr) += static_cast<PDFValueType>(
1026                        m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );
1027
1028         if( this->m_UseExplicitPDFDerivatives ) {
1029           // Compute the cubicBSplineDerivative for later repeated use.
1030           double cubicBSplineDerivativeValue =
1031             m_CubicBSplineDerivativeKernel->Evaluate(
1032               movingImageParzenWindowArg );
1033
1034           // Compute PDF derivative contribution.
1035           this->ComputePDFDerivatives( nFixedImageSamples,
1036                                        pdfMovingIndex,
1037                                        movingImageGradientValue,
1038                                        cubicBSplineDerivativeValue );
1039         }
1040
1041       }  //end parzen windowing for loop
1042
1043     } //end if-block check sampleOk
1044
1045     ++nFixedImageSamples;
1046
1047   } // end iterating over fixed image spatial sample container for loop
1048
1049   itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
1050                  << nSamples << " / " << m_NumberOfSpatialSamples
1051                  << std::endl );
1052
1053   if( nSamples < m_NumberOfSpatialSamples / 16 ) {
1054     itkExceptionMacro( "Too many samples map outside moving image buffer: "
1055                        << nSamples << " / " << m_NumberOfSpatialSamples
1056                        << std::endl );
1057   }
1058
1059   this->m_NumberOfPixelsCounted = nSamples;
1060
1061   /**
1062    * Normalize the PDFs, compute moving image marginal PDF
1063    *
1064    */
1065   typedef ImageRegionIterator<JointPDFType> JointPDFIteratorType;
1066   JointPDFIteratorType jointPDFIterator ( m_JointPDF,
1067                                           m_JointPDF->GetBufferedRegion() );
1068
1069   jointPDFIterator.GoToBegin();
1070
1071
1072   // Compute joint PDF normalization factor
1073   // (to ensure joint PDF sum adds to 1.0)
1074   double jointPDFSum = 0.0;
1075
1076   while( !jointPDFIterator.IsAtEnd() ) {
1077     jointPDFSum += jointPDFIterator.Get();
1078     ++jointPDFIterator;
1079   }
1080
1081   if ( jointPDFSum == 0.0 ) {
1082     itkExceptionMacro( "Joint PDF summed to zero" );
1083   }
1084
1085
1086   // Normalize the PDF bins
1087   jointPDFIterator.GoToEnd();
1088   while( !jointPDFIterator.IsAtBegin() ) {
1089     --jointPDFIterator;
1090     jointPDFIterator.Value() /= static_cast<PDFValueType>( jointPDFSum );
1091   }
1092
1093
1094   // Normalize the fixed image marginal PDF
1095   double fixedPDFSum = 0.0;
1096   for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) {
1097     fixedPDFSum += m_FixedImageMarginalPDF[bin];
1098   }
1099
1100   if ( fixedPDFSum == 0.0 ) {
1101     itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
1102   }
1103
1104   for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
1105     m_FixedImageMarginalPDF[bin] /= static_cast<PDFValueType>( fixedPDFSum );
1106   }
1107
1108
1109   // Compute moving image marginal PDF by summing over fixed image bins.
1110   typedef ImageLinearIteratorWithIndex<JointPDFType> JointPDFLinearIterator;
1111   JointPDFLinearIterator linearIter(
1112     m_JointPDF, m_JointPDF->GetBufferedRegion() );
1113
1114   linearIter.SetDirection( 1 );
1115   linearIter.GoToBegin();
1116   unsigned int movingIndex1 = 0;
1117
1118   while( !linearIter.IsAtEnd() ) {
1119
1120     double sum = 0.0;
1121
1122     while( !linearIter.IsAtEndOfLine() ) {
1123       sum += linearIter.Get();
1124       ++linearIter;
1125     }
1126
1127     m_MovingImageMarginalPDF[movingIndex1] = static_cast<PDFValueType>(sum);
1128
1129     linearIter.NextLine();
1130     ++movingIndex1;
1131
1132   }
1133
1134   double nFactor = 1.0 / ( m_MovingImageBinSize
1135                            * static_cast<double>( nSamples ) );
1136
1137   if( this->m_UseExplicitPDFDerivatives ) {
1138     // Normalize the joint PDF derivatives by the test image binsize and nSamples
1139     typedef ImageRegionIterator<JointPDFDerivativesType>
1140     JointPDFDerivativesIteratorType;
1141     JointPDFDerivativesIteratorType jointPDFDerivativesIterator (
1142       m_JointPDFDerivatives,
1143       m_JointPDFDerivatives->GetBufferedRegion()
1144     );
1145
1146     jointPDFDerivativesIterator.GoToBegin();
1147
1148     while( !jointPDFDerivativesIterator.IsAtEnd() ) {
1149       jointPDFDerivativesIterator.Value() *= nFactor;
1150       ++jointPDFDerivativesIterator;
1151     }
1152   }
1153
1154   /**
1155    * Compute the metric by double summation over histogram.
1156    */
1157
1158   // Setup pointer to point to the first bin
1159   JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
1160
1161   // Initialize sum to zero
1162   double sum = 0.0;
1163
1164   for( unsigned int fixedIndex = 0;
1165        fixedIndex < m_NumberOfHistogramBins;
1166        ++fixedIndex ) {
1167     double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
1168
1169     for( unsigned int movingIndex = 0; movingIndex < m_NumberOfHistogramBins;
1170          ++movingIndex, jointPDFPtr++ ) {
1171       double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
1172       double jointPDFValue = *(jointPDFPtr);
1173
1174       // check for non-zero bin contribution
1175       if( jointPDFValue > 1e-16 &&  movingImagePDFValue > 1e-16 ) {
1176
1177         double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
1178
1179         if( fixedImagePDFValue > 1e-16) {
1180           sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
1181         }
1182
1183         if( this->m_UseExplicitPDFDerivatives ) {
1184           // move joint pdf derivative pointer to the right position
1185           JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1186                                          + ( fixedIndex  * m_JointPDFDerivatives->GetOffsetTable()[2] )
1187                                          + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1188
1189           for( unsigned int parameter=0; parameter < m_NumberOfParameters; ++parameter, derivPtr++ ) {
1190             // Ref: eqn 23 of Thevenaz & Unser paper [3]
1191             derivative[parameter] -= (*derivPtr) * pRatio;
1192           }  // end for-loop over parameters
1193         } else {
1194           this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor;
1195         }
1196       }  // end if-block to check non-zero bin contribution
1197     }  // end for-loop over moving index
1198   }  // end for-loop over fixed index
1199
1200   if( !(this->m_UseExplicitPDFDerivatives ) ) {
1201     // Second pass: This one is done for accumulating the contributions
1202     //              to the derivative array.
1203
1204     nFixedImageSamples = 0;
1205
1206     for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
1207
1208       // Get moving image value
1209       MovingImagePointType mappedPoint;
1210       bool sampleOk;
1211       double movingImageValue;
1212
1213       this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
1214                             sampleOk, movingImageValue );
1215
1216       if( sampleOk ) {
1217         // Get moving image derivative at the mapped position
1218         ImageDerivativesType movingImageGradientValue;
1219         this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue );
1220
1221
1222         /**
1223          * Compute this sample's contribution to the marginal
1224          *   and joint distributions.
1225          *
1226          */
1227
1228         // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
1229         double movingImageParzenWindowTerm =
1230           movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
1231         unsigned int movingImageParzenWindowIndex =
1232           static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
1233
1234         // Make sure the extreme values are in valid bins
1235         if ( movingImageParzenWindowIndex < 2 ) {
1236           movingImageParzenWindowIndex = 2;
1237         } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
1238           movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
1239         }
1240
1241
1242         // Move the pointer to the fist affected bin
1243         int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
1244
1245         for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
1246              + 2;
1247              pdfMovingIndex++ ) {
1248
1249           // Update PDF for the current intensity pair
1250           double movingImageParzenWindowArg =
1251             static_cast<double>( pdfMovingIndex ) -
1252             static_cast<double>(movingImageParzenWindowTerm);
1253
1254           // Compute the cubicBSplineDerivative for later repeated use.
1255           double cubicBSplineDerivativeValue =
1256             m_CubicBSplineDerivativeKernel->Evaluate(
1257               movingImageParzenWindowArg );
1258
1259           // Compute PDF derivative contribution.
1260           this->ComputePDFDerivatives( nFixedImageSamples,
1261                                        pdfMovingIndex,
1262                                        movingImageGradientValue,
1263                                        cubicBSplineDerivativeValue );
1264
1265
1266         }  //end parzen windowing for loop
1267
1268       } //end if-block check sampleOk
1269
1270       ++nFixedImageSamples;
1271
1272     } // end iterating over fixed image spatial sample container for loop
1273
1274
1275     derivative = this->m_MetricDerivative;
1276   }
1277
1278   value = static_cast<MeasureType>( -1.0 * sum );
1279 }
1280
1281
1282 /**
1283  * Get the match measure derivative
1284  */
1285 template < class TFixedImage, class TMovingImage  >
1286 void
1287 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1288 ::GetDerivative( const ParametersType& parameters,
1289                  DerivativeType & derivative ) const
1290 {
1291   MeasureType value;
1292   // call the combined version
1293   this->GetValueAndDerivative( parameters, value, derivative );
1294 }
1295
1296
1297 /**
1298  * Compute image derivatives using a central difference function
1299  * if we are not using a BSplineInterpolator, which includes
1300  * derivatives.
1301  */
1302 template < class TFixedImage, class TMovingImage >
1303 void
1304 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1305 ::ComputeImageDerivatives(
1306   const MovingImagePointType& mappedPoint,
1307   ImageDerivativesType& gradient ) const
1308 {
1309
1310   if( m_InterpolatorIsBSpline ) {
1311     // Computed moving image gradient using derivative BSpline kernel.
1312     gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint );
1313   } else {
1314     // For all generic interpolator use central differencing.
1315     gradient = m_DerivativeCalculator->Evaluate( mappedPoint );
1316   }
1317
1318 }
1319
1320
1321 /**
1322  * Transform a point from FixedImage domain to MovingImage domain.
1323  * This function also checks if mapped point is within support region.
1324  */
1325 template < class TFixedImage, class TMovingImage >
1326 void
1327 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1328 ::TransformPoint(
1329   unsigned int sampleNumber,
1330   const ParametersType& parameters,
1331   MovingImagePointType& mappedPoint,
1332   bool& sampleOk,
1333   double& movingImageValue ) const
1334 {
1335
1336   if ( !m_TransformIsBSpline ) {
1337     // Use generic transform to compute mapped position
1338     mappedPoint = this->m_Transform->TransformPoint(
1339                     m_FixedImageSamples[sampleNumber].FixedImagePointValue );
1340
1341     // Check if mapped point inside image buffer
1342     sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint );
1343   } else {
1344
1345     if( this->m_UseCachingOfBSplineWeights ) {
1346       // If the transform is BSplineDeformable, we can use the precomputed
1347       // weights and indices to obtained the mapped position
1348       //
1349       const WeightsValueType * weights =
1350         m_BSplineTransformWeightsArray[sampleNumber];
1351       const IndexValueType   * indices =
1352         m_BSplineTransformIndicesArray[sampleNumber];
1353       mappedPoint.Fill( 0.0 );
1354
1355       if ( m_WithinSupportRegionArray[sampleNumber] ) {
1356         for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ ) {
1357           for ( unsigned int j = 0; j < FixedImageDimension; j++ ) {
1358             mappedPoint[j] += weights[k] *
1359                               parameters[ indices[k] + m_ParametersOffset[j] ];
1360           }
1361         }
1362       }
1363
1364       for( unsigned int j = 0; j < FixedImageDimension; j++ ) {
1365         mappedPoint[j] += m_PreTransformPointsArray[sampleNumber][j];
1366       }
1367
1368       // Check if mapped point inside image buffer
1369       sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint );
1370
1371       // Check if mapped point is within the support region of a grid point.
1372       // This is neccessary for computing the metric gradient
1373       sampleOk = sampleOk && m_WithinSupportRegionArray[sampleNumber];
1374     } else {
1375       // If not caching values, we invoke the Transform to recompute the
1376       // mapping of the point.
1377       this->m_BSplineTransform->TransformPoint(
1378         this->m_FixedImageSamples[sampleNumber].FixedImagePointValue,
1379         mappedPoint, this->m_Weights, this->m_Indices, sampleOk);
1380
1381       // Check if mapped point inside image buffer
1382       sampleOk = sampleOk && this->m_Interpolator->IsInsideBuffer( mappedPoint );
1383     }
1384
1385   }
1386
1387   // If user provided a mask over the Moving image
1388   if ( this->m_MovingImageMask ) {
1389     // Check if mapped point is within the support region of the moving image
1390     // mask
1391     sampleOk = sampleOk && this->m_MovingImageMask->IsInside( mappedPoint );
1392   }
1393
1394
1395   if ( sampleOk ) {
1396     movingImageValue = this->m_Interpolator->Evaluate( mappedPoint );
1397
1398     if ( movingImageValue < m_MovingImageTrueMin ||
1399          movingImageValue > m_MovingImageTrueMax ) {
1400       // need to throw out this sample as it will not fall into a valid bin
1401       sampleOk = false;
1402     }
1403   }
1404 }
1405
1406
1407 /**
1408  * Compute PDF derivatives contribution for each parameter
1409  */
1410 template < class TFixedImage, class TMovingImage >
1411 void
1412 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1413 ::ComputePDFDerivatives(
1414   unsigned int sampleNumber,
1415   int pdfMovingIndex,
1416   const ImageDerivativesType& movingImageGradientValue,
1417   double cubicBSplineDerivativeValue ) const
1418 {
1419
1420   const int pdfFixedIndex =
1421     m_FixedImageSamples[sampleNumber].FixedImageParzenWindowIndex;
1422
1423   JointPDFValueType * derivPtr = NULL;
1424   double precomputedWeight = 0.0;
1425
1426   if( this->m_UseExplicitPDFDerivatives ) {
1427     // Update bins in the PDF derivatives for the current intensity pair
1428     derivPtr = m_JointPDFDerivatives->GetBufferPointer() +
1429                ( pdfFixedIndex  * m_JointPDFDerivatives->GetOffsetTable()[2] ) +
1430                ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1431   } else {
1432     // Recover the precomputed weight for this specific PDF bin
1433     precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex];
1434   }
1435
1436   if( !m_TransformIsBSpline ) {
1437
1438     /**
1439      * Generic version which works for all transforms.
1440      */
1441
1442     // Compute the transform Jacobian.
1443     typedef typename TransformType::JacobianType JacobianType;
1444 #if ITK_VERSION_MAJOR >= 4
1445       JacobianType jacobian;
1446       this->m_Transform->ComputeJacobianWithRespectToParameters( m_FixedImageSamples[sampleNumber].FixedImagePointValue, jacobian );
1447 #else
1448       const JacobianType & jacobian =
1449         this->m_Transform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue );
1450 #endif
1451
1452     for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) {
1453       double innerProduct = 0.0;
1454       for ( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1455         innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1456       }
1457
1458       const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1459
1460       if( this->m_UseExplicitPDFDerivatives ) {
1461         *(derivPtr) -= derivativeContribution;
1462         ++derivPtr;
1463       } else {
1464         this->m_MetricDerivative[mu] += precomputedWeight * derivativeContribution;
1465       }
1466     }
1467
1468   } else {
1469     const WeightsValueType * weights = NULL;
1470     const IndexValueType   * indices = NULL;
1471
1472     if( this->m_UseCachingOfBSplineWeights ) {
1473       /**
1474        * If the transform is of type BSplineDeformableTransform,
1475        * we can obtain a speed up by only processing the affected parameters.
1476        * Note that these pointers are just pointing to pre-allocated rows
1477        * of the caching arrays. There is therefore, no need to free this
1478        * memory.
1479        */
1480       weights = m_BSplineTransformWeightsArray[sampleNumber];
1481       indices = m_BSplineTransformIndicesArray[sampleNumber];
1482     } else {
1483 #if ITK_VERSION_MAJOR >= 4
1484       m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition(
1485         m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
1486 #else
1487       m_BSplineTransform->GetJacobian(
1488         m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
1489 #endif
1490     }
1491
1492     for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1493
1494       double innerProduct;
1495       int parameterIndex;
1496
1497       for( unsigned int mu = 0; mu < m_NumBSplineWeights; mu++ ) {
1498
1499         /* The array weights contains the Jacobian values in a 1-D array
1500          * (because for each parameter the Jacobian is non-zero in only 1 of the
1501          * possible dimensions) which is multiplied by the moving image
1502          * gradient. */
1503         if( this->m_UseCachingOfBSplineWeights ) {
1504           innerProduct = movingImageGradientValue[dim] * weights[mu];
1505           parameterIndex = indices[mu] + m_ParametersOffset[dim];
1506         } else {
1507           innerProduct = movingImageGradientValue[dim] * this->m_Weights[mu];
1508           parameterIndex = this->m_Indices[mu] + this->m_ParametersOffset[dim];
1509         }
1510
1511         const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1512
1513         if( this->m_UseExplicitPDFDerivatives ) {
1514           JointPDFValueType * ptr = derivPtr + parameterIndex;
1515           *(ptr) -= derivativeContribution;
1516         } else {
1517           this->m_MetricDerivative[parameterIndex] += precomputedWeight * derivativeContribution;
1518         }
1519
1520       } //end mu for loop
1521     } //end dim for loop
1522
1523   } // end if-block transform is BSpline
1524
1525 }
1526
1527
1528 // Method to reinitialize the seed of the random number generator
1529 template < class TFixedImage, class TMovingImage  > void
1530 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1531 ::ReinitializeSeed()
1532 {
1533   Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
1534 }
1535
1536 // Method to reinitialize the seed of the random number generator
1537 template < class TFixedImage, class TMovingImage  > void
1538 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1539 ::ReinitializeSeed(int seed)
1540 {
1541   Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(
1542     seed);
1543 }
1544
1545
1546 /**
1547  * Cache pre-transformed points, weights and indices.
1548  * This method is only called if the flag UseCachingOfBSplineWeights is ON.
1549  */
1550 template < class TFixedImage, class TMovingImage >
1551 void
1552 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1553 ::PreComputeTransformValues()
1554 {
1555   // Create all zero dummy transform parameters
1556   ParametersType dummyParameters( this->m_Transform->GetNumberOfParameters() );
1557   dummyParameters.Fill( 0.0 );
1558   this->m_Transform->SetParameters( dummyParameters );
1559
1560   // Cycle through each sampled fixed image point
1561   BSplineTransformWeightsType weights( m_NumBSplineWeights );
1562   BSplineTransformIndexArrayType indices( m_NumBSplineWeights );
1563   bool valid;
1564   MovingImagePointType mappedPoint;
1565
1566   // Declare iterators for iteration over the sample container
1567   typename FixedImageSpatialSampleContainer::const_iterator fiter;
1568   typename FixedImageSpatialSampleContainer::const_iterator fend =
1569     m_FixedImageSamples.end();
1570   unsigned long counter = 0;
1571
1572   for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ ) {
1573     m_BSplineTransform->TransformPoint(
1574       m_FixedImageSamples[counter].FixedImagePointValue,
1575       mappedPoint, weights, indices, valid );
1576
1577     for( unsigned long k = 0; k < m_NumBSplineWeights; k++ ) {
1578       m_BSplineTransformWeightsArray[counter][k] = weights[k];
1579       m_BSplineTransformIndicesArray[counter][k] = indices[k];
1580     }
1581
1582     m_PreTransformPointsArray[counter]      = mappedPoint;
1583     m_WithinSupportRegionArray[counter]     = valid;
1584
1585   }
1586
1587 }
1588
1589
1590 } // end namespace itk
1591
1592
1593 #endif
1594
1595 #endif