X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FitkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx;h=929d53452645837eb6cd70058b45de2b57248386;hb=HEAD;hp=1b38cae012c10569ca6d8add8bcf8ab490d6fb0c;hpb=c4de479fec231c7d53555dcd21d308f06aad17ec;p=clitk.git diff --git a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx index 1b38cae..929d534 100644 --- a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx +++ b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx @@ -41,1559 +41,5 @@ // 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