File: itkAvantsMutualInformationRegistrationFunction.h

package info (click to toggle)
ants 1.9.2%2Bsvn680.dfsg-4
  • links: PTS, VCS
  • area: main
  • in suites: wheezy
  • size: 12,136 kB
  • sloc: cpp: 41,966; sh: 2,545; perl: 216; makefile: 43
file content (703 lines) | stat: -rw-r--r-- 26,944 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
/*=========================================================================

  Program:   Advanced Normalization Tools
  Module:    $RCSfile: itkAvantsMutualInformationRegistrationFunction.h,v $
  Language:  C++
  Date:      $Date: 2009/01/08 15:14:48 $
  Version:   $Revision: 1.20 $

  Copyright (c) ConsortiumOfANTS. All rights reserved.
  See accompanying COPYING.txt or 
 http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkAvantsMutualInformationRegistrationFunction_h
#define __itkAvantsMutualInformationRegistrationFunction_h

#include "itkImageFileWriter.h"
#include "itkImageToImageMetric.h"
#include "itkAvantsPDEDeformableRegistrationFunction.h"
#include "itkCovariantVector.h"
#include "itkPoint.h"
#include "itkIndex.h"
#include "itkBSplineKernelFunction.h"
#include "itkBSplineDerivativeKernelFunction.h"
#include "itkCentralDifferenceImageFunction.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineDeformableTransform.h"
#include "itkTranslationTransform.h"
#include "itkArray2D.h"
#include "itkImageBase.h"
#include "itkTransform.h"
#include "itkInterpolateImageFunction.h"
#include "itkSingleValuedCostFunction.h"
#include "itkExceptionObject.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkSpatialObject.h"
#include "itkConstNeighborhoodIterator.h"

namespace itk
{

/** \class AvantsMutualInformationRegistrationFunction
 * \brief Computes the mutual information between two images to be 
 * registered using the method of Avants et al.
 *
 * AvantsMutualInformationRegistrationFunction computes the mutual 
 * information between a fixed and moving image to be registered.
 *
 * This class is templated over the FixedImage type and the MovingImage 
 * type.
 *
 * The fixed and moving images are set via methods SetFixedImage() and
 * SetMovingImage(). This metric makes use of user specified Transform and
 * Interpolator. The Transform is used to map points from the fixed image to
 * the moving image domain. The Interpolator is used to evaluate the image
 * intensity at user specified geometric points in the moving image.
 * The Transform and Interpolator are set via methods SetTransform() and
 * SetInterpolator().
 *
 * If a BSplineInterpolationFunction is used, this class obtain
 * image derivatives from the BSpline interpolator. Otherwise, 
 * image derivatives are computed using central differencing.
 *
 * \warning This metric assumes that the moving image has already been
 * connected to the interpolator outside of this class. 
 *
 * The method GetValue() computes of the mutual information
 * while method GetValueAndDerivative() computes
 * both the mutual information and its derivatives with respect to the
 * transform parameters.
 *
 * The calculations are based on the method of Avants et al [1,2]
 * where the probability density distribution are estimated using
 * Parzen histograms. Since the fixed image PDF does not contribute
 * to the derivatives, it does not need to be smooth. Hence, 
 * a zero order (box car) BSpline kernel is used
 * for the fixed image intensity PDF. On the other hand, to ensure
 * smoothness a third order BSpline kernel is used for the 
 * moving image intensity PDF.
 *
 * On Initialize(), the FixedImage is uniformly sampled within
 * the FixedImageRegion. The number of samples used can be set
 * via SetNumberOfSpatialSamples(). Typically, the number of
 * spatial samples used should increase with the image size.
 *
 * During each call of GetValue(), GetDerivatives(),
 * GetValueAndDerivatives(), marginal and joint intensity PDF's
 * values are estimated at discrete position or bins. 
 * The number of bins used can be set via SetNumberOfHistogramBins().
 * To handle data with arbitray magnitude and dynamic range, 
 * the image intensity is scale such that any contribution to the
 * histogram will fall into a valid bin.
 *
 * One the PDF's have been contructed, the mutual information
 * is obtained by doubling summing over the discrete PDF values.
 *
 *
 * Notes: 
 * 1. This class returns the negative mutual information value.
 * 2. This class in not thread safe due the private data structures
 *     used to the store the sampled points and the marginal and joint pdfs.
 *
 * References:
 * [1] "Nonrigid multimodality image registration"
 *      D. Avants, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
 *      Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620.
 * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations"
 *      D. Avants, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank
 *      IEEE Transactions in Medical Imaging. Vol.22, No.1, 
        January 2003. pp.120-128.
 * [3] "Optimization of Mutual Information for MultiResolution Image
 *      Registration"
 *      P. Thevenaz and M. Unser
 *      IEEE Transactions in Image Processing, 9(12) December 2000.
 *
 * \ingroup RegistrationMetrics
 * \ingroup ThreadUnSafe
 */
  template <class TFixedImage,class TMovingImage , class TDeformationField>
