1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
19 /*=========================================================================
21 Program: Insight Segmentation & Registration Toolkit
22 Module: $RCSfile: itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx,v $
24 Date: $Date: 2010/06/14 17:32:07 $
25 Version: $Revision: 1.1 $
27 Copyright (c) Insight Software Consortium. All rights reserved.
28 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
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.
34 =========================================================================*/
35 #ifndef __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
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"
47 #include "vnl/vnl_vector.txx"
48 #include "vnl/vnl_c_vector.txx"
56 template < class TFixedImage, class TMovingImage >
57 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
58 ::MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
60 m_NumberOfHistogramBins = 50;
62 this->SetComputeGradient(false); // don't use the default gradient for now
64 // Initialize PDFs to NULL
66 m_JointPDFDerivatives = NULL;
69 m_FixedImageMarginalPDF = NULL;
70 m_MovingImageMarginalPDF = NULL;
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;
79 m_CubicBSplineDerivativeKernel = NULL;
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;
91 this->m_ThreaderMetricDerivative = NULL;
92 this->m_UseExplicitPDFDerivatives = true;
93 this->m_ImplicitDerivativesSecondPass = false;
96 template < class TFixedImage, class TMovingImage >
97 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
98 ::~MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
101 if(m_FixedImageMarginalPDF != NULL) {
102 delete [] m_FixedImageMarginalPDF;
104 m_FixedImageMarginalPDF = NULL;
105 if(m_MovingImageMarginalPDF != NULL) {
106 delete [] m_MovingImageMarginalPDF;
108 m_MovingImageMarginalPDF = NULL;
109 if(m_ThreaderJointPDF != NULL) {
110 delete [] m_ThreaderJointPDF;
112 m_ThreaderJointPDF = NULL;
114 if(m_ThreaderJointPDFDerivatives != NULL) {
115 delete [] m_ThreaderJointPDFDerivatives;
117 m_ThreaderJointPDFDerivatives = NULL;
119 if(m_ThreaderFixedImageMarginalPDF != NULL) {
120 delete [] m_ThreaderFixedImageMarginalPDF;
122 m_ThreaderFixedImageMarginalPDF = NULL;
124 if(m_ThreaderJointPDFStartBin != NULL) {
125 delete [] m_ThreaderJointPDFStartBin;
127 m_ThreaderJointPDFStartBin = NULL;
129 if(m_ThreaderJointPDFEndBin != NULL) {
130 delete [] m_ThreaderJointPDFEndBin;
132 m_ThreaderJointPDFEndBin = NULL;
134 if(m_ThreaderJointPDFSum != NULL) {
135 delete [] m_ThreaderJointPDFSum;
137 m_ThreaderJointPDFSum = NULL;
139 if( this->m_ThreaderMetricDerivative != NULL ) {
140 delete [] this->m_ThreaderMetricDerivative;
142 this->m_ThreaderMetricDerivative = NULL;
146 * Print out internal information about this class
148 template < class TFixedImage, class TMovingImage >
150 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
151 ::PrintSelf(std::ostream& os, Indent indent) const
154 Superclass::PrintSelf(os, indent);
156 os << indent << "NumberOfHistogramBins: ";
157 os << this->m_NumberOfHistogramBins << std::endl;
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;
183 template <class TFixedImage, class TMovingImage>
185 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
186 ::Initialize(void) throw ( ExceptionObject )
188 this->Superclass::Initialize();
189 this->Superclass::MultiThreadingInitialize();
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();
198 m_FixedImageTrueMin = fixedImageStats->GetMinimum();
199 m_FixedImageTrueMax = fixedImageStats->GetMaximum();
200 double fixedImageMin = m_FixedImageTrueMin;
201 double fixedImageMax = m_FixedImageTrueMax;
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();
211 m_MovingImageTrueMin = movingImageStats->GetMinimum();
212 m_MovingImageTrueMax = movingImageStats->GetMaximum();
213 double movingImageMin = m_MovingImageTrueMin;
214 double movingImageMax = m_MovingImageTrueMax;
216 itkDebugMacro( " FixedImageMin: " << fixedImageMin <<
217 " FixedImageMax: " << fixedImageMax << std::endl );
218 itkDebugMacro( " MovingImageMin: " << movingImageMin <<
219 " MovingImageMax: " << movingImageMax << std::endl );
223 * Compute binsize for the histograms.
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.
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
238 const int padding = 2; // this will pad by 2 bins
240 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin )
241 / static_cast<double>( m_NumberOfHistogramBins
243 m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize
244 - static_cast<double>( padding );
246 m_MovingImageBinSize = ( movingImageMax - movingImageMin )
247 / static_cast<double>( m_NumberOfHistogramBins
249 m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize
250 - static_cast<double>( padding );
252 itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin );
253 itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin );
254 itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize );
255 itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize );
259 * Allocate memory for the marginal PDF and initialize values
260 * to zero. The marginal PDFs are stored as std::vector.
262 if(m_FixedImageMarginalPDF != NULL) {
263 delete [] m_FixedImageMarginalPDF;
265 m_FixedImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
266 if(m_MovingImageMarginalPDF != NULL) {
267 delete [] m_MovingImageMarginalPDF;
269 m_MovingImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
272 * Allocate memory for the joint PDF and joint PDF derivatives.
273 * The joint PDF and joint PDF derivatives are store as itk::Image.
275 m_JointPDF = JointPDFType::New();
276 m_JointPDFDerivatives = JointPDFDerivativesType::New();
278 // Instantiate a region, index, size
279 JointPDFRegionType jointPDFRegion;
280 JointPDFIndexType jointPDFIndex;
281 JointPDFSizeType jointPDFSize;
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 );
289 JointPDFDerivativesRegionType jointPDFDerivativesRegion;
292 // Now allocate memory according to the user-selected method.
294 if( this->m_UseExplicitPDFDerivatives ) {
295 this->m_JointPDFDerivatives = JointPDFDerivativesType::New();
297 JointPDFDerivativesIndexType jointPDFDerivativesIndex;
298 JointPDFDerivativesSizeType jointPDFDerivativesSize;
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,
305 jointPDFDerivativesIndex.Fill( 0 );
306 jointPDFDerivativesSize[0] = this->m_NumberOfParameters;
307 jointPDFDerivativesSize[1] = this->m_NumberOfHistogramBins;
308 jointPDFDerivativesSize[2] = this->m_NumberOfHistogramBins;
310 jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex );
311 jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize );
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);
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.
326 this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins );
327 this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() );
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 );
337 jointPDFRegion.SetIndex( jointPDFIndex );
338 jointPDFRegion.SetSize( jointPDFSize );
340 // Set the regions and allocate
341 m_JointPDF->SetRegions( jointPDFRegion );
342 m_JointPDF->Allocate();
344 m_JointPDFBufferSize = jointPDFSize[0] * jointPDFSize[1] * sizeof(PDFValueType);
348 * Setup the kernels used for the Parzen windows.
350 m_CubicBSplineKernel = CubicBSplineFunctionType::New();
351 m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New();
354 * Pre-compute the fixed image parzen window index for
355 * each point of the fixed image sample points list.
357 this->ComputeFixedImageParzenWindowIndices( this->m_FixedImageSamples );
360 if(m_ThreaderFixedImageMarginalPDF != NULL) {
361 delete [] m_ThreaderFixedImageMarginalPDF;
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];
368 if(m_ThreaderJointPDF != NULL) {
369 delete [] m_ThreaderJointPDF;
371 m_ThreaderJointPDF = new typename
372 JointPDFType::Pointer[this->m_NumberOfThreads-1];
374 if(m_ThreaderJointPDFStartBin != NULL) {
375 delete [] m_ThreaderJointPDFStartBin;
377 m_ThreaderJointPDFStartBin = new int[this->m_NumberOfThreads];
379 if(m_ThreaderJointPDFEndBin != NULL) {
380 delete [] m_ThreaderJointPDFEndBin;
382 m_ThreaderJointPDFEndBin = new int[this->m_NumberOfThreads];
384 if(m_ThreaderJointPDFSum != NULL) {
385 delete [] m_ThreaderJointPDFSum;
387 m_ThreaderJointPDFSum = new double[this->m_NumberOfThreads];
389 unsigned int threadID;
391 int binRange = m_NumberOfHistogramBins / this->m_NumberOfThreads;
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();
398 m_ThreaderJointPDFStartBin[threadID] = threadID * binRange;
399 m_ThreaderJointPDFEndBin[threadID] = (threadID + 1) * binRange - 1;
402 m_ThreaderJointPDFStartBin[this->m_NumberOfThreads-1] =
403 (this->m_NumberOfThreads - 1 ) * binRange;
405 m_ThreaderJointPDFEndBin[this->m_NumberOfThreads-1] = m_NumberOfHistogramBins - 1;
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;
413 m_ThreaderJointPDFDerivatives = NULL;
415 if(m_ThreaderMetricDerivative != NULL) {
416 delete [] m_ThreaderMetricDerivative;
418 m_ThreaderMetricDerivative = NULL;
421 if( this->m_UseExplicitPDFDerivatives ) {
422 m_ThreaderJointPDFDerivatives = new typename
423 JointPDFDerivativesType::Pointer[this->m_NumberOfThreads-1];
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();
432 m_ThreaderMetricDerivative = new DerivativeType[this->m_NumberOfThreads-1];
434 for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
435 this->m_ThreaderMetricDerivative[threadID] = DerivativeType( this->GetNumberOfParameters() );
441 * Uniformly sample the fixed image domain using a random walk
443 template < class TFixedImage, class TMovingImage >
445 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
446 ::ComputeFixedImageParzenWindowIndices(
447 FixedImageSampleContainer& samples )
450 typename FixedImageSampleContainer::iterator iter;
451 typename FixedImageSampleContainer::const_iterator end=samples.end();
453 for( iter=samples.begin(); iter != end; ++iter ) {
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 );
461 // Make sure the extreme values are in valid bins
464 } else if ( pindex > (m_NumberOfHistogramBins - 3) ) {
465 pindex = m_NumberOfHistogramBins - 3;
468 (*iter).valueIndex = pindex;
473 template < class TFixedImage, class TMovingImage >
475 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
476 ::GetValueThreadPreProcess( unsigned int threadID,
477 bool withinSampleThread ) const
480 this->Superclass::GetValueThreadPreProcess( threadID, withinSampleThread );
483 memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
485 m_JointPDFBufferSize );
486 memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
487 *m_NumberOfHistogramBins]),
489 m_NumberOfHistogramBins*sizeof(PDFValueType) );
491 // zero-th thread uses the variables directly
492 memset( m_JointPDF->GetBufferPointer(),
494 m_JointPDFBufferSize );
495 memset( m_FixedImageMarginalPDF,
497 m_NumberOfHistogramBins*sizeof(PDFValueType) );
501 template < class TFixedImage, class TMovingImage >
503 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
504 ::GetValueThreadProcessSample( unsigned int threadID,
505 unsigned long fixedImageSample,
506 const MovingImagePointType & itkNotUsed(mappedPoint),
507 double movingImageValue) const
510 * Compute this sample's contribution to the marginal and
511 * joint distributions.
515 if(movingImageValue < m_MovingImageTrueMin) {
517 } else if(movingImageValue > m_MovingImageTrueMax) {
521 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
522 double movingImageParzenWindowTerm = movingImageValue
523 / m_MovingImageBinSize
524 - m_MovingImageNormalizedMin;
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;
534 unsigned int fixedImageParzenWindowIndex =
535 this->m_FixedImageSamples[fixedImageSample].valueIndex;
537 m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
538 + fixedImageParzenWindowIndex] += 1;
540 m_FixedImageMarginalPDF[fixedImageParzenWindowIndex] += 1;
543 // Pointer to affected bin to be updated
544 JointPDFValueType *pdfPtr;
546 pdfPtr = m_ThreaderJointPDF[threadID-1]->GetBufferPointer() +
547 ( fixedImageParzenWindowIndex
548 * m_ThreaderJointPDF[threadID-1]
549 ->GetOffsetTable()[1] );
551 pdfPtr = m_JointPDF->GetBufferPointer() +
552 ( fixedImageParzenWindowIndex
553 * m_JointPDF->GetOffsetTable()[1] );
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;
561 double movingImageParzenWindowArg =
562 static_cast<double>( pdfMovingIndex )
563 - movingImageParzenWindowTerm;
565 while( pdfMovingIndex <= pdfMovingIndexMax ) {
566 *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
568 movingImageParzenWindowArg ) );
569 movingImageParzenWindowArg += 1;
576 template < class TFixedImage, class TMovingImage >
578 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
579 ::GetValueThreadPostProcess( unsigned int threadID,
580 bool itkNotUsed(withinSampleThread) ) const
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++);
606 for(i = m_ThreaderJointPDFStartBin[threadID];
607 i <= m_ThreaderJointPDFEndBin[threadID];
609 m_FixedImageMarginalPDF[i] += m_ThreaderFixedImageMarginalPDF[
610 (t*m_NumberOfHistogramBins) + i];
613 double jointPDFSum = 0.0;
614 pdfPtr = pdfPtrStart;
615 for(i = 0; i < maxI; i++) {
616 jointPDFSum += *(pdfPtr++);
619 m_ThreaderJointPDFSum[threadID-1] = jointPDFSum;
621 m_JointPDFSum = jointPDFSum;
625 template < class TFixedImage, class TMovingImage >
626 typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
628 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
629 ::GetValue( const ParametersType & parameters ) const
631 // Set up the parameters in the transform
632 this->m_Transform->SetParameters( parameters );
633 this->m_Parameters = parameters;
635 // MUST BE CALLED TO INITIATE PROCESSING
636 this->GetValueMultiThreadedInitiate();
638 // MUST BE CALLED TO INITIATE PROCESSING
639 this->GetValueMultiThreadedPostProcessInitiate();
641 for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
642 m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
644 if ( m_JointPDFSum == 0.0 ) {
645 itkExceptionMacro( "Joint PDF summed to zero" );
648 memset( m_MovingImageMarginalPDF,
650 m_NumberOfHistogramBins*sizeof(PDFValueType) );
652 JointPDFValueType * pdfPtr;
653 PDFValueType * movingMarginalPtr;
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++);
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
675 // Normalize the fixed image marginal PDF
676 if ( fixedPDFSum == 0.0 ) {
677 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
679 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
680 m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
684 * Compute the metric by double summation over histogram.
687 // Setup pointer to point to the first bin
688 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
692 for( unsigned int fixedIndex = 0;
693 fixedIndex < m_NumberOfHistogramBins;
695 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
697 for( unsigned int movingIndex = 0;
698 movingIndex < m_NumberOfHistogramBins;
699 ++movingIndex, jointPDFPtr++ ) {
700 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
701 double jointPDFValue = *(jointPDFPtr);
703 // check for non-zero bin contribution
704 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
706 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
707 if( fixedImagePDFValue > 1e-16) {
708 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
711 } // end if-block to check non-zero bin contribution
712 } // end for-loop over moving index
713 } // end for-loop over fixed index
715 return static_cast<MeasureType>( -1.0 * sum );
720 template < class TFixedImage, class TMovingImage >
722 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
723 ::GetValueAndDerivativeThreadPreProcess( unsigned int threadID,
724 bool itkNotUsed(withinSampleThread) ) const
727 memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
729 m_JointPDFBufferSize );
731 memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
732 * m_NumberOfHistogramBins]),
734 m_NumberOfHistogramBins*sizeof(PDFValueType) );
736 if( this->m_UseExplicitPDFDerivatives ) {
737 memset( m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer(),
739 m_JointPDFDerivativesBufferSize );
742 memset( m_JointPDF->GetBufferPointer(),
744 m_JointPDFBufferSize );
745 memset( m_FixedImageMarginalPDF,
747 m_NumberOfHistogramBins*sizeof(PDFValueType) );
749 if( this->m_UseExplicitPDFDerivatives ) {
750 memset( m_JointPDFDerivatives->GetBufferPointer(),
752 m_JointPDFDerivativesBufferSize );
757 template < class TFixedImage, class TMovingImage >
759 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
760 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
761 unsigned long fixedImageSample,
762 const MovingImagePointType & itkNotUsed(mappedPoint),
763 double movingImageValue,
764 const ImageDerivativesType &
765 movingImageGradientValue) const
768 * Compute this sample's contribution to the marginal
769 * and joint distributions.
772 if(movingImageValue < m_MovingImageTrueMin) {
774 } else if(movingImageValue > m_MovingImageTrueMax) {
778 unsigned int fixedImageParzenWindowIndex =
779 this->m_FixedImageSamples[fixedImageSample].valueIndex;
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 );
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;
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.
799 ++m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
800 + fixedImageParzenWindowIndex];
802 ++m_FixedImageMarginalPDF[fixedImageParzenWindowIndex];
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.
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.
818 // Pointer to affected bin to be updated
819 JointPDFValueType *pdfPtr;
821 pdfPtr = m_ThreaderJointPDF[threadID-1]
822 ->GetBufferPointer() +
823 ( fixedImageParzenWindowIndex
824 * m_NumberOfHistogramBins );
826 pdfPtr = m_JointPDF->GetBufferPointer() +
827 ( fixedImageParzenWindowIndex
828 * m_NumberOfHistogramBins );
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;
836 double movingImageParzenWindowArg = static_cast<double>( pdfMovingIndex )
837 - static_cast<double>( movingImageParzenWindowTerm );
839 while( pdfMovingIndex <= pdfMovingIndexMax ) {
840 *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
842 movingImageParzenWindowArg ) );
844 if( this->m_UseExplicitPDFDerivatives || this->m_ImplicitDerivativesSecondPass ) {
845 // Compute the cubicBSplineDerivative for later repeated use.
846 double cubicBSplineDerivativeValue =
847 m_CubicBSplineDerivativeKernel->Evaluate( movingImageParzenWindowArg );
849 // Compute PDF derivative contribution.
850 this->ComputePDFDerivatives( threadID,
853 movingImageGradientValue,
854 cubicBSplineDerivativeValue );
857 movingImageParzenWindowArg += 1;
864 template < class TFixedImage, class TMovingImage >
866 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
867 ::GetValueAndDerivativeThreadPostProcess( unsigned int threadID,
868 bool withinSampleThread ) const
870 this->GetValueThreadPostProcess( threadID, withinSampleThread );
872 if( this->m_UseExplicitPDFDerivatives ) {
873 const unsigned int rowSize = this->m_NumberOfParameters * m_NumberOfHistogramBins;
875 const unsigned int maxI =
876 rowSize * ( m_ThreaderJointPDFEndBin[threadID]
877 - m_ThreaderJointPDFStartBin[threadID] + 1 );
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()
891 tPdfDPtrEnd = tPdfDPtr + maxI;
892 // for(i = 0; i < maxI; i++)
893 while(tPdfDPtr < tPdfDPtrEnd) {
894 *(pdfDPtr++) += *(tPdfDPtr++);
898 double nFactor = 1.0 / (m_MovingImageBinSize
899 * this->m_NumberOfPixelsCounted);
901 pdfDPtr = pdfDPtrStart;
902 tPdfDPtrEnd = pdfDPtrStart + maxI;
903 //for(int i = 0; i < maxI; i++)
904 while(pdfDPtr < tPdfDPtrEnd) {
905 *(pdfDPtr++) *= nFactor;
911 * Get the both Value and Derivative Measure
913 template < class TFixedImage, class TMovingImage >
915 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
916 ::GetValueAndDerivative( const ParametersType & parameters,
918 DerivativeType & derivative) const
920 // Set output values to zero
921 value = NumericTraits< MeasureType >::Zero;
923 if( this->m_UseExplicitPDFDerivatives ) {
924 // Set output values to zero
925 if(derivative.GetSize() != this->m_NumberOfParameters) {
926 derivative = DerivativeType( this->m_NumberOfParameters );
928 memset( derivative.data_block(),
930 this->m_NumberOfParameters * sizeof(double) );
932 this->m_PRatioArray.Fill( 0.0 );
933 this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero );
934 for(unsigned int threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++ ) {
935 this->m_ThreaderMetricDerivative[threadID].Fill( NumericTraits< MeasureType >::Zero );
937 this->m_ImplicitDerivativesSecondPass = false;
940 // Set up the parameters in the transform
941 this->m_Transform->SetParameters( parameters );
942 this->m_Parameters = parameters;
944 // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
945 this->GetValueAndDerivativeMultiThreadedInitiate();
947 // CALL IF DOING THREADED POST PROCESSING
948 this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
950 for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
951 m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
953 if ( m_JointPDFSum == 0.0 ) {
954 itkExceptionMacro( "Joint PDF summed to zero" );
957 memset( m_MovingImageMarginalPDF,
959 m_NumberOfHistogramBins*sizeof(PDFValueType) );
961 JointPDFValueType * pdfPtr;
962 PDFValueType * movingMarginalPtr;
964 double fixedPDFSum = 0.0;
965 const double normalizationFactor = 1.0 / m_JointPDFSum;
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++);
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
985 // Normalize the fixed image marginal PDF
986 if ( fixedPDFSum == 0.0 ) {
987 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
989 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
990 m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
994 * Compute the metric by double summation over histogram.
997 // Setup pointer to point to the first bin
998 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
1000 // Initialize sum to zero
1003 const double nFactor = 1.0 / (m_MovingImageBinSize
1004 * this->m_NumberOfPixelsCounted);
1006 for( unsigned int fixedIndex = 0;
1007 fixedIndex < m_NumberOfHistogramBins;
1009 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
1011 for( unsigned int movingIndex = 0;
1012 movingIndex < m_NumberOfHistogramBins;
1013 ++movingIndex, jointPDFPtr++ ) {
1014 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
1015 double jointPDFValue = *(jointPDFPtr);
1017 // check for non-zero bin contribution
1018 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
1020 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
1022 if( fixedImagePDFValue > 1e-16) {
1023 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
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] );
1032 for( unsigned int parameter=0; parameter < this->m_NumberOfParameters; ++parameter, derivPtr++ ) {
1034 // Ref: eqn 23 of Thevenaz & Unser paper [3]
1035 derivative[parameter] -= (*derivPtr) * pRatio;
1037 } // end for-loop over parameters
1039 this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor;
1041 } // end if-block to check non-zero bin contribution
1042 } // end for-loop over moving index
1043 } // end for-loop over fixed index
1045 if( !(this->m_UseExplicitPDFDerivatives ) ) {
1046 // Second pass: This one is done for accumulating the contributions
1047 // to the derivative array.
1049 this->m_ImplicitDerivativesSecondPass = true;
1051 // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
1052 this->GetValueAndDerivativeMultiThreadedInitiate();
1054 // CALL IF DOING THREADED POST PROCESSING
1055 this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
1057 // Consolidate the contributions from each one of the threads to the total
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];
1066 derivative = this->m_MetricDerivative;
1069 value = static_cast<MeasureType>( -1.0 * sum );
1074 * Get the match measure derivative
1076 template < class TFixedImage, class TMovingImage >
1078 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1079 ::GetDerivative( const ParametersType & parameters,
1080 DerivativeType & derivative ) const
1083 // call the combined version
1084 this->GetValueAndDerivative( parameters, value, derivative );
1089 * Compute PDF derivatives contribution for each parameter
1091 template < class TFixedImage, class TMovingImage >
1093 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1094 ::ComputePDFDerivatives( unsigned int threadID,
1095 unsigned int sampleNumber,
1097 const ImageDerivativesType & movingImageGradientValue,
1098 double cubicBSplineDerivativeValue ) const
1100 // Update bins in the PDF derivatives for the current intensity pair
1101 // Could pre-compute
1102 JointPDFDerivativesValueType * derivPtr;
1104 double precomputedWeight = 0.0;
1106 const int pdfFixedIndex = this->m_FixedImageSamples[sampleNumber].valueIndex;
1108 DerivativeType * derivativeHelperArray = NULL;
1110 if( this->m_UseExplicitPDFDerivatives ) {
1112 derivPtr = m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer()
1113 + ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1114 + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1116 derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1117 + ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1118 + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1122 // Recover the precomputed weight for this specific PDF bin
1123 precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex];
1125 derivativeHelperArray = &(this->m_ThreaderMetricDerivative[threadID-1]);
1127 derivativeHelperArray = &(this->m_MetricDerivative);
1132 if( !this->m_TransformIsBSpline ) {
1134 * Generic version which works for all transforms.
1137 // Compute the transform Jacobian.
1138 // Should pre-compute
1139 typedef typename TransformType::JacobianType JacobianType;
1141 // Need to use one of the threader transforms if we're
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;
1150 transform = this->m_ThreaderTransform[threadID - 1];
1152 transform = this->m_Transform;
1155 const JacobianType& jacobian =
1156 transform->GetJacobian( this->m_FixedImageSamples[sampleNumber].point );
1158 // for ( unsigned int mu = 0; mu < this->m_NumberOfParameters; mu++ )
1160 // double innerProduct = 0.0;
1161 // for ( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ )
1163 // innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1166 // const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1168 // if( this->m_UseExplicitPDFDerivatives )
1170 // *(derivPtr) -= derivativeContribution;
1175 // (*derivativeHelperArray)[mu] += precomputedWeight * derivativeContribution;
1180 unsigned int mu, dim;
1181 double innerProduct,derivativeContribution;
1182 for ( mu = 0; mu < this->m_NumberOfParameters; mu+=3 ) {
1183 // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
1184 if (jacobian[0][mu])
1185 for ( dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1186 innerProduct = jacobian[dim][mu+dim] * movingImageGradientValue[dim];
1187 derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1189 if( this->m_UseExplicitPDFDerivatives ) {
1190 *(derivPtr) -= derivativeContribution;
1193 (*derivativeHelperArray)[mu+dim] += precomputedWeight * derivativeContribution;
1196 else if( this->m_UseExplicitPDFDerivatives ) derivPtr+=3;
1199 const WeightsValueType * weights = NULL;
1200 const IndexValueType * indices = NULL;
1203 BSplineTransformWeightsType * weightsHelper = NULL;
1204 BSplineTransformIndexArrayType * indicesHelper = NULL;
1206 if( this->m_UseCachingOfBSplineWeights ) {
1208 // If the transform is of type BSplineDeformableTransform, we can obtain
1209 // a speed up by only processing the affected parameters. Note that
1210 // these pointers are just pointing to pre-allocated rows of the caching
1211 // arrays. There is therefore, no need to free this memory.
1213 weights = this->m_BSplineTransformWeightsArray[sampleNumber];
1214 indices = this->m_BSplineTransformIndicesArray[sampleNumber];
1216 if( threadID > 0 ) {
1217 weightsHelper = &(this->m_ThreaderBSplineTransformWeights[threadID-1]);
1218 indicesHelper = &(this->m_ThreaderBSplineTransformIndices[threadID-1]);
1220 weightsHelper = &(this->m_BSplineTransformWeights);
1221 indicesHelper = &(this->m_BSplineTransformIndices);
1224 this->m_BSplineTransform->GetJacobian(
1225 this->m_FixedImageSamples[sampleNumber].point,
1226 *weightsHelper, *indicesHelper );
1229 for( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1231 double innerProduct;
1234 for( unsigned int mu = 0; mu < this->m_NumBSplineWeights; mu++ ) {
1236 /* The array weights contains the Jacobian values in a 1-D array
1237 * (because for each parameter the Jacobian is non-zero in only 1 of the
1238 * possible dimensions) which is multiplied by the moving image
1240 if( this->m_UseCachingOfBSplineWeights ) {
1241 innerProduct = movingImageGradientValue[dim] * weights[mu];
1242 parameterIndex = indices[mu] + this->m_BSplineParametersOffset[dim];
1244 innerProduct = movingImageGradientValue[dim] * (*weightsHelper)[mu];
1245 parameterIndex = (*indicesHelper)[mu] + this->m_BSplineParametersOffset[dim];
1248 const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1250 if( this->m_UseExplicitPDFDerivatives ) {
1251 JointPDFValueType * ptr = derivPtr + parameterIndex;
1252 *(ptr) -= derivativeContribution;
1254 (*derivativeHelperArray)[parameterIndex] += precomputedWeight * derivativeContribution;
1258 } //end dim for loop
1260 } // end if-block transform is BSpline
1265 } // end namespace itk