DGtal
1.2.0

Part of the Geometry package.
The Fast Marching Method (FMM for short) is a numerical technique for tracking the evolution of an expanding front within the levelset framework (see [95] for a review).
If the front evolves in the normal direction at unit speed, the arrival time at a point \( x \) corresponds to the distance of \( x \) to the initial front.
That is why FMM can be used to compute a (signed) distance field from an initial point set (which may be located around an interface or not).
However, FMM can be used to other tasks: e.g. the extrapolation of a scalar field out of an interface, such that it is constant on rays normal to the interface [2].
Related examples are exampleFMM2D.cpp and exampleFMM3D.cpp.
Let us assume that we compute in the fundamental lattice \( \mathbb{Z}^n \) the distance field from a given lattice point set. The distance value of each lattice point is computed by marching out from the initial set of lattice points for which the distance values are known. This set is an initialization of the socalled accepted point set, which incrementally grows until all lattice points of the computation domain have been accepted or until a stopping criterion has been met.
Each lattice point adjacent to one of the accepted points is put into the socalled candidate point set.
The accepted point set and the candidate point set are depicted below. See Data structures for details about how these sets are stored.
For each candidate point, a tentative value is computed for its distance, using only the values of the accepted points lying in its neighborhood. This task is delegated to an instance of a model of CPointFunctor, which is defined as L2FirstOrderLocalDistance by default. However, you are free to use L2SecondOrderLocalDistance, for second order convergence, L1LocalDistance and LInfLocalDistance for other norms. See Computing distances for further details.
Then, the point of smallest tentative value is iteratively added to the set of accepted points. The tentative values of the candidates adjacent to the newly added point are updated using the distance value of the newly added point.
The class FMM, which implements the fast marching method, is based on two external data structures and one internal data structure:
In practice, you can choose either ImageContainerBySTLVector and DigitalSetBySTLSet if you work on the whole domain:
or ImageContainerBySTLMap and DigitalSetFromMap if you work on a small fraction of the domain:
a data structure for the candidate point set. Three operations are required:
This data structure is usually a minheap data structure with back pointers from the lattice to the minheap (see [95]). The memory cost of such solution is however high. That is why, we implemented the candidate point set as a STL set of pairs <point, tentative value>. Instead of updating the tentative values, we insert a new pair <point, tentative value>. This solution is less memory consumming and experimentally (nearly) as efficient as the former one.
In a numerical point of view, the key point of the algorithm is to accurately approximate the distance value of a given point from the distance value of its neighbors.
Four models of CPointFunctor can be used for this task:
We detail below how L2FirstOrderLocalDistance, which is used by default, approximate the Euclidean distance at some point \( x \), from the available distance values of the lattice points lying in the 1neighborhood of \( x \).
The distance value is computed such that the upwind gradient of the distance image is one. Let us denote by \( d \) the distance value at \( x \). Moreover, for a given quadrant, let us denote by \( d_i \) the distance value of the unique 1neighbor of \( x \) along the \( i \) axis. The distance value \( d \) is taken as the minimum solution of the following quadratic equation
over all the quadrants:
\( \sum_{i = 1 \ldots n } ( d  d_i )^2 = 1 \)
Obviously, if the distance value \( d' \) of only one 1neighbor is available, then \( d = d' \pm 1 \) ( \( + \) for positive distance values, \(  \) for negative ones).
If the distance values of two (or more than two in nD) 1neighbors are available, the previous quadratic equation should have two solutions and we take the larger one for \( d \).
The numerical approach implemented in L2FirstOrderLocalDistance::operator() to compute \( d \), even with poor initial data or with numerical errors, is detailed in [79] [p. 7273]. Note that if \( h \) is the lattice spacing, the approximation error tends to zero in \( O(h) \) (firstorder convergence).
In DGtal, the class FMM is parametrized by the following types:
After the inclusion of the following file
you may define the FMM type as follows:
Since an instance of FMM keeps a reference to the accepted point set, a reference to the distance image, and a constant reference to the point predicate, you must be sure that these three objects have been constructed and exist before constructing a FMM instance.
In this example, we start from one point whose distance value is zero:
We now construct a FMM instance as follows:
After the call of the method FMM::compute(), the distance value of all the accepted points is available in the distance image.
Here is the computed distance field:
The whole example may be found in exampleFMM2D.cpp.
The computation ends either if all the points represented by the point predicate have been accepted, or if one of the following stopping criteria have been met:
In order to use these stopping criteria, you must provide the corresponding thresholds at construction as follows:
With the above code, we stop the computation as soon as the number of accepted points equals to maximalNumberOfPoints or the distance value of the next candidate point exceed maximalDistance.
Note that you can also perform the computation step by step. This approach may be useful to better control the computation and its end. You must use FMM::computeOneStep(), which takes as input and returns the last accepted point and its distance value, as follows:
As reported in [95], the Fast Marching Method has numerous applications.
A signed distance field can be used to
In exampleFMM3D.cpp, we show how to compute a signed distance field to a digital surface.
First, we define the FMM type as follows:
Then, we construct and initialize the external data structures as follows:
In the above snippet, we insert all the 3D lattice points incident to the digital surface in the initial point set and we assign to them an initial distance value (equal to \( 0.5 \) for outer points and \( 0.5 \) for inner points) by merely calling the static method FMM::initFromBelsRange().
Finally, we instanciate the FMM class and we perform the whole computation:
The computation is limited to a narrow band around the digital surface: each point of the band are located at a distance less than maximalDistance to the digital surface.
The result is displayed below.