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://oncora1.lyon.fnclcc.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 #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     const JacobianType& jacobian =
1445       this->m_Transform->GetJacobian(
1446         m_FixedImageSamples[sampleNumber].FixedImagePointValue );
1447
1448     for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) {
1449       double innerProduct = 0.0;
1450       for ( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1451         innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1452       }
1453
1454       const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1455
1456       if( this->m_UseExplicitPDFDerivatives ) {
1457         *(derivPtr) -= derivativeContribution;
1458         ++derivPtr;
1459       } else {
1460         this->m_MetricDerivative[mu] += precomputedWeight * derivativeContribution;
1461       }
1462     }
1463
1464   } else {
1465     const WeightsValueType * weights = NULL;
1466     const IndexValueType   * indices = NULL;
1467
1468     if( this->m_UseCachingOfBSplineWeights ) {
1469       /**
1470        * If the transform is of type BSplineDeformableTransform,
1471        * we can obtain a speed up by only processing the affected parameters.
1472        * Note that these pointers are just pointing to pre-allocated rows
1473        * of the caching arrays. There is therefore, no need to free this
1474        * memory.
1475        */
1476       weights = m_BSplineTransformWeightsArray[sampleNumber];
1477       indices = m_BSplineTransformIndicesArray[sampleNumber];
1478     } else {
1479       m_BSplineTransform->GetJacobian(
1480         m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
1481     }
1482
1483     for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1484
1485       double innerProduct;
1486       int parameterIndex;
1487
1488       for( unsigned int mu = 0; mu < m_NumBSplineWeights; mu++ ) {
1489
1490         /* The array weights contains the Jacobian values in a 1-D array
1491          * (because for each parameter the Jacobian is non-zero in only 1 of the
1492          * possible dimensions) which is multiplied by the moving image
1493          * gradient. */
1494         if( this->m_UseCachingOfBSplineWeights ) {
1495           innerProduct = movingImageGradientValue[dim] * weights[mu];
1496           parameterIndex = indices[mu] + m_ParametersOffset[dim];
1497         } else {
1498           innerProduct = movingImageGradientValue[dim] * this->m_Weights[mu];
1499           parameterIndex = this->m_Indices[mu] + this->m_ParametersOffset[dim];
1500         }
1501
1502         const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1503
1504         if( this->m_UseExplicitPDFDerivatives ) {
1505           JointPDFValueType * ptr = derivPtr + parameterIndex;
1506           *(ptr) -= derivativeContribution;
1507         } else {
1508           this->m_MetricDerivative[parameterIndex] += precomputedWeight * derivativeContribution;
1509         }
1510
1511       } //end mu for loop
1512     } //end dim for loop
1513
1514   } // end if-block transform is BSpline
1515
1516 }
1517
1518
1519 // Method to reinitialize the seed of the random number generator
1520 template < class TFixedImage, class TMovingImage  > void
1521 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1522 ::ReinitializeSeed()
1523 {
1524   Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
1525 }
1526
1527 // Method to reinitialize the seed of the random number generator
1528 template < class TFixedImage, class TMovingImage  > void
1529 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1530 ::ReinitializeSeed(int seed)
1531 {
1532   Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(
1533     seed);
1534 }
1535
1536
1537 /**
1538  * Cache pre-transformed points, weights and indices.
1539  * This method is only called if the flag UseCachingOfBSplineWeights is ON.
1540  */
1541 template < class TFixedImage, class TMovingImage >
1542 void
1543 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1544 ::PreComputeTransformValues()
1545 {
1546   // Create all zero dummy transform parameters
1547   ParametersType dummyParameters( this->m_Transform->GetNumberOfParameters() );
1548   dummyParameters.Fill( 0.0 );
1549   this->m_Transform->SetParameters( dummyParameters );
1550
1551   // Cycle through each sampled fixed image point
1552   BSplineTransformWeightsType weights( m_NumBSplineWeights );
1553   BSplineTransformIndexArrayType indices( m_NumBSplineWeights );
1554   bool valid;
1555   MovingImagePointType mappedPoint;
1556
1557   // Declare iterators for iteration over the sample container
1558   typename FixedImageSpatialSampleContainer::const_iterator fiter;
1559   typename FixedImageSpatialSampleContainer::const_iterator fend =
1560     m_FixedImageSamples.end();
1561   unsigned long counter = 0;
1562
1563   for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ ) {
1564     m_BSplineTransform->TransformPoint(
1565       m_FixedImageSamples[counter].FixedImagePointValue,
1566       mappedPoint, weights, indices, valid );
1567
1568     for( unsigned long k = 0; k < m_NumBSplineWeights; k++ ) {
1569       m_BSplineTransformWeightsArray[counter][k] = weights[k];
1570       m_BSplineTransformIndicesArray[counter][k] = indices[k];
1571     }
1572
1573     m_PreTransformPointsArray[counter]      = mappedPoint;
1574     m_WithinSupportRegionArray[counter]     = valid;
1575
1576   }
1577
1578 }
1579
1580
1581 } // end namespace itk
1582
1583
1584 #endif
1585
1586 #endif