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