#ifndef _itkESMDemonsRegistrationFunction_h_ #define _itkESMDemonsRegistrationFunction_h_ #include "itkPDEDeformableRegistrationFunction.h" #include "itkPoint.h" #include "itkCovariantVector.h" #include "itkInterpolateImageFunction.h" #include "itkLinearInterpolateImageFunction.h" #include "itkCentralDifferenceImageFunction.h" #include "itkWarpImageFilter.h" namespace itk { /** * \class ESMDemonsRegistrationFunction * * \brief Fast implementation of the symmetric demons registration force * * This class provides a substantially faster implementation of the * symmetric demons registration force. Speed is improved by keeping * a deformed copy of the moving image for gradient evaluation. * * The moving image should not be saturated. We indeed use * NumericTraits::Max() as a special value. * * \author Tom Vercauteren, INRIA & Mauna Kea Technologies * * \sa SymmetricForcesDemonsRegistrationFunction * \sa SymmetricForcesDemonsRegistrationFilter * \sa DemonsRegistrationFilter * \sa DemonsRegistrationFunction * \ingroup FiniteDifferenceFunctions * */ template class ITK_EXPORT ESMDemonsRegistrationFunction : public PDEDeformableRegistrationFunction< TFixedImage, TMovingImage, TDeformationField> { public: /** Standard class typedefs. */ typedef ESMDemonsRegistrationFunction Self; typedef PDEDeformableRegistrationFunction< TFixedImage, TMovingImage, TDeformationField > Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Method for creation through the object factory. */ itkNewMacro(Self); /** Run-time type information (and related methods). */ itkTypeMacro( ESMDemonsRegistrationFunction, PDEDeformableRegistrationFunction ); /** MovingImage image type. */ typedef typename Superclass::MovingImageType MovingImageType; typedef typename Superclass::MovingImagePointer MovingImagePointer; /** FixedImage image type. */ typedef typename Superclass::FixedImageType FixedImageType; typedef typename Superclass::FixedImagePointer FixedImagePointer; typedef typename FixedImageType::IndexType IndexType; typedef typename FixedImageType::SizeType SizeType; typedef typename FixedImageType::SpacingType SpacingType; /** Deformation field type. */ typedef typename Superclass::DeformationFieldType DeformationFieldType; typedef typename Superclass::DeformationFieldTypePointer DeformationFieldTypePointer; /** Inherit some enums from the superclass. */ itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension); /** Inherit some enums from the superclass. */ typedef typename Superclass::PixelType PixelType; typedef typename Superclass::RadiusType RadiusType; typedef typename Superclass::NeighborhoodType NeighborhoodType; typedef typename Superclass::FloatOffsetType FloatOffsetType; typedef typename Superclass::TimeStepType TimeStepType; /** Interpolator type. */ typedef double CoordRepType; typedef InterpolateImageFunction InterpolatorType; typedef typename InterpolatorType::Pointer InterpolatorPointer; typedef typename InterpolatorType::PointType PointType; typedef LinearInterpolateImageFunction DefaultInterpolatorType; /** Warper type */ typedef WarpImageFilter WarperType; typedef typename WarperType::Pointer WarperPointer; /** Covariant vector type. */ typedef CovariantVector CovariantVectorType; /** Gradient calculator types. */ typedef CentralDifferenceImageFunction GradientCalculatorType; typedef typename GradientCalculatorType::Pointer GradientCalculatorPointer; /** Set the moving image interpolator. */ void SetMovingImageInterpolator( InterpolatorType * ptr ) { m_MovingImageInterpolator = ptr; m_MovingImageWarper->SetInterpolator( ptr ); } /** Get the moving image interpolator. */ InterpolatorType * GetMovingImageInterpolator(void) { return m_MovingImageInterpolator; } /** This class uses a constant timestep of 1. */ virtual TimeStepType ComputeGlobalTimeStep(void * itkNotUsed(GlobalData)) const { return m_TimeStep; } /** Return a pointer to a global data structure that is passed to * this object from the solver at each calculation. */ virtual void *GetGlobalDataPointer() const { GlobalDataStruct *global = new GlobalDataStruct(); global->m_SumOfSquaredDifference = 0.0; global->m_NumberOfPixelsProcessed = 0L; global->m_SumOfSquaredChange = 0; return global; } /** Release memory for global data structure. */ virtual void ReleaseGlobalDataPointer( void *GlobalData ) const; /** Set the object's state before each iteration. */ virtual void InitializeIteration(); /** This method is called by a finite difference solver image filter at * each pixel that does not lie on a data set boundary */ virtual PixelType ComputeUpdate(const NeighborhoodType &neighborhood, void *globalData, const FloatOffsetType &offset = FloatOffsetType(0.0)); /** Get the metric value. The metric value is the mean square difference * in intensity between the fixed image and transforming moving image * computed over the the overlapping region between the two images. */ virtual double GetMetric() const { return m_Metric; } /** Get the rms change in deformation field. */ virtual const double &GetRMSChange() const { return m_RMSChange; } /** Set/Get the threshold below which the absolute difference of * intensity yields a match. When the intensities match between a * moving and fixed image pixel, the update vector (for that * iteration) will be the zero vector. Default is 0.001. */ virtual void SetIntensityDifferenceThreshold(double); virtual double GetIntensityDifferenceThreshold() const; /** Set/Get the maximum update step length. In Thirion this is 0.5. * Setting it to 0 implies no restriction (beware of numerical * instability in this case. */ virtual void SetMaximumUpdateStepLength(double sm){ this->m_MaximumUpdateStepLength =sm;}; virtual double GetMaximumUpdateStepLength() const { return this->m_MaximumUpdateStepLength;}; enum GradientType { Symmetric=0, Fixed, Moving }; virtual void SetUseGradientType( GradientType gtype ) { m_UseGradientType = gtype; } virtual GradientType GetUseGradientType() const { return m_UseGradientType; } protected: ESMDemonsRegistrationFunction(); ~ESMDemonsRegistrationFunction() {} void PrintSelf(std::ostream& os, Indent indent) const; /** FixedImage image neighborhood iterator type. */ typedef ConstNeighborhoodIterator FixedImageNeighborhoodIteratorType; /** A global data type for this class of equation. Used to store * iterators for the fixed image. */ struct GlobalDataStruct { double m_SumOfSquaredDifference; unsigned long m_NumberOfPixelsProcessed; double m_SumOfSquaredChange; }; private: ESMDemonsRegistrationFunction(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented /** Cache fixed image information. */ SpacingType m_FixedImageSpacing; PointType m_FixedImageOrigin; double m_Normalizer; /** Function to compute derivatives of the fixed image. */ GradientCalculatorPointer m_FixedImageGradientCalculator; GradientType m_UseGradientType; /** Function to interpolate the moving image. */ InterpolatorPointer m_MovingImageInterpolator; /** Filter to warp moving image for fast gradient computation. */ WarperPointer m_MovingImageWarper; /** The global timestep. */ TimeStepType m_TimeStep; /** Threshold below which the denominator term is considered zero. */ double m_DenominatorThreshold; /** Threshold below which two intensity value are assumed to match. */ double m_IntensityDifferenceThreshold; /** Maximum update step length in pixels (default is 0.5 as in Thirion). */ double m_MaximumUpdateStepLength; /** The metric value is the mean square difference in intensity between * the fixed image and transforming moving image computed over the * the overlapping region between the two images. */ mutable double m_Metric; mutable double m_SumOfSquaredDifference; mutable unsigned long m_NumberOfPixelsProcessed; mutable double m_RMSChange; mutable double m_SumOfSquaredChange; /** Mutex lock to protect modification to metric. */ mutable SimpleFastMutexLock m_MetricCalculationLock; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkESMDemonsRegistrationFunction.txx" #endif #endif