X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FitkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h;h=5342fea88784ff665a45c30be9abfb53ac662994;hb=HEAD;hp=091e36cb57f5eac7a56e58a381475806e97a88e2;hpb=cff1bb0cd9b0850379e20c0594ccd511a3a38133;p=clitk.git diff --git a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h index 091e36c..5342fea 100644 --- a/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h +++ b/registration/itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.h @@ -39,497 +39,5 @@ // 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.h" -#else - -#include "itkImageToImageMetric.h" -#include "itkCovariantVector.h" -#include "itkPoint.h" -#include "itkIndex.h" -#include "itkBSplineKernelFunction.h" -#include "itkBSplineDerivativeKernelFunction.h" -#include "itkCentralDifferenceImageFunction.h" -#include "itkBSplineInterpolateImageFunction.h" -#include "itkBSplineDeformableTransform.h" -#include "itkArray2D.h" - -namespace itk -{ - -/** \class MattesMutualInformationImageToImageMetricFor3DBLUTFFD - * \brief Computes the mutual information between two images to be - * registered using the method of Mattes et al. - * - * MattesMutualInformationImageToImageMetricFor3DBLUTFFD computes the mutual - * information between a fixed and moving image to be registered. - * - * This class is templated over the FixedImage type and the MovingImage - * type. - * - * The fixed and moving images are set via methods SetFixedImage() and - * SetMovingImage(). This metric makes use of user specified Transform and - * Interpolator. The Transform is used to map points from the fixed image to - * the moving image domain. The Interpolator is used to evaluate the image - * intensity at user specified geometric points in the moving image. - * The Transform and Interpolator are set via methods SetTransform() and - * SetInterpolator(). - * - * If a BSplineInterpolationFunction is used, this class obtain - * image derivatives from the BSpline interpolator. Otherwise, - * image derivatives are computed using central differencing. - * - * \warning This metric assumes that the moving image has already been - * connected to the interpolator outside of this class. - * - * The method GetValue() computes of the mutual information - * while method GetValueAndDerivative() computes - * both the mutual information and its derivatives with respect to the - * transform parameters. - * - * The calculations are based on the method of Mattes et al [1,2] - * where the probability density distribution are estimated using - * Parzen histograms. Since the fixed image PDF does not contribute - * to the derivatives, it does not need to be smooth. Hence, - * a zero order (box car) BSpline kernel is used - * for the fixed image intensity PDF. On the other hand, to ensure - * smoothness a third order BSpline kernel is used for the - * moving image intensity PDF. - * - * On Initialize(), the FixedImage is uniformly sampled within - * the FixedImageRegion. The number of samples used can be set - * via SetNumberOfSpatialSamples(). Typically, the number of - * spatial samples used should increase with the image size. - * - * The option UseAllPixelOn() disables the random sampling and uses - * all the pixels of the FixedImageRegion in order to estimate the - * joint intensity PDF. - * - * During each call of GetValue(), GetDerivatives(), - * GetValueAndDerivatives(), marginal and joint intensity PDF's - * values are estimated at discrete position or bins. - * The number of bins used can be set via SetNumberOfHistogramBins(). - * To handle data with arbitray magnitude and dynamic range, - * the image intensity is scale such that any contribution to the - * histogram will fall into a valid bin. - * - * One the PDF's have been contructed, the mutual information - * is obtained by doubling summing over the discrete PDF values. - * - * - * Notes: - * 1. This class returns the negative mutual information value. - * 2. This class in not thread safe due the private data structures - * used to the store the sampled points and the marginal and joint pdfs. - * - * References: - * [1] "Nonrigid multimodality image registration" - * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank - * Medical Imaging 2001: Image Processing, 2001, pp. 1609-1620. - * [2] "PET-CT Image Registration in the Chest Using Free-form Deformations" - * D. Mattes, D. R. Haynor, H. Vesselle, T. Lewellen and W. Eubank - * IEEE Transactions in Medical Imaging. Vol.22, No.1, - January 2003. pp.120-128. - * [3] "Optimization of Mutual Information for MultiResolution Image - * Registration" - * P. Thevenaz and M. Unser - * IEEE Transactions in Image Processing, 9(12) December 2000. - * - * \ingroup RegistrationMetrics - * \ingroup ThreadUnSafe - */ -template -class ITK_EXPORT MattesMutualInformationImageToImageMetricFor3DBLUTFFD : - public ImageToImageMetric< TFixedImage, TMovingImage > -{ -public: - - /** Standard class typedefs. */ - typedef MattesMutualInformationImageToImageMetricFor3DBLUTFFD Self; - typedef ImageToImageMetric< TFixedImage, TMovingImage > Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(MattesMutualInformationImageToImageMetricFor3DBLUTFFD, ImageToImageMetric); - - /** Types inherited from Superclass. */ - typedef typename Superclass::TransformType TransformType; - typedef typename Superclass::TransformPointer TransformPointer; - typedef typename Superclass::TransformJacobianType TransformJacobianType; - typedef typename Superclass::InterpolatorType InterpolatorType; - typedef typename Superclass::MeasureType MeasureType; - typedef typename Superclass::DerivativeType DerivativeType; - typedef typename Superclass::ParametersType ParametersType; - typedef typename Superclass::FixedImageType FixedImageType; - typedef typename Superclass::MovingImageType MovingImageType; - typedef typename Superclass::FixedImageConstPointer FixedImageConstPointer; - typedef typename Superclass::MovingImageConstPointer MovingImageCosntPointer; - typedef typename Superclass::InputPointType InputPointType; - typedef typename Superclass::OutputPointType OutputPointType; - - typedef typename Superclass::CoordinateRepresentationType - CoordinateRepresentationType; - - /** Index and Point typedef support. */ - typedef typename FixedImageType::IndexType FixedImageIndexType; - typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType; - typedef typename MovingImageType::IndexType MovingImageIndexType; - typedef typename TransformType::InputPointType FixedImagePointType; - typedef typename TransformType::OutputPointType MovingImagePointType; - - /** The moving image dimension. */ - itkStaticConstMacro( MovingImageDimension, unsigned int, - MovingImageType::ImageDimension ); - - /** - * Initialize the Metric by - * (1) making sure that all the components are present and plugged - * together correctly, - * (2) uniformly select NumberOfSpatialSamples within - * the FixedImageRegion, and - * (3) allocate memory for pdf data structures. */ - virtual void Initialize(void) throw ( ExceptionObject ); - - /** Get the derivatives of the match measure. */ - void GetDerivative( const ParametersType& parameters, - DerivativeType & Derivative ) const; - - /** Get the value. */ - MeasureType GetValue( const ParametersType& parameters ) const; - - /** Get the value and derivatives for single valued optimizers. */ - void GetValueAndDerivative( const ParametersType& parameters, - MeasureType& Value, - DerivativeType& Derivative ) const; - - /** Number of spatial samples to used to compute metric */ - itkSetClampMacro( NumberOfSpatialSamples, unsigned long, - 1, NumericTraits::max() ); - itkGetConstReferenceMacro( NumberOfSpatialSamples, unsigned long); - - /** Number of bins to used in the histogram. Typical value is 50. */ - itkSetClampMacro( NumberOfHistogramBins, unsigned long, - 1, NumericTraits::max() ); - itkGetConstReferenceMacro( NumberOfHistogramBins, unsigned long); - - /** Reinitialize the seed of the random number generator that selects the - * sample of pixels used for estimating the image histograms and the joint - * histogram. By nature, this metric is not deterministic, since at each run - * it may select a different set of pixels. By initializing the random number - * generator seed to the same value you can restore determinism. On the other - * hand, calling the method ReinitializeSeed() without arguments will use the - * clock from your machine in order to have a very random initialization of - * the seed. This will indeed increase the non-deterministic behavior of the - * metric. */ - void ReinitializeSeed(); - void ReinitializeSeed(int); - - /** Select whether the metric will be computed using all the pixels on the - * fixed image region, or only using a set of randomly selected pixels. */ - itkSetMacro(UseAllPixels,bool); - itkGetConstReferenceMacro(UseAllPixels,bool); - itkBooleanMacro(UseAllPixels); - - /** This variable selects the method to be used for computing the Metric - * derivatives with respect to the Transform parameters. Two modes of - * computation are available. The choice between one and the other is a - * trade-off between computation speed and memory allocations. The two modes - * are described in detail below: - * - * UseExplicitPDFDerivatives = True - * will compute the Metric derivative by first calculating the derivatives of - * each one of the Joint PDF bins with respect to each one of the Transform - * parameters and then accumulating these contributions in the final metric - * derivative array by using a bin-specific weight. The memory required for - * storing the intermediate derivatives is a 3D array of doubles with size - * equals to the product of (number of histogram bins)^2 times number of - * transform parameters. This method is well suited for Transform with a small - * number of parameters. - * - * UseExplicitPDFDerivatives = False will compute the Metric derivative by - * first computing the weights for each one of the Joint PDF bins and caching - * them into an array. Then it will revisit each one of the PDF bins for - * computing its weighted contribution to the full derivative array. In this - * method an extra 2D array is used for storing the weights of each one of - * the PDF bins. This is an array of doubles with size equals to (number of - * histogram bins)^2. This method is well suited for Transforms with a large - * number of parameters, such as, BSplineDeformableTransforms. */ - itkSetMacro(UseExplicitPDFDerivatives,bool); - itkGetConstReferenceMacro(UseExplicitPDFDerivatives,bool); - itkBooleanMacro(UseExplicitPDFDerivatives); - - /** This boolean flag is only relevant when this metric is used along - * with a BSplineDeformableTransform. The flag enables/disables the - * caching of values computed when a physical point is mapped through - * the BSplineDeformableTransform. In particular it will cache the - * values of the BSpline weights for that points, and the indexes - * indicating what BSpline-grid nodes are relevant for that specific - * point. This caching is made optional due to the fact that the - * memory arrays used for the caching can reach large sizes even for - * moderate image size problems. For example, for a 3D image of - * 256^3, using 20% of pixels, these arrays will take about 1 - * Gigabyte of RAM for storage. The ratio of computing time between - * using the cache or not using the cache can reach 1:5, meaning that - * using the caching can provide a five times speed up. It is - * therefore, interesting to enable the caching, if enough memory is - * available for it. The caching is enabled by default, in order to - * preserve backward compatibility with previous versions of ITK. */ - itkSetMacro(UseCachingOfBSplineWeights,bool); - itkGetConstReferenceMacro(UseCachingOfBSplineWeights,bool); - itkBooleanMacro(UseCachingOfBSplineWeights); - -protected: - - MattesMutualInformationImageToImageMetricFor3DBLUTFFD(); - virtual ~MattesMutualInformationImageToImageMetricFor3DBLUTFFD() {}; - void PrintSelf(std::ostream& os, Indent indent) const; - - /** \class FixedImageSpatialSample - * A fixed image spatial sample consists of the fixed domain point - * and the fixed image value at that point. */ - /// @cond - class FixedImageSpatialSample - { - public: - FixedImageSpatialSample():FixedImageValue(0.0) { - FixedImagePointValue.Fill(0.0); - } - ~FixedImageSpatialSample() {}; - - FixedImagePointType FixedImagePointValue; - double FixedImageValue; - unsigned int FixedImageParzenWindowIndex; - }; - /// @endcond - - /** FixedImageSpatialSample typedef support. */ - typedef std::vector - FixedImageSpatialSampleContainer; - - /** Container to store a set of points and fixed image values. */ - FixedImageSpatialSampleContainer m_FixedImageSamples; - - /** Uniformly select a sample set from the fixed image domain. */ - virtual void SampleFixedImageDomain( - FixedImageSpatialSampleContainer& samples); - - /** Gather all the pixels from the fixed image domain. */ - virtual void SampleFullFixedImageDomain( - FixedImageSpatialSampleContainer& samples); - - /** Transform a point from FixedImage domain to MovingImage domain. - * This function also checks if mapped point is within support region. */ - virtual void TransformPoint( unsigned int sampleNumber, - const ParametersType& parameters, - MovingImagePointType& mappedPoint, - bool& sampleWithinSupportRegion, - double& movingImageValue ) const; - -private: - - //purposely not implemented - MattesMutualInformationImageToImageMetricFor3DBLUTFFD(const Self&); - //purposely not implemented - void operator=(const Self&); - - - /** The marginal PDFs are stored as std::vector. */ - typedef float PDFValueType; - typedef std::vector MarginalPDFType; - - /** The fixed image marginal PDF. */ - mutable MarginalPDFType m_FixedImageMarginalPDF; - - /** The moving image marginal PDF. */ - mutable MarginalPDFType m_MovingImageMarginalPDF; - - /** Helper array for storing the values of the JointPDF ratios. */ - typedef double PRatioType; - typedef Array2D< PRatioType > PRatioArrayType; - mutable PRatioArrayType m_PRatioArray; - - /** Helper variable for accumulating the derivative of the metric. */ - mutable DerivativeType m_MetricDerivative; - - /** Typedef for the joint PDF and PDF derivatives are stored as ITK Images. */ - typedef Image JointPDFType; - typedef JointPDFType::IndexType JointPDFIndexType; - typedef JointPDFType::PixelType JointPDFValueType; - typedef JointPDFType::RegionType JointPDFRegionType; - typedef JointPDFType::SizeType JointPDFSizeType; - typedef Image JointPDFDerivativesType; - typedef JointPDFDerivativesType::IndexType JointPDFDerivativesIndexType; - typedef JointPDFDerivativesType::PixelType JointPDFDerivativesValueType; - typedef JointPDFDerivativesType::RegionType JointPDFDerivativesRegionType; - typedef JointPDFDerivativesType::SizeType JointPDFDerivativesSizeType; - - /** The joint PDF and PDF derivatives. */ - typename JointPDFType::Pointer m_JointPDF; - typename JointPDFDerivativesType::Pointer m_JointPDFDerivatives; - - unsigned long m_NumberOfSpatialSamples; - unsigned long m_NumberOfParameters; - - /** Variables to define the marginal and joint histograms. */ - unsigned long m_NumberOfHistogramBins; - double m_MovingImageNormalizedMin; - double m_FixedImageNormalizedMin; - double m_MovingImageTrueMin; - double m_MovingImageTrueMax; - double m_FixedImageBinSize; - double m_MovingImageBinSize; - - /** Typedefs for BSpline kernel and derivative functions. */ - typedef BSplineKernelFunction<3> CubicBSplineFunctionType; - typedef BSplineDerivativeKernelFunction<3> CubicBSplineDerivativeFunctionType; - - /** Cubic BSpline kernel for computing Parzen histograms. */ - typename CubicBSplineFunctionType::Pointer m_CubicBSplineKernel; - typename CubicBSplineDerivativeFunctionType::Pointer - m_CubicBSplineDerivativeKernel; - - /** Precompute fixed image parzen window indices. */ - virtual void ComputeFixedImageParzenWindowIndices( - FixedImageSpatialSampleContainer& samples ); - - /** - * Types and variables related to image derivative calculations. - * If a BSplineInterpolationFunction is used, this class obtain - * image derivatives from the BSpline interpolator. Otherwise, - * image derivatives are computed using central differencing. - */ - typedef CovariantVector< double, - itkGetStaticConstMacro(MovingImageDimension) > - ImageDerivativesType; - - /** Compute image derivatives at a point. */ - virtual void ComputeImageDerivatives( const MovingImagePointType& mappedPoint, - ImageDerivativesType& gradient ) const; - - /** Boolean to indicate if the interpolator BSpline. */ - bool m_InterpolatorIsBSpline; - - /** Typedefs for using BSpline interpolator. */ - typedef - BSplineInterpolateImageFunction - BSplineInterpolatorType; - - /** Pointer to BSplineInterpolator. */ - typename BSplineInterpolatorType::Pointer m_BSplineInterpolator; - - /** Typedefs for using central difference calculator. */ - typedef CentralDifferenceImageFunction - DerivativeFunctionType; - - /** Pointer to central difference calculator. */ - typename DerivativeFunctionType::Pointer m_DerivativeCalculator; - - - /** Compute PDF derivative contribution for each parameter. */ - virtual void ComputePDFDerivatives( unsigned int sampleNumber, - int movingImageParzenWindowIndex, - const ImageDerivativesType& - movingImageGradientValue, - double cubicBSplineDerivativeValue - ) const; - - /** - * Types and variables related to BSpline deformable transforms. - * If the transform is of type third order BSplineDeformableTransform, - * then we can speed up the metric derivative calculation by - * only inspecting the parameters within the support region - * of a mapped point. - */ - - /** Boolean to indicate if the transform is BSpline deformable. */ - bool m_TransformIsBSpline; - - /** The number of BSpline parameters per image dimension. */ - long m_NumParametersPerDim; - - /** - * The number of BSpline transform weights is the number of - * of parameter in the support region (per dimension ). */ - unsigned long m_NumBSplineWeights; - - /** The fixed image dimension. */ - itkStaticConstMacro( FixedImageDimension, unsigned int, - FixedImageType::ImageDimension ); - - /** - * Enum of the deformabtion field spline order. - */ - enum { DeformationSplineOrder = 3 }; - - /** - * Typedefs for the BSplineDeformableTransform. - */ - typedef BSplineDeformableTransform< - CoordinateRepresentationType, - ::itk::GetImageDimension::ImageDimension, - DeformationSplineOrder> BSplineTransformType; - typedef typename BSplineTransformType::WeightsType - BSplineTransformWeightsType; - typedef typename BSplineTransformType::ParameterIndexArrayType - BSplineTransformIndexArrayType; - - /** - * Variables used when transform is of type BSpline deformable. - */ - typename BSplineTransformType::Pointer m_BSplineTransform; - - /** - * Cache pre-transformed points, weights, indices and - * within support region flag. - */ - typedef typename BSplineTransformWeightsType::ValueType WeightsValueType; - typedef Array2D BSplineTransformWeightsArrayType; - typedef typename BSplineTransformIndexArrayType::ValueType IndexValueType; - typedef Array2D BSplineTransformIndicesArrayType; - typedef std::vector MovingImagePointArrayType; - typedef std::vector BooleanArrayType; - - BSplineTransformWeightsArrayType m_BSplineTransformWeightsArray; - BSplineTransformIndicesArrayType m_BSplineTransformIndicesArray; - MovingImagePointArrayType m_PreTransformPointsArray; - BooleanArrayType m_WithinSupportRegionArray; - - typedef FixedArray::ImageDimension> - ParametersOffsetType; - ParametersOffsetType m_ParametersOffset; - - bool m_UseAllPixels; - - virtual void PreComputeTransformValues(); - - bool m_ReseedIterator; - int m_RandomSeed; - - // Selection of explicit or implicit computation of PDF derivatives - // with respect to Transform parameters. - bool m_UseExplicitPDFDerivatives; - - // Variables needed for optionally caching values when using a BSpline transform. - bool m_UseCachingOfBSplineWeights; - mutable BSplineTransformWeightsType m_Weights; - mutable BSplineTransformIndexArrayType m_Indices; - -}; - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkMattesMutualInformationImageToImageMetricFor3DBLUTFFD.txx" -#endif - -#endif - #endif