/*========================================================================= * * Copyright NumFOCUS * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * https://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef itkFiniteDifferenceFunction_h #define itkFiniteDifferenceFunction_h #include "itkLightObject.h" #include "itkConstNeighborhoodIterator.h" #include "itkVector.h" namespace itk { /** * \class FiniteDifferenceFunction * * This class is a component object of the finite difference solver hierarchy * (see FiniteDifferenceImageFilter). It defines a generic interface for a * function object that computes a single scalar value from a neighborhood of * values. Examples of the application of this class are the various flavors * of AnisotropicDiffusionFunction and LevelSetFunction objects. * * These functions calculate the incremental change at a pixel in the solution * image from one iteration of the Partial Differential Equation (PDE) solver to * the next. * * \par * Subclasses of FiniteDifferenceImageFilter (solvers) call the * ComputeUpdate() method of this class to compute \f$ \Delta u^n_{\mathbf{i}} * \f$ at each \f$ \mathbf{i} \f$ in \f$ u \f$. Because the size of the time * step for each iteration of the p.d.e. solution depends on the particular * calculations done, this function object is also responsible for computing * that time step (see ComputeGlobalTimeStep() ). * * \par How to use this class * FiniteDifferenceFunction must be subclassed to add functionality for * ComputeUpdate, ComputeGlobalTimeStep, and Get/ReleaseGlobalDataPointer. * * \par A note on thread safety. * The ComputeUpdate() methods of this filter are declared as const to enforce * thread-safety during execution of FiniteDifferenceImageFilter solver * algorithms. The InitializeIteration() method is intended to provide a safe * way to modify the state of the object between threaded calculations of * solvers. * * \todo Possibly subclass this object from Function. Stumbling blocks here * are the specialized api of FiniteDifferenceFunction. * * \ingroup Functions * \ingroup ITKFiniteDifference */ template class ITK_TEMPLATE_EXPORT FiniteDifferenceFunction : public LightObject { public: ITK_DISALLOW_COPY_AND_MOVE(FiniteDifferenceFunction); /** Standard class type aliases. */ using Self = FiniteDifferenceFunction; using Superclass = LightObject; using Pointer = SmartPointer; using ConstPointer = SmartPointer; /** \see LightObject::GetNameOfClass() */ itkOverrideGetNameOfClassMacro(FiniteDifferenceFunction); /** Extract some parameters from the image type */ using ImageType = TImageType; using PixelType = typename ImageType::PixelType; using PixelRealType = double; /** Save image dimension. */ static constexpr unsigned int ImageDimension = ImageType::ImageDimension; /** Define the TimeStepType to always be double. */ using TimeStepType = double; /** The default boundary condition for finite difference * functions that is used unless overridden in the Evaluate() method. */ using DefaultBoundaryConditionType = ZeroFluxNeumannBoundaryCondition; /** Neighborhood radius type */ using RadiusType = typename ConstNeighborhoodIterator::RadiusType; /** The type of data structure that is passed to this function object * to evaluate at a pixel that does not lie on a data set boundary. */ using NeighborhoodType = ConstNeighborhoodIterator; /** The type of data structure that holds the scales with which the * neighborhood is weighted to properly account for spacing and neighborhood radius. */ using NeighborhoodScalesType = Vector; /** A floating point offset from an image grid location. Used for * interpolation among grid values in a neighborhood. */ using FloatOffsetType = Vector; /** This method allows the function to set its state before each iteration * of the finite difference solver (image filter) that uses it. This is * a thread-safe time to manipulate the object's state. * * An example of how this can be used: the AnisotropicDiffusionFunction class * of FiniteDifferenceFunction use this method to pre-calculate an average * gradient magnitude across the entire image region. This value is set in * the function object and used by the ComputeUpdate() methods that are called * at each pixel as a constant. */ 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. The width of the * data set boundary is defined by the width of the neighborhood being * evaluated. * * The neighborhood argument is the neighborhood data. * The globalData argument is a pointer to a data structure that holds * values that need to be persistent across calls to this function, but * cannot be stored in the function object itself for thread-safety reasons. * Examples are values needed to compute the time-step for an iteration * of the solver. * \sa InitializeIteration * \sa ComputeGlobalTimeStep */ virtual PixelType ComputeUpdate(const NeighborhoodType & neighborhood, void * globalData, const FloatOffsetType & offset = FloatOffsetType(0.0)) = 0; /** Sets the radius of the neighborhood this FiniteDifferenceFunction * needs to perform its calculations. */ void SetRadius(const RadiusType & r); /** Returns the radius of the neighborhood this FiniteDifferenceFunction * needs to perform its calculations. */ const RadiusType & GetRadius() const; /** Set the ScaleCoefficients for the difference * operators. The defaults a 1.0. These can be set to take the image * spacing into account. */ void SetScaleCoefficients(const PixelRealType vals[ImageDimension]); /** Returns the current scale coefficients. */ void GetScaleCoefficients(PixelRealType vals[ImageDimension]) const; /** Compute the scales that weight the neighborhood during difference operations * to properly account for spacing and neighborhood radius */ const NeighborhoodScalesType ComputeNeighborhoodScales() const; /** Computes the time step for an update given a global data structure. * The data used in the computation may take different forms depending on * the nature of the equations. This global data cannot be kept in the * instance of the equation object itself since the equation object must * remain stateless for thread safety. The global data is therefore managed * for each thread by the finite difference solver filters. */ virtual TimeStepType ComputeGlobalTimeStep(void * GlobalData) const = 0; /** Returns a pointer to a global data structure that is passed to this * object from the solver at each calculation. The idea is that the * solver holds the state of any global values needed to calculate the * time step, while the equation object performs the actual calculations. * * The global data should also be initialized in this method. * */ virtual void * GetGlobalDataPointer() const = 0; /** When the finite difference solver filter has finished using a global * data pointer, it passes it to this method, which frees the memory. * * The solver cannot free the memory because it does not know the type * to which the pointer points. */ virtual void ReleaseGlobalDataPointer(void * GlobalData) const = 0; protected: FiniteDifferenceFunction(); ~FiniteDifferenceFunction() override = default; void PrintSelf(std::ostream & os, Indent indent) const override; RadiusType m_Radius{}; PixelRealType m_ScaleCoefficients[ImageDimension]{}; }; } // end namespace itk #ifndef ITK_MANUAL_INSTANTIATION # include "itkFiniteDifferenceFunction.hxx" #endif #endif