class ITK_EXPORT AvantsMutualInformationRegistrationFunction :
  public AvantsPDEDeformableRegistrationFunction< TFixedImage, TMovingImage , TDeformationField>
{
public:
  /** Standard class typedefs. */
  typedef AvantsMutualInformationRegistrationFunction    Self;
  typedef AvantsPDEDeformableRegistrationFunction< TFixedImage,
    TMovingImage, TDeformationField >    Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro( AvantsMutualInformationRegistrationFunction, 
    AvantsPDEDeformableRegistrationFunction );

  /** MovingImage image type. */
  typedef typename Superclass::MovingImageType     MovingImageType;
  typedef typename Superclass::MovingImagePointer  MovingImagePointer;

  /** FixedImage image type. */
  typedef typename Superclass::FixedImageType     FixedImageType;
  typedef typename Superclass::FixedImagePointer  FixedImagePointer;
  typedef typename FixedImageType::IndexType      IndexType;
  typedef typename FixedImageType::SizeType       SizeType;
  typedef typename FixedImageType::SpacingType    SpacingType;
  
  /** Deformation field type. */
  typedef typename Superclass::VectorType VectorType;
  typedef typename Superclass::DeformationFieldType    DeformationFieldType;
  typedef typename Superclass::DeformationFieldTypePointer   
    DeformationFieldTypePointer;

  /** Inherit some enums from the superclass. */
  itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);

  /** Inherit some enums from the superclass. */
  typedef typename Superclass::PixelType     PixelType;
  typedef typename Superclass::RadiusType    RadiusType;
  typedef typename Superclass::NeighborhoodType    NeighborhoodType;
  typedef typename Superclass::FloatOffsetType  FloatOffsetType;
  typedef typename Superclass::TimeStepType TimeStepType;


  /** Interpolator type. */
  typedef double CoordRepType;
  typedef    //       //    LinearInterpolateImageFunction<MovingImageType,CoordRepType>
    BSplineInterpolateImageFunction<MovingImageType,CoordRepType> 
    InterpolatorType;
  typedef typename InterpolatorType::Pointer         InterpolatorPointer;
  typedef typename InterpolatorType::PointType       PointType;
  typedef InterpolatorType DefaultInterpolatorType;
  //  typedef LinearInterpolateImageFunction<MovingImageType,CoordRepType>
  //DefaultInterpolatorType;

  /** Covariant vector type. */
  typedef CovariantVector<double,itkGetStaticConstMacro(ImageDimension)> CovariantVectorType;

  /** Gradient calculator type. */
  typedef CentralDifferenceImageFunction<FixedImageType> GradientCalculatorType;
  typedef typename GradientCalculatorType::Pointer   GradientCalculatorPointer;

  /** Set the moving image interpolator. */
  void SetMovingImageInterpolator( InterpolatorType * ptr )
  { m_MovingImageInterpolator = ptr; }
  
  /** Get the moving image interpolator. */
  InterpolatorType * GetMovingImageInterpolator(void)
    { return m_MovingImageInterpolator; }
  
  /** This class uses a constant timestep of 1. */
  virtual TimeStepType ComputeGlobalTimeStep(void *itkNotUsed(GlobalData)) const
    { return 1; }

  /** Return a pointer to a global data structure that is passed to
   * this object from the solver at each calculation.  */
  virtual void *GetGlobalDataPointer() const
  {
        GlobalDataStruct *global = new GlobalDataStruct();
    //    global->m_SumOfSquaredDifference  = 0.0;
    /// global->m_NumberOfPixelsProcessed = 0L;
    // global->m_SumOfSquaredChange      = 0;
       return global;
  }

  
  /** Release memory for global data structure. */
  virtual void ReleaseGlobalDataPointer( void *GlobalData ) const
  { 
     delete (GlobalDataStruct *) GlobalData;  
  }

  /** Set the object's state before each iteration. */
  virtual void InitializeIteration();


  typedef double CoordinateRepresentationType;

  /** Types inherited from Superclass. */
  typedef TranslationTransform<CoordinateRepresentationType, 
    //                    itkGetStaticConstMacro(ImageDimension),
                    itkGetStaticConstMacro(ImageDimension)> TransformType;

  typedef ImageToImageMetric< TFixedImage, TMovingImage > Metricclass;

  typedef typename TransformType::Pointer            TransformPointer;
  typedef typename Metricclass::TransformJacobianType    TransformJacobianType;
  //  typedef typename Metricclass::InterpolatorType         InterpolatorType;
  typedef typename Metricclass::MeasureType              MeasureType;
  typedef typename Metricclass::DerivativeType           DerivativeType;
  typedef typename TransformType::ParametersType           ParametersType;
  typedef typename Metricclass::FixedImageConstPointer   FixedImageConstPointer;
  typedef typename Metricclass::MovingImageConstPointer  MovingImageCosntPointer;
  // typedef typename TransformType::CoordinateRepresentationType  CoordinateRepresentationType;

  /** Index and Point typedef support. */
  typedef typename FixedImageType::IndexType            FixedImageIndexType;
  typedef typename FixedImageIndexType::IndexValueType  FixedImageIndexValueType;
  typedef typename MovingImageType::IndexType           MovingImageIndexType;
  typedef typename TransformType::InputPointType        FixedImagePointType;
  typedef typename TransformType::OutputPointType       MovingImagePointType;



  /** The marginal PDFs are stored as std::vector. */
  typedef float PDFValueType;
  //  typedef std::vector<PDFValueType> MarginalPDFType;
  typedef Image<PDFValueType,1> MarginalPDFType;
  typedef typename MarginalPDFType::IndexType MarginalPDFIndexType;
  /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */
  typedef Image<PDFValueType,2> JointPDFType;
  typedef Image<PDFValueType,3> JointPDFDerivativesType;
  typedef JointPDFType::IndexType                JointPDFIndexType;
  typedef JointPDFType::PixelType                JointPDFValueType;
  typedef JointPDFType::RegionType              JointPDFRegionType;
  typedef JointPDFType::SizeType                JointPDFSizeType;
  typedef JointPDFDerivativesType::IndexType    JointPDFDerivativesIndexType;
  typedef JointPDFDerivativesType::PixelType    JointPDFDerivativesValueType;
  typedef JointPDFDerivativesType::RegionType    JointPDFDerivativesRegionType;
  typedef JointPDFDerivativesType::SizeType      JointPDFDerivativesSizeType;



  /**  Get the value and derivatives for single valued optimizers. */
  double GetValueAndDerivative( IndexType index, 
                              MeasureType& Value, DerivativeType& Derivative1 , DerivativeType& Derivative2 ) ;


  /**  Get the value and derivatives for single valued optimizers. */
  double GetValueAndDerivativeInv( IndexType index, 
                              MeasureType& Value, DerivativeType& Derivative1 , DerivativeType& Derivative2 ) ;


  /** Number of spatial samples to used to compute metric */
  //  itkSetClampMacro( NumberOfSpatialSamples, unsigned long,
  //                1, NumericTraits<unsigned long>::max() );
  //itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long); 

  /** Number of bins to used in the histogram. Typical value is 50. */
  //  itkSetClampMacro( NumberOfHistogramBins, unsigned long,
  //                1, NumericTraits<unsigned long>::max() );
  // itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long);   
  void SetNumberOfHistogramBins(unsigned long nhb) { m_NumberOfHistogramBins=nhb+2*this->m_Padding;}
  unsigned long GetNumberOfHistogramBins() {return m_NumberOfHistogramBins;}

  void SetTransform(TransformPointer t){m_Transform=t;}
  TransformPointer GetTransform(){return m_Transform;}
  void SetInterpolator(InterpolatorPointer t){m_Interpolator=t;}
  InterpolatorPointer GetInterpolator(){ return m_Interpolator; }

  void GetProbabilities();


  double ComputeMutualInformation()
    {
     /** Ray Razlighi changes : May 3, 2010: improves convergence.
1- The padding is lowered to 2.
2- The MI energy is changed to actual MI, please note the value of this MI is in the range of [0  min(H(x),H(y))]. in case you need it to be normalized.
3- The Natural logarithm is changed to log2.
4- In the ComputeMutualInformation() the iterator range has been changed to cover the entire PDF.
5- In the FitIndexInBins() the truncation by floor has been modified to round.
6- The normalization is done based on NomberOfHistogramBins-1 instead of NomberOfHistogramBins. */
      float px,py,pxy;
      double mival = 0;
      double mi;
      unsigned long ct = 0;
      typename JointPDFType::IndexType index;
    
	  for (unsigned int ii=0; ii<m_NumberOfHistogramBins; ii++)
	  {
	  MarginalPDFIndexType mind;
	  mind[0]=ii;
	  px = m_FixedImageMarginalPDF->GetPixel(mind);
	  for (unsigned int jj=0; jj<m_NumberOfHistogramBins; jj++)
	    {
	      mind[0]=jj;	      
	      py = m_MovingImageMarginalPDF->GetPixel(mind);
	      float denom = px*py;
	      index[0]=ii;
	      index[1]=jj;
	      //pxy=m_JointPDF->GetPixel(index);
	      
	      JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
		( ii* m_NumberOfHistogramBins );
	      // Move the pointer to the first affected bin
	      int pdfMovingIndex = static_cast<int>( jj );
	      pdfPtr += pdfMovingIndex;
	      pxy=*(pdfPtr);
      
	      mi=0;
	      if (fabs(denom) > 0 ) 
		{
		  if (pxy/denom > 0) 
		    {
		      //true mi		      
		      mi = pxy*log(pxy/denom);
		      //test mi		      
		      //mi = 1.0 + log(pxy/denom);
		      ct++;
		    }
		    
		}
	      
	      mival += mi;
	    }
	//	  std::cout << " II " << ii << " JJ " << ii << " pxy " << pxy << " px " << px << std::endl;  

	}
	this->m_Energy = -1.0*mival/log(2);
	return this->m_Energy;
    }



  virtual VectorType ComputeUpdateInv(const NeighborhoodType &neighborhood,
                                   void *globalData,
                                   const FloatOffsetType &offset = FloatOffsetType(0.0))
  {
    VectorType update;
    update.Fill(0.0);   
    IndexType oindex = neighborhood.GetIndex();
    
    FixedImageType* img =const_cast<FixedImageType *>(this->Superclass::m_FixedImage.GetPointer()); 
    if (!img) return update;
    typename FixedImageType::SpacingType spacing=img->GetSpacing();
    typename FixedImageType::SizeType imagesize=img->GetLargestPossibleRegion().GetSize();
    
    for (unsigned int dd=0; dd<ImageDimension; dd++)
      {
	if ( oindex[dd] < 1 || 
	     oindex[dd] >= static_cast<typename IndexType::IndexValueType>(imagesize[dd]-1) ) 
	  return update;
      }    
    
    CovariantVectorType fixedGradient;
    double loce=0.0;
    ParametersType fdvec1(ImageDimension);
    ParametersType fdvec2(ImageDimension);    
    fdvec1.Fill(0);
    fdvec2.Fill(0);
    fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex( oindex ); 
    double nccm1=0;
    loce=this->GetValueAndDerivativeInv(oindex,nccm1,fdvec1,fdvec2); 
    float eps=10;
    if ( loce > eps ) loce=eps;
    if ( loce < eps*(-1.0) ) loce=eps*(-1.0);
    for (int imd=0; imd<ImageDimension; imd++) update[imd]=loce*fixedGradient[imd]*spacing[imd]*(1);
    if (this->m_MetricImage) this->m_MetricImage->SetPixel(oindex,loce);
    return update;

  }


  /*   Normalizing the image to the range of [0 1] */
  double GetMovingParzenTerm( double intensity )
  {
    double windowTerm = static_cast<double>( intensity ) -  this->m_MovingImageTrueMin;
    windowTerm = windowTerm / ( this->m_MovingImageTrueMax - this->m_MovingImageTrueMin   )  ;
    return windowTerm ; 
  }

  double GetFixedParzenTerm( double intensity )
  {
    double windowTerm = static_cast<double>( intensity ) -  this->m_FixedImageTrueMin;
    windowTerm = windowTerm / ( this->m_FixedImageTrueMax - this->m_FixedImageTrueMin   )  ;
    return  windowTerm ; 
  }

	
	
  /* find the image index in the number of histogram bins */	
	unsigned int FitIndexInBins( double windowTerm ) 
	{
		unsigned int movingImageParzenWindowIndex  = 
			static_cast<unsigned int>( this->m_Padding + round( windowTerm * (float)(this->m_NumberOfHistogramBins - 1 - this->m_Padding))) ;

		// Make sure the extreme values are in valid bins
		if ( movingImageParzenWindowIndex < this->m_Padding )
        {
			movingImageParzenWindowIndex = this->m_Padding -1 ;
        }
		else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - this->m_Padding ) )
        {
			movingImageParzenWindowIndex = m_NumberOfHistogramBins - this->m_Padding  - 1;
        }

		return movingImageParzenWindowIndex;
	}


  double FitContIndexInBins( double windowTerm ) {
    return ( (double) this->m_Padding + windowTerm*(float)(this->m_NumberOfHistogramBins-this->m_Padding));
  }


  virtual VectorType ComputeUpdate(const NeighborhoodType &neighborhood,
                                   void *globalData,
                                   const FloatOffsetType &offset = FloatOffsetType(0.0))
  {
    VectorType update;
    update.Fill(0.0);   
    IndexType oindex = neighborhood.GetIndex();
    
    FixedImageType* img =const_cast<FixedImageType *>(this->Superclass::m_MovingImage.GetPointer()); 
    if (!img) return update;
    typename FixedImageType::SpacingType spacing=img->GetSpacing();
    typename FixedImageType::SizeType imagesize=img->GetLargestPossibleRegion().GetSize();
    
    for (unsigned int dd=0; dd<ImageDimension; dd++)
      {
	if ( oindex[dd] < 1 || 
	     oindex[dd] >= static_cast<typename IndexType::IndexValueType>(imagesize[dd]-1) ) 
	  return update;
      }    
    
    CovariantVectorType movingGradient;
    double loce=0.0;
    ParametersType fdvec1(ImageDimension);
    ParametersType fdvec2(ImageDimension);    
    fdvec1.Fill(0);
    fdvec2.Fill(0);
    movingGradient = m_MovingImageGradientCalculator->EvaluateAtIndex( oindex ); 
	
    double nccm1=0;
    loce=this->GetValueAndDerivative(oindex,nccm1,fdvec1,fdvec2); 

    float eps=10;
    if ( loce > eps ) loce=eps;
    if ( loce < eps*(-1.0) ) loce=eps*(-1.0);

    for (int imd=0; imd<ImageDimension; imd++) update[imd]=loce*movingGradient[imd]*spacing[imd]*(1);
    return update;
  }


  void WriteImages()
  {
    
    if (this->m_MetricImage)
    {
      typedef ImageFileWriter<FixedImageType> writertype;
      typename writertype::Pointer w= writertype::New();
      w->SetInput(  this->m_MetricImage);
      w->SetFileName("ZZmetric.nii.gz");
      w->Write();
    }   

  }
  

  void SetOpticalFlow(bool b){ m_OpticalFlow = b; }

  typename JointPDFType::Pointer GetJointPDF() { return m_JointPDF; }
  typename JointPDFType::Pointer GetJointHist() { return m_JointHist; }


  void SetFixedImageMask( FixedImageType* img) {m_FixedImageMask=img; }



 /** FixedImage image neighborhood iterator type. */
  typedef ConstNeighborhoodIterator<FixedImageType> FixedImageNeighborhoodIteratorType;
  
  /** A global data type for this class of equation. Used to store
   * iterators for the fixed image. */
  struct GlobalDataStruct
   {
     FixedImageNeighborhoodIteratorType   m_FixedImageIterator;
   };



  /** The fixed image marginal PDF. */
  mutable MarginalPDFType::Pointer m_FixedImageMarginalPDF;

  /** The moving image marginal PDF. */
  mutable MarginalPDFType::Pointer m_MovingImageMarginalPDF;

 /** The joint PDF and PDF derivatives. */
  typename JointPDFType::Pointer m_JointPDF;

  unsigned long m_NumberOfSpatialSamples;
  unsigned long m_NumberOfParameters;

  /** Variables to define the marginal and joint histograms. */
  unsigned long m_NumberOfHistogramBins;
  double m_MovingImageNormalizedMin;
  double m_FixedImageNormalizedMin;
  double m_MovingImageTrueMin;
  double m_MovingImageTrueMax;
  double m_FixedImageTrueMin;
  double m_FixedImageTrueMax;
  double m_FixedImageBinSize;
  double m_MovingImageBinSize;

