Fast Marching Method
In the family of boundary value problem solving methods,
Fast Marching method is another numerical method. It solves the Eikonal
equation
Introduction
Fast marching method was introduced by Sethian. It was
originally designed for forward only propagation of the boundary with positive
values for speed, but has later been revised
for including both positive and negative values for speed and this is
now called Generalized Fast Marching
Method.
Motivation
This method allows us to convert the problem of continuously forward or backward moving front
into a stationary formulation where the arrival surface gives the front at zero
–level. Furthermore, using a grid with unit-distance and numerical techniques, this method is very
fast, efficient and highly stable.
Terminology
Arrival surface: The arrival surface results from the travel
time function that gives the arrival of the front at a given point
Narrow band: This is the set of active data points being
considered for the calculation of the next level of the front. The travel time
of each of the data points is computed based on the gradient between the
previous point and the current as well as the gradient between current and the
next point.
Initial Data
The speed at a given data point is known. Each data point
can have a different speed
Algorithm
The fast marching method solves boundary value problems that describe how a
closed curve evolves as a function of time t with a speed f(x) in the normal
direction at a point x on the curve. The speed is given but the time it takes
for the contour to pass the point x is calculated by the equation. Different
points may have different speeds. The region grows outward from the seed of the
area. The more the points considered, the better the contour. The fast
marching method makes use of a stationary surface that gives the arrival
information of the moving boundary. From the initial curve, build a little bit
of the surface upwards, with each iteration always starting from the initial
position. It's called fast because it uses a fast heap sort algorithm
that can tell the proper data point to update and hence to build the surface.
For each of the N points at any given height h, the next point is determined
and always in the order from the priority queue. The surface grows with height
h in a monotonically increasing order because previously evaluated grid points
are not revisited. This is because of non-negative speed and because the heap
guarantees that the grid point with the minimum distance is selected first out
of the data points in the sweep of N points at any level. Hence the entire
operation of determining the next contour takes NlogN time.
This is a numerical method because it tries to view the curve as a set of
discrete markers. The fast marching method has applications in image processing
such as for region growing.
In the fast march method, let us say there are N points on the zero-level of
the surface we track in a N*N*N unit distance grid. The curve is pushed
outwards in unit-distances. Then, each step involves the following:
1) The point with the smallest distance is removed from the priority queue
and its value is frozen so that we don't have to revisit it again and our
direction is strictly outward.
2) Points are added to the priority queue to maintain unit-thickness
3) The distance of the neighbors of the removed point are recomputed. The
finite difference is a multivariable quadratic equation.
we did not discuss the partial differential equation
(PDE) or how to solve it numerically. Specifically, we did not discuss the
distance function (time cost) based on which the next data point is selected.
Here we mention it with an example.
We start from a data point whose travel time is frozen i.e. we don't change
it again and proceed outwards to the neighboring points from it. We put these
neighboring points in a priority queue where the points are sorted based on the
minimum travel time. Then the top point in this priority queue is removed which
gives the minimum travel time. This point is then frozen. The neighboring
points are updated with the Euclidean norm of the gradients. The gradients are
along the x,y,z axes and are computed as derivatives (inverse of speed) defined
as the difference in the travel time between the new point and the old point on
the same axes over a unit distance. We take the maximum of this
gradient. The gradients are also taken between the next point and the new
point in the normal direction. We take the minimum of this gradient. The
Eucledian norm is nothing but a sum of squares just like we measure a
hypotenuse. So we take the squares of the maximum and minimum respectively for
each of the three axes. This gives us a quadratic in travel time t at the
new point and speed v that has a well-defined solution as one of either
the travel time at the previous point + (unit-distance / speed) or the next
point + (unit-distance/speed). In the fast marching method, we simplify this
computation to take forward only direction and at zero-level plane so only x
and y co-ordinates.
We repeat this for all other neighboring points around the one we removed
from the priority queue.
It might so happen that one or more of the neighboring points are in the
active set i.e. in the priority queue. If we take a new point for
which we want to compute the travel time as the center, then there are four
neighbours around it and one of them is the smallest of all the trial values.
Then there are four cases we need to consider. 1) where none of the
neighbors are in the active set, 2) one of the neighbors is in the active set
3) two of the neighbors are in the active set and 4) all three of
these neighbors are in the active set.
In the case where none of the neighbors are in the active set, the travel
time of the new point is the sum of that of A and the inverse speed function f
(f is based on the unit-distance).
In the case where one of the neighboring points b is already in the active
set, then we know that its travel time is more than the one we just removed and
also that the travel time of the new point will be lesser than the sum of the
travel time for b and the inverse speed function f.
If two of the neighbors are in the active set and one of them is frozen,
then this point is in the forward direction and will be the contributor to the
inverse speed function since the frozen value will not be changed and this
degenerates to the case 1)
If three of the neighbors are in the active set, the smallest values in each
co-ordinate direction is taken and this degenerates to the cases discussed
above.
Thus we obtain the travel times of the neighboring points as well and add
them to the priority queue
The neighboring points are typically referred to as the narrow band and we
locate the grid point with the minimum travel time.
Computational Complexity
The entire operation of determining the next contour takes NlogN time.
Error Analysis
There are two limitations of this approach.
This
method does not apply to continuous boundaries
There
is a potential large error in the first step of the numerical calculation.
Courtesy : Al Khalifa, S. Fomel : Fast Marching Implementation, Sethian : Fast Marching Method.