1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
|
/*=========================================================================
*
* 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 <typename TImageType>
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<Self>;
using ConstPointer = SmartPointer<const Self>;
/** \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<ImageType>;
/** Neighborhood radius type */
using RadiusType = typename ConstNeighborhoodIterator<TImageType>::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<TImageType, DefaultBoundaryConditionType>;
/** 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<PixelRealType, Self::ImageDimension>;
/** A floating point offset from an image grid location. Used for
* interpolation among grid values in a neighborhood. */
using FloatOffsetType = Vector<float, Self::ImageDimension>;
/** 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
|
