/*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv Authors belong to: - University of LYON http://www.universite-lyon.fr/ - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the copyright notices for more information. It is distributed under dual licence - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ===========================================================================**/ /*========================================================================= Program: Insight Segmentation & Registration Toolkit Module: $RCSfile: itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx,v $ Language: C++ Date: $Date: 2010/06/14 17:32:07 $ Version: $Revision: 1.1 $ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx #define __itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD_txx // First make sure that the configuration is available. // This line can be removed once the optimized versions // gets integrated into the main directories. #include "itkConfigure.h" #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS #include "itkOptMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx" #else #include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h" #include "itkBSplineInterpolateImageFunction.h" #include "itkCovariantVector.h" #include "itkImageRandomConstIteratorWithIndex.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionIterator.h" #include "itkImageIterator.h" #include "vnl/vnl_math.h" #if ITK_VERSION_MAJOR >= 4 #include "itkBSplineTransform.h" #else #include "itkBSplineDeformableTransform.h" #endif namespace itk { /** * Constructor */ template < class TFixedImage, class TMovingImage > MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::MattesMutualInformationImageToImageMetricFor3DBLUTFFD() { m_NumberOfSpatialSamples = 500; m_NumberOfHistogramBins = 50; this->SetComputeGradient(false); // don't use the default gradient for now m_InterpolatorIsBSpline = false; m_TransformIsBSpline = false; // Initialize PDFs to NULL m_JointPDF = NULL; m_JointPDFDerivatives = NULL; m_UseExplicitPDFDerivatives = true; typename BSplineTransformType::Pointer transformer = BSplineTransformType::New(); this->SetTransform (transformer); typename BSplineInterpolatorType::Pointer interpolator = BSplineInterpolatorType::New(); this->SetInterpolator (interpolator); // Initialize memory m_MovingImageNormalizedMin = 0.0; m_FixedImageNormalizedMin = 0.0; m_MovingImageTrueMin = 0.0; m_MovingImageTrueMax = 0.0; m_FixedImageBinSize = 0.0; m_MovingImageBinSize = 0.0; m_CubicBSplineDerivativeKernel = NULL; m_BSplineInterpolator = NULL; m_DerivativeCalculator = NULL; m_NumParametersPerDim = 0; m_NumBSplineWeights = 0; m_BSplineTransform = NULL; m_NumberOfParameters = 0; m_UseAllPixels = false; m_ReseedIterator = false; m_RandomSeed = -1; m_UseCachingOfBSplineWeights = true; } /** * Print out internal information about this class */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "NumberOfSpatialSamples: "; os << m_NumberOfSpatialSamples << std::endl; os << indent << "NumberOfHistogramBins: "; os << m_NumberOfHistogramBins << std::endl; os << indent << "UseAllPixels: "; os << m_UseAllPixels << std::endl; // Debugging information os << indent << "NumberOfParameters: "; os << m_NumberOfParameters << std::endl; os << indent << "FixedImageNormalizedMin: "; os << m_FixedImageNormalizedMin << std::endl; os << indent << "MovingImageNormalizedMin: "; os << m_MovingImageNormalizedMin << std::endl; os << indent << "MovingImageTrueMin: "; os << m_MovingImageTrueMin << std::endl; os << indent << "MovingImageTrueMax: "; os << m_MovingImageTrueMax << std::endl; os << indent << "FixedImageBinSize: "; os << m_FixedImageBinSize << std::endl; os << indent << "MovingImageBinSize: "; os << m_MovingImageBinSize << std::endl; os << indent << "InterpolatorIsBSpline: "; os << m_InterpolatorIsBSpline << std::endl; os << indent << "TransformIsBSpline: "; os << m_TransformIsBSpline << std::endl; os << indent << "UseCachingOfBSplineWeights: "; os << m_UseCachingOfBSplineWeights << std::endl; os << indent << "UseExplicitPDFDerivatives: "; os << m_UseExplicitPDFDerivatives << std::endl; } /** * Initialize */ template void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::Initialize(void) throw ( ExceptionObject ) { this->Superclass::Initialize(); // Cache the number of transformation parameters m_NumberOfParameters = this->m_Transform->GetNumberOfParameters(); /** * Compute the minimum and maximum for the FixedImage over * the FixedImageRegion. * * NB: We can't use StatisticsImageFilter to do this because * the filter computes the min/max for the largest possible region */ double fixedImageMin = NumericTraits::max(); double fixedImageMax = NumericTraits::NonpositiveMin(); typedef ImageRegionConstIterator FixedIteratorType; FixedIteratorType fixedImageIterator( this->m_FixedImage, this->GetFixedImageRegion() ); for ( fixedImageIterator.GoToBegin(); !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) { double sample = static_cast( fixedImageIterator.Get() ); if ( sample < fixedImageMin ) { fixedImageMin = sample; } if ( sample > fixedImageMax ) { fixedImageMax = sample; } } /** * Compute the minimum and maximum for the entire moving image * in the buffer. */ double movingImageMin = NumericTraits::max(); double movingImageMax = NumericTraits::NonpositiveMin(); typedef ImageRegionConstIterator MovingIteratorType; MovingIteratorType movingImageIterator( this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() ); for ( movingImageIterator.GoToBegin(); !movingImageIterator.IsAtEnd(); ++movingImageIterator) { double sample = static_cast( movingImageIterator.Get() ); if ( sample < movingImageMin ) { movingImageMin = sample; } if ( sample > movingImageMax ) { movingImageMax = sample; } } m_MovingImageTrueMin = movingImageMin; m_MovingImageTrueMax = movingImageMax; itkDebugMacro( " FixedImageMin: " << fixedImageMin << " FixedImageMax: " << fixedImageMax << std::endl ); itkDebugMacro( " MovingImageMin: " << movingImageMin << " MovingImageMax: " << movingImageMax << std::endl ); /** * Compute binsize for the histograms. * * The binsize for the image intensities needs to be adjusted so that * we can avoid dealing with boundary conditions using the cubic * spline as the Parzen window. We do this by increasing the size * of the bins so that the joint histogram becomes "padded" at the * borders. Because we are changing the binsize, * we also need to shift the minimum by the padded amount in order to * avoid minimum values filling in our padded region. * * Note that there can still be non-zero bin values in the padded region, * it's just that these bins will never be a central bin for the Parzen * window. * */ const int padding = 2; // this will pad by 2 bins m_FixedImageBinSize = ( fixedImageMax - fixedImageMin ) / static_cast( m_NumberOfHistogramBins - 2 * padding ); m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize - static_cast( padding ); m_MovingImageBinSize = ( movingImageMax - movingImageMin ) / static_cast( m_NumberOfHistogramBins - 2 * padding ); m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize - static_cast( padding ); itkDebugMacro( "FixedImageNormalizedMin: " << m_FixedImageNormalizedMin ); itkDebugMacro( "MovingImageNormalizedMin: " << m_MovingImageNormalizedMin ); itkDebugMacro( "FixedImageBinSize: " << m_FixedImageBinSize ); itkDebugMacro( "MovingImageBinSize; " << m_MovingImageBinSize ); if( m_UseAllPixels ) { m_NumberOfSpatialSamples = this->GetFixedImageRegion().GetNumberOfPixels(); } /** * Allocate memory for the fixed image sample container. */ m_FixedImageSamples.resize( m_NumberOfSpatialSamples ); /** * Allocate memory for the marginal PDF and initialize values * to zero. The marginal PDFs are stored as std::vector. */ m_FixedImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 ); m_MovingImageMarginalPDF.resize( m_NumberOfHistogramBins, 0.0 ); /** * Allocate memory for the joint PDF and joint PDF derivatives. * The joint PDF and joint PDF derivatives are store as itk::Image. */ m_JointPDF = JointPDFType::New(); // Instantiate a region, index, size JointPDFRegionType jointPDFRegion; JointPDFIndexType jointPDFIndex; JointPDFSizeType jointPDFSize; // Deallocate the memory that may have been allocated for // previous runs of the metric. this->m_JointPDFDerivatives = NULL; // by destroying the dynamic array this->m_PRatioArray.SetSize( 1, 1 ); // and by allocating very small the static ones this->m_MetricDerivative = DerivativeType( 1 ); // // Now allocate memory according to the user-selected method. // if( this->m_UseExplicitPDFDerivatives ) { this->m_JointPDFDerivatives = JointPDFDerivativesType::New(); JointPDFDerivativesRegionType jointPDFDerivativesRegion; JointPDFDerivativesIndexType jointPDFDerivativesIndex; JointPDFDerivativesSizeType jointPDFDerivativesSize; // For the derivatives of the joint PDF define a region starting from {0,0,0} // with size {m_NumberOfParameters,m_NumberOfHistogramBins, // m_NumberOfHistogramBins}. The dimension represents transform parameters, // fixed image parzen window index and moving image parzen window index, // respectively. jointPDFDerivativesIndex.Fill( 0 ); jointPDFDerivativesSize[0] = m_NumberOfParameters; jointPDFDerivativesSize[1] = m_NumberOfHistogramBins; jointPDFDerivativesSize[2] = m_NumberOfHistogramBins; jointPDFDerivativesRegion.SetIndex( jointPDFDerivativesIndex ); jointPDFDerivativesRegion.SetSize( jointPDFDerivativesSize ); // Set the regions and allocate m_JointPDFDerivatives->SetRegions( jointPDFDerivativesRegion ); m_JointPDFDerivatives->Allocate(); } else { /** Allocate memory for helper array that will contain the pRatios * for each bin of the joint histogram. This is part of the effort * for flattening the computation of the PDF Jacobians. */ this->m_PRatioArray.SetSize( this->m_NumberOfHistogramBins, this->m_NumberOfHistogramBins ); this->m_MetricDerivative = DerivativeType( this->GetNumberOfParameters() ); } // For the joint PDF define a region starting from {0,0} // with size {m_NumberOfHistogramBins, m_NumberOfHistogramBins}. // The dimension represents fixed image parzen window index // and moving image parzen window index, respectively. jointPDFIndex.Fill( 0 ); jointPDFSize.Fill( m_NumberOfHistogramBins ); jointPDFRegion.SetIndex( jointPDFIndex ); jointPDFRegion.SetSize( jointPDFSize ); // Set the regions and allocate m_JointPDF->SetRegions( jointPDFRegion ); m_JointPDF->Allocate(); /** * Setup the kernels used for the Parzen windows. */ m_CubicBSplineKernel = CubicBSplineFunctionType::New(); m_CubicBSplineDerivativeKernel = CubicBSplineDerivativeFunctionType::New(); if( m_UseAllPixels ) { /** * Take all the pixels within the fixed image region) * to create the sample points list. */ this->SampleFullFixedImageDomain( m_FixedImageSamples ); } else { /** * Uniformly sample the fixed image (within the fixed image region) * to create the sample points list. */ this->SampleFixedImageDomain( m_FixedImageSamples ); } /** * Pre-compute the fixed image parzen window index for * each point of the fixed image sample points list. */ this->ComputeFixedImageParzenWindowIndices( m_FixedImageSamples ); /** * Check if the interpolator is of type BSplineInterpolateImageFunction. * If so, we can make use of its EvaluateDerivatives method. * Otherwise, we instantiate an external central difference * derivative calculator. * * TODO: Also add it the possibility of using the default gradient * provided by the superclass. * */ m_InterpolatorIsBSpline = true; BSplineInterpolatorType * testPtr = dynamic_cast( this->m_Interpolator.GetPointer() ); if ( !testPtr ) { m_InterpolatorIsBSpline = false; m_DerivativeCalculator = DerivativeFunctionType::New(); #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION m_DerivativeCalculator->UseImageDirectionOn(); #endif m_DerivativeCalculator->SetInputImage( this->m_MovingImage ); m_BSplineInterpolator = NULL; itkDebugMacro( "Interpolator is not BSpline" ); } else { m_BSplineInterpolator = testPtr; #ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION m_BSplineInterpolator->UseImageDirectionOn(); #endif m_DerivativeCalculator = NULL; itkDebugMacro( "Interpolator is BSpline" ); } /** * Check if the transform is of type BSplineDeformableTransform. * * If so, several speed up features are implemented. * [1] Precomputing the results of bulk transform for each sample point. * [2] Precomputing the BSpline weights for each sample point, * to be used later to directly compute the deformation vector * [3] Precomputing the indices of the parameters within the * the support region of each sample point. */ m_TransformIsBSpline = true; BSplineTransformType * testPtr2 = dynamic_cast( this->m_Transform.GetPointer() ); if( !testPtr2 ) { m_TransformIsBSpline = false; m_BSplineTransform = NULL; itkDebugMacro( "Transform is not BSplineDeformable" ); } else { m_BSplineTransform = testPtr2; m_NumParametersPerDim = m_BSplineTransform->GetNumberOfParametersPerDimension(); m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights(); itkDebugMacro( "Transform is BSplineDeformable" ); } if( this->m_TransformIsBSpline ) { // First, deallocate memory that may have been used from previous run of the Metric this->m_BSplineTransformWeightsArray.SetSize( 1, 1 ); this->m_BSplineTransformIndicesArray.SetSize( 1, 1 ); this->m_PreTransformPointsArray.resize( 1 ); this->m_WithinSupportRegionArray.resize( 1 ); this->m_Weights.SetSize( 1 ); this->m_Indices.SetSize( 1 ); if( this->m_UseCachingOfBSplineWeights ) { m_BSplineTransformWeightsArray.SetSize( m_NumberOfSpatialSamples, m_NumBSplineWeights ); m_BSplineTransformIndicesArray.SetSize( m_NumberOfSpatialSamples, m_NumBSplineWeights ); m_PreTransformPointsArray.resize( m_NumberOfSpatialSamples ); m_WithinSupportRegionArray.resize( m_NumberOfSpatialSamples ); this->PreComputeTransformValues(); } else { this->m_Weights.SetSize( this->m_NumBSplineWeights ); this->m_Indices.SetSize( this->m_NumBSplineWeights ); } for ( unsigned int j = 0; j < FixedImageDimension; j++ ) { m_ParametersOffset[j] = j * m_BSplineTransform->GetNumberOfParametersPerDimension(); } } } /** * Uniformly sample the fixed image domain using a random walk */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples ) { // Set up a random interator within the user specified fixed image region. typedef ImageRandomConstIteratorWithIndex RandomIterator; RandomIterator randIter( this->m_FixedImage, this->GetFixedImageRegion() ); randIter.SetNumberOfSamples( m_NumberOfSpatialSamples ); randIter.GoToBegin(); typename FixedImageSpatialSampleContainer::iterator iter; typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); if( this->m_FixedImageMask ) { InputPointType inputPoint; iter=samples.begin(); int count = 0; int samples_found = 0; int maxcount = m_NumberOfSpatialSamples * 10; while( iter != end ) { if ( count > maxcount ) { #if 0 itkExceptionMacro( "Drew too many samples from the mask (is it too small?): " << maxcount << std::endl ); #else samples.resize(samples_found); // this->SetNumberOfSpatialSamples(sample_found); break; #endif } count++; // Get sampled index FixedImageIndexType index = randIter.GetIndex(); // Check if the Index is inside the mask, translate index to point this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); // If not inside the mask, ignore the point if( !this->m_FixedImageMask->IsInside( inputPoint ) ) { ++randIter; // jump to another random position continue; } // Get sampled fixed image value (*iter).FixedImageValue = randIter.Get(); // Translate index to point (*iter).FixedImagePointValue = inputPoint; samples_found++; // Jump to random position ++randIter; ++iter; } } else { for( iter=samples.begin(); iter != end; ++iter ) { // Get sampled index FixedImageIndexType index = randIter.GetIndex(); // Get sampled fixed image value (*iter).FixedImageValue = randIter.Get(); // Translate index to point this->m_FixedImage->TransformIndexToPhysicalPoint( index, (*iter).FixedImagePointValue ); // Jump to random position ++randIter; } } } /** * Sample the fixed image domain using all pixels in the Fixed image region */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples ) { // Set up a region interator within the user specified fixed image region. typedef ImageRegionConstIteratorWithIndex RegionIterator; RegionIterator regionIter( this->m_FixedImage, this->GetFixedImageRegion() ); regionIter.GoToBegin(); typename FixedImageSpatialSampleContainer::iterator iter; typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); if( this->m_FixedImageMask ) { InputPointType inputPoint; iter=samples.begin(); unsigned long nSamplesPicked = 0; while( iter != end && !regionIter.IsAtEnd() ) { // Get sampled index FixedImageIndexType index = regionIter.GetIndex(); // Check if the Index is inside the mask, translate index to point this->m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint ); // If not inside the mask, ignore the point if( !this->m_FixedImageMask->IsInside( inputPoint ) ) { ++regionIter; // jump to next pixel continue; } // Get sampled fixed image value (*iter).FixedImageValue = regionIter.Get(); // Translate index to point (*iter).FixedImagePointValue = inputPoint; ++regionIter; ++iter; ++nSamplesPicked; } // If we picked fewer samples than the desired number, // resize the container if (nSamplesPicked != this->m_NumberOfSpatialSamples) { this->m_NumberOfSpatialSamples = nSamplesPicked; samples.resize(this->m_NumberOfSpatialSamples); } } else { // not restricting sample throwing to a mask // cannot sample more than the number of pixels in the image region if ( this->m_NumberOfSpatialSamples > this->GetFixedImageRegion().GetNumberOfPixels()) { this->m_NumberOfSpatialSamples = this->GetFixedImageRegion().GetNumberOfPixels(); samples.resize(this->m_NumberOfSpatialSamples); } for( iter=samples.begin(); iter != end; ++iter ) { // Get sampled index FixedImageIndexType index = regionIter.GetIndex(); // Get sampled fixed image value (*iter).FixedImageValue = regionIter.Get(); // Translate index to point this->m_FixedImage->TransformIndexToPhysicalPoint( index, (*iter).FixedImagePointValue ); ++regionIter; } } } /** * Uniformly sample the fixed image domain using a random walk */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::ComputeFixedImageParzenWindowIndices( FixedImageSpatialSampleContainer& samples ) { typename FixedImageSpatialSampleContainer::iterator iter; typename FixedImageSpatialSampleContainer::const_iterator end=samples.end(); for( iter=samples.begin(); iter != end; ++iter ) { // Determine parzen window arguments (see eqn 6 of Mattes paper [2]). double windowTerm = static_cast( (*iter).FixedImageValue ) / m_FixedImageBinSize - m_FixedImageNormalizedMin; unsigned int pindex = static_cast( vcl_floor(windowTerm ) ); // Make sure the extreme values are in valid bins if ( pindex < 2 ) { pindex = 2; } else if ( pindex > (m_NumberOfHistogramBins - 3) ) { pindex = m_NumberOfHistogramBins - 3; } (*iter).FixedImageParzenWindowIndex = pindex; } } /** * Get the match Measure */ template < class TFixedImage, class TMovingImage > typename MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::MeasureType MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::GetValue( const ParametersType& parameters ) const { // Reset marginal pdf to all zeros. // Assumed the size has already been set to NumberOfHistogramBins // in Initialize(). for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) { m_FixedImageMarginalPDF[j] = 0.0; m_MovingImageMarginalPDF[j] = 0.0; } // Reset the joint pdfs to zero m_JointPDF->FillBuffer( 0.0 ); // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); // Declare iterators for iteration over the sample container typename FixedImageSpatialSampleContainer::const_iterator fiter; typename FixedImageSpatialSampleContainer::const_iterator fend = m_FixedImageSamples.end(); unsigned long nSamples=0; unsigned long nFixedImageSamples=0; for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { // Get moving image value MovingImagePointType mappedPoint; bool sampleOk; double movingImageValue; this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, sampleOk, movingImageValue ); ++nFixedImageSamples; if( sampleOk ) { ++nSamples; /** * Compute this sample's contribution to the marginal and * joint distributions. * */ // Determine parzen window arguments (see eqn 6 of Mattes paper [2]) double movingImageParzenWindowTerm = movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; unsigned int movingImageParzenWindowIndex = static_cast( vcl_floor(movingImageParzenWindowTerm ) ); // Make sure the extreme values are in valid bins if ( movingImageParzenWindowIndex < 2 ) { movingImageParzenWindowIndex = 2; } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; } // Since a zero-order BSpline (box car) kernel is used for // the fixed image marginal pdf, we need only increment the // fixedImageParzenWindowIndex by value of 1.0. m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] += static_cast( 1 ); /** * The region of support of the parzen window determines which bins * of the joint PDF are effected by the pair of image values. * Since we are using a cubic spline for the moving image parzen * window, four bins are affected. The fixed image parzen window is * a zero-order spline (box car) and thus effects only one bin. * * The PDF is arranged so that moving image bins corresponds to the * zero-th (column) dimension and the fixed image bins corresponds * to the first (row) dimension. * */ // Pointer to affected bin to be updated JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() + ( (*fiter).FixedImageParzenWindowIndex * m_JointPDF->GetOffsetTable()[1] ); // Move the pointer to the first affected bin int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; pdfPtr += pdfMovingIndex; for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) + 2; pdfMovingIndex++, pdfPtr++ ) { // Update PDF for the current intensity pair double movingImageParzenWindowArg = static_cast( pdfMovingIndex ) - static_cast( movingImageParzenWindowTerm ); *(pdfPtr) += static_cast( m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) ); } //end parzen windowing for loop } //end if-block check sampleOk } // end iterating over fixed image spatial sample container for loop itkDebugMacro( "Ratio of voxels mapping into moving image buffer: " << nSamples << " / " << m_NumberOfSpatialSamples << std::endl ); if( nSamples < m_NumberOfSpatialSamples / 16 ) { itkExceptionMacro( "Too many samples map outside moving image buffer: " << nSamples << " / " << m_NumberOfSpatialSamples << std::endl ); } this->m_NumberOfPixelsCounted = nSamples; /** * Normalize the PDFs, compute moving image marginal PDF * */ typedef ImageRegionIterator JointPDFIteratorType; JointPDFIteratorType jointPDFIterator ( m_JointPDF, m_JointPDF->GetBufferedRegion() ); jointPDFIterator.GoToBegin(); // Compute joint PDF normalization factor // (to ensure joint PDF sum adds to 1.0) double jointPDFSum = 0.0; while( !jointPDFIterator.IsAtEnd() ) { jointPDFSum += jointPDFIterator.Get(); ++jointPDFIterator; } if ( jointPDFSum == 0.0 ) { itkExceptionMacro( "Joint PDF summed to zero" ); } // Normalize the PDF bins jointPDFIterator.GoToEnd(); while( !jointPDFIterator.IsAtBegin() ) { --jointPDFIterator; jointPDFIterator.Value() /= static_cast( jointPDFSum ); } // Normalize the fixed image marginal PDF double fixedPDFSum = 0.0; for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) { fixedPDFSum += m_FixedImageMarginalPDF[bin]; } if ( fixedPDFSum == 0.0 ) { itkExceptionMacro( "Fixed image marginal PDF summed to zero" ); } for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) { m_FixedImageMarginalPDF[bin] /= static_cast( fixedPDFSum ); } // Compute moving image marginal PDF by summing over fixed image bins. typedef ImageLinearIteratorWithIndex JointPDFLinearIterator; JointPDFLinearIterator linearIter( m_JointPDF, m_JointPDF->GetBufferedRegion() ); linearIter.SetDirection( 1 ); linearIter.GoToBegin(); unsigned int movingIndex1 = 0; while( !linearIter.IsAtEnd() ) { double sum = 0.0; while( !linearIter.IsAtEndOfLine() ) { sum += linearIter.Get(); ++linearIter; } m_MovingImageMarginalPDF[movingIndex1] = static_cast(sum); linearIter.NextLine(); ++movingIndex1; } /** * Compute the metric by double summation over histogram. */ // Setup pointer to point to the first bin JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer(); // Initialze sum to zero double sum = 0.0; for( unsigned int fixedIndex = 0; fixedIndex < m_NumberOfHistogramBins; ++fixedIndex ) { double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex]; for( unsigned int movingIndex = 0; movingIndex < m_NumberOfHistogramBins; ++movingIndex, jointPDFPtr++ ) { double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex]; double jointPDFValue = *(jointPDFPtr); // check for non-zero bin contribution if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) { double pRatio = vcl_log(jointPDFValue / movingImagePDFValue ); if( fixedImagePDFValue > 1e-16) { sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) ); } } // end if-block to check non-zero bin contribution } // end for-loop over moving index } // end for-loop over fixed index return static_cast( -1.0 * sum ); } /** * Get the both Value and Derivative Measure */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::GetValueAndDerivative( const ParametersType& parameters, MeasureType& value, DerivativeType& derivative) const { // Set output values to zero value = NumericTraits< MeasureType >::Zero; if( this->m_UseExplicitPDFDerivatives ) { m_JointPDFDerivatives->FillBuffer( 0.0 ); derivative = DerivativeType( this->GetNumberOfParameters() ); derivative.Fill( NumericTraits< MeasureType >::Zero ); } else { this->m_MetricDerivative.Fill( NumericTraits< MeasureType >::Zero ); this->m_PRatioArray.Fill( 0.0 ); } // Reset marginal pdf to all zeros. // Assumed the size has already been set to NumberOfHistogramBins // in Initialize(). for ( unsigned int j = 0; j < m_NumberOfHistogramBins; j++ ) { m_FixedImageMarginalPDF[j] = 0.0; m_MovingImageMarginalPDF[j] = 0.0; } // Reset the joint pdfs to zero m_JointPDF->FillBuffer( 0.0 ); // Set up the parameters in the transform this->m_Transform->SetParameters( parameters ); // Declare iterators for iteration over the sample container typename FixedImageSpatialSampleContainer::const_iterator fiter; typename FixedImageSpatialSampleContainer::const_iterator fend = m_FixedImageSamples.end(); unsigned long nSamples=0; unsigned long nFixedImageSamples=0; for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { // Get moving image value MovingImagePointType mappedPoint; bool sampleOk; double movingImageValue; this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, sampleOk, movingImageValue ); if( sampleOk ) { ++nSamples; // Get moving image derivative at the mapped position ImageDerivativesType movingImageGradientValue; this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue ); /** * Compute this sample's contribution to the marginal * and joint distributions. * */ // Determine parzen window arguments (see eqn 6 of Mattes paper [2]) double movingImageParzenWindowTerm = movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; unsigned int movingImageParzenWindowIndex = static_cast( vcl_floor(movingImageParzenWindowTerm ) ); // Make sure the extreme values are in valid bins if ( movingImageParzenWindowIndex < 2 ) { movingImageParzenWindowIndex = 2; } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; } // Since a zero-order BSpline (box car) kernel is used for // the fixed image marginal pdf, we need only increment the // fixedImageParzenWindowIndex by value of 1.0. m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] += static_cast( 1 ); /** * The region of support of the parzen window determines which bins * of the joint PDF are effected by the pair of image values. * Since we are using a cubic spline for the moving image parzen * window, four bins are effected. The fixed image parzen window is * a zero-order spline (box car) and thus effects only one bin. * * The PDF is arranged so that moving image bins corresponds to the * zero-th (column) dimension and the fixed image bins corresponds * to the first (row) dimension. * */ // Pointer to affected bin to be updated JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() + ( (*fiter).FixedImageParzenWindowIndex * m_NumberOfHistogramBins ); // Move the pointer to the fist affected bin int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; pdfPtr += pdfMovingIndex; for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) + 2; pdfMovingIndex++, pdfPtr++ ) { // Update PDF for the current intensity pair double movingImageParzenWindowArg = static_cast( pdfMovingIndex ) - static_cast(movingImageParzenWindowTerm); *(pdfPtr) += static_cast( m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) ); if( this->m_UseExplicitPDFDerivatives ) { // Compute the cubicBSplineDerivative for later repeated use. double cubicBSplineDerivativeValue = m_CubicBSplineDerivativeKernel->Evaluate( movingImageParzenWindowArg ); // Compute PDF derivative contribution. this->ComputePDFDerivatives( nFixedImageSamples, pdfMovingIndex, movingImageGradientValue, cubicBSplineDerivativeValue ); } } //end parzen windowing for loop } //end if-block check sampleOk ++nFixedImageSamples; } // end iterating over fixed image spatial sample container for loop itkDebugMacro( "Ratio of voxels mapping into moving image buffer: " << nSamples << " / " << m_NumberOfSpatialSamples << std::endl ); if( nSamples < m_NumberOfSpatialSamples / 16 ) { itkExceptionMacro( "Too many samples map outside moving image buffer: " << nSamples << " / " << m_NumberOfSpatialSamples << std::endl ); } this->m_NumberOfPixelsCounted = nSamples; /** * Normalize the PDFs, compute moving image marginal PDF * */ typedef ImageRegionIterator JointPDFIteratorType; JointPDFIteratorType jointPDFIterator ( m_JointPDF, m_JointPDF->GetBufferedRegion() ); jointPDFIterator.GoToBegin(); // Compute joint PDF normalization factor // (to ensure joint PDF sum adds to 1.0) double jointPDFSum = 0.0; while( !jointPDFIterator.IsAtEnd() ) { jointPDFSum += jointPDFIterator.Get(); ++jointPDFIterator; } if ( jointPDFSum == 0.0 ) { itkExceptionMacro( "Joint PDF summed to zero" ); } // Normalize the PDF bins jointPDFIterator.GoToEnd(); while( !jointPDFIterator.IsAtBegin() ) { --jointPDFIterator; jointPDFIterator.Value() /= static_cast( jointPDFSum ); } // Normalize the fixed image marginal PDF double fixedPDFSum = 0.0; for( unsigned int bin = 0; bin < m_NumberOfHistogramBins; bin++ ) { fixedPDFSum += m_FixedImageMarginalPDF[bin]; } if ( fixedPDFSum == 0.0 ) { itkExceptionMacro( "Fixed image marginal PDF summed to zero" ); } for( unsigned int bin=0; bin < m_NumberOfHistogramBins; bin++ ) { m_FixedImageMarginalPDF[bin] /= static_cast( fixedPDFSum ); } // Compute moving image marginal PDF by summing over fixed image bins. typedef ImageLinearIteratorWithIndex JointPDFLinearIterator; JointPDFLinearIterator linearIter( m_JointPDF, m_JointPDF->GetBufferedRegion() ); linearIter.SetDirection( 1 ); linearIter.GoToBegin(); unsigned int movingIndex1 = 0; while( !linearIter.IsAtEnd() ) { double sum = 0.0; while( !linearIter.IsAtEndOfLine() ) { sum += linearIter.Get(); ++linearIter; } m_MovingImageMarginalPDF[movingIndex1] = static_cast(sum); linearIter.NextLine(); ++movingIndex1; } double nFactor = 1.0 / ( m_MovingImageBinSize * static_cast( nSamples ) ); if( this->m_UseExplicitPDFDerivatives ) { // Normalize the joint PDF derivatives by the test image binsize and nSamples typedef ImageRegionIterator JointPDFDerivativesIteratorType; JointPDFDerivativesIteratorType jointPDFDerivativesIterator ( m_JointPDFDerivatives, m_JointPDFDerivatives->GetBufferedRegion() ); jointPDFDerivativesIterator.GoToBegin(); while( !jointPDFDerivativesIterator.IsAtEnd() ) { jointPDFDerivativesIterator.Value() *= nFactor; ++jointPDFDerivativesIterator; } } /** * Compute the metric by double summation over histogram. */ // Setup pointer to point to the first bin JointPDFValueType * jointPDFPtr = m_JointPDF->GetBufferPointer(); // Initialize sum to zero double sum = 0.0; for( unsigned int fixedIndex = 0; fixedIndex < m_NumberOfHistogramBins; ++fixedIndex ) { double fixedImagePDFValue = m_FixedImageMarginalPDF[fixedIndex]; for( unsigned int movingIndex = 0; movingIndex < m_NumberOfHistogramBins; ++movingIndex, jointPDFPtr++ ) { double movingImagePDFValue = m_MovingImageMarginalPDF[movingIndex]; double jointPDFValue = *(jointPDFPtr); // check for non-zero bin contribution if( jointPDFValue > 1e-16 && movingImagePDFValue > 1e-16 ) { double pRatio = vcl_log(jointPDFValue / movingImagePDFValue ); if( fixedImagePDFValue > 1e-16) { sum += jointPDFValue * ( pRatio - vcl_log(fixedImagePDFValue ) ); } if( this->m_UseExplicitPDFDerivatives ) { // move joint pdf derivative pointer to the right position JointPDFValueType * derivPtr = m_JointPDFDerivatives->GetBufferPointer() + ( fixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] ) + ( movingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] ); for( unsigned int parameter=0; parameter < m_NumberOfParameters; ++parameter, derivPtr++ ) { // Ref: eqn 23 of Thevenaz & Unser paper [3] derivative[parameter] -= (*derivPtr) * pRatio; } // end for-loop over parameters } else { this->m_PRatioArray[fixedIndex][movingIndex] = pRatio * nFactor; } } // end if-block to check non-zero bin contribution } // end for-loop over moving index } // end for-loop over fixed index if( !(this->m_UseExplicitPDFDerivatives ) ) { // Second pass: This one is done for accumulating the contributions // to the derivative array. nFixedImageSamples = 0; for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter ) { // Get moving image value MovingImagePointType mappedPoint; bool sampleOk; double movingImageValue; this->TransformPoint( nFixedImageSamples, parameters, mappedPoint, sampleOk, movingImageValue ); if( sampleOk ) { // Get moving image derivative at the mapped position ImageDerivativesType movingImageGradientValue; this->ComputeImageDerivatives( mappedPoint, movingImageGradientValue ); /** * Compute this sample's contribution to the marginal * and joint distributions. * */ // Determine parzen window arguments (see eqn 6 of Mattes paper [2]). double movingImageParzenWindowTerm = movingImageValue / m_MovingImageBinSize - m_MovingImageNormalizedMin; unsigned int movingImageParzenWindowIndex = static_cast( vcl_floor(movingImageParzenWindowTerm ) ); // Make sure the extreme values are in valid bins if ( movingImageParzenWindowIndex < 2 ) { movingImageParzenWindowIndex = 2; } else if ( movingImageParzenWindowIndex > (m_NumberOfHistogramBins - 3) ) { movingImageParzenWindowIndex = m_NumberOfHistogramBins - 3; } // Move the pointer to the fist affected bin int pdfMovingIndex = static_cast( movingImageParzenWindowIndex ) - 1; for (; pdfMovingIndex <= static_cast( movingImageParzenWindowIndex ) + 2; pdfMovingIndex++ ) { // Update PDF for the current intensity pair double movingImageParzenWindowArg = static_cast( pdfMovingIndex ) - static_cast(movingImageParzenWindowTerm); // Compute the cubicBSplineDerivative for later repeated use. double cubicBSplineDerivativeValue = m_CubicBSplineDerivativeKernel->Evaluate( movingImageParzenWindowArg ); // Compute PDF derivative contribution. this->ComputePDFDerivatives( nFixedImageSamples, pdfMovingIndex, movingImageGradientValue, cubicBSplineDerivativeValue ); } //end parzen windowing for loop } //end if-block check sampleOk ++nFixedImageSamples; } // end iterating over fixed image spatial sample container for loop derivative = this->m_MetricDerivative; } value = static_cast( -1.0 * sum ); } /** * Get the match measure derivative */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::GetDerivative( const ParametersType& parameters, DerivativeType & derivative ) const { MeasureType value; // call the combined version this->GetValueAndDerivative( parameters, value, derivative ); } /** * Compute image derivatives using a central difference function * if we are not using a BSplineInterpolator, which includes * derivatives. */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::ComputeImageDerivatives( const MovingImagePointType& mappedPoint, ImageDerivativesType& gradient ) const { if( m_InterpolatorIsBSpline ) { // Computed moving image gradient using derivative BSpline kernel. gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint ); } else { // For all generic interpolator use central differencing. gradient = m_DerivativeCalculator->Evaluate( mappedPoint ); } } /** * Transform a point from FixedImage domain to MovingImage domain. * This function also checks if mapped point is within support region. */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::TransformPoint( unsigned int sampleNumber, const ParametersType& parameters, MovingImagePointType& mappedPoint, bool& sampleOk, double& movingImageValue ) const { if ( !m_TransformIsBSpline ) { // Use generic transform to compute mapped position mappedPoint = this->m_Transform->TransformPoint( m_FixedImageSamples[sampleNumber].FixedImagePointValue ); // Check if mapped point inside image buffer sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint ); } else { if( this->m_UseCachingOfBSplineWeights ) { // If the transform is BSplineDeformable, we can use the precomputed // weights and indices to obtained the mapped position // const WeightsValueType * weights = m_BSplineTransformWeightsArray[sampleNumber]; const IndexValueType * indices = m_BSplineTransformIndicesArray[sampleNumber]; mappedPoint.Fill( 0.0 ); if ( m_WithinSupportRegionArray[sampleNumber] ) { for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ ) { for ( unsigned int j = 0; j < FixedImageDimension; j++ ) { mappedPoint[j] += weights[k] * parameters[ indices[k] + m_ParametersOffset[j] ]; } } } for( unsigned int j = 0; j < FixedImageDimension; j++ ) { mappedPoint[j] += m_PreTransformPointsArray[sampleNumber][j]; } // Check if mapped point inside image buffer sampleOk = this->m_Interpolator->IsInsideBuffer( mappedPoint ); // Check if mapped point is within the support region of a grid point. // This is neccessary for computing the metric gradient sampleOk = sampleOk && m_WithinSupportRegionArray[sampleNumber]; } else { // If not caching values, we invoke the Transform to recompute the // mapping of the point. this->m_BSplineTransform->TransformPoint( this->m_FixedImageSamples[sampleNumber].FixedImagePointValue, mappedPoint, this->m_Weights, this->m_Indices, sampleOk); // Check if mapped point inside image buffer sampleOk = sampleOk && this->m_Interpolator->IsInsideBuffer( mappedPoint ); } } // If user provided a mask over the Moving image if ( this->m_MovingImageMask ) { // Check if mapped point is within the support region of the moving image // mask sampleOk = sampleOk && this->m_MovingImageMask->IsInside( mappedPoint ); } if ( sampleOk ) { movingImageValue = this->m_Interpolator->Evaluate( mappedPoint ); if ( movingImageValue < m_MovingImageTrueMin || movingImageValue > m_MovingImageTrueMax ) { // need to throw out this sample as it will not fall into a valid bin sampleOk = false; } } } /** * Compute PDF derivatives contribution for each parameter */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::ComputePDFDerivatives( unsigned int sampleNumber, int pdfMovingIndex, const ImageDerivativesType& movingImageGradientValue, double cubicBSplineDerivativeValue ) const { const int pdfFixedIndex = m_FixedImageSamples[sampleNumber].FixedImageParzenWindowIndex; JointPDFValueType * derivPtr = NULL; double precomputedWeight = 0.0; if( this->m_UseExplicitPDFDerivatives ) { // Update bins in the PDF derivatives for the current intensity pair derivPtr = m_JointPDFDerivatives->GetBufferPointer() + ( pdfFixedIndex * m_JointPDFDerivatives->GetOffsetTable()[2] ) + ( pdfMovingIndex * m_JointPDFDerivatives->GetOffsetTable()[1] ); } else { // Recover the precomputed weight for this specific PDF bin precomputedWeight = this->m_PRatioArray[pdfFixedIndex][pdfMovingIndex]; } if( !m_TransformIsBSpline ) { /** * Generic version which works for all transforms. */ // Compute the transform Jacobian. typedef typename TransformType::JacobianType JacobianType; #if ITK_VERSION_MAJOR >= 4 JacobianType jacobian; this->m_Transform->ComputeJacobianWithRespectToParameters( m_FixedImageSamples[sampleNumber].FixedImagePointValue, jacobian ); #else const JacobianType & jacobian = this->m_Transform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue ); #endif for ( unsigned int mu = 0; mu < m_NumberOfParameters; mu++ ) { double innerProduct = 0.0; for ( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) { innerProduct += jacobian[dim][mu] * movingImageGradientValue[dim]; } const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue; if( this->m_UseExplicitPDFDerivatives ) { *(derivPtr) -= derivativeContribution; ++derivPtr; } else { this->m_MetricDerivative[mu] += precomputedWeight * derivativeContribution; } } } else { const WeightsValueType * weights = NULL; const IndexValueType * indices = NULL; if( this->m_UseCachingOfBSplineWeights ) { /** * If the transform is of type BSplineDeformableTransform, * we can obtain a speed up by only processing the affected parameters. * Note that these pointers are just pointing to pre-allocated rows * of the caching arrays. There is therefore, no need to free this * memory. */ weights = m_BSplineTransformWeightsArray[sampleNumber]; indices = m_BSplineTransformIndicesArray[sampleNumber]; } else { #if ITK_VERSION_MAJOR >= 4 m_BSplineTransform->ComputeJacobianFromBSplineWeightsWithRespectToPosition( m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices ); #else m_BSplineTransform->GetJacobian( m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices ); #endif } for( unsigned int dim = 0; dim < FixedImageDimension; dim++ ) { double innerProduct; int parameterIndex; for( unsigned int mu = 0; mu < m_NumBSplineWeights; mu++ ) { /* The array weights contains the Jacobian values in a 1-D array * (because for each parameter the Jacobian is non-zero in only 1 of the * possible dimensions) which is multiplied by the moving image * gradient. */ if( this->m_UseCachingOfBSplineWeights ) { innerProduct = movingImageGradientValue[dim] * weights[mu]; parameterIndex = indices[mu] + m_ParametersOffset[dim]; } else { innerProduct = movingImageGradientValue[dim] * this->m_Weights[mu]; parameterIndex = this->m_Indices[mu] + this->m_ParametersOffset[dim]; } const double derivativeContribution = innerProduct * cubicBSplineDerivativeValue; if( this->m_UseExplicitPDFDerivatives ) { JointPDFValueType * ptr = derivPtr + parameterIndex; *(ptr) -= derivativeContribution; } else { this->m_MetricDerivative[parameterIndex] += precomputedWeight * derivativeContribution; } } //end mu for loop } //end dim for loop } // end if-block transform is BSpline } // Method to reinitialize the seed of the random number generator template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::ReinitializeSeed() { Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(); } // Method to reinitialize the seed of the random number generator template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::ReinitializeSeed(int seed) { Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed( seed); } /** * Cache pre-transformed points, weights and indices. * This method is only called if the flag UseCachingOfBSplineWeights is ON. */ template < class TFixedImage, class TMovingImage > void MattesMutualInformationImageToImageMetricFor3DBLUTFFD ::PreComputeTransformValues() { // Create all zero dummy transform parameters ParametersType dummyParameters( this->m_Transform->GetNumberOfParameters() ); dummyParameters.Fill( 0.0 ); this->m_Transform->SetParameters( dummyParameters ); // Cycle through each sampled fixed image point BSplineTransformWeightsType weights( m_NumBSplineWeights ); BSplineTransformIndexArrayType indices( m_NumBSplineWeights ); bool valid; MovingImagePointType mappedPoint; // Declare iterators for iteration over the sample container typename FixedImageSpatialSampleContainer::const_iterator fiter; typename FixedImageSpatialSampleContainer::const_iterator fend = m_FixedImageSamples.end(); unsigned long counter = 0; for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ ) { m_BSplineTransform->TransformPoint( m_FixedImageSamples[counter].FixedImagePointValue, mappedPoint, weights, indices, valid ); for( unsigned long k = 0; k < m_NumBSplineWeights; k++ ) { m_BSplineTransformWeightsArray[counter][k] = weights[k]; m_BSplineTransformIndicesArray[counter][k] = indices[k]; } m_PreTransformPointsArray[counter] = mappedPoint; m_WithinSupportRegionArray[counter] = valid; } } } // end namespace itk #endif #endif