protected:

  AvantsMutualInformationRegistrationFunction();
  virtual ~AvantsMutualInformationRegistrationFunction() {};
  void PrintSelf(std::ostream& os, Indent indent) const;

 

private:

  AvantsMutualInformationRegistrationFunction(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented


   typename JointPDFType::Pointer m_JointHist;
  typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives;


  /** Typedefs for BSpline kernel and derivative functions. */
  typedef BSplineKernelFunction<3> CubicBSplineFunctionType;
  typedef BSplineDerivativeKernelFunction<3> 
  CubicBSplineDerivativeFunctionType;

  /** Cubic BSpline kernel for computing Parzen histograms. */
  typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel;
  typename CubicBSplineDerivativeFunctionType::Pointer 
  m_CubicBSplineDerivativeKernel;


  /**
   * Types and variables related to image derivative calculations.
   * If a BSplineInterpolationFunction is used, this class obtain
   * image derivatives from the BSpline interpolator. Otherwise, 
   * image derivatives are computed using central differencing.
   */
  typedef CovariantVector< double,
                           itkGetStaticConstMacro(ImageDimension) > ImageDerivativesType;


  /** Boolean to indicate if the interpolator BSpline. */
  bool m_InterpolatorIsBSpline;

  // boolean to determine if we use mono-modality assumption
  bool m_OpticalFlow;

  /** Typedefs for using BSpline interpolator. */
  typedef
  BSplineInterpolateImageFunction<MovingImageType,
                                  CoordinateRepresentationType> BSplineInterpolatorType;

  /** Pointer to BSplineInterpolator. */
  typename BSplineInterpolatorType::Pointer m_BSplineInterpolator;

  /** Typedefs for using central difference calculator. */
  typedef CentralDifferenceImageFunction<MovingImageType,
                                         CoordinateRepresentationType> DerivativeFunctionType;

  /** Pointer to central difference calculator. */
  typename DerivativeFunctionType::Pointer m_DerivativeCalculator;


  /**
   * Types and variables related to BSpline deformable transforms.
   * If the transform is of type third order BSplineDeformableTransform,
   * then we can speed up the metric derivative calculation by
   * only inspecting the parameters within the support region
   * of a mapped point.
   */

  /** Boolean to indicate if the transform is BSpline deformable. */
  bool m_TransformIsBSpline;

  /** The number of BSpline parameters per image dimension. */
  long m_NumParametersPerDim;

  /** 
   * The number of BSpline transform weights is the number of
   * of parameter in the support region (per dimension ). */   
  unsigned long m_NumBSplineWeights;


  /** 
   * Enum of the deformabtion field spline order. 
   */
  enum { DeformationSplineOrder = 3 };

  /**
   * Typedefs for the BSplineDeformableTransform.
   */
  typedef BSplineDeformableTransform<
    CoordinateRepresentationType,
    ::itk::GetImageDimension<FixedImageType>::ImageDimension,
    DeformationSplineOrder> BSplineTransformType;
  typedef typename BSplineTransformType::WeightsType 
  BSplineTransformWeightsType;
  typedef typename BSplineTransformType::ParameterIndexArrayType 
  BSplineTransformIndexArrayType;

  /**
   * Variables used when transform is of type BSpline deformable.
   */
  typename BSplineTransformType::Pointer m_BSplineTransform;

  /**
   * Cache pre-transformed points, weights, indices and 
   * within support region flag.
   */
  typedef typename BSplineTransformWeightsType::ValueType WeightsValueType;
  typedef          Array2D<WeightsValueType> BSplineTransformWeightsArrayType;
  typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType;
  typedef          Array2D<IndexValueType> BSplineTransformIndicesArrayType;
  typedef          std::vector<MovingImagePointType> MovingImagePointArrayType;
  typedef          std::vector<bool> BooleanArrayType;

  BSplineTransformWeightsArrayType      m_BSplineTransformWeightsArray;
  BSplineTransformIndicesArrayType      m_BSplineTransformIndicesArray;
  MovingImagePointArrayType             m_PreTransformPointsArray;
  BooleanArrayType                      m_WithinSupportRegionArray;
  
  typename TFixedImage::SpacingType                  m_FixedImageSpacing;
  typename TFixedImage::PointType                  m_FixedImageOrigin;

  typedef FixedArray<unsigned long, 
    ::itk::GetImageDimension<FixedImageType>::ImageDimension> ParametersOffsetType;
  ParametersOffsetType                  m_ParametersOffset;

  mutable TransformPointer    m_Transform;
  InterpolatorPointer         m_Interpolator;
  InterpolatorPointer             m_FixedImageInterpolator;
  InterpolatorPointer             m_MovingImageInterpolator;

  GradientCalculatorPointer       m_FixedImageGradientCalculator;
  GradientCalculatorPointer       m_MovingImageGradientCalculator;

  FixedImagePointer m_FixedImageMask;
  MovingImagePointer m_MovingImageMask;

  double m_NormalizeMetric;
  float m_Normalizer;

  typedef BSplineInterpolateImageFunction<JointPDFType,double> pdfintType;
  typename pdfintType::Pointer pdfinterpolator;

  typedef BSplineInterpolateImageFunction<JointPDFDerivativesType,double> dpdfintType;
  typename dpdfintType::Pointer dpdfinterpolator;

  typedef BSplineInterpolateImageFunction<MarginalPDFType,double> pdfintType2;
  typename pdfintType2::Pointer pdfinterpolator2;
  typename pdfintType2::Pointer pdfinterpolator3;

  unsigned int m_Padding;

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAvantsMutualInformationRegistrationFunction.cxx"
#endif

#endif