The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method
The present post follows A numerical solution to the 1D eikonal equation and a Python implementation, extends its contents and deals with the description of the numerical solution of the two-dimensional eikonal equation.
The two-dimensional eikonal equation is the following:
where U is the eikonal function, D is the domain in which we want to numerically solve the equation, Γ is the target set, f is the reciprocal of the speed and
represents the coordinates of a generic point in 2D. Of course, the eikonal equation takes the above form also for more than two dimensions, but we will stick here on the two-dimensional case only.
We assume the following two-dimensional computational domain D
and we introduce the following discretization
with
In other words, we use the same discretization step both along x and along y.
When we are far from the boundaries of the computational domain, the one-dimensional update scheme of A numerical solution to the 1D eikonal equation and a Python implementation is generalized into
where
and
At the left boundary, m=0 and we can use only forward differences to discretize the derivative, so that:
Similarly, at the right boundary, m=M and we can use only backward differences to discretize the derivative, so that:
The same applies to the top and bottom boundaries.
In the scheme above, finite differences in the sense of A numerical solution to the 1D eikonal equation and a Python implementation have been used for the evaluation of the derivatives at the points in the interior of D and the one sided differences at the boundary of the computational domain. The above discretization is an upwind discretization scheme in the sense of A numerical solution to the 1D eikonal equation and a Python implementation and is consistent, namely, progressive refinements of the discretizations will lead the numerical solution to converge to that of the continuous problem.
Let us suppose, without loss of generality, that
The case when
can be dealt with analogously.
Let us also set
A solution to the numerically discretized eikonal equation can be found as follows, distinguishing between two cases.
Case #1
Case #1 refers to the condition in which
Let us indeed consider the tentative solution
Note that, since
then
Case #1 assumption implies that
which, in turn, entails that
Accordingly and according also to the ReLU operator, the numerical discretization of the eikonal equation becomes
and the tentative solution is also the solution, namely
Case #2
Case #2 refers to the condition in which
Under such an assumption,
then, the solution should be searched for as the root of the following quadratic equation
The two possible solutions of such a quadratic equation are
The two solutions can be rewritten as
To be consistent with the solution of the above quadratic equation, then we must have
Accordingly, since
then we must choose the solution to the quadratic equation with positive sign. Again to guarantee
we must enforce the following condition to the argument of the square root
which means
consistently with Case #2 assumption.
As already seen for the one-dimensional case, the unknown should be initialized at a value larger than the largest value of f multiplied by the diameter Dmax of the computational domain. In other words, the unknown should be initialized at a value larger than
Above, we have discussed about the numerical update of the solution of the eikonal equation at a certain grid point mn. To fully solve the equation at all points of the computational domain, a full scanning would be needed to iteratively update the solution until convergence. A such approach was proposed by Rouy and Tourin in the paper entitled A viscosity solutions approach to shape-from-shading, published in 1992 on the SIAM Journal on Numerical Analysis. However, for a two-dimensional domain of NxN points, the computational complexity would be O(N²). Here, we discuss an alternative approach known as fast sweeping approach.
In A numerical solution to the 1D eikonal equation and a Python implementation, we have seen that two sweeps were enough to reach convergence in the one-dimensional case since the minimum-time curves could develop only from left to right and from right to left. In the two-dimensional case, the minimum-time curves may develop along arbitrary directions, so that the question arises about which sweeps are now necessary.
The fast sweeping algorithm performs four sweeps:
The above discussed numerical scheme is increasing with the iteration number and converges monotonically to the solution of the numerical discretization of the eikonal equation.
Moreover,
Starting from the target set Γ, the information reaches all the points of D when
where the underlining denotes the vector of all the samples of U in the computational domain. When this occurs, the iterations can be terminated. Even if the solution has not reached true convergence, it is already as accurate as that corresponding to convergence.
In the case when the target set Γ is a solitary point, namely,
then, after 4 sweeps,
where d is the distance between two points and
is the solution to numerical discretization of the eikonal equation with a discretization step equal to h.
In the case when the target set Γ is a collection of points, namely,
then, on letting
after 4 sweeps, we have
Finally, in the case when Γ is a curve or a surface, on denoting by
the exact solution to the discretized problem, then
A Python implementation for the fast sweeping approach in two dimensions is available at the GitHub page.
The code uses the eikonal_routines Python library which includes the update_eikonal_2D function which performs the update of the eikonal according to the above indicated rules.
