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: itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.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 __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
36 #define __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx
39 // First make sure that the configuration is available.
40 // This line can be removed once the optimized versions
41 // gets integrated into the main directories.
42 #include "itkConfigure.h"
44 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
45 #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx"
49 #include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h"
50 #include "itkBSplineInterpolateImageFunction.h"
51 #include "itkCovariantVector.h"
52 #include "itkImageRandomConstIteratorWithIndex.h"
53 #include "itkImageRegionConstIterator.h"
54 #include "itkImageRegionIterator.h"
55 #include "itkImageIterator.h"
56 #include "vnl/vnl_math.h"
57 #if ITK_VERSION_MAJOR >= 4
58 #include "itkBSplineTransform.h"
60 #include "itkBSplineDeformableTransform.h"
70 template < class TFixedImage, class TMovingImage >
71 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
72 ::MattesMutualInformationImageToImageMetricFor3DBLUTFFD()
75 m_NumberOfSpatialSamples = 500;
76 m_NumberOfHistogramBins = 50;
78 this->SetComputeGradient(false); // don't use the default gradient for now
80 m_InterpolatorIsBSpline = false;
81 m_TransformIsBSpline = false;
83 // Initialize PDFs to NULL
85 m_JointPDFDerivatives = NULL;
87 m_UseExplicitPDFDerivatives = true;
89 typename BSplineTransformType::Pointer transformer =
90 BSplineTransformType::New();
91 this->SetTransform (transformer);
93 typename BSplineInterpolatorType::Pointer interpolator =
94 BSplineInterpolatorType::New();
95 this->SetInterpolator (interpolator);
98 m_MovingImageNormalizedMin = 0.0;
99 m_FixedImageNormalizedMin = 0.0;
100 m_MovingImageTrueMin = 0.0;
101 m_MovingImageTrueMax = 0.0;
102 m_FixedImageBinSize = 0.0;
103 m_MovingImageBinSize = 0.0;
104 m_CubicBSplineDerivativeKernel = NULL;
105 m_BSplineInterpolator = NULL;
106 m_DerivativeCalculator = NULL;
107 m_NumParametersPerDim = 0;
108 m_NumBSplineWeights = 0;
109 m_BSplineTransform = NULL;
110 m_NumberOfParameters = 0;
111 m_UseAllPixels = false;
112 m_ReseedIterator = false;
114 m_UseCachingOfBSplineWeights = true;
119 * Print out internal information about this class
121 template < class TFixedImage, class TMovingImage >
123 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
124 ::PrintSelf(std::ostream& os, Indent indent) const
127 Superclass::PrintSelf(os, indent);
129 os << indent << "NumberOfSpatialSamples: ";
130 os << m_NumberOfSpatialSamples << std::endl;
131 os << indent << "NumberOfHistogramBins: ";
132 os << m_NumberOfHistogramBins << std::endl;
133 os << indent << "UseAllPixels: ";
134 os << m_UseAllPixels << std::endl;
136 // Debugging information
137 os << indent << "NumberOfParameters: ";
138 os << m_NumberOfParameters << std::endl;
139 os << indent << "FixedImageNormalizedMin: ";
140 os << m_FixedImageNormalizedMin << std::endl;
141 os << indent << "MovingImageNormalizedMin: ";
142 os << m_MovingImageNormalizedMin << std::endl;
143 os << indent << "MovingImageTrueMin: ";
144 os << m_MovingImageTrueMin << std::endl;
145 os << indent << "MovingImageTrueMax: ";
146 os << m_MovingImageTrueMax << std::endl;
147 os << indent << "FixedImageBinSize: ";
148 os << m_FixedImageBinSize << std::endl;
149 os << indent << "MovingImageBinSize: ";
150 os << m_MovingImageBinSize << std::endl;
151 os << indent << "InterpolatorIsBSpline: ";
152 os << m_InterpolatorIsBSpline << std::endl;
153 os << indent << "TransformIsBSpline: ";
154 os << m_TransformIsBSpline << std::endl;
155 os << indent << "UseCachingOfBSplineWeights: ";
156 os << m_UseCachingOfBSplineWeights << std::endl;
157 os << indent << "UseExplicitPDFDerivatives: ";
158 os << m_UseExplicitPDFDerivatives << std::endl;
166 template <class TFixedImage, class TMovingImage>
168 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
169 ::Initialize(void) throw ( ExceptionObject )
171 this->Superclass::Initialize();
173 // Cache the number of transformation parameters
174 m_NumberOfParameters = this->m_Transform->GetNumberOfParameters();
177 * Compute the minimum and maximum for the FixedImage over
178 * the FixedImageRegion.
180 * NB: We can't use StatisticsImageFilter to do this because
181 * the filter computes the min/max for the largest possible region
183 double fixedImageMin = NumericTraits<double>::max();
184 double fixedImageMax = NumericTraits<double>::NonpositiveMin();
186 typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
187 FixedIteratorType fixedImageIterator(
188 this->m_FixedImage, this->GetFixedImageRegion() );
190 for ( fixedImageIterator.GoToBegin();
191 !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
193 double sample = static_cast<double>( fixedImageIterator.Get() );
195 if ( sample < fixedImageMin ) {
196 fixedImageMin = sample;
199 if ( sample > fixedImageMax ) {
200 fixedImageMax = sample;
206 * Compute the minimum and maximum for the entire moving image
209 double movingImageMin = NumericTraits<double>::max();
210 double movingImageMax = NumericTraits<double>::NonpositiveMin();
212 typedef ImageRegionConstIterator<MovingImageType> MovingIteratorType;
213 MovingIteratorType movingImageIterator(
214 this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() );
216 for ( movingImageIterator.GoToBegin();
217 !movingImageIterator.IsAtEnd(); ++movingImageIterator) {
218 double sample = static_cast<double>( movingImageIterator.Get() );
220 if ( sample < movingImageMin ) {
221 movingImageMin = sample;
224 if ( sample > movingImageMax ) {
225 movingImageMax = sample;
229 m_MovingImageTrueMin = movingImageMin;
230 m_MovingImageTrueMax = movingImageMax;
232 itkDebugMacro( " FixedImageMin: " << fixedImageMin <<
233 " FixedImageMax: " << fixedImageMax << std::endl );
234 itkDebugMacro( " MovingImageMin: " << movingImageMin <<
235 " MovingImageMax: " << movingImageMax << std::endl );
239 * Compute binsize for the histograms.
241 * The binsize for the image intensities needs to be adjusted so that
242 * we can avoid dealing with boundary conditions using the cubic
243 * spline as the Parzen window. We do this by increasing the size
244 * of the bins so that the joint histogram becomes "padded" at the
245 * borders. Because we are changing the binsize,
246 * we also need to shift the minimum by the padded amount in order to
247 * avoid minimum values filling in our padded region.
249 * Note that there can still be non-zero bin values in the padded region,
250 * it's just that these bins will never be a central bin for the Parzen
254 const int padding = 2; // this will pad by 2 bins
256 m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) /
257 static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
258 m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize -
259 static_cast<double>( padding );
261 m_MovingImageBinSize = ( movingImageMax - movingImageMin ) /
262 static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
263 m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize -
264 static_cast<double>( padding );
267 itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin );
268 itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin );
269 itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize );
270 itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize );
272 if( m_UseAllPixels ) {
273 m_NumberOfSpatialSamples =
274 this->GetFixedImageRegion().GetNumberOfPixels();
278 * Allocate memory for the fixed image sample container.
280 m_FixedImageSamples.resize( m_NumberOfSpatialSamples );
284 * Allocate memory for the marginal PDF and initialize values
285 * to zero. The marginal PDFs are stored as std::vector.
287 m_FixedImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 );
288 m_MovingImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 );
291 * Allocate memory for the joint PDF and joint PDF derivatives.
292 * The joint PDF and joint PDF derivatives are store as itk::Image.
294 m_JointPDF = JointPDFType::New();
296 // Instantiate a region, index, size
297 JointPDFRegionType jointPDFRegion;
298 JointPDFIndexType jointPDFIndex;
299 JointPDFSizeType jointPDFSize;
301 // Deallocate the memory that may have been allocated for
302 // previous runs of the metric.
303 this->m_JointPDFDerivatives = NULL; // by destroying the dynamic array
304 this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones
305 this->m_MetricDerivative = DerivativeType( 1 );
308 // Now allocate memory according to the user-selected method.
310 if( this->m_UseExplicitPDFDerivatives ) {
311 this->m_JointPDFDerivatives = JointPDFDerivativesType::New();
312 JointPDFDerivativesRegionType jointPDFDerivativesRegion;
313 JointPDFDerivativesIndexType jointPDFDerivativesIndex;
314 JointPDFDerivativesSizeType jointPDFDerivativesSize;
316 // For the derivatives of the joint PDF define a region starting from {0,0,0}
317 // with size {m_NumberOfParameters,m_NumberOfHistogramBins,
318 // m_NumberOfHistogramBins}. The dimension represents transform parameters,
319 // fixed image parzen window index and moving image parzen window index,
321 jointPDFDerivativesIndex.Fill( 0 );
322 jointPDFDerivativesSize[0] = m_NumberOfParameters;
323 jointPDFDerivativesSize[1] = m_NumberOfHistogramBins;
324 jointPDFDerivativesSize[2] = m_NumberOfHistogramBins;
326 jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex );
327 jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize );
329 // Set the regions and allocate
330 m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion );
331 m_JointPDFDerivatives->Allocate();
333 /** Allocate memory for helper array that will contain the pRatios
334 * for each bin of the joint histogram. This is part of the effort
335 * for flattening the computation of the PDF Jacobians.
337 this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins );
338 this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() );
341 // For the joint PDF define a region starting from {0,0}
342 // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}.
343 // The dimension represents fixed image parzen window index
344 // and moving image parzen window index, respectively.
345 jointPDFIndex.Fill( 0 );
346 jointPDFSize.Fill( m_NumberOfHistogramBins );
348 jointPDFRegion.SetIndex( jointPDFIndex );
349 jointPDFRegion.SetSize( jointPDFSize );
351 // Set the regions and allocate
352 m_JointPDF->SetRegions( jointPDFRegion );
353 m_JointPDF->Allocate();
357 * Setup the kernels used for the Parzen windows.
359 m_CubicBSplineKernel = CubicBSplineFunctionType::New();
360 m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New();
363 if( m_UseAllPixels ) {
365 * Take all the pixels within the fixed image region)
366 * to create the sample points list.
368 this->SampleFullFixedImageDomain( m_FixedImageSamples );
371 * Uniformly sample the fixed image (within the fixed image region)
372 * to create the sample points list.
374 this->SampleFixedImageDomain( m_FixedImageSamples );
378 * Pre-compute the fixed image parzen window index for
379 * each point of the fixed image sample points list.
381 this->ComputeFixedImageParzenWindowIndices( m_FixedImageSamples );
384 * Check if the interpolator is of type BSplineInterpolateImageFunction.
385 * If so, we can make use of its EvaluateDerivatives method.
386 * Otherwise, we instantiate an external central difference
387 * derivative calculator.
389 * TODO: Also add it the possibility of using the default gradient
390 * provided by the superclass.
393 m_InterpolatorIsBSpline = true;
395 BSplineInterpolatorType * testPtr = dynamic_cast<BSplineInterpolatorType *>(
396 this->m_Interpolator.GetPointer() );
398 m_InterpolatorIsBSpline = false;
400 m_DerivativeCalculator = DerivativeFunctionType::New();
402 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
403 m_DerivativeCalculator->UseImageDirectionOn();
406 m_DerivativeCalculator->SetInputImage( this->m_MovingImage );
408 m_BSplineInterpolator = NULL;
409 itkDebugMacro( "Interpolator is not BSpline" );
411 m_BSplineInterpolator = testPtr;
413 #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
414 m_BSplineInterpolator->UseImageDirectionOn();
417 m_DerivativeCalculator = NULL;
418 itkDebugMacro( "Interpolator is BSpline" );
422 * Check if the transform is of type BSplineDeformableTransform.
424 * If so, several speed up features are implemented.
425 * [1] Precomputing the results of bulk transform for each sample point.
426 * [2] Precomputing the BSpline weights for each sample point,
427 * to be used later to directly compute the deformation vector
428 * [3] Precomputing the indices of the parameters within the
429 * the support region of each sample point.
431 m_TransformIsBSpline = true;
433 BSplineTransformType * testPtr2 = dynamic_cast<BSplineTransformType *>(
434 this->m_Transform.GetPointer() );
436 m_TransformIsBSpline = false;
437 m_BSplineTransform = NULL;
438 itkDebugMacro( "Transform is not BSplineDeformable" );
440 m_BSplineTransform = testPtr2;
441 m_NumParametersPerDim =
442 m_BSplineTransform->GetNumberOfParametersPerDimension();
443 m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights();
444 itkDebugMacro( "Transform is BSplineDeformable" );
447 if( this->m_TransformIsBSpline ) {
448 // First, deallocate memory that may have been used from previous run of the Metric
449 this->m_BSplineTransformWeightsArray.SetSize( 1, 1 );
450 this->m_BSplineTransformIndicesArray.SetSize( 1, 1 );
451 this->m_PreTransformPointsArray.resize( 1 );
452 this->m_WithinSupportRegionArray.resize( 1 );
453 this->m_Weights.SetSize( 1 );
454 this->m_Indices.SetSize( 1 );
457 if( this->m_UseCachingOfBSplineWeights ) {
458 m_BSplineTransformWeightsArray.SetSize(
459 m_NumberOfSpatialSamples, m_NumBSplineWeights );
460 m_BSplineTransformIndicesArray.SetSize(
461 m_NumberOfSpatialSamples, m_NumBSplineWeights );
462 m_PreTransformPointsArray.resize( m_NumberOfSpatialSamples );
463 m_WithinSupportRegionArray.resize( m_NumberOfSpatialSamples );
465 this->PreComputeTransformValues();
467 this->m_Weights.SetSize( this->m_NumBSplineWeights );
468 this->m_Indices.SetSize( this->m_NumBSplineWeights );
471 for ( unsigned int j = 0; j < FixedImageDimension; j++ ) {
472 m_ParametersOffset[j] = j *
473 m_BSplineTransform->GetNumberOfParametersPerDimension();
481 * Uniformly sample the fixed image domain using a random walk
483 template < class TFixedImage, class TMovingImage >
485 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
486 ::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples )
489 // Set up a random interator within the user specified fixed image region.
490 typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
491 RandomIterator randIter( this->m_FixedImage, this->GetFixedImageRegion() );
493 randIter.SetNumberOfSamples( m_NumberOfSpatialSamples );
494 randIter.GoToBegin();
496 typename FixedImageSpatialSampleContainer::iterator iter;
497 typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
499 if( this->m_FixedImageMask ) {
501 InputPointType inputPoint;
503 iter=samples.begin();
505 int samples_found = 0;
506 int maxcount = m_NumberOfSpatialSamples * 10;
507 while( iter != end ) {
509 if ( count > maxcount ) {
512 "Drew too many samples from the mask (is it too small?): "
513 << maxcount << std::endl );
515 samples.resize(samples_found);
516 // this->SetNumberOfSpatialSamples(sample_found);
523 FixedImageIndexType index = randIter.GetIndex();
524 // Check if the Index is inside the mask, translate index to point
525 this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
527 // If not inside the mask, ignore the point
528 if( !this->m_FixedImageMask->IsInside( inputPoint ) ) {
529 ++randIter; // jump to another random position
533 // Get sampled fixed image value
534 (*iter).FixedImageValue = randIter.Get();
535 // Translate index to point
536 (*iter).FixedImagePointValue = inputPoint;
538 // Jump to random position
543 for( iter=samples.begin(); iter != end; ++iter ) {
545 FixedImageIndexType index = randIter.GetIndex();
546 // Get sampled fixed image value
547 (*iter).FixedImageValue = randIter.Get();
548 // Translate index to point
549 this->m_FixedImage->TransformIndexToPhysicalPoint( index,
550 (*iter).FixedImagePointValue );
551 // Jump to random position
559 * Sample the fixed image domain using all pixels in the Fixed image region
561 template < class TFixedImage, class TMovingImage >
563 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
564 ::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples )
567 // Set up a region interator within the user specified fixed image region.
568 typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
569 RegionIterator regionIter( this->m_FixedImage, this->GetFixedImageRegion() );
571 regionIter.GoToBegin();
573 typename FixedImageSpatialSampleContainer::iterator iter;
574 typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
576 if( this->m_FixedImageMask ) {
577 InputPointType inputPoint;
579 iter=samples.begin();
580 unsigned long nSamplesPicked = 0;
582 while( iter != end && !regionIter.IsAtEnd() ) {
584 FixedImageIndexType index = regionIter.GetIndex();
585 // Check if the Index is inside the mask, translate index to point
586 this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
588 // If not inside the mask, ignore the point
589 if( !this->m_FixedImageMask->IsInside( inputPoint ) ) {
590 ++regionIter; // jump to next pixel
594 // Get sampled fixed image value
595 (*iter).FixedImageValue = regionIter.Get();
596 // Translate index to point
597 (*iter).FixedImagePointValue = inputPoint;
604 // If we picked fewer samples than the desired number,
605 // resize the container
606 if (nSamplesPicked != this->m_NumberOfSpatialSamples) {
607 this->m_NumberOfSpatialSamples = nSamplesPicked;
608 samples.resize(this->m_NumberOfSpatialSamples);
610 } else { // not restricting sample throwing to a mask
612 // cannot sample more than the number of pixels in the image region
613 if ( this->m_NumberOfSpatialSamples
614 > this->GetFixedImageRegion().GetNumberOfPixels()) {
615 this->m_NumberOfSpatialSamples
616 = this->GetFixedImageRegion().GetNumberOfPixels();
617 samples.resize(this->m_NumberOfSpatialSamples);
620 for( iter=samples.begin(); iter != end; ++iter ) {
622 FixedImageIndexType index = regionIter.GetIndex();
623 // Get sampled fixed image value
624 (*iter).FixedImageValue = regionIter.Get();
625 // Translate index to point
626 this->m_FixedImage->TransformIndexToPhysicalPoint( index,
627 (*iter).FixedImagePointValue );
634 * Uniformly sample the fixed image domain using a random walk
636 template < class TFixedImage, class TMovingImage >
638 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
639 ::ComputeFixedImageParzenWindowIndices(
640 FixedImageSpatialSampleContainer& samples )
643 typename FixedImageSpatialSampleContainer::iterator iter;
644 typename FixedImageSpatialSampleContainer::const_iterator end=samples.end();
646 for( iter=samples.begin(); iter != end; ++iter ) {
648 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
650 static_cast<double>( (*iter).FixedImageValue ) / m_FixedImageBinSize -
651 m_FixedImageNormalizedMin;
652 unsigned int pindex = static_cast<unsigned int>( vcl_floor(windowTerm ) );
654 // Make sure the extreme values are in valid bins
657 } else if ( pindex > (m_NumberOfHistogramBins - 3) ) {
658 pindex = m_NumberOfHistogramBins - 3;
661 (*iter).FixedImageParzenWindowIndex = pindex;
668 * Get the match Measure
670 template < class TFixedImage, class TMovingImage >
671 typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
673 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
674 ::GetValue( const ParametersType& parameters ) const
677 // Reset marginal pdf to all zeros.
678 // Assumed the size has already been set to NumberOfHistogramBins
680 for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) {
681 m_FixedImageMarginalPDF[j] = 0.0;
682 m_MovingImageMarginalPDF[j] = 0.0;
685 // Reset the joint pdfs to zero
686 m_JointPDF->FillBuffer( 0.0 );
688 // Set up the parameters in the transform
689 this->m_Transform->SetParameters( parameters );
692 // Declare iterators for iteration over the sample container
693 typename FixedImageSpatialSampleContainer::const_iterator fiter;
694 typename FixedImageSpatialSampleContainer::const_iterator fend =
695 m_FixedImageSamples.end();
697 unsigned long nSamples=0;
698 unsigned long nFixedImageSamples=0;
701 for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
703 // Get moving image value
704 MovingImagePointType mappedPoint;
706 double movingImageValue;
708 this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
709 sampleOk, movingImageValue );
711 ++nFixedImageSamples;
718 * Compute this sample's contribution to the marginal and
719 * joint distributions.
723 // Determine parzen window arguments (see eqn 6 of Mattes paper [2])
724 double movingImageParzenWindowTerm =
725 movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
726 unsigned int movingImageParzenWindowIndex =
727 static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
729 // Make sure the extreme values are in valid bins
730 if ( movingImageParzenWindowIndex < 2 ) {
731 movingImageParzenWindowIndex = 2;
732 } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
733 movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
737 // Since a zero-order BSpline (box car) kernel is used for
738 // the fixed image marginal pdf, we need only increment the
739 // fixedImageParzenWindowIndex by value of 1.0.
740 m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] +=
741 static_cast<PDFValueType>( 1 );
744 * The region of support of the parzen window determines which bins
745 * of the joint PDF are effected by the pair of image values.
746 * Since we are using a cubic spline for the moving image parzen
747 * window, four bins are affected. The fixed image parzen window is
748 * a zero-order spline (box car) and thus effects only one bin.
750 * The PDF is arranged so that moving image bins corresponds to the
751 * zero-th (column) dimension and the fixed image bins corresponds
752 * to the first (row) dimension.
756 // Pointer to affected bin to be updated
757 JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
758 ( (*fiter).FixedImageParzenWindowIndex
759 * m_JointPDF->GetOffsetTable()[1] );
761 // Move the pointer to the first affected bin
762 int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
763 pdfPtr += pdfMovingIndex;
765 for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
767 pdfMovingIndex++, pdfPtr++ ) {
769 // Update PDF for the current intensity pair
770 double movingImageParzenWindowArg =
771 static_cast<double>( pdfMovingIndex ) -
772 static_cast<double>( movingImageParzenWindowTerm );
774 *(pdfPtr) += static_cast<PDFValueType>(
775 m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );
777 } //end parzen windowing for loop
779 } //end if-block check sampleOk
780 } // end iterating over fixed image spatial sample container for loop
782 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
783 << nSamples << " / " << m_NumberOfSpatialSamples
786 if( nSamples < m_NumberOfSpatialSamples / 16 ) {
787 itkExceptionMacro( "Too many samples map outside moving image buffer: "
788 << nSamples << " / " << m_NumberOfSpatialSamples
792 this->m_NumberOfPixelsCounted = nSamples;
796 * Normalize the PDFs, compute moving image marginal PDF
799 typedef ImageRegionIterator<JointPDFType> JointPDFIteratorType;
800 JointPDFIteratorType jointPDFIterator ( m_JointPDF,
801 m_JointPDF->GetBufferedRegion() );
803 jointPDFIterator.GoToBegin();
805 // Compute joint PDF normalization factor
806 // (to ensure joint PDF sum adds to 1.0)
807 double jointPDFSum = 0.0;
809 while( !jointPDFIterator.IsAtEnd() ) {
810 jointPDFSum += jointPDFIterator.Get();
814 if ( jointPDFSum == 0.0 ) {
815 itkExceptionMacro( "Joint PDF summed to zero" );
819 // Normalize the PDF bins
820 jointPDFIterator.GoToEnd();
821 while( !jointPDFIterator.IsAtBegin() ) {
823 jointPDFIterator.Value() /= static_cast<PDFValueType>( jointPDFSum );
827 // Normalize the fixed image marginal PDF
828 double fixedPDFSum = 0.0;
829 for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) {
830 fixedPDFSum += m_FixedImageMarginalPDF[bin];
833 if ( fixedPDFSum == 0.0 ) {
834 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
837 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
838 m_FixedImageMarginalPDF[bin] /= static_cast<PDFValueType>( fixedPDFSum );
842 // Compute moving image marginal PDF by summing over fixed image bins.
843 typedef ImageLinearIteratorWithIndex<JointPDFType> JointPDFLinearIterator;
844 JointPDFLinearIterator linearIter( m_JointPDF,
845 m_JointPDF->GetBufferedRegion() );
847 linearIter.SetDirection( 1 );
848 linearIter.GoToBegin();
849 unsigned int movingIndex1 = 0;
851 while( !linearIter.IsAtEnd() ) {
855 while( !linearIter.IsAtEndOfLine() ) {
856 sum += linearIter.Get();
860 m_MovingImageMarginalPDF[movingIndex1] = static_cast<PDFValueType>(sum);
862 linearIter.NextLine();
868 * Compute the metric by double summation over histogram.
871 // Setup pointer to point to the first bin
872 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
874 // Initialze sum to zero
877 for( unsigned int fixedIndex = 0;
878 fixedIndex < m_NumberOfHistogramBins;
881 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
883 for( unsigned int movingIndex = 0;
884 movingIndex < m_NumberOfHistogramBins;
885 ++movingIndex, jointPDFPtr++ ) {
887 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
888 double jointPDFValue = *(jointPDFPtr);
890 // check for non-zero bin contribution
891 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
893 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
894 if( fixedImagePDFValue > 1e-16) {
895 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
898 } // end if-block to check non-zero bin contribution
899 } // end for-loop over moving index
900 } // end for-loop over fixed index
902 return static_cast<MeasureType>( -1.0 * sum );
908 * Get the both Value and Derivative Measure
910 template < class TFixedImage, class TMovingImage >
912 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
913 ::GetValueAndDerivative(
914 const ParametersType& parameters,
916 DerivativeType& derivative) const
919 // Set output values to zero
920 value = NumericTraits< MeasureType >::Zero;
922 if( this->m_UseExplicitPDFDerivatives ) {
923 m_JointPDFDerivatives->FillBuffer( 0.0 );
924 derivative = DerivativeType( this->GetNumberOfParameters() );
925 derivative.Fill( NumericTraits< MeasureType >::Zero );
927 this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero );
928 this->m_PRatioArray.Fill( 0.0 );
931 // Reset marginal pdf to all zeros.
932 // Assumed the size has already been set to NumberOfHistogramBins
934 for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) {
935 m_FixedImageMarginalPDF[j] = 0.0;
936 m_MovingImageMarginalPDF[j] = 0.0;
939 // Reset the joint pdfs to zero
940 m_JointPDF->FillBuffer( 0.0 );
943 // Set up the parameters in the transform
944 this->m_Transform->SetParameters( parameters );
947 // Declare iterators for iteration over the sample container
948 typename FixedImageSpatialSampleContainer::const_iterator fiter;
949 typename FixedImageSpatialSampleContainer::const_iterator fend =
950 m_FixedImageSamples.end();
952 unsigned long nSamples=0;
953 unsigned long nFixedImageSamples=0;
955 for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
957 // Get moving image value
958 MovingImagePointType mappedPoint;
960 double movingImageValue;
963 this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
964 sampleOk, movingImageValue );
969 // Get moving image derivative at the mapped position
970 ImageDerivativesType movingImageGradientValue;
971 this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue );
975 * Compute this sample's contribution to the marginal
976 * and joint distributions.
980 // Determine parzen window arguments (see eqn 6 of Mattes paper [2])
981 double movingImageParzenWindowTerm =
982 movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
983 unsigned int movingImageParzenWindowIndex =
984 static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
986 // Make sure the extreme values are in valid bins
987 if ( movingImageParzenWindowIndex < 2 ) {
988 movingImageParzenWindowIndex = 2;
989 } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
990 movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
994 // Since a zero-order BSpline (box car) kernel is used for
995 // the fixed image marginal pdf, we need only increment the
996 // fixedImageParzenWindowIndex by value of 1.0.
997 m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] +=
998 static_cast<PDFValueType>( 1 );
1001 * The region of support of the parzen window determines which bins
1002 * of the joint PDF are effected by the pair of image values.
1003 * Since we are using a cubic spline for the moving image parzen
1004 * window, four bins are effected. The fixed image parzen window is
1005 * a zero-order spline (box car) and thus effects only one bin.
1007 * The PDF is arranged so that moving image bins corresponds to the
1008 * zero-th (column) dimension and the fixed image bins corresponds
1009 * to the first (row) dimension.
1013 // Pointer to affected bin to be updated
1014 JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
1015 ( (*fiter).FixedImageParzenWindowIndex * m_NumberOfHistogramBins );
1017 // Move the pointer to the fist affected bin
1018 int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
1019 pdfPtr += pdfMovingIndex;
1021 for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
1023 pdfMovingIndex++, pdfPtr++ ) {
1024 // Update PDF for the current intensity pair
1025 double movingImageParzenWindowArg =
1026 static_cast<double>( pdfMovingIndex ) -
1027 static_cast<double>(movingImageParzenWindowTerm);
1029 *(pdfPtr) += static_cast<PDFValueType>(
1030 m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );
1032 if( this->m_UseExplicitPDFDerivatives ) {
1033 // Compute the cubicBSplineDerivative for later repeated use.
1034 double cubicBSplineDerivativeValue =
1035 m_CubicBSplineDerivativeKernel->Evaluate(
1036 movingImageParzenWindowArg );
1038 // Compute PDF derivative contribution.
1039 this->ComputePDFDerivatives( nFixedImageSamples,
1041 movingImageGradientValue,
1042 cubicBSplineDerivativeValue );
1045 } //end parzen windowing for loop
1047 } //end if-block check sampleOk
1049 ++nFixedImageSamples;
1051 } // end iterating over fixed image spatial sample container for loop
1053 itkDebugMacro( "Ratio of voxels mapping into moving image buffer: "
1054 << nSamples << " / " << m_NumberOfSpatialSamples
1057 if( nSamples < m_NumberOfSpatialSamples / 16 ) {
1058 itkExceptionMacro( "Too many samples map outside moving image buffer: "
1059 << nSamples << " / " << m_NumberOfSpatialSamples
1063 this->m_NumberOfPixelsCounted = nSamples;
1066 * Normalize the PDFs, compute moving image marginal PDF
1069 typedef ImageRegionIterator<JointPDFType> JointPDFIteratorType;
1070 JointPDFIteratorType jointPDFIterator ( m_JointPDF,
1071 m_JointPDF->GetBufferedRegion() );
1073 jointPDFIterator.GoToBegin();
1076 // Compute joint PDF normalization factor
1077 // (to ensure joint PDF sum adds to 1.0)
1078 double jointPDFSum = 0.0;
1080 while( !jointPDFIterator.IsAtEnd() ) {
1081 jointPDFSum += jointPDFIterator.Get();
1085 if ( jointPDFSum == 0.0 ) {
1086 itkExceptionMacro( "Joint PDF summed to zero" );
1090 // Normalize the PDF bins
1091 jointPDFIterator.GoToEnd();
1092 while( !jointPDFIterator.IsAtBegin() ) {
1094 jointPDFIterator.Value() /= static_cast<PDFValueType>( jointPDFSum );
1098 // Normalize the fixed image marginal PDF
1099 double fixedPDFSum = 0.0;
1100 for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) {
1101 fixedPDFSum += m_FixedImageMarginalPDF[bin];
1104 if ( fixedPDFSum == 0.0 ) {
1105 itkExceptionMacro( "Fixed image marginal PDF summed to zero" );
1108 for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) {
1109 m_FixedImageMarginalPDF[bin] /= static_cast<PDFValueType>( fixedPDFSum );
1113 // Compute moving image marginal PDF by summing over fixed image bins.
1114 typedef ImageLinearIteratorWithIndex<JointPDFType> JointPDFLinearIterator;
1115 JointPDFLinearIterator linearIter(
1116 m_JointPDF, m_JointPDF->GetBufferedRegion() );
1118 linearIter.SetDirection( 1 );
1119 linearIter.GoToBegin();
1120 unsigned int movingIndex1 = 0;
1122 while( !linearIter.IsAtEnd() ) {
1126 while( !linearIter.IsAtEndOfLine() ) {
1127 sum += linearIter.Get();
1131 m_MovingImageMarginalPDF[movingIndex1] = static_cast<PDFValueType>(sum);
1133 linearIter.NextLine();
1138 double nFactor = 1.0 / ( m_MovingImageBinSize
1139 * static_cast<double>( nSamples ) );
1141 if( this->m_UseExplicitPDFDerivatives ) {
1142 // Normalize the joint PDF derivatives by the test image binsize and nSamples
1143 typedef ImageRegionIterator<JointPDFDerivativesType>
1144 JointPDFDerivativesIteratorType;
1145 JointPDFDerivativesIteratorType jointPDFDerivativesIterator (
1146 m_JointPDFDerivatives,
1147 m_JointPDFDerivatives->GetBufferedRegion()
1150 jointPDFDerivativesIterator.GoToBegin();
1152 while( !jointPDFDerivativesIterator.IsAtEnd() ) {
1153 jointPDFDerivativesIterator.Value() *= nFactor;
1154 ++jointPDFDerivativesIterator;
1159 * Compute the metric by double summation over histogram.
1162 // Setup pointer to point to the first bin
1163 JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer();
1165 // Initialize sum to zero
1168 for( unsigned int fixedIndex = 0;
1169 fixedIndex < m_NumberOfHistogramBins;
1171 double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex];
1173 for( unsigned int movingIndex = 0; movingIndex < m_NumberOfHistogramBins;
1174 ++movingIndex, jointPDFPtr++ ) {
1175 double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex];
1176 double jointPDFValue = *(jointPDFPtr);
1178 // check for non-zero bin contribution
1179 if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) {
1181 double pRatio = vcl_log(jointPDFValue / movingImagePDFValue );
1183 if( fixedImagePDFValue > 1e-16) {
1184 sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) );
1187 if( this->m_UseExplicitPDFDerivatives ) {
1188 // move joint pdf derivative pointer to the right position
1189 JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer()
1190 + ( fixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] )
1191 + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1193 for( unsigned int parameter=0; parameter < m_NumberOfParameters; ++parameter, derivPtr++ ) {
1194 // Ref: eqn 23 of Thevenaz & Unser paper [3]
1195 derivative[parameter] -= (*derivPtr) * pRatio;
1196 } // end for-loop over parameters
1198 this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor;
1200 } // end if-block to check non-zero bin contribution
1201 } // end for-loop over moving index
1202 } // end for-loop over fixed index
1204 if( !(this->m_UseExplicitPDFDerivatives ) ) {
1205 // Second pass: This one is done for accumulating the contributions
1206 // to the derivative array.
1208 nFixedImageSamples = 0;
1210 for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) {
1212 // Get moving image value
1213 MovingImagePointType mappedPoint;
1215 double movingImageValue;
1217 this->TransformPoint( nFixedImageSamples, parameters, mappedPoint,
1218 sampleOk, movingImageValue );
1221 // Get moving image derivative at the mapped position
1222 ImageDerivativesType movingImageGradientValue;
1223 this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue );
1227 * Compute this sample's contribution to the marginal
1228 * and joint distributions.
1232 // Determine parzen window arguments (see eqn 6 of Mattes paper [2]).
1233 double movingImageParzenWindowTerm =
1234 movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin;
1235 unsigned int movingImageParzenWindowIndex =
1236 static_cast<unsigned int>( vcl_floor(movingImageParzenWindowTerm ) );
1238 // Make sure the extreme values are in valid bins
1239 if ( movingImageParzenWindowIndex < 2 ) {
1240 movingImageParzenWindowIndex = 2;
1241 } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) {
1242 movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3;
1246 // Move the pointer to the fist affected bin
1247 int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
1249 for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
1251 pdfMovingIndex++ ) {
1253 // Update PDF for the current intensity pair
1254 double movingImageParzenWindowArg =
1255 static_cast<double>( pdfMovingIndex ) -
1256 static_cast<double>(movingImageParzenWindowTerm);
1258 // Compute the cubicBSplineDerivative for later repeated use.
1259 double cubicBSplineDerivativeValue =
1260 m_CubicBSplineDerivativeKernel->Evaluate(
1261 movingImageParzenWindowArg );
1263 // Compute PDF derivative contribution.
1264 this->ComputePDFDerivatives( nFixedImageSamples,
1266 movingImageGradientValue,
1267 cubicBSplineDerivativeValue );
1270 } //end parzen windowing for loop
1272 } //end if-block check sampleOk
1274 ++nFixedImageSamples;
1276 } // end iterating over fixed image spatial sample container for loop
1279 derivative = this->m_MetricDerivative;
1282 value = static_cast<MeasureType>( -1.0 * sum );
1287 * Get the match measure derivative
1289 template < class TFixedImage, class TMovingImage >
1291 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1292 ::GetDerivative( const ParametersType& parameters,
1293 DerivativeType & derivative ) const
1296 // call the combined version
1297 this->GetValueAndDerivative( parameters, value, derivative );
1302 * Compute image derivatives using a central difference function
1303 * if we are not using a BSplineInterpolator, which includes
1306 template < class TFixedImage, class TMovingImage >
1308 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1309 ::ComputeImageDerivatives(
1310 const MovingImagePointType& mappedPoint,
1311 ImageDerivativesType& gradient ) const
1314 if( m_InterpolatorIsBSpline ) {
1315 // Computed moving image gradient using derivative BSpline kernel.
1316 gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint );
1318 // For all generic interpolator use central differencing.
1319 gradient = m_DerivativeCalculator->Evaluate( mappedPoint );
1326 * Transform a point from FixedImage domain to MovingImage domain.
1327 * This function also checks if mapped point is within support region.
1329 template < class TFixedImage, class TMovingImage >
1331 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1333 unsigned int sampleNumber,
1334 const ParametersType& parameters,
1335 MovingImagePointType& mappedPoint,
1337 double& movingImageValue ) const
1340 if ( !m_TransformIsBSpline ) {
1341 // Use generic transform to compute mapped position
1342 mappedPoint = this->m_Transform->TransformPoint(
1343 m_FixedImageSamples[sampleNumber].FixedImagePointValue );
1345 // Check if mapped point inside image buffer
1346 sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint );
1349 if( this->m_UseCachingOfBSplineWeights ) {
1350 // If the transform is BSplineDeformable, we can use the precomputed
1351 // weights and indices to obtained the mapped position
1353 const WeightsValueType * weights =
1354 m_BSplineTransformWeightsArray[sampleNumber];
1355 const IndexValueType * indices =
1356 m_BSplineTransformIndicesArray[sampleNumber];
1357 mappedPoint.Fill( 0.0 );
1359 if ( m_WithinSupportRegionArray[sampleNumber] ) {
1360 for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ ) {
1361 for ( unsigned int j = 0; j < FixedImageDimension; j++ ) {
1362 mappedPoint[j] += weights[k] *
1363 parameters[ indices[k] + m_ParametersOffset[j] ];
1368 for( unsigned int j = 0; j < FixedImageDimension; j++ ) {
1369 mappedPoint[j] += m_PreTransformPointsArray[sampleNumber][j];
1372 // Check if mapped point inside image buffer
1373 sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint );
1375 // Check if mapped point is within the support region of a grid point.
1376 // This is neccessary for computing the metric gradient
1377 sampleOk = sampleOk && m_WithinSupportRegionArray[sampleNumber];
1379 // If not caching values, we invoke the Transform to recompute the
1380 // mapping of the point.
1381 this->m_BSplineTransform->TransformPoint(
1382 this->m_FixedImageSamples[sampleNumber].FixedImagePointValue,
1383 mappedPoint, this->m_Weights, this->m_Indices, sampleOk);
1385 // Check if mapped point inside image buffer
1386 sampleOk = sampleOk && this->m_Interpolator->IsInsideBuffer( mappedPoint );
1391 // If user provided a mask over the Moving image
1392 if ( this->m_MovingImageMask ) {
1393 // Check if mapped point is within the support region of the moving image
1395 sampleOk = sampleOk && this->m_MovingImageMask->IsInside( mappedPoint );
1400 movingImageValue = this->m_Interpolator->Evaluate( mappedPoint );
1402 if ( movingImageValue < m_MovingImageTrueMin ||
1403 movingImageValue > m_MovingImageTrueMax ) {
1404 // need to throw out this sample as it will not fall into a valid bin
1412 * Compute PDF derivatives contribution for each parameter
1414 template < class TFixedImage, class TMovingImage >
1416 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1417 ::ComputePDFDerivatives(
1418 unsigned int sampleNumber,
1420 const ImageDerivativesType& movingImageGradientValue,
1421 double cubicBSplineDerivativeValue ) const
1424 const int pdfFixedIndex =
1425 m_FixedImageSamples[sampleNumber].FixedImageParzenWindowIndex;
1427 JointPDFValueType * derivPtr = NULL;
1428 double precomputedWeight = 0.0;
1430 if( this->m_UseExplicitPDFDerivatives ) {
1431 // Update bins in the PDF derivatives for the current intensity pair
1432 derivPtr = m_JointPDFDerivatives->GetBufferPointer() +
1433 ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] ) +
1434 ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] );
1436 // Recover the precomputed weight for this specific PDF bin
1437 precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex];
1440 if( !m_TransformIsBSpline ) {
1443 * Generic version which works for all transforms.
1446 // Compute the transform Jacobian.
1447 typedef typename TransformType::JacobianType JacobianType;
1448 #if ITK_VERSION_MAJOR >= 4
1449 JacobianType jacobian;
1450 this->m_Transform->ComputeJacobianWithRespectToParameters( m_FixedImageSamples[sampleNumber].FixedImagePointValue, jacobian );
1452 const JacobianType & jacobian =
1453 this->m_Transform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue );
1456 for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) {
1457 double innerProduct = 0.0;
1458 for ( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1459 innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim];
1462 const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1464 if( this->m_UseExplicitPDFDerivatives ) {
1465 *(derivPtr) -= derivativeContribution;
1468 this->m_MetricDerivative[mu] += precomputedWeight * derivativeContribution;
1473 const WeightsValueType * weights = NULL;
1474 const IndexValueType * indices = NULL;
1476 if( this->m_UseCachingOfBSplineWeights ) {
1478 * If the transform is of type BSplineDeformableTransform,
1479 * we can obtain a speed up by only processing the affected parameters.
1480 * Note that these pointers are just pointing to pre-allocated rows
1481 * of the caching arrays. There is therefore, no need to free this
1484 weights = m_BSplineTransformWeightsArray[sampleNumber];
1485 indices = m_BSplineTransformIndicesArray[sampleNumber];
1487 #if ITK_VERSION_MAJOR >= 4
1488 m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition(
1489 m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
1491 m_BSplineTransform->GetJacobian(
1492 m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
1496 for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) {
1498 double innerProduct;
1501 for( unsigned int mu = 0; mu < m_NumBSplineWeights; mu++ ) {
1503 /* The array weights contains the Jacobian values in a 1-D array
1504 * (because for each parameter the Jacobian is non-zero in only 1 of the
1505 * possible dimensions) which is multiplied by the moving image
1507 if( this->m_UseCachingOfBSplineWeights ) {
1508 innerProduct = movingImageGradientValue[dim] * weights[mu];
1509 parameterIndex = indices[mu] + m_ParametersOffset[dim];
1511 innerProduct = movingImageGradientValue[dim] * this->m_Weights[mu];
1512 parameterIndex = this->m_Indices[mu] + this->m_ParametersOffset[dim];
1515 const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue;
1517 if( this->m_UseExplicitPDFDerivatives ) {
1518 JointPDFValueType * ptr = derivPtr + parameterIndex;
1519 *(ptr) -= derivativeContribution;
1521 this->m_MetricDerivative[parameterIndex] += precomputedWeight * derivativeContribution;
1525 } //end dim for loop
1527 } // end if-block transform is BSpline
1532 // Method to reinitialize the seed of the random number generator
1533 template < class TFixedImage, class TMovingImage > void
1534 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1535 ::ReinitializeSeed()
1537 Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
1540 // Method to reinitialize the seed of the random number generator
1541 template < class TFixedImage, class TMovingImage > void
1542 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1543 ::ReinitializeSeed(int seed)
1545 Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(
1551 * Cache pre-transformed points, weights and indices.
1552 * This method is only called if the flag UseCachingOfBSplineWeights is ON.
1554 template < class TFixedImage, class TMovingImage >
1556 MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
1557 ::PreComputeTransformValues()
1559 // Create all zero dummy transform parameters
1560 ParametersType dummyParameters( this->m_Transform->GetNumberOfParameters() );
1561 dummyParameters.Fill( 0.0 );
1562 this->m_Transform->SetParameters( dummyParameters );
1564 // Cycle through each sampled fixed image point
1565 BSplineTransformWeightsType weights( m_NumBSplineWeights );
1566 BSplineTransformIndexArrayType indices( m_NumBSplineWeights );
1568 MovingImagePointType mappedPoint;
1570 // Declare iterators for iteration over the sample container
1571 typename FixedImageSpatialSampleContainer::const_iterator fiter;
1572 typename FixedImageSpatialSampleContainer::const_iterator fend =
1573 m_FixedImageSamples.end();
1574 unsigned long counter = 0;
1576 for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ ) {
1577 m_BSplineTransform->TransformPoint(
1578 m_FixedImageSamples[counter].FixedImagePointValue,
1579 mappedPoint, weights, indices, valid );
1581 for( unsigned long k = 0; k < m_NumBSplineWeights; k++ ) {
1582 m_BSplineTransformWeightsArray[counter][k] = weights[k];
1583 m_BSplineTransformIndicesArray[counter][k] = indices[k];
1586 m_PreTransformPointsArray[counter] = mappedPoint;
1587 m_WithinSupportRegionArray[counter] = valid;
1594 } // end namespace itk