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://www.centreleonberard.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.h"
48 #include "vnl/vnl_c_vector.h"
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 #if ( ( ITK_VERSION_MAJOR == 4 ) && ( ITK_VERSION_MINOR > 12 ) || ( ITK_VERSION_MAJOR > 4 ))
189 ::Initialize(void) throw ( ExceptionObject )
192 this->Superclass::Initialize();
193 this->Superclass::MultiThreadingInitialize();
195 typedef StatisticsImageFilter<FixedImageType> FixedImageStatisticsFilterType;
196 typename FixedImageStatisticsFilterType::Pointer fixedImageStats =
197 FixedImageStatisticsFilterType::New();
198 fixedImageStats->SetInput( this->m_FixedImage );
199 #if ITK_VERSION_MAJOR <= 4
200 fixedImageStats->SetNumberOfThreads( this->m_NumberOfThreads );
202 fixedImageStats->SetNumberOfWorkUnits( this->m_NumberOfWorkUnits );
204 fixedImageStats->Update();
206 m_FixedImageTrueMin = fixedImageStats->GetMinimum();
207 m_FixedImageTrueMax = fixedImageStats->GetMaximum();
208 double fixedImageMin = m_FixedImageTrueMin;
209 double fixedImageMax = m_FixedImageTrueMax;
211 typedef StatisticsImageFilter<MovingImageType>
212 MovingImageStatisticsFilterType;
213 typename MovingImageStatisticsFilterType::Pointer movingImageStats =
214 MovingImageStatisticsFilterType::New();
215 movingImageStats->SetInput( this->m_MovingImage );
216 #if ITK_VERSION_MAJOR <= 4
217 movingImageStats->SetNumberOfThreads( this->m_NumberOfThreads );
219 movingImageStats->SetNumberOfWorkUnits( this->m_NumberOfWorkUnits );
221 movingImageStats->Update();
223 m_MovingImageTrueMin = movingImageStats->GetMinimum();
224 m_MovingImageTrueMax = movingImageStats->GetMaximum();
225 double movingImageMin = m_MovingImageTrueMin;
226 double movingImageMax = m_MovingImageTrueMax;
228 itkDebugMacro( " FixedImageMin: " << fixedImageMin <<
229 " FixedImageMax: " << fixedImageMax << std::endl );
230 itkDebugMacro( " MovingImageMin: " << movingImageMin <<
231 " MovingImageMax: " << movingImageMax << std::endl );
235 * Compute binsize for the histograms.
237 * The binsize for the image intensities needs to be adjusted so that
238 * we can avoid dealing with boundary conditions using the cubic
239 * spline as the Parzen window. We do this by increasing the size
240 * of the bins so that the joint histogram becomes "padded" at the
241 * borders. Because we are changing the binsize,
242 * we also need to shift the minimum by the padded amount in order to
243 * avoid minimum values filling in our padded region.
245 * Note that there can still be non-zero bin values in the padded region,
246 * it's just that these bins will never be a central bin for the Parzen
250 const int padding = 2; // this will pad by 2 bins
252 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin )
253 / static_cast<double>( m_NumberOfHistogramBins
255 m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize
256 - static_cast<double>( padding );
258 m_MovingImageBinSize = ( movingImageMax - movingImageMin )
259 / static_cast<double>( m_NumberOfHistogramBins
261 m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize
262 - static_cast<double>( padding );
264 itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin );
265 itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin );
266 itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize );
267 itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize );
271 * Allocate memory for the marginal PDF and initialize values
272 * to zero. The marginal PDFs are stored as std::vector.
274 if(m_FixedImageMarginalPDF != NULL) {
275 delete [] m_FixedImageMarginalPDF;
277 m_FixedImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
278 if(m_MovingImageMarginalPDF != NULL) {
279 delete [] m_MovingImageMarginalPDF;
281 m_MovingImageMarginalPDF = new PDFValueType[m_NumberOfHistogramBins];
284 * Allocate memory for the joint PDF and joint PDF derivatives.
285 * The joint PDF and joint PDF derivatives are store as itk::Image.
287 m_JointPDF = JointPDFType::New();
288 m_JointPDFDerivatives = JointPDFDerivativesType::New();
290 // Instantiate a region, index, size
291 JointPDFRegionType jointPDFRegion;
292 JointPDFIndexType jointPDFIndex;
293 JointPDFSizeType jointPDFSize;
295 // Deallocate the memory that may have been allocated for
296 // previous runs of the metric.
297 this->m_JointPDFDerivatives = NULL; // by destroying the dynamic array
298 this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones
299 this->m_MetricDerivative = DerivativeType( 1 );
301 JointPDFDerivativesRegionType jointPDFDerivativesRegion;
304 // Now allocate memory according to the user-selected method.
306 if( this->m_UseExplicitPDFDerivatives ) {
307 this->m_JointPDFDerivatives = JointPDFDerivativesType::New();
309 JointPDFDerivativesIndexType jointPDFDerivativesIndex;
310 JointPDFDerivativesSizeType jointPDFDerivativesSize;
312 // For the derivatives of the joint PDF define a region starting from {0,0,0}
313 // with size {m_NumberOfParameters,m_NumberOfHistogramBins,
314 // m_NumberOfHistogramBins}. The dimension represents transform parameters,
315 // fixed image parzen window index and moving image parzen window index,
317 jointPDFDerivativesIndex.Fill( 0 );
318 jointPDFDerivativesSize[0] = this->m_NumberOfParameters;
319 jointPDFDerivativesSize[1] = this->m_NumberOfHistogramBins;
320 jointPDFDerivativesSize[2] = this->m_NumberOfHistogramBins;
322 jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex );
323 jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize );
325 // Set the regions and allocate
326 m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion );
327 m_JointPDFDerivatives->Allocate();
328 m_JointPDFDerivativesBufferSize = jointPDFDerivativesSize[0]
329 * jointPDFDerivativesSize[1]
330 * jointPDFDerivativesSize[2]
331 * sizeof(JointPDFDerivativesValueType);
334 /** Allocate memory for helper array that will contain the pRatios
335 * for each bin of the joint histogram. This is part of the effort
336 * for flattening the computation of the PDF Jacobians.
338 this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins );
339 this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() );
342 // For the joint PDF define a region starting from {0,0}
343 // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}.
344 // The dimension represents fixed image parzen window index
345 // and moving image parzen window index, respectively.
346 jointPDFIndex.Fill( 0 );
347 jointPDFSize.Fill( m_NumberOfHistogramBins );
349 jointPDFRegion.SetIndex( jointPDFIndex );
350 jointPDFRegion.SetSize( jointPDFSize );
352 // Set the regions and allocate
353 m_JointPDF->SetRegions( jointPDFRegion );
354 m_JointPDF->Allocate();
356 m_JointPDFBufferSize = jointPDFSize[0] * jointPDFSize[1] * sizeof(PDFValueType);
360 * Setup the kernels used for the Parzen windows.
362 m_CubicBSplineKernel = CubicBSplineFunctionType::New();
363 m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New();
366 * Pre-compute the fixed image parzen window index for
367 * each point of the fixed image sample points list.
369 this->ComputeFixedImageParzenWindowIndices( this->m_FixedImageSamples );
372 if(m_ThreaderFixedImageMarginalPDF != NULL) {
373 delete [] m_ThreaderFixedImageMarginalPDF;
375 // Assumes number of threads doesn't change between calls to Initialize
376 m_ThreaderFixedImageMarginalPDF = new
377 #if ITK_VERSION_MAJOR <= 4
378 PDFValueType[(this->m_NumberOfThreads-1) * m_NumberOfHistogramBins];
380 PDFValueType[(this->m_NumberOfWorkUnits-1) * m_NumberOfHistogramBins];
383 if(m_ThreaderJointPDF != NULL) {
384 delete [] m_ThreaderJointPDF;
386 m_ThreaderJointPDF = new typename
387 #if ITK_VERSION_MAJOR <= 4
388 JointPDFType::Pointer[this->m_NumberOfThreads-1];
390 JointPDFType::Pointer[this->m_NumberOfWorkUnits-1];
393 if(m_ThreaderJointPDFStartBin != NULL) {
394 delete [] m_ThreaderJointPDFStartBin;
396 #if ITK_VERSION_MAJOR <= 4
397 m_ThreaderJointPDFStartBin = new int[this->m_NumberOfThreads];
399 m_ThreaderJointPDFStartBin = new int[this->m_NumberOfWorkUnits];
402 if(m_ThreaderJointPDFEndBin != NULL) {
403 delete [] m_ThreaderJointPDFEndBin;
405 #if ITK_VERSION_MAJOR <= 4
406 m_ThreaderJointPDFEndBin = new int[this->m_NumberOfThreads];
408 m_ThreaderJointPDFEndBin = new int[this->m_NumberOfWorkUnits];
411 if(m_ThreaderJointPDFSum != NULL) {
412 delete [] m_ThreaderJointPDFSum;
414 #if ITK_VERSION_MAJOR <= 4
415 m_ThreaderJointPDFSum = new double[this->m_NumberOfThreads];
417 m_ThreaderJointPDFSum = new double[this->m_NumberOfWorkUnits];
420 unsigned int threadID;
422 #if ITK_VERSION_MAJOR <= 4
423 int binRange = m_NumberOfHistogramBins / this->m_NumberOfThreads;
425 int binRange = m_NumberOfHistogramBins / this->m_NumberOfWorkUnits;
428 #if ITK_VERSION_MAJOR <= 4
429 for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
431 for(threadID = 0; threadID < this->m_NumberOfWorkUnits-1; threadID++) {
433 m_ThreaderJointPDF[threadID] = JointPDFType::New();
434 m_ThreaderJointPDF[threadID]->SetRegions( jointPDFRegion );
435 m_ThreaderJointPDF[threadID]->Allocate();
437 m_ThreaderJointPDFStartBin[threadID] = threadID * binRange;
438 m_ThreaderJointPDFEndBin[threadID] = (threadID + 1) * binRange - 1;
441 #if ITK_VERSION_MAJOR <= 4
442 m_ThreaderJointPDFStartBin[this->m_NumberOfThreads-1] = (this->m_NumberOfThreads - 1 ) * binRange;
443 m_ThreaderJointPDFEndBin[this->m_NumberOfThreads-1] = m_NumberOfHistogramBins - 1;
445 m_ThreaderJointPDFStartBin[this->m_NumberOfWorkUnits-1] = (this->m_NumberOfWorkUnits - 1 ) * binRange;
446 m_ThreaderJointPDFEndBin[this->m_NumberOfWorkUnits-1] = m_NumberOfHistogramBins - 1;
449 // Release memory of arrays that may have been used for
450 // previous executions of this metric with different settings
451 // of the memory caching flags.
452 if(m_ThreaderJointPDFDerivatives != NULL) {
453 delete [] m_ThreaderJointPDFDerivatives;
455 m_ThreaderJointPDFDerivatives = NULL;
457 if(m_ThreaderMetricDerivative != NULL) {
458 delete [] m_ThreaderMetricDerivative;
460 m_ThreaderMetricDerivative = NULL;
463 if( this->m_UseExplicitPDFDerivatives ) {
464 m_ThreaderJointPDFDerivatives = new typename
465 #if ITK_VERSION_MAJOR <= 4
466 JointPDFDerivativesType::Pointer[this->m_NumberOfThreads-1];
468 JointPDFDerivativesType::Pointer[this->m_NumberOfWorkUnits-1];
471 #if ITK_VERSION_MAJOR <= 4
472 for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
474 for(threadID = 0; threadID < this->m_NumberOfWorkUnits-1; threadID++) {
476 m_ThreaderJointPDFDerivatives[threadID] = JointPDFDerivativesType::New();
477 m_ThreaderJointPDFDerivatives[threadID]->SetRegions(
478 jointPDFDerivativesRegion );
479 m_ThreaderJointPDFDerivatives[threadID]->Allocate();
482 #if ITK_VERSION_MAJOR <= 4
483 m_ThreaderMetricDerivative = new DerivativeType[this->m_NumberOfThreads-1];
485 m_ThreaderMetricDerivative = new DerivativeType[this->m_NumberOfWorkUnits-1];
488 #if ITK_VERSION_MAJOR <= 4
489 for(threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++) {
491 for(threadID = 0; threadID < this->m_NumberOfWorkUnits-1; threadID++) {
493 this->m_ThreaderMetricDerivative[threadID] = DerivativeType( this->GetNumberOfParameters() );
499 * Uniformly sample the fixed image domain using a random walk
501 template < class TFixedImage, class TMovingImage >
503 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
504 ::ComputeFixedImageParzenWindowIndices(
505 FixedImageSampleContainer& samples )
508 typename FixedImageSampleContainer::iterator iter;
509 typename FixedImageSampleContainer::const_iterator end=samples.end();
511 for( iter=samples.begin(); iter != end; ++iter ) {
513 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
514 double windowTerm = static_cast<double>( (*iter).value )
515 / m_FixedImageBinSize
516 - m_FixedImageNormalizedMin;
517 unsigned int pindex = static_cast<unsigned int>( windowTerm );
519 // Make sure the extreme values are in valid bins
522 } else if ( pindex > (m_NumberOfHistogramBins - 3) ) {
523 pindex = m_NumberOfHistogramBins - 3;
526 (*iter).valueIndex = pindex;
531 template < class TFixedImage, class TMovingImage >
533 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
534 ::GetValueThreadPreProcess( unsigned int threadID,
535 bool withinSampleThread ) const
538 this->Superclass::GetValueThreadPreProcess( threadID, withinSampleThread );
541 memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
543 m_JointPDFBufferSize );
544 memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
545 *m_NumberOfHistogramBins]),
547 m_NumberOfHistogramBins*sizeof(PDFValueType) );
549 // zero-th thread uses the variables directly
550 memset( m_JointPDF->GetBufferPointer(),
552 m_JointPDFBufferSize );
553 memset( m_FixedImageMarginalPDF,
555 m_NumberOfHistogramBins*sizeof(PDFValueType) );
559 template < class TFixedImage, class TMovingImage >
561 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
562 ::GetValueThreadProcessSample( unsigned int threadID,
563 unsigned long fixedImageSample,
564 const MovingImagePointType & itkNotUsed(mappedPoint),
565 double movingImageValue) const
568 * Compute this sample's contribution to the marginal and
569 * joint distributions.
573 if(movingImageValue < m_MovingImageTrueMin) {
575 } else if(movingImageValue > m_MovingImageTrueMax) {
579 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
580 double movingImageParzenWindowTerm = movingImageValue
581 / m_MovingImageBinSize
582 - m_MovingImageNormalizedMin;
584 unsigned int movingImageParzenWindowIndex =
585 static_cast<unsigned int>( movingImageParzenWindowTerm );
586 if( movingImageParzenWindowIndex < 2 ) {
587 movingImageParzenWindowIndex = 2;
588 } else if( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
589 movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
592 unsigned int fixedImageParzenWindowIndex =
593 this->m_FixedImageSamples[fixedImageSample].valueIndex;
595 m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
596 + fixedImageParzenWindowIndex] += 1;
598 m_FixedImageMarginalPDF[fixedImageParzenWindowIndex] += 1;
601 // Pointer to affected bin to be updated
602 JointPDFValueType *pdfPtr;
604 pdfPtr = m_ThreaderJointPDF[threadID-1]->GetBufferPointer() +
605 ( fixedImageParzenWindowIndex
606 * m_ThreaderJointPDF[threadID-1]
607 ->GetOffsetTable()[1] );
609 pdfPtr = m_JointPDF->GetBufferPointer() +
610 ( fixedImageParzenWindowIndex
611 * m_JointPDF->GetOffsetTable()[1] );
614 // Move the pointer to the first affected bin
615 int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
616 pdfPtr += pdfMovingIndex;
617 int pdfMovingIndexMax = static_cast<int>(movingImageParzenWindowIndex) + 2;
619 double movingImageParzenWindowArg =
620 static_cast<double>( pdfMovingIndex )
621 - movingImageParzenWindowTerm;
623 while( pdfMovingIndex <= pdfMovingIndexMax ) {
624 *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
626 movingImageParzenWindowArg ) );
627 movingImageParzenWindowArg += 1;
634 template < class TFixedImage, class TMovingImage >
636 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
637 ::GetValueThreadPostProcess( unsigned int threadID,
638 bool itkNotUsed(withinSampleThread) ) const
643 maxI = m_NumberOfHistogramBins
644 * ( m_ThreaderJointPDFEndBin[threadID]
645 - m_ThreaderJointPDFStartBin[threadID] + 1);
646 JointPDFValueType *pdfPtr;
647 JointPDFValueType *pdfPtrStart;
648 pdfPtrStart = m_JointPDF->GetBufferPointer()
649 + ( m_ThreaderJointPDFStartBin[threadID]
650 * m_JointPDF->GetOffsetTable()[1] );
651 JointPDFValueType *tPdfPtr;
652 JointPDFValueType *tPdfPtrEnd;
653 unsigned int tPdfPtrOffset;
654 tPdfPtrOffset = ( m_ThreaderJointPDFStartBin[threadID]
655 * m_JointPDF->GetOffsetTable()[1] );
656 #if ITK_VERSION_MAJOR <= 4
657 for(t=0; t<this->m_NumberOfThreads-1; t++) {
659 for(t=0; t<this->m_NumberOfWorkUnits-1; t++) {
661 pdfPtr = pdfPtrStart;
662 tPdfPtr = m_ThreaderJointPDF[t]->GetBufferPointer() + tPdfPtrOffset;
663 tPdfPtrEnd = tPdfPtr + maxI;
664 //for(i=0; i < maxI; i++)
665 while(tPdfPtr < tPdfPtrEnd) {
666 *(pdfPtr++) += *(tPdfPtr++);
668 for(i = m_ThreaderJointPDFStartBin[threadID];
669 i <= m_ThreaderJointPDFEndBin[threadID];
671 m_FixedImageMarginalPDF[i] += m_ThreaderFixedImageMarginalPDF[
672 (t*m_NumberOfHistogramBins) + i];
675 double jointPDFSum = 0.0;
676 pdfPtr = pdfPtrStart;
677 for(i = 0; i < maxI; i++) {
678 jointPDFSum += *(pdfPtr++);
681 m_ThreaderJointPDFSum[threadID-1] = jointPDFSum;
683 m_JointPDFSum = jointPDFSum;
687 template < class TFixedImage, class TMovingImage >
688 typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
690 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
691 ::GetValue( const ParametersType & parameters ) const
693 // Set up the parameters in the transform
694 this->m_Transform->SetParameters( parameters );
696 // MUST BE CALLED TO INITIATE PROCESSING
697 this->GetValueMultiThreadedInitiate();
699 // MUST BE CALLED TO INITIATE PROCESSING
700 this->GetValueMultiThreadedPostProcessInitiate();
702 #if ITK_VERSION_MAJOR <= 4
703 for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
705 for(unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits-1; threadID++) {
707 m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
709 if ( m_JointPDFSum == 0.0 ) {
710 itkExceptionMacro( "Joint PDF summed to zero" );
713 memset( m_MovingImageMarginalPDF,
715 m_NumberOfHistogramBins*sizeof(PDFValueType) );
717 JointPDFValueType * pdfPtr;
718 PDFValueType * movingMarginalPtr;
720 double fixedPDFSum = 0.0;
721 double nFactor = 1.0 / m_JointPDFSum;
722 pdfPtr = m_JointPDF->GetBufferPointer();
723 for(i=0; i<m_NumberOfHistogramBins; i++) {
724 fixedPDFSum += m_FixedImageMarginalPDF[i];
725 movingMarginalPtr = m_MovingImageMarginalPDF;
726 for(j=0; j<m_NumberOfHistogramBins; j++) {
727 *(pdfPtr) *= nFactor;
728 *(movingMarginalPtr++) += *(pdfPtr++);
732 if( this->m_NumberOfPixelsCounted <
733 this->m_NumberOfFixedImageSamples / 16 ) {
734 itkExceptionMacro( "Too many samples map outside moving image buffer: "
735 << this->m_NumberOfPixelsCounted << " / "
736 << this->m_NumberOfFixedImageSamples
740 // Normalize the fixed image marginal PDF
741 if ( fixedPDFSum == 0.0 ) {
742 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
744 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
745 m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
749 * Compute the metric by double summation over histogram.
752 // Setup pointer to point to the first bin
753 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
757 for( unsigned int fixedIndex = 0;
758 fixedIndex < m_NumberOfHistogramBins;
760 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
762 for( unsigned int movingIndex = 0;
763 movingIndex < m_NumberOfHistogramBins;
764 ++movingIndex, jointPDFPtr++ ) {
765 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
766 double jointPDFValue = *(jointPDFPtr);
768 // check for non-zero bin contribution
769 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
771 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
772 if( fixedImagePDFValue > 1e-16) {
773 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
776 } // end if-block to check non-zero bin contribution
777 } // end for-loop over moving index
778 } // end for-loop over fixed index
780 return static_cast<MeasureType>( -1.0 * sum );
785 template < class TFixedImage, class TMovingImage >
787 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
788 ::GetValueAndDerivativeThreadPreProcess( unsigned int threadID,
789 bool itkNotUsed(withinSampleThread) ) const
792 memset( m_ThreaderJointPDF[threadID-1]->GetBufferPointer(),
794 m_JointPDFBufferSize );
796 memset( &(m_ThreaderFixedImageMarginalPDF[(threadID-1)
797 * m_NumberOfHistogramBins]),
799 m_NumberOfHistogramBins*sizeof(PDFValueType) );
801 if( this->m_UseExplicitPDFDerivatives ) {
802 memset( m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer(),
804 m_JointPDFDerivativesBufferSize );
807 memset( m_JointPDF->GetBufferPointer(),
809 m_JointPDFBufferSize );
810 memset( m_FixedImageMarginalPDF,
812 m_NumberOfHistogramBins*sizeof(PDFValueType) );
814 if( this->m_UseExplicitPDFDerivatives ) {
815 memset( m_JointPDFDerivatives->GetBufferPointer(),
817 m_JointPDFDerivativesBufferSize );
822 template < class TFixedImage, class TMovingImage >
824 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
825 ::GetValueAndDerivativeThreadProcessSample( unsigned int threadID,
826 unsigned long fixedImageSample,
827 const MovingImagePointType & itkNotUsed(mappedPoint),
828 double movingImageValue,
829 const ImageDerivativesType &
830 movingImageGradientValue) const
833 * Compute this sample's contribution to the marginal
834 * and joint distributions.
837 if(movingImageValue < m_MovingImageTrueMin) {
839 } else if(movingImageValue > m_MovingImageTrueMax) {
843 unsigned int fixedImageParzenWindowIndex =
844 this->m_FixedImageSamples[fixedImageSample].valueIndex;
846 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
847 double movingImageParzenWindowTerm = movingImageValue
848 / m_MovingImageBinSize
849 - m_MovingImageNormalizedMin;
850 unsigned int movingImageParzenWindowIndex =
851 static_cast<unsigned int>( movingImageParzenWindowTerm );
853 // Make sure the extreme values are in valid bins
854 if ( movingImageParzenWindowIndex < 2 ) {
855 movingImageParzenWindowIndex = 2;
856 } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
857 movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
860 // Since a zero-order BSpline (box car) kernel is used for
861 // the fixed image marginal pdf, we need only increment the
862 // fixedImageParzenWindowIndex by value of 1.0.
864 ++m_ThreaderFixedImageMarginalPDF[(threadID-1)*m_NumberOfHistogramBins
865 + fixedImageParzenWindowIndex];
867 ++m_FixedImageMarginalPDF[fixedImageParzenWindowIndex];
871 * The region of support of the parzen window determines which bins
872 * of the joint PDF are effected by the pair of image values.
873 * Since we are using a cubic spline for the moving image parzen
874 * window, four bins are effected. The fixed image parzen window is
875 * a zero-order spline (box car) and thus effects only one bin.
877 * The PDF is arranged so that moving image bins corresponds to the
878 * zero-th (column) dimension and the fixed image bins corresponds
879 * to the first (row) dimension.
883 // Pointer to affected bin to be updated
884 JointPDFValueType *pdfPtr;
886 pdfPtr = m_ThreaderJointPDF[threadID-1]
887 ->GetBufferPointer() +
888 ( fixedImageParzenWindowIndex
889 * m_NumberOfHistogramBins );
891 pdfPtr = m_JointPDF->GetBufferPointer() +
892 ( fixedImageParzenWindowIndex
893 * m_NumberOfHistogramBins );
896 // Move the pointer to the fist affected bin
897 int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
898 pdfPtr += pdfMovingIndex;
899 int pdfMovingIndexMax = static_cast<int>(movingImageParzenWindowIndex) + 2;
901 double movingImageParzenWindowArg = static_cast<double>( pdfMovingIndex )
902 - static_cast<double>( movingImageParzenWindowTerm );
904 while( pdfMovingIndex <= pdfMovingIndexMax ) {
905 *(pdfPtr++) += static_cast<PDFValueType>( m_CubicBSplineKernel
907 movingImageParzenWindowArg ) );
909 if( this->m_UseExplicitPDFDerivatives || this->m_ImplicitDerivativesSecondPass ) {
910 // Compute the cubicBSplineDerivative for later repeated use.
911 double cubicBSplineDerivativeValue =
912 m_CubicBSplineDerivativeKernel->Evaluate( movingImageParzenWindowArg );
914 // Compute PDF derivative contribution.
915 this->ComputePDFDerivatives( threadID,
918 movingImageGradientValue,
919 cubicBSplineDerivativeValue );
922 movingImageParzenWindowArg += 1;
929 template < class TFixedImage, class TMovingImage >
931 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
932 ::GetValueAndDerivativeThreadPostProcess( unsigned int threadID,
933 bool withinSampleThread ) const
935 this->GetValueThreadPostProcess( threadID, withinSampleThread );
937 if( this->m_UseExplicitPDFDerivatives ) {
938 const unsigned int rowSize = this->m_NumberOfParameters * m_NumberOfHistogramBins;
940 const unsigned int maxI =
941 rowSize * ( m_ThreaderJointPDFEndBin[threadID]
942 - m_ThreaderJointPDFStartBin[threadID] + 1 );
944 JointPDFDerivativesValueType *pdfDPtr;
945 JointPDFDerivativesValueType *pdfDPtrStart;
946 pdfDPtrStart = m_JointPDFDerivatives->GetBufferPointer()
947 + ( m_ThreaderJointPDFStartBin[threadID] * rowSize );
948 JointPDFDerivativesValueType *tPdfDPtr;
949 JointPDFDerivativesValueType *tPdfDPtrEnd;
950 unsigned int tPdfDPtrOffset;
951 tPdfDPtrOffset = m_ThreaderJointPDFStartBin[threadID] * rowSize;
952 #if ITK_VERSION_MAJOR <= 4
953 for(unsigned int t=0; t<this->m_NumberOfThreads-1; t++) {
955 for(unsigned int t=0; t<this->m_NumberOfWorkUnits-1; t++) {
957 pdfDPtr = pdfDPtrStart;
958 tPdfDPtr = m_ThreaderJointPDFDerivatives[t]->GetBufferPointer()
960 tPdfDPtrEnd = tPdfDPtr + maxI;
961 // for(i = 0; i < maxI; i++)
962 while(tPdfDPtr < tPdfDPtrEnd) {
963 *(pdfDPtr++) += *(tPdfDPtr++);
967 double nFactor = 1.0 / (m_MovingImageBinSize
968 * this->m_NumberOfPixelsCounted);
970 pdfDPtr = pdfDPtrStart;
971 tPdfDPtrEnd = pdfDPtrStart + maxI;
972 //for(int i = 0; i < maxI; i++)
973 while(pdfDPtr < tPdfDPtrEnd) {
974 *(pdfDPtr++) *= nFactor;
980 * Get the both Value and Derivative Measure
982 template < class TFixedImage, class TMovingImage >
984 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
985 ::GetValueAndDerivative( const ParametersType & parameters,
987 DerivativeType & derivative) const
989 // Set output values to zero
990 value = NumericTraits< MeasureType >::Zero;
992 if( this->m_UseExplicitPDFDerivatives ) {
993 // Set output values to zero
994 if(derivative.GetSize() != this->m_NumberOfParameters) {
995 derivative = DerivativeType( this->m_NumberOfParameters );
997 memset( derivative.data_block(),
999 this->m_NumberOfParameters * sizeof(double) );
1001 this->m_PRatioArray.Fill( 0.0 );
1002 this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero );
1003 #if ITK_VERSION_MAJOR <= 4
1004 for(unsigned int threadID = 0; threadID < this->m_NumberOfThreads-1; threadID++ ) {
1006 for(unsigned int threadID = 0; threadID < this->m_NumberOfWorkUnits-1; threadID++ ) {
1008 this->m_ThreaderMetricDerivative[threadID].Fill( NumericTraits< MeasureType >::Zero );
1010 this->m_ImplicitDerivativesSecondPass = false;
1013 // Set up the parameters in the transform
1014 this->m_Transform->SetParameters( parameters );
1016 // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
1017 this->GetValueAndDerivativeMultiThreadedInitiate();
1019 // CALL IF DOING THREADED POST PROCESSING
1020 this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
1022 #if ITK_VERSION_MAJOR <= 4
1023 for(unsigned int threadID = 0; threadID<this->m_NumberOfThreads-1; threadID++) {
1025 for(unsigned int threadID = 0; threadID<this->m_NumberOfWorkUnits-1; threadID++) {
1027 m_JointPDFSum += m_ThreaderJointPDFSum[threadID];
1029 if ( m_JointPDFSum == 0.0 ) {
1030 itkExceptionMacro( "Joint PDF summed to zero" );
1033 memset( m_MovingImageMarginalPDF,
1035 m_NumberOfHistogramBins*sizeof(PDFValueType) );
1037 JointPDFValueType * pdfPtr;
1038 PDFValueType * movingMarginalPtr;
1040 double fixedPDFSum = 0.0;
1041 const double normalizationFactor = 1.0 / m_JointPDFSum;
1043 pdfPtr = m_JointPDF->GetBufferPointer();
1044 for(i=0; i<m_NumberOfHistogramBins; i++) {
1045 fixedPDFSum += m_FixedImageMarginalPDF[i];
1046 movingMarginalPtr = m_MovingImageMarginalPDF;
1047 for(j=0; j<m_NumberOfHistogramBins; j++) {
1048 *(pdfPtr) *= normalizationFactor;
1049 *(movingMarginalPtr++) += *(pdfPtr++);
1053 if( this->m_NumberOfPixelsCounted <
1054 this->m_NumberOfFixedImageSamples / 16 ) {
1055 itkExceptionMacro( "Too many samples map outside moving image buffer: "
1056 << this->m_NumberOfPixelsCounted << " / "
1057 << this->m_NumberOfFixedImageSamples
1061 // Normalize the fixed image marginal PDF
1062 if ( fixedPDFSum == 0.0 ) {
1063 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
1065 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
1066 m_FixedImageMarginalPDF[bin] /= fixedPDFSum;
1070 * Compute the metric by double summation over histogram.
1073 // Setup pointer to point to the first bin
1074 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
1076 // Initialize sum to zero
1079 const double nFactor = 1.0 / (m_MovingImageBinSize
1080 * this->m_NumberOfPixelsCounted);
1082 for( unsigned int fixedIndex = 0;
1083 fixedIndex < m_NumberOfHistogramBins;
1085 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
1087 for( unsigned int movingIndex = 0;
1088 movingIndex < m_NumberOfHistogramBins;
1089 ++movingIndex, jointPDFPtr++ ) {
1090 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
1091 double jointPDFValue = *(jointPDFPtr);
1093 // check for non-zero bin contribution
1094 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
1096 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
1098 if( fixedImagePDFValue > 1e-16) {
1099 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
1102 if( this->m_UseExplicitPDFDerivatives ) {
1103 // move joint pdf derivative pointer to the right position
1104 JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1105 + ( fixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1106 + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1108 for( unsigned int parameter=0; parameter < this->m_NumberOfParameters; ++parameter, derivPtr++ ) {
1110 // Ref: eqn 23 of Thevenaz & Unser paper [3]
1111 derivative[parameter] -= (*derivPtr) * pRatio;
1113 } // end for-loop over parameters
1115 this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor;
1117 } // end if-block to check non-zero bin contribution
1118 } // end for-loop over moving index
1119 } // end for-loop over fixed index
1121 if( !(this->m_UseExplicitPDFDerivatives ) ) {
1122 // Second pass: This one is done for accumulating the contributions
1123 // to the derivative array.
1125 this->m_ImplicitDerivativesSecondPass = true;
1127 // MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
1128 this->GetValueAndDerivativeMultiThreadedInitiate();
1130 // CALL IF DOING THREADED POST PROCESSING
1131 this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
1133 // Consolidate the contributions from each one of the threads to the total
1135 #if ITK_VERSION_MAJOR <= 4
1136 for(unsigned int t = 0; t < this->m_NumberOfThreads-1; t++ ) {
1138 for(unsigned int t = 0; t < this->m_NumberOfWorkUnits-1; t++ ) {
1140 DerivativeType * source = &(this->m_ThreaderMetricDerivative[t]);
1141 for(unsigned int pp=0; pp < this->m_NumberOfParameters; pp++ ) {
1142 this->m_MetricDerivative[pp] += (*source)[pp];
1146 derivative = this->m_MetricDerivative;
1149 value = static_cast<MeasureType>( -1.0 * sum );
1154 * Get the match measure derivative
1156 template < class TFixedImage, class TMovingImage >
1158 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1159 ::GetDerivative( const ParametersType & parameters,
1160 DerivativeType & derivative ) const
1163 // call the combined version
1164 this->GetValueAndDerivative( parameters, value, derivative );
1169 * Compute PDF derivatives contribution for each parameter
1171 template < class TFixedImage, class TMovingImage >
1173 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1174 ::ComputePDFDerivatives( unsigned int threadID,
1175 unsigned int sampleNumber,
1177 const ImageDerivativesType & movingImageGradientValue,
1178 double cubicBSplineDerivativeValue ) const
1180 // Update bins in the PDF derivatives for the current intensity pair
1181 // Could pre-compute
1182 JointPDFDerivativesValueType * derivPtr;
1184 double precomputedWeight = 0.0;
1186 const int pdfFixedIndex = this->m_FixedImageSamples[sampleNumber].valueIndex;
1188 DerivativeType * derivativeHelperArray = NULL;
1190 if( this->m_UseExplicitPDFDerivatives ) {
1192 derivPtr = m_ThreaderJointPDFDerivatives[threadID-1]->GetBufferPointer()
1193 + ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1194 + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1196 derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1197 + ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1198 + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1202 // Recover the precomputed weight for this specific PDF bin
1203 precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex];
1205 derivativeHelperArray = &(this->m_ThreaderMetricDerivative[threadID-1]);
1207 derivativeHelperArray = &(this->m_MetricDerivative);
1212 if( !this->m_TransformIsBSpline ) {
1214 * Generic version which works for all transforms.
1217 // Compute the transform Jacobian.
1218 // Should pre-compute
1219 typedef typename TransformType::JacobianType JacobianType;
1221 // Need to use one of the threader transforms if we're
1224 // Use a raw pointer here to avoid the overhead of smart pointers.
1225 // For instance, Register and UnRegister have mutex locks around
1226 // the reference counts.
1227 TransformType* transform;
1230 transform = this->m_ThreaderTransform[threadID - 1];
1232 transform = this->m_Transform;
1235 JacobianType jacobian;
1236 transform->ComputeJacobianWithRespectToParameters(this->m_FixedImageSamples[sampleNumber].point, jacobian);
1238 // for ( unsigned int mu = 0; mu < this->m_NumberOfParameters; mu++ )
1240 // double innerProduct = 0.0;
1241 // for ( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ )
1243 // innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1246 // const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1248 // if( this->m_UseExplicitPDFDerivatives )
1250 // *(derivPtr) -= derivativeContribution;
1255 // (*derivativeHelperArray)[mu] += precomputedWeight * derivativeContribution;
1260 unsigned int mu, dim;
1261 double innerProduct,derivativeContribution;
1262 for ( mu = 0; mu < this->m_NumberOfParameters; mu+=3 ) {
1263 // JV only for 3D Space BLUT FFD: if J(0, par)=0, then J(1,par+1)=0 & ...
1264 if (jacobian[0][mu])
1265 for ( dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1266 innerProduct = jacobian[dim][mu+dim] * movingImageGradientValue[dim];
1267 derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1269 if( this->m_UseExplicitPDFDerivatives ) {
1270 *(derivPtr) -= derivativeContribution;
1273 (*derivativeHelperArray)[mu+dim] += precomputedWeight * derivativeContribution;
1276 else if( this->m_UseExplicitPDFDerivatives ) derivPtr+=3;
1279 const WeightsValueType * weights = NULL;
1280 const IndexValueType * indices = NULL;
1283 BSplineTransformWeightsType * weightsHelper = NULL;
1284 BSplineTransformIndexArrayType * indicesHelper = NULL;
1286 if( this->m_UseCachingOfBSplineWeights ) {
1288 // If the transform is of type BSplineDeformableTransform, we can obtain
1289 // a speed up by only processing the affected parameters. Note that
1290 // these pointers are just pointing to pre-allocated rows of the caching
1291 // arrays. There is therefore, no need to free this memory.
1293 weights = this->m_BSplineTransformWeightsArray[sampleNumber];
1294 indices = this->m_BSplineTransformIndicesArray[sampleNumber];
1296 if( threadID > 0 ) {
1297 weightsHelper = &(this->m_ThreaderBSplineTransformWeights[threadID-1]);
1298 indicesHelper = &(this->m_ThreaderBSplineTransformIndices[threadID-1]);
1300 weightsHelper = &(this->m_BSplineTransformWeights);
1301 indicesHelper = &(this->m_BSplineTransformIndices);
1304 this->m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition(
1305 this->m_FixedImageSamples[sampleNumber].point,
1306 *weightsHelper, *indicesHelper );
1309 for( unsigned int dim = 0; dim < Superclass::FixedImageDimension; dim++ ) {
1311 double innerProduct;
1314 for( unsigned int mu = 0; mu < this->m_NumBSplineWeights; mu++ ) {
1316 /* The array weights contains the Jacobian values in a 1-D array
1317 * (because for each parameter the Jacobian is non-zero in only 1 of the
1318 * possible dimensions) which is multiplied by the moving image
1320 if( this->m_UseCachingOfBSplineWeights ) {
1321 innerProduct = movingImageGradientValue[dim] * weights[mu];
1322 parameterIndex = indices[mu] + this->m_BSplineParametersOffset[dim];
1324 innerProduct = movingImageGradientValue[dim] * (*weightsHelper)[mu];
1325 parameterIndex = (*indicesHelper)[mu] + this->m_BSplineParametersOffset[dim];
1328 const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1330 if( this->m_UseExplicitPDFDerivatives ) {
1331 JointPDFValueType * ptr = derivPtr + parameterIndex;
1332 *(ptr) -= derivativeContribution;
1334 (*derivativeHelperArray)[parameterIndex] += precomputedWeight * derivativeContribution;
1338 } //end dim for loop
1340 } // end if-block transform is BSpline
1345 } // end namespace itk