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