+/*=========================================================================
+ 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://oncora1.lyon.fnclcc.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"
+#include "itkBSplineDeformableTransform.h"
+
+namespace itk
+{
+
+
+/**
+ * Constructor
+ */
+template < class TFixedImage, class TMovingImage >
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<TFixedImage,TMovingImage>
+::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 <class TFixedImage, class TMovingImage>
+void
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<double>::max();
+ double fixedImageMax = NumericTraits<double>::NonpositiveMin();
+
+ typedef ImageRegionConstIterator<FixedImageType> FixedIteratorType;
+ FixedIteratorType fixedImageIterator(
+ this->m_FixedImage, this->GetFixedImageRegion() );
+
+ for ( fixedImageIterator.GoToBegin();
+ !fixedImageIterator.IsAtEnd(); ++fixedImageIterator ) {
+
+ double sample = static_cast<double>( 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<double>::max();
+ double movingImageMax = NumericTraits<double>::NonpositiveMin();
+
+ typedef ImageRegionConstIterator<MovingImageType> MovingIteratorType;
+ MovingIteratorType movingImageIterator(
+ this->m_MovingImage, this->m_MovingImage->GetBufferedRegion() );
+
+ for ( movingImageIterator.GoToBegin();
+ !movingImageIterator.IsAtEnd(); ++movingImageIterator) {
+ double sample = static_cast<double>( 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<double>( m_NumberOfHistogramBins - 2 * padding );
+ m_FixedImageNormalizedMin = fixedImageMin / m_FixedImageBinSize -
+ static_cast<double>( padding );
+
+ m_MovingImageBinSize = ( movingImageMax - movingImageMin ) /
+ static_cast<double>( m_NumberOfHistogramBins - 2 * padding );
+ m_MovingImageNormalizedMin = movingImageMin / m_MovingImageBinSize -
+ static_cast<double>( 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<BSplineInterpolatorType *>(
+ 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<BSplineTransformType *>(
+ 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<TFixedImage,TMovingImage>
+::SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples )
+{
+
+ // Set up a random interator within the user specified fixed image region.
+ typedef ImageRandomConstIteratorWithIndex<FixedImageType> 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<TFixedImage,TMovingImage>
+::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer& samples )
+{
+
+ // Set up a region interator within the user specified fixed image region.
+ typedef ImageRegionConstIteratorWithIndex<FixedImageType> 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<TFixedImage,TMovingImage>
+::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<double>( (*iter).FixedImageValue ) / m_FixedImageBinSize -
+ m_FixedImageNormalizedMin;
+ unsigned int pindex = static_cast<unsigned int>( 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<TFixedImage,TMovingImage>
+::MeasureType
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<unsigned int>( 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<PDFValueType>( 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<int>( movingImageParzenWindowIndex ) - 1;
+ pdfPtr += pdfMovingIndex;
+
+ for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
+ + 2;
+ pdfMovingIndex++, pdfPtr++ ) {
+
+ // Update PDF for the current intensity pair
+ double movingImageParzenWindowArg =
+ static_cast<double>( pdfMovingIndex ) -
+ static_cast<double>( movingImageParzenWindowTerm );
+
+ *(pdfPtr) += static_cast<PDFValueType>(
+ 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<JointPDFType> 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<PDFValueType>( 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<PDFValueType>( fixedPDFSum );
+ }
+
+
+ // Compute moving image marginal PDF by summing over fixed image bins.
+ typedef ImageLinearIteratorWithIndex<JointPDFType> 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<PDFValueType>(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<MeasureType>( -1.0 * sum );
+
+}
+
+
+/**
+ * Get the both Value and Derivative Measure
+ */
+template < class TFixedImage, class TMovingImage >
+void
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<unsigned int>( 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<PDFValueType>( 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<int>( movingImageParzenWindowIndex ) - 1;
+ pdfPtr += pdfMovingIndex;
+
+ for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
+ + 2;
+ pdfMovingIndex++, pdfPtr++ ) {
+ // Update PDF for the current intensity pair
+ double movingImageParzenWindowArg =
+ static_cast<double>( pdfMovingIndex ) -
+ static_cast<double>(movingImageParzenWindowTerm);
+
+ *(pdfPtr) += static_cast<PDFValueType>(
+ 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<JointPDFType> 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<PDFValueType>( 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<PDFValueType>( fixedPDFSum );
+ }
+
+
+ // Compute moving image marginal PDF by summing over fixed image bins.
+ typedef ImageLinearIteratorWithIndex<JointPDFType> 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<PDFValueType>(sum);
+
+ linearIter.NextLine();
+ ++movingIndex1;
+
+ }
+
+ double nFactor = 1.0 / ( m_MovingImageBinSize
+ * static_cast<double>( nSamples ) );
+
+ if( this->m_UseExplicitPDFDerivatives ) {
+ // Normalize the joint PDF derivatives by the test image binsize and nSamples
+ typedef ImageRegionIterator<JointPDFDerivativesType>
+ 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<unsigned int>( 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<int>( movingImageParzenWindowIndex ) - 1;
+
+ for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
+ + 2;
+ pdfMovingIndex++ ) {
+
+ // Update PDF for the current intensity pair
+ double movingImageParzenWindowArg =
+ static_cast<double>( pdfMovingIndex ) -
+ static_cast<double>(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<MeasureType>( -1.0 * sum );
+}
+
+
+/**
+ * Get the match measure derivative
+ */
+template < class TFixedImage, class TMovingImage >
+void
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<TFixedImage,TMovingImage>
+::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<TFixedImage,TMovingImage>
+::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<TFixedImage,TMovingImage>
+::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;
+ const JacobianType& jacobian =
+ this->m_Transform->GetJacobian(
+ m_FixedImageSamples[sampleNumber].FixedImagePointValue );
+
+ 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 {
+ m_BSplineTransform->GetJacobian(
+ m_FixedImageSamples[sampleNumber].FixedImagePointValue, m_Weights, m_Indices );
+ }
+
+ 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<TFixedImage,TMovingImage>
+::ReinitializeSeed()
+{
+ Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
+}
+
+// Method to reinitialize the seed of the random number generator
+template < class TFixedImage, class TMovingImage > void
+MattesMutualInformationImageToImageMetricFor3DBLUTFFD<TFixedImage,TMovingImage>
+::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<TFixedImage,TMovingImage>
+::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