Two-dimensional Fast Marching
The present post follows a series of posts: Escaping the maze: the Eikonal Equation unveiled, A numerical solution to the 1D eikonal equation and a Python implementation and The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method. Preliminarly reading the above posts is recommended.
Let us consider the classical eikonal equation
in which we assume
and that the computational domain is
and C is the unit circle of countour Γ.
It is easy to convince oneself that the solution to the eikonal equation in this case is
As shown in the figure above, the level curves, namely, the curves on which U is constant, are concentric circles while, for each point in the plane outside C, the paths allowing, for the point at hand, to reach Γ in the minimum time, namely, the characteristic curves, are radial lines. In this case, the solution to the eikonal equation is differentiable at each point. It should be noticed that the characteristic curves are always orthogonal to the level curves, as it can be evinced from the last part of the post Escaping the maze: the Eikonal Equation unveiled. Again in this case, the characteristic curves can be constructed as curves orthogonal to Γ which, by definition, is a level curve itself.
Let us now consider the case when the surface Γ is concave. In this case, we could think of building up the characteristic curves as curves orthogonal to Γ and, consequently, to construct the level curves starting from the characteristic ones.
As an alterative, starting from each point on Γ, we could trace circles with a certain radius a centered at that point, set up the level curves as the evenlope of the circles and finally trace the characteristic curves by exploiting the orthogonality to the level ones.
Nevertheless, we would notice how, doing so, the characteristic curves would intersect, giving rise to regions of non-differentiability in the level curves. The intersection points represent ambiguities in the determination of the minimum paths. Indeed, starting from the intersection points, the direction to take to reach the target set Γ would not be univocal. The non-differentiability of the level curves would be the counterpart of the non-differentiability of the eikonal function.
Just as in the one-dimensional case reported in A numerical solution to the 1D eikonal equation and a Python implementation, to restore the differentiability of the solution, the viscuous solution to the differential equation
is searched for with a vanishing ε. The fast marching method determines a solution to the viscuous equation for vanishing ε.
As seen in A numerical solution to the 1D eikonal equation and a Python implementation and in The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method, the solution is constructed starting from the boundary on which the initial conditions are assigned and stepping away from it. Moreover, let us recall the numerical discretization of the eikonal equation reported in The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method:
where
and
As can be seen, the new value of the eikonal should be larger than the neighboring values a₁ and a₂. Therefore, the information propagates starting from the least eikonal values towards points with larger eikonal. The Fast Marching Method exploits such a kind of causality to set up an efficient sweeping of the computational grid and to significantly reduce the computational complexity as compared to the and Tourin scheme.
Due to the need of setting up an iterative process, at the first step, only the points around Γ can be updated, when points around Γ means points that are not spaced to Γ more than one sample, either vertically or horizontally. All the other points of the computational domain will be downwind in the sense that they are located far away from Γ and not yet reached by the information carried on by the initial conditions.
Setting up a narrowband, namely, the set of all the points around Γ becomes natural. The idea behind the narrowband is to have a band of points separating the upwind region, which has been already traversed by the information necessary to the update of the solution, from the downwind region, not yet reached by the information coming from Γ. The equation numerically discretizing the eikonal equation from the only narroband point at each iterative step is the only region of the computational domain in which the solution evolves at each iteration.
At each iteration step, the minimum eikonals in the narrowband are labelled as set since the information flows from lower to higher eikonal points and so the minima cannot be influenced by the neighbors having a higher value of U. The set points are removed from the narrowband and dealt with as new initial conditions for the downwind region by adding its neighboors to the narrowband itself.
Let us now consider a two-dimensional problem in which the initial condition is assigned by the value of the eikonal at a single point. The algorithm initializes all the grid points to an “infinite” value, except for the point on which the initial condition is set as vanishing.
The point at which the initial condition is assigned is labelled as set since it has the least eikonal value and removed from the narrowband. The neighboring points and the point representing Γ, as mentioned, form the narrowband.
Among the four remaining points, the one with the least value of U is chosen. Obviously, any of the points in the narrowband can be chosen since we are at the first iteration step and the symmetry of the problem is spherical. Such a point is considered as set while the narrowband is advanced by including its neighboring points.
The numerical solution at the narrowband points is updated and the procedure is restarted: the narrowband node with minimum eikonal is determined, considered as set, the narrowband advanced etc.
To guarantee the efficiency of the approach, an algorithm to fastly determine the minimum of the narrowband is necessary. The list of coordinates of the narrowband points and of their respective eikonal values is then managed by a priority queue. It is based on a binary tree for which the eikonal value associated to each node is less then or equal to the eikonal values at the children nodes. In this way, the minimum of the narrowband eikonal is always associated to the root node and the cost for the search of the minimum is O(1). The sorted insertion of an element in the queue costs O(MlogM), where M is the dimension of the priority queue. Assuming M=N², the overall Fast Marching algorithm complexity, including the management of the priority queue, becomes O(N²logN) which is more convenient than the O(N³) of the Rouy and Tourin approach.
The script implementing the described approach is available at the link.
The first steps of the algorithms are the same as those described in the post: The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method.
The priority queue is managed by a proper library:
The initialization is similar to that of the Fast Sweeping Method. The only difference is the definition of an isSet matrix appointed to record if the eikonal at the corresponding node has been set or not. At the initialization, the points of Γ are appended to the narrowband (active list).
Following the initialization, the iterations are looped until the active list is not empty.
The code evolves along the following steps:
- the node on top of the priority queue, namely, that corresponding to the least narrowband eikonal value, is withdrawn;
- the withdrawn node is labelled as set;
- the eikonals at the neighbooring nodes to the withdrawn node are updated;
- each neighboor not already belonging to the narrowband is appended to the narrowband.
