registration/itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx

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