/*========================================================================= Program: Insight Segmentation & Registration Toolkit Module: $RCSfile: itkAnytimeImageToImageMetric2.h,v $ Language: C++ Date: $Date: 2007/09/28 21:38:43 $ Version: $Revision: 1.18 $ 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 __itkAnytimeImageToImageMetric2_h #define __itkAnytimeImageToImageMetric2_h #include "itkAnytimeImageToImageMetric.h" #include "itkAnytimeImageToImageMetric2Base.h" #include "itkArray2D.h" #include "itkContinuousIndex.h" #include "itkTransformBase.h" #include "itkCentralDifferenceImageFunction.h" #include "itkBSplineInterpolateImageFunction.h" #include "itkResampleImageFilter.h" namespace itk { /** \class AnytimeImageToImageMetric2 * \brief Computes similarity between regions of two images. * * This Class is templated over the type of the two input images. * It expects a Transform and an Interpolator to be plugged in. * This particular class is the base class for a hierarchy of * similarity metrics. * * This class computes a value that measures the similarity * between the Fixed image and the transformed Moving image. * The Interpolator is used to compute intensity values on * non-grid positions resulting from mapping points through * the Transform. * * * \ingroup RegistrationMetrics * */ template class ITK_EXPORT AnytimeImageToImageMetric2 : virtual public AnytimeImageToImageMetric2Base, public AnytimeImageToImageMetric { public: /** Standard class typedefs. */ typedef AnytimeImageToImageMetric2 Self; typedef AnytimeImageToImageMetric Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Type used for representing point components */ typedef typename Superclass::ParametersValueType CoordinateRepresentationType; typedef itk::Image PriorityImageType; /** Run-time type information (and related methods). */ itkTypeMacro(AnytimeImageToImageMetric2, AnytimeImageToImageMetric); /** Have to explicitly identify the methods from the base class that i want to overload by name but not signature. sigh... C++... sheesh */ using Superclass::GetValueAndDerivative; using Superclass::GetValue; using Superclass::GetDerivative; unsigned int m_NumberOfChannels; mutable unsigned long m_TotalPixelsProcessed; unsigned long m_NumberOfSpatialSamples; unsigned long m_NumberOfAvailableSamples; int m_RandomSeed; PerformanceProfileTable m_PerformanceProfile; PriorityImageType * m_PriorityImage; itkSetObjectMacro(PriorityImage,PriorityImageType); itkSetMacro(NumberOfChannels,unsigned int); itkSetMacro(TotalPixelsProcessed,unsigned long); itkGetMacro(TotalPixelsProcessed,unsigned long); itkGetMacro(NumberOfAvailableSamples,unsigned long); itkSetMacro(RandomSeed,int); itkSetMacro(NumberOfEvaluations,unsigned long); itkGetMacro(NumberOfEvaluations,unsigned long); itkSetMacro(NumberOfValueEvaluations,unsigned long); itkGetMacro(NumberOfValueEvaluations,unsigned long); itkSetMacro(NumberOfDerivativeEvaluations,unsigned long); itkGetMacro(NumberOfDerivativeEvaluations,unsigned long); itkSetMacro(NumberOfHessianEvaluations,unsigned long); itkGetMacro(NumberOfHessianEvaluations,unsigned long); itkGetConstReferenceMacro(NumberOfChannels,unsigned int); /** Number of spatial samples to used to compute metric */ itkSetClampMacro( NumberOfSpatialSamples, unsigned long, 1, NumericTraits::max() ); itkGetConstMacro( NumberOfSpatialSamples, unsigned long); enum GradientCalculationMethod{GradientImage,BSpline,Calculator,Inverse,ESM,InverseCompositional}; GradientCalculationMethod m_GradientMethod; itkGetMacro(GradientMethod,GradientCalculationMethod); itkSetMacro(GradientMethod,GradientCalculationMethod); enum HessianCalculationMethod{Fast,Diagonal,Full,True}; HessianCalculationMethod m_HessianMethod; itkGetMacro(HessianMethod,HessianCalculationMethod); itkSetMacro(HessianMethod,HessianCalculationMethod); typedef typename Superclass::FixedImageType FixedImageType; typedef typename Superclass::MovingImageType MovingImageType; typedef typename Superclass::TransformType TransformType; typedef typename Superclass::TransformPointer TransformPointer; typedef typename Superclass::InputPointType InputPointType; typedef typename Superclass::OutputPointType OutputPointType; typedef typename Superclass::TransformParametersType TransformParametersType; typedef typename Superclass::TransformJacobianType TransformJacobianType; typedef typename Superclass::InterpolatorType InterpolatorType; typedef typename Superclass::RealType RealType; typedef typename Superclass::GradientImageType GradientImageType; typedef typename Superclass::GradientPixelType GradientPixelType; typedef typename Superclass::GradientImagePointer GradientImagePointer; typedef typename Superclass::GradientImageFilterType GradientImageFilterType; typedef typename Superclass::GradientImageFilterPointer GradientImageFilterPointer; typedef typename Superclass::InterpolatorPointer InterpolatorPointer; typedef typename Superclass::FixedImageMaskType FixedImageMaskType; typedef typename Superclass::FixedImageMaskPointer FixedImageMaskPointer; typedef typename Superclass::MovingImageMaskType MovingImageMaskType; typedef typename Superclass::MovingImageMaskPointer MovingImageMaskPointer; typedef typename Superclass::MeasureType MeasureType; typedef typename Superclass::DerivativeType DerivativeType; typedef typename Superclass::ParametersType ParametersType; /** Index and Point typedef support. */ typedef typename FixedImageType::IndexType FixedImageIndexType; typedef typename FixedImageType::PixelType FixedImagePixelType; typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType; typedef typename MovingImageType::IndexType MovingImageIndexType; typedef typename TransformType::InputPointType FixedImagePointType; typedef typename TransformType::OutputPointType MovingImagePointType; //typedef typename Superclass::OutputPointType OutputPointType; typedef typename OutputPointType::CoordRepType CoordRepType; typedef ContinuousIndex MovingImageContinuousIndexType; typedef typename Superclass::RandomIterator RandomIterator; typedef typename Superclass::OutputPixelImageType OutputPixelImageType; typedef typename Superclass::PixelImageWriterType PixelImageWriterType; typedef Array2D HessianType; // for dpsidphi typedef typename MatrixOffsetTransformBase::MatrixType MatrixType; typedef typename MatrixOffsetTransformBase::OffsetType VectorType; typedef typename MatrixOffsetTransformBase::InputPointType PointType; public: /** Typedefs for using central difference calculator. */ typedef CentralDifferenceImageFunction DerivativeFunctionType; /** Pointer to central difference calculator. */ typename DerivativeFunctionType::Pointer m_DerivativeCalculator; virtual void GetValueAndDerivativeAndHessian_calc(bool calcvalue,bool calcderivative,bool calchessian, bool fixedmask,bool movingmask,const ParametersType & parameters, MeasureType & value, DerivativeType & derivative, HessianType & H, int channel) const; virtual MeasureType GetValue( const ParametersType& parameters, const int channel ) const ; /** Get the derivatives of the match measure. */ virtual void GetDerivative( const ParametersType & parameters, DerivativeType & derivative, int channel ) const ; /** Get value and derivatives for multiple valued optimizers. */ virtual void GetValueAndDerivative( const ParametersType & parameters, MeasureType& Value, DerivativeType& Derivative, int channel ) const ; virtual void GetHessian( const ParametersType& parameters, HessianType& Hessian, const int channel ) const ; virtual void GetValueAndHessian( const ParametersType& parameters, MeasureType& Value, HessianType& Hessian, const int channel ) const ; /** Get the derivatives of the match measure. */ virtual void GetDerivativeAndHessian( const ParametersType & parameters, DerivativeType & derivative, HessianType& Hessian, const int channel ) const; /** Get value and derivatives for multiple valued optimizers. */ virtual void GetValueAndDerivativeAndHessian( const ParametersType & parameters, MeasureType& Value, DerivativeType& Derivative, HessianType& Hessian, const int channel ) const; virtual void GetHessian( const ParametersType& parameters, HessianType& Hessian) const; virtual void GetValueAndHessian( const ParametersType& parameters, MeasureType& Value, HessianType& Hessian ) const ; /** Get the derivatives of the match measure. */ virtual void GetDerivativeAndHessian( const ParametersType & parameters, DerivativeType & derivative, HessianType& Hessian ) const ; /** Get value and derivatives for multiple valued optimizers. */ virtual void GetValueAndDerivativeAndHessian( const ParametersType & parameters, MeasureType& Value, DerivativeType& Derivative, HessianType& Hessian) const ; virtual void GetHessian( const ParametersType& parameters, HessianType& Hessian, const int channel, const double level ) const ; virtual void GetValueAndHessian( const ParametersType& parameters, MeasureType& Value,HessianType& Hessian, const int channel, const double level ) const ; /** Get the derivatives of the match measure. */ virtual void GetDerivativeAndHessian( const ParametersType & parameters, DerivativeType & derivative, HessianType& Hessian, int channel, double level ) const; /** Get value and derivatives for multiple valued optimizers. */ virtual void GetValueAndDerivativeAndHessian( const ParametersType & parameters, MeasureType& Value, DerivativeType& Derivative, HessianType& Hessian, int channel,double level ) const ; virtual void GetHessian( const ParametersType& parameters, HessianType& Hessian, const double level ) const ; virtual void GetValueAndHessian( const ParametersType& parameters, MeasureType& Value, HessianType& Hessian, const double level ) const ; /** Get the derivatives of the match measure. */ virtual void GetDerivativeAndHessian( const ParametersType & parameters, DerivativeType & derivative, HessianType& Hessian, double level ) const; /** Get value and derivatives for multiple valued optimizers. */ virtual void GetValueAndDerivativeAndHessian( const ParametersType & parameters, MeasureType& Value, DerivativeType& Derivative, HessianType& Hessian, double level ) const ; void Initialize(void) throw ( ExceptionObject ); /** 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; protected: AnytimeImageToImageMetric2(){ m_PriorityImage=NULL; m_RandomSeed=0; m_GradientMethod=GradientImage; m_HessianMethod=Fast; m_DerivativeCalculator = NULL; m_InterpolatorIsBSpline = false; m_FixedGradientImage= NULL; m_IdentityTransform=NULL; m_NumberOfEvaluations=0; m_NumberOfValueEvaluations=0; m_NumberOfDerivativeEvaluations=0; m_NumberOfHessianEvaluations=0; m_DPsiDPhiComputed=false; }; virtual void ComputeDpsidphi() const; typedef Array2D ConversionMatrixType; mutable ConversionMatrixType dpsidphi; mutable HessianType* m_HessianCache[NUMCHANNELS]; mutable HessianType* m_HessianSumCache[NUMCHANNELS]; mutable HessianType* m_HessianSumCacheESM[NUMCHANNELS]; mutable HessianType m_ICHessian; mutable bool m_ESMMode; mutable bool m_ICHessianComputed; mutable bool m_DPsiDPhiComputed; TransformPointer m_IdentityTransform; GradientImagePointer m_FixedGradientImage; virtual void ComputeFixedGradient() ; mutable unsigned long m_NumberOfEvaluations; mutable unsigned long m_NumberOfValueEvaluations; mutable unsigned long m_NumberOfDerivativeEvaluations; mutable unsigned long m_NumberOfHessianEvaluations; virtual ~AnytimeImageToImageMetric2() {}; /** * A fixed image spatial sample consists of the fixed domain point * and the fixed image value at that point. */ class FixedImageSpatialSample { public: FixedImageSpatialSample():FixedImageValue(NumericTraits::Zero) { FixedImagePointValue.Fill(0.0); // FixedImageJacobian.Fill(0.0); FixedImageIndex.Fill(0); //FixedImageParameterGradient.Fill(0.0); } ~FixedImageSpatialSample() {}; ParametersType FixedImageParameterGradient; //typename TransformType::JacobianType FixedImageJacobian; FixedImageIndexType FixedImageIndex; FixedImagePointType FixedImagePointValue; FixedImagePixelType FixedImageValue; }; virtual TransformBase::Pointer GetTransformCopy()const; /** FixedImageSpatialSample typedef support. */ typedef std::vector FixedImageSpatialSampleContainer; /** Container to store a set of points and fixed image values. */ FixedImageSpatialSampleContainer m_FixedImageSamples; /** Compute image derivatives at a point. */ // virtual void ComputeImageDerivatives( const FixedImageIndexType& fixedIndex, // const FixedImagePointType& fixedPoint, // const MovingImagePointType& mappedPoint, // ParametersType& parametergradient, // const MovingImageContinuousIndexType& tempindex) const; virtual void ComputeImageDerivatives( typename FixedImageSpatialSampleContainer::const_iterator fiter, const MovingImagePointType& mappedPoint, ParametersType& parametergradient, const MovingImageContinuousIndexType& tempindex) const; virtual void ConvertDerivative( GradientCalculationMethod gradmethod,ParametersType& derivative,ParametersType& derivativeESM ) const; virtual void ConvertHessian( GradientCalculationMethod gradmethod,HessianType& Hessian,HessianType& HessianESM ) const; virtual void ConvertDerivativeAndHessian( GradientCalculationMethod gradmethod,ParametersType& derivative,HessianType& Hessian,ParametersType& derivativeESM,HessianType& HessianESM ) const; virtual void SampleJacobians( FixedImageSpatialSampleContainer& samples ); /** Vars which control the computation **/ mutable typename FixedImageSpatialSampleContainer::const_iterator m_Index[NUMCHANNELS]; /** Cache the parameters **/ /** Uniformly select a sample set from the fixed image domain. */ virtual void SampleFixedImageDomain( FixedImageSpatialSampleContainer& samples); private: AnytimeImageToImageMetric2(const Self&) {}; //purposely not implemented void operator=(const Self&); //purposely not implemented }; class MetricEvaluationEvent : public AnyEvent { public: typedef MetricEvaluationEvent Self; typedef AnyEvent Superclass; typedef double MeasureType; // i'd rather define it from the metric class, but its typedef Array ParametersType; // templated MeasureType Value; ParametersType Parameters; bool ValueEvaluated; bool DerivativeEvaluated; bool HessianEvaluated; MetricEvaluationEvent(const MeasureType& ValueIn, const ParametersType& ParametersIn,const bool & veIn,const bool & deIn,const bool & heIn) { Value=ValueIn; Parameters=ParametersIn; ValueEvaluated=veIn; DerivativeEvaluated=deIn; HessianEvaluated=heIn; } MetricEvaluationEvent() { ValueEvaluated=false; DerivativeEvaluated=false; HessianEvaluated=false; } virtual ~MetricEvaluationEvent() {} virtual const char * GetEventName() const { return "MetricEvaluationEvent"; } virtual bool CheckEvent(const ::itk::EventObject* e) const { return dynamic_cast(e); } virtual ::itk::EventObject* MakeObject() const { return new Self; } MetricEvaluationEvent(const Self&s) :AnyEvent(s){}; private: void operator=(const Self&); }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkAnytimeImageToImageMetric2.txx" #endif #endif