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