Corner cases, singularities, and dynamic factoring
JJournal of Scientific Computing manuscript No. (will be inserted by the editor)
Corner cases, singularities, and dynamic factoring
Dongping Qi · Alexander Vladimirsky
Received: date / Accepted: date
Abstract
In Eikonal equations, rarefaction is a common phenomenon known to de-grade the rate of convergence of numerical methods. The “factoring” approach al-leviates this difficulty by deriving a PDE for a new (locally smooth) variable whilecapturing the rarefaction-related singularity in a known (non-smooth) “factor”. Pre-viously this technique was successfully used to address rarefaction fans arising atpoint sources. In this paper we show how similar ideas can be used to factor the 2Drarefactions arising due to nonsmoothness of domain boundaries or discontinuities inPDE coefficients. Locations and orientations of such rarefaction fans are not knownin advance and we construct a “just-in-time factoring” method that identifies them dy-namically. The resulting algorithm is a generalization of the Fast Marching Methodoriginally introduced for the regular (unfactored) Eikonal equations. We show thatour approach restores the first-order convergence and illustrate it using a range ofmaze navigation examples with non-permeable and “slowly permeable” obstacles.
Keywords
Eikonal · rarefaction fans · factoring · Fast Marching · robotic navigation Mathematics Subject Classification (2000) · · · · The Eikonal equation arises in many application domains including the geometricoptics, optimal control theory, and robotic navigation. It is a first-order non-linear
The second author’s work is supported in part by the National Science Foundation grant DMS-1738010and the Simons Foundation Fellowship.D. QiCenter for Applied Mathematics, Cornell University, Ithaca, NY 14853E-mail: [email protected]. VladimirskyDept. of Mathematics and Center for Applied Mathematics, Cornell University, Ithaca, NY 14853E-mail: [email protected] a r X i v : . [ m a t h . NA ] D ec Dongping Qi, Alexander Vladimirsky partial differential equation of the form (cid:40) | ∇ u ( x ) | F ( x ) = , u ( x ) = , x ∈ Q . (1)In control-theoretic context, the function u ( x ) can be interpreted as the minimumtime needed to reach the exit set Q starting from x , with F specifying the speed ofmotion. The characteristics of this PDE coincide with the gradient lines of u and yieldmin-time-optimal trajectories to Q . This equation typically does not have a classicalsolution, making it necessary to select a weak Lipschitz-continuous solution that isphysically meaningful. The theory of viscosity solutions [7] accomplishes this andguarantees the existence and uniqueness for a very broad set of problems. Viscositysolutions exhibit both shocklines (where characteristics run into each other) and rar-efaction fans (where many characteristics emanate from the same point). Typically,the former receive most of the attention in numerical treatment since they lead to adiscontinuity in ∇ u . But in this paper we focus on the latter, which result in a blowup in second derivatives of u on parts of the domain where ∇ u remains continuousand bounded. (See Figure 1.) (A) (B)Fig. 1: A simple example with F = Q consisting of a single “point source” x = ( . , ) . (Left) Level sets of u ( x ) , the distance to x on the domain with an obstacle. Dashed lines show the opti-mal/shortest paths, some of which follow a part of the obstacle boundary. The “shockline” (where ∇ u isundefined) consists of all points starting from which the shortest path is not unique. The point source and3 out of 4 obstacle corners generate “rarefaction fans” of characteristics. (Right) The level curves of thelog-Frobenius-norm of Hessian, illustrating the blow up of the second derivatives of u due to rarefactionfans. Since the local truncation error of numerical discretizations is proportional tohigher derivatives, it is not surprising that rarefaction fans degrade the rate of con-vergence of standard numerical methods. Special “factoring techniques” have beenintroduced to alleviate this problem for rarefaction fans arising from point sourceboundary conditions [10, 15, 16, 17, 20, 27]. The idea is to represent u as a product ynamic Factoring 3 [10] or a sum [16] of two functions: one capturing the right type of singularity at thepoint source and another with bounded derivatives except at shocklines. The formeris known based on F restricted to Q ; the latter is a priori unknown and recovered bysolving a modified PDE by a version of one of the efficient methods (e.g., [21, 29] or[32]) originally developed for an “un-factored” Eikonal equation (1). We review thisprior work in section 2 and show that, despite its better rate of convergence, factoringcan also have detrimental effects on parts of the domain far from the point source.We also consider a localized version of factoring, which often improves the accuracyand is more suitable for problems with multiple point sources.However, point sources are not the only origin of rarefaction fans. They can alsoarise due to non-smoothness of the boundary; e.g., at some of the “obstacle corners”in 2D maze navigation problems (see Figure 1). The degradation of accuracy leadsto numerical artifacts in computed trajectories passing near such corners. Unlike inthe point source case, these rarefaction fans are not radially symmetric; moreover,their locations and geometry have to be determined dynamically. We handle this bydeveloping a “just-in-time localized factoring” method and verifying its rate of con-vergence numerically (section 3). Even more complicated fans can arise at corners of“slowly permeable obstacles” (i.e., problems with piecewise-continuous speed func-tion F ). The extension of our dynamic factoring to this case is covered in section4. Throughout the paper, we present our approach using a Cartesian grid discretiza-tion of the Eikonal equation with grid aligned obstacles (and discontinuities of F ).However, the ideas presented here have a broader applicability. Arbitrary polygonalobstacles can be treated in exactly the same way if the discretization is posed onan obstacle-fitted triangulated mesh. Dynamic factoring for more general Hamilton-Jacobi-Bellman PDEs would work very similarly, though the factoring will need toaccount for the anisotropy in rarefaction fans. We conclude by discussing these andother future extensions as well as the limitations of our approach in section 5. We begin by examining the classical Eikonal equation (1) in 2D with a “single point-source” boundary condition: Q = { x } and u ( x ) =
0. Throughout this paper, wewill assume that the controlled process is restricted to a closed set Ω ⊂ R . I.e., u ( x ) is the minimal-time from x ∈ Ω to Q ⊂ Ω without leaving Ω though possiblytraveling along parts of ∂ Ω ; see the trajectories traveling along the obstacle boundaryin Figure 1(A). This makes u an Ω -constrained viscosity solution of (1) (see Chapter4.5 in [2]), but for the purposes of numerical implementation we can simply imposethe boundary condition u = + ∞ on R \ Ω .A common approach for discretizing (1) on a uniform Cartesian grid is to useupwind finite differences:max { D − xi , j u , − D + xi , j u , } + max { D − yi , j u , − D + yi , j u , } = F i , j , (2) Dongping Qi, Alexander Vladimirsky using the standard finite difference notation D − xi , j u = u i , j − u i − , j h , D + xi , j u = u i + , j − u i , j h , D − yi , j u = u i , j − u i , j − h , D + yi , j u = u i , j + − u i , j h . (3)The discretized system of equations (2) is coupled and non-linear. We postpone thediscussion of fast algorithms used to solve it until section 2.1 and for now focus onthe rate of convergence of the numerical solution to the viscosity solution of PDE(1). Since this discretization is globally first-order accurate, the local truncation erroris proportional to the second derivatives of u , which blow up as we approach x .Because of this blow up, even if we assume that local truncation errors accumulatelinearly, the global error would decreases as O ( h log h ) instead of the expected O ( h ) . To illustrate the convergence rate, we will consider a simple example with aknown analytic solution [16] on a square domain : F ( x ) = s + v · ( x − x ) = ⇒ u ( x ) = | v | acosh (cid:32) + s | v | | x − x | F ( x ) (cid:33) . (4)Figure 2 shows the solution and the convergence plots for the parameter values x =( , ) , s = , v = ( . , ) . The thick black line in Figure 2(B) corresponds to solving(2) on Ω and clearly shows the sublinear convergence. (Throughout this paper, alllogarithms of errors and grid sizes are reported and plotted in base e . )To alleviate this issue, one could simply “enlarge the exit set” by choosing some( h -independent) constant radius r > , initializing u = B = B r ( x ) = { x ∈ Ω | | x − x | ≤ r } , and solving (2) on the rest of the grid. This avoids the rarefac-tion fan (since only one characteristic stems from each point on ∂ B ), but introducesa O ( r ) difference compared to the solution of the original point-source-based prob-lem. One can also use a better approximation of u on B ; e.g., using T ( x ) = | x − x | F ( x ) ,we already reduce this additional error to O ( r ) . The latter approach is based on as-suming that F (rather than u ) is constant on B , in which case the characteristics arestraight lines. Of course, one could also take a truly Lagrangian approach, employingray-tracing to compute a highly accurate value of u at all gridpoints on B , but thisbecomes increasingly expensive as h →
0. Instead, we evaluate the feasibility of an“approximate Lagrangian” initialization technique, where the characteristics are stillassumed to be straight on B , but u ( x ) is approximated more accurately by integrating1 / F on the line segment from x to x . In all the figures of this section, we slightlyabuse the notation and call this approach Lagrangian, reporting the results for r = . F , the errordue to this “characteristics are straight” assumption prevents the overall first orderconvergence on the entire grid. This formula for u ( x ) is derived on the unbounded domain Ω ∞ = { x ∈ R | F ( x ) > } but remainsvalid on Ω = [ , ] × [ , ] as long as Ω ∞ -optimal trajectories from every x ∈ Ω to x stay entirely inside Ω , which is the case for all examples considered in this section. The linearity of F ( x ) can be used to showthat all optimal paths are circular arcs, whose radii are monotone decreasing in | v | . (When v = , theseradii are infinite; i.e., all optimal paths are straight lines and u ( x ) = s | x − x | . )ynamic Factoring 5(A) level sets of u ( x ) Log(1/h) -7.5-7-6.5-6-5.5-5-4.5-4-3.5
Log ( e rr o r) L-inf error unit slopeoriginalglobal factoringlocal factoringLagrangian (B) L ∞ errorFig. 2: (Left) Level curves of the value function computed from the formula (4) with x = ( , ) , s = , v =( . , ) . (Right) The corresponding convergence plot (based on the L ∞ error) for several discretizationapproaches, with the thin black line of slope ( − ) included to aid the visual comparison. For this example,the convergence plot based on the L error is very similar and thus omitted. A factored Eikonal equation was proposed in [10] as a method for dealing withpoint-source rarefaction fans without introducing any special approximations on B and recovering the first order of accuracy on the entire domain. The main idea isto split the original value function u ( x ) into two functions: one of them ( T , definedabove) encodes the right type of singularity at the point source, while the other ( τ ( x ) ,our new unknown) will be smooth near x . In addition to point-sources, [10] also usedfactoring to treat “plane-wave sources” (i.e., constant Dirichlet boundary conditionsspecified on a straight line in 2D). In the latter case, there is no singularities in thesolutions at the boundary (so, the numerics for the original/unfactored version is stillfirst-order accurate), but a factored version still yields lower errors in many examples.The original factoring in [10] used an ansatz u ( x ) = T ( x ) τ ( x ) to derive a newfactored PDE for τ , which was then solved with the boundary condition τ ( x ) = . In this paper we rely on a slightly simpler additive splitting introduced in [16] :assuming that u ( x ) = T ( x ) + τ ( x ) , we find τ by solving | ∇ T ( x ) + ∇ τ ( x ) | F ( x ) = , (5)with the boundary condition τ ( x ) = . Similarly to (1) this can be discretized usingupwind finite differences and then solved on Ω (see section 2.2). In all our con-vergence figures we refer to this approach as the “global factoring” (shown by asolid blue line). In Figures 2-4 it is clear that global factoring has a clean linearconvergence. Moreover, in Figure 2 it exhibits smaller errors than either the origi-nal/“unfactored” discretization or even the Lagrangian-initialized version regardlessof the grid resolution. However, as examples in Figures 3 and 4 show, on relativelycoarse grids it can be actually less accurate than the unfactored discretization – partic-ularly when the characteristics are far from straight lines. This phenomenon has not Throughout the paper we refer to this approach as “additive factoring” to stay consistent with theterminology used in prior literature. Dongping Qi, Alexander Vladimirsky been examined in prior literature, but it is hardly surprising: farther from the pointsource, u and T can be quite different, and if the second derivatives of u are smaller,this will result in larger local truncation errors when computing τ . (A) level sets Log(1/h) -8.5-8-7.5-7-6.5-6-5.5-5
Log ( e rr o r) L-inf error unit slopeoriginalglobal factoringlocal factoringLagrangian (B) L ∞ error Log(1/h) -10-9.5-9-8.5-8-7.5-7-6.5-6-5.5
Log ( e rr o r) L-1 error unit slopeoriginalglobal factoringlocal factoringLagrangian (C) L errorFig. 3: (Left) Level curves of the value function u ( x ) computed from the formula (4) with x = ( , ) , s = . , v = ( , ) . (Center & Right) The corresponding convergence plot for several discretization ap-proaches (based on the L ∞ & L errors respectively) . We further examine a “localized factoring” version of this idea. For small r , thiscould be posed as a 2-stage process: first solve (5) on B and then switch to solving(1) on Ω \ B . However, we have found that another interpretation is more suitable,particularly when characteristics are highly curved: solve (5) on the entire Ω butdefining T ( x ) = Ω \ B . We note that several recent papers have already consid-ered such hybrid/localized factoring motivated by decreasing the computational cost(since ∇ T = T when pursuing a higher-order accurate discretization [17]. Here, however,we show that the localized factoring (shown by a blue dashed lines on convergenceplots) has some accuracy advantages even with the first-order upwind discretization.In our first example (Figure 2), u ≈ T remains true on the whole Ω and characteris-tics are fairly close to straight lines; so, the global factoring is more accurate. But inFigure 3 this is no longer the case, and the localized factoring is clearly preferable.Localized factoring is also often advantageous (and more natural) when dealingwith multiple point sources. Consider, for example, the speed function specified in(4) with s = . , v = ( , ) and two point sources: x = ( , ) and x = ( , . ) . (See Figure 4.) If Q = { x } or Q = { x } , the respective value functions u and u are specified by formula (4). For the two point-sources case, Q = { x , x } , the valuefunction u ( x ) = min ( u ( x ) , u ( x )) is no longer smooth: ∇ u is undefined at the pointsfrom which the optimal paths to x and x are equally good. Since the characteristicsrun into the shockline rather than originate from it, this does not degrade the rate ofconvergence, but the rarefaction fans remain a challenge. In global factoring, there isa number of choices to capture in T the singularities at both point sources. We use T ( x ) = | x − x | F ( x ) , T ( x ) = | x − x | F ( x ) , T ( x ) = min (cid:0) T ( x ) , T ( x ) (cid:1) , ynamic Factoring 7 but the convergence results based on T = T + T are qualitatively similar (thoughthe L error becomes significantly larger). For the localized factoring, we simplyset T ( x ) = T i ( x ) on B r ( x i ) and T ( x ) = L ∞ and L norms. (A) level sets Log(1/h) -9-8.5-8-7.5-7-6.5-6-5.5-5
Log ( e rr o r) L-inf error unit slopeoriginalglobal factoringlocal factoringLagrangian (B) L ∞ error Log(1/h) -10.5-10-9.5-9-8.5-8-7.5-7
Log ( e rr o r) L-1 error unit slopeoriginalglobal factoringlocal factoringLagrangian (C) L errorFig. 4: (Left) Level curves of the value function u ( x ) for s = . , v = ( , ) and two point sources at x =( , ) and x = ( . , ) . (Center & Right) The corresponding convergence plot for several discretizationapproaches (based on the L ∞ & L errors respectively) . U i starting from each node x i is posed using the min-time-to-exit values ( U j ’s)at its neighboring nodes ( x j ’s). Dijkstra’s method [9] is perhaps the most famousof label-setting algorithms for graphs with positive edge-weights. It is based on theidea of monotone causality : an optimal path from x i starts with a transition to someneighboring node x j ∗ and, since all edge weights are positive, this implies U i > U j ∗ . Thus, even if we don’t know x j ∗ a priori, U i can be still computed based on the setof all smaller neighboring values. Dijkstra’s method exploits this observation to dy-namically uncover the correct node-ordering, decoupling the system of equations for U i ’s. The nodes are split into three sets: Accepted (with permanently fixed U values),Considered (with tentatively computed U values) and Far (with U values assumedto be + ∞ ). At each stage of the algorithm, the smallest of Considered U values isdeclared Accepted and the values at its immediate not-yet-Accepted neighbors arerecomputed. On a graph with M nodes and bounded node connectivity, this yieldsthe overall complexity of O ( M log M ) due to the use of a heap to sort the Consid-ered values. An additional useful feature of this approach is that the U values arefixed/accepted after a small number of updates on an (incrementally growing) part of Dongping Qi, Alexander Vladimirsky the graph. Many label-correcting algorithms aim to mimic this property, but withoutusing expensive data structures to sort any of the values. Unlike in label-setting al-gorithms, they cannot provide an a priori upper bound on the number of times each U i might be updated. As a result, their worst-case complexity is O ( M ) , but on manytypes of graphs their average-case behavior has been observed to be at least as goodas that of label-setting techniques [4].For Eikonal PDEs discretized on Cartesian grids, the Dijkstra-like approach wasfirst introduced in Tsitsiklis’s Algorithm [29] and Sethian’s Fast Marching Method(FMM) [21]. The latter was further extended to simplicial mesh discretizations in R n and on manifolds [12, 24], to higher-order accurate numerical schemes [22, 24], andto Hamilton-Jacobi-Bellman equations, which can be viewed as anisotropic general-izations of the Eikonal [1, 8, 18, 25, 26]. The major difficulty in applying Dijkstra-likeideas in continuous setting is that, unlike in problems on graphs, a value of u at a grid-point x i , j will depend on several neighboring values used to approximate ∇ u ( x i , j ) . Toobtain the same monotone causality, this u ( x i , j ) has to be always larger than all ofsuch contributing u values adjacent to x i , j . This property is enjoyed by only somediscretiations of (1), including the first-order accurate upwind scheme (2): a directverification shows that u x and u y are never approximated using any neighbors largerthat u i , j . Another popular class of efficient Eikonal solvers is Fast Sweeping Methods(FSM) [32], which solve the system (2) by Gauss-Seidel iterations, but changing theorder in which the gridpoints are updated from iteration to iteration. When the di-rection of the current sweep is aligned with the general direction of characteristics,many gridpoints will receive correct values in a single “sweep”. If marching-typemethods attempt to uncover the correct gridpoint-ordering dynamically, in FSM theidea is to alternate through a number of geometrically motivated orderings. (In 2Dproblems: from northeast, from northwest, from southwest, and from southeast.) Theresulting Eikonal solvers have O ( M ) complexity, but with a hidden constant factor,which depends on F and the grid orientation, and cannot be bound a priori. Thesetechniques have also been extended to anisotropic (and even non-convex) Hamilton-Jacobi equations [11, 28], as well as higher-order finite-difference (e.g., [31]) anddiscontinuous Galerkin (e.g., [30]) discretizations. Hybrid two-scale methods, com-bining the best features of marching and sweeping, were more recently introduced in[5, 6]. We also refer to [5] for a comprehensive review of other fast solvers inspiredby label-correcting algorithms.2.2 Modified Fast Marching for factored EikonalWe start by simplifying the original upwind discretization scheme (2) for the un-factored Eikonal. Focusing on a gridpoint x i , j we define its smallest horizontal andvertical neighboring values: u H = min ( u i − , j , u i + , j ) and u V = min ( u i , j − , u i , j + ) . Ifboth of these values are needed to compute u i , j , then (2) becomes equivalent to aquadratic equation ( u i , j − u H ) + ( u i , j − u V ) = h / F i , j . We are interested in its small-est real root satisfying an upwinding condition u i , j ≥ max ( u H , u V ) . If there is no rootsatisfying it, then u i , j should instead be computed from a one-sided-update u i , j = ynamic Factoring 9 min ( u H , u V )+ h / F i , j , which corresponds to the case where either max { D − xi , j u , − D + xi , j u , } or max { D − yi , j u , − D + yi , j u , } evaluates to zero. This procedure is monotone causal byconstruction and its equivalence to (2) was demonstrated in [21], making a Dijkstra-like computational approach suitable.Recalling the “additive factoring” ansatz u ( x ) = T ( x ) + τ ( x ) , we now definethe upwind vertical and horizontal neighboring values for the new unknown τ i , j butbasing the comparison on u rather than on τ itself and using the flags k H , k V ∈ {− , } to identify the selected neighbors. More specifically, (cid:40) τ H = τ i − , j and k H = , if ( T i − , j + τ i − , j ) < ( T i + , j + τ i + , j ) ; τ H = τ i + , j and k H = − , otherwise;with ( τ V , k V ) similarly defined based on the vertical neighbors. Since the partialderivatives of T are known, the corresponding quadratic equation is (cid:18) k H ∂ T i , j ∂ x + τ i , j − τ H h (cid:19) + (cid:18) k V ∂ T i , j ∂ y + τ i , j − τ v h (cid:19) = F i , j . (6)We are interested in its smallest real root satisfying a similarly modified upwindingconditionT i , j + τ i , j ≥ max (cid:18) min k = ± (cid:8) T i + k , j + τ i + k , j (cid:9) , min k = ± (cid:8) T i , j + k + τ i , j + k (cid:9)(cid:19) . (7)If there is no such root, then τ i , j should instead be computed from a one-sided-update as the smaller of the two values corresponding to k H ∂ T i , j ∂ x + τ i , j − τ H h = F i , j , and k V ∂ T i , j ∂ y + τ i , j − τ V h = F i , j . In other words, τ i , j = min (cid:26) τ H − hk H ∂ T i , j ∂ x , τ V − hk V ∂ T i , j ∂ y (cid:27) + hF i , j . (8)The above is a full recipe for computing τ i , j if all of the neighboring grid values arealready known. But since τ is a priori known only on Q , this yields a large coupledsystem of discretized equations (one per each gridpoint in Ω \ Q ). This system wasfirst treated iteratively via Fast Sweeping [16], but the monotone causality makes aDijkstra-like approach applicable as well. As in the original Fast Marching, the Con-sidered gridpoints are sorted using a binary heap with a O ( M log M ) computationalcomplexity; however, the sorting criterion is based on u rather than on τ values. Theresulting method is summarized in Algorithm 1. It is very similar to a modified FMMrecently introduced for the case of “multiplicative factoring” in [27]. Algorithm 1:
Modified Fast Marching Method
Input: source point x , speed function F ( x ) Initialize τ ( x ) : = τ ( x ) : = + ∞ for all gridpoints x (cid:54) = x ;Initialize Considered : = { x } and Accepted : = /0; while Considered (cid:54) = /0 do Find a Considered gridpoint ˆ x whose ( T + τ ) value is the smallest;Move ˆ x to Accepted list; for all not-yet-Accepted neighbors x i , j of ˆ x do Find τ H , τ V and solve the quadratic equation (6) for τ newi , j ; if τ newi , j does not satisfy the upwinding condition (7) then Use a one-sided update formula (8) to (re-)compute τ newi , j ; endif τ newi , j < τ i , j then Set τ i , j : = τ newi , j ; endendend Even though all prior work on factored Eikonal equation was focused on isolatedpoint sources, there are other well-known situations where rarefaction fans can arise.As a simple example in Figure 1 shows, they can easily develop at the corners ofobstacles (which are viewed as a part of R \ Ω ) or, more generally at any points on ∂ Ω where the boundary is non-smooth and the interior angle is larger than π . Thesenon-point-source rarefaction fans result in a similar degradation of convergence ratefor standard numerical methods and also lead to unpleasant artifacts in optimal trajec-tory approximations obtained by following ( − ∇ u ) to the target set Q . Figure 5 showsseveral such trajectories in a “maze navigation” problem. All of these trajectoriesshould be piecewise-linear, with their directions only changing at obstacle corners. Azoomed version in Figure 5(B) clearly shows that they often approach an obstacle tooearly, following its boundary to the corner and yielding longer paths. Similar artifactsare common in determining parts of the domain visible by an observer [23] and inmultiobjective path-planning [13, 19]. A natural question (and the focus of this paper)is whether factoring techniques can be used to alleviate this problem. In section 3.1we demonstrate experimentally that the “global factoring” is not suitable for this task.On the other hand, the localized factoring works, but adopting it to corner-inducedrarefaction fans presents two new challenges. First of all, not all obstacle cornersproduce this effect; e.g., see the lower left corner in Figure 1(A). Definition 3.1.
An obstacle corner ˜ x is “regular” if the characteristic leading to itfrom Q points into that obstacle. (I.e., if an optimal trajectory starts from ˜ x in thedirection a , then ( − a ) should point into the obstacle.) An obstacle corner is “rarefy-ing” if it is not regular. ynamic Factoring 11(A) multiple obstacles (B) distorted trajectoriesFig. 5: A maze navigation example: non-permiable obstacles with F = Ω . (Left) The levelsets of the value function u computed by the Fast Marching Method on a 240 ×
240 grid and approximateoptimal trajectories to the origin from 12 starting locations. (Right) A zoomed version to highlight theincorrect direction of “optimal” trajectories in the rarefaction fans at obstacle corners.
So, even though the rarefying corners are not known in advance, we can identifythem dynamically, checking the above condition when the corresponding corner ˜ x becomes Accepted in Fast Marching Method and approximating the optimal a = − ∇ u ( ˜ x ) | ∇ u ( ˜ x ) | using ˜ x ’s previously Accepted upwind neighbors. The resulting “just-in-timelocalized factoring” method is detailed in Algorithm 2. It maintains a list of identifiedrarefaction fans, with each entry containing the center of the fan (either a point sourceor a rarefying corner) and the corresponding localized function T . The algorithm isformulated in terms of u , with the corresponding τ values computed on the fly oncethe appropriate localized T is selected. (A) level sets of u (B) domain splitting (C) level sets of T Fig. 6: A simple example with one rarefying corner. (Left) level curves of the domain-restricted distanceto a point source. (Center) a dynamic domain splitting based on the rarefaction fan. (Right) the level setsof a “cone+plane” function T capturing the correct rarefaction behavior. The second difficulty is to define a suitable T that will be used for factoring whenupdating all not-yet-Accepted gridpoints in B r ( ˜ x ) . Intuitively, it might seem that a
Algorithm 2:
Just-in-time Localized Factoring
Input: source point x , speed function F ( x ) , fixed radius r Initialize u ( x ) : = u ( x ) : = + ∞ for all gridpoints x (cid:54) = x ;Initialize Considered : = { x } and Accepted : = /0;Initialize FanList : = { ( x , T x = | x − x | / F ( x )) } ; while Considered (cid:54) = /0 do Find a Considered gridpoint ˆ x whose u value is the smallest;Move ˆ x to Accepted list; if ˆ x is a rarefying corner then Build a suitable T ˆ x using formula (9);Add an entry (cid:0) ˆ x , T ˆ x (cid:1) to FanList; endfor all not-yet-Accepted neighbors x i , j of ˆ x do Check if x i , j is within distance r from any ˜ x on FanListand identify the appropriate T = T ˜ x (use T = u values at x i , j & its neighbors,define their τ values as τ = ( u − T ) based on T selected for x i , j ;Find τ H , τ V and solve the quadratic equation (6) for τ newi , j ; if τ newi , j does not satisfy the upwinding condition (7) then Use a one-sided update formula (8) to (re-)compute τ newi , j ; end u newi , j : = τ newi , j + T i , j ; if u new < u i , j then Set u i , j : = u newi , j ; endendend cone-like T = | x − ˜ x | F ( ˜ x ) is the right choice, similarly to our handling of point sources.However, as we show in section 3.1, this choice does not yield the desired rate ofconvergence. This is due to the fact that such corner-born rarefaction fans are notradially symmetric. They only exist for u > u ( ˜ x ) in the sector between a part ofobstacle boundary and the characteristic passing through ˜ x ; see Figure 1 and an evensimpler example in Figure 6. Note that, outside of that rarefaction sector, the secondderivatives of u are bounded. But using the ansatz u = T + τ with the above cone-like T would introduce unbounded second derivatives in τ on non-rarefying partsof B r ( ˜ x ) , thus degrading the rate of convergence. Therefore, we need to construct T which is cone-like only in the correct sector and remains smooth on the entire not-yet-Accepted portion of the domain. This yields a “cone+plane” version of T shownin Figure 6(C). Assuming that c is a unit vector bisecting the obstacle corner at ˜ x , wecan now split the plane into two sets S = { x ∈ Ω | ( x − ˜ x ) is between c and ( − a ) } ; S = Ω \ S . ynamic Factoring 13 Analytically, we can define T as follows T ( x ) = | x − ˜ x | F ( ˜ x ) , x ∈ S − a · ( x − ˜ x ) F ( ˜ x ) , x ∈ S . (9)The resulting T is not continuous along c , but this will not matter since the discon-tinuity is hidden within the obstacle. The gradient of T is also continuous wherever T is, though the second derivatives are bounded but discontinuous along ( − a ) . Ournumerical results show that this T fully recovers the first-order convergence of thenumerical solution. Remark 1.
The bisector of obstacle corner is not the only choice for c . Any direc-tions falling inside the obstacle will work just as well since the idea is to “hide” thediscontinuity line of T . For rectangular obstacles, it might feel more natural to choose c orthogonal to the characteristic direction a . However, we prefer the bisector simplybecause it is a safe choice for arbitrary polygonal obstacles, which can be handledby a version of Algorithm 2 on (obstacle-fitted) triangulated meshes.Another possibility is to hide the T ’s discontinuity in the part of Ω alreadyAccepted (e.g., along a ) by the time this rarefying corner is identified. This is theapproach we use in section 4, when dealing with “slowly permeable obstacles”. Remark 2.
Our approach uses local factoring with a continuous T on the entireB r ( ˜ x ) ∩ Ω . A variant of the same idea is to employ local factoring with a standardcone-like T = | x − ˜ x | F ( ˜ x ) but only when updating gridpoints from a subset B r ( ˜ x ) ∩ Ω ∩ S .(The difference is that one would use T = when updating gridpoints in B r ( ˜ x ) ∩ Ω ∩ S .) While we do not formally include this variant in our convergence studies insection 3.1, its performance appears to be quite similar. E.g., for the above examplefrom Figure 6, the alternative version also shows the first-order convergence, but withL ∞ errors ≈ larger than those resulting from the “cone+plane” formula (9) . We close this section by discussing a subtle property implicitly used in our approach.The above construction relies on having a sufficiently accurate representation of thecharacteristic direction a at each rarefying corner. This might seem unreasonable: ifour numerical approximation of u is only O ( h ) accurate (as is the case in formulas(2) and (6-8)), then one could expect the resulting finite difference approximationof ∇ u to be completely inaccurate. The same argument would imply that optimaltrajectories also cannot be reliably approximated based on any first-order accuraterepresentation of the value function. However, there is ample experimental evidencethat such trajectory approximation works in practice. See, for example, the optimaltrajectories in Figure 5,8,9,13 and in many optimal control and seismic imaging pub-lications with and without factoring. The fact that this gradient approximation is infact O ( h ) accurate is also confirmed by a numerical study in [3] and is instrumen-tal for constructing other (compact stencil, second-order) schemes for the Eikonal[3, 24]. A plausible explanation for this “superconvergence” phenomenon is that theerror in u -approximation is sufficiently “smooth”, resulting in a convergent ∇ u –approximation despite the use of divided differences. To the best of our knowledge,this property has not been proven for general Eikonal PDEs, though it has been rig-orously demonstrated for the distance-to-a-point computations and for constant co-efficient linear advection equations[14, Appendix B]. In our current context, we usethe same idea to conjecture that the a -dependent approximation of the localized T issufficiently accurate to recover the full first-order accuracy in u computations withdynamic factoring. This conjecture appears to be fully confirmed by the convergencerates observed in our numerical experiments throughout this paper.3.1 Numerical ExamplesAs a first numerical test, we consider a simple example from Figure 6(A): a domain-constrained distance u to the origin on Ω = [ , ] × [ , ] \ Ω ob , with a single rectan-gular obstacle Ω ob = ( , . ) × ( . , ) . Since F = , all optimal paths are piecewise-linear. According to Definition 3.1, we find that the corner at ˜ x = ( . , . ) is rarefy-ing. As a result, the problem contains two rarefaction fans: one at this corner and theother at the point source x = ( , ) . We test the accuracy of several methods:1. original:
Original (un-factored) Eikonal solved on the entire Ω with the originalFast Marching Method.2. global cone: Global factoring using T ( x ) = | x − x | / F ( x ) on the entire Ω withAlgorithm 1.3. global 2 cones: Global factoring using T ( x ) = | x − x | / F ( x ) + | x − ˜ x | / F ( ˜ x ) onthe entire Ω with Algorithm 1.4. switching cones: Global factoring using T ( x ) = | x − x | / F ( x ) until ˜ x is ac-cepted and then switch to global factoring using T ( x ) = | x − ˜ x | / F ( ˜ x ) on the restof Ω . localized cone only: Just-in-time localized factoring Algorithm 2 with T ( x ) = | x − x | / F ( x ) on B r ( x ) and another cone-like T ( x ) = | x − ˜ x | / F ( ˜ x ) on B r ( ˜ x ) .6. localized cone+plane: Just-in-time localized factoring Algorithm 2 with T ( x ) = | x − x | / F ( x ) on B r ( x ) and a dynamically defined “cone+plane” T ( x ) specifiedby formula (9) on B r ( ˜ x ) .The first four of these are included to show that the corner-induced rarefaction fans doindeed degrade the rate of convergence and the issue cannot be addressed by globalfactoring. Accuracy of all methods is tested using a range of gridsizes ( h = − k ,where k = , . . . , r = .
18. As Figure 7(B)clearly shows, only the last method actually exhibits the first-order of convergence.Even though the usual global factoring (method 2) starts out with smaller errors oncoarser meshes, it becomes worse than our preferred approach (method 6) for smallervalues of h . The fact that method 5 has a similar performance degradation proves theimportance of choosing the correct localized factoring function T . ynamic Factoring 15(A) level sets Log(1/h) -7.5-7-6.5-6-5.5-5-4.5-4-3.5
Log ( e rr o r) unit slopeoriginalglobal coneglobal 2 conesswitching coneslocal cone onlylocal cone+plane (B) L ∞ errorFig. 7: Convergence of several methods for a simple obstacle case of Figure 6(A). Since we do not rely on knowing ahead of time which corners are rarefying,the just-in-time localized factoring algorithm is excellent for “maze-navigation prob-lems” with numerous obstacles and possibly inhomogeneous speed function F . Be-fore running the algorithm, we create a binary array to identify all gridpoints con-tained inside obstacles. This array is then used to identify obstacle corners in Al-gorithm 2, and whenever a corner is found to be rarefying by Definition 3.1, a fac-toring procedure is locally applied with an appropriately chosen “additive factor” T ( x ) . Below we show two examples based on a “maze” from Figure 5. In Figure8 we explore the version with F = . In Figure 9 we use an inhomogeneous speed F ( x , y ) = + . ( π x ) sin ( π y ) . In both cases, the white lines are used to indicatethe (local) boundaries of corner-induced rarefaction fans. Twelve sample trajectoriesare shown to demonstrate that the trajectory distortions near the corners are avoidedby just-in-time factoring. The convergence is tested using the gridsizes h = − k ,where k = , . . . , , and the “ground truth” is computed on a much finer grid with h = / . Figure 10 shows that, unlike the Fast Marching Method for the originalEikonal, our approach is globally first-order accurate in both examples.
Rarefaction fans can also arise due to discontinuities in F . Here we consider a simplesubclass of such problems: a generalization of maze navigation examples from theprevious section to account for “slowly permeable obstacles”. We will assume thatobstacles are described by an open set Ω ob ⊂ ( Ω \ Q ) and the speed F is lower insidethem. In the simplest setting, F is piecewise constant with a discontinuity on ∂ Ω ob and 0 < F ob < F f ree . We will use ϒ = F f ree / F ob to measure the severity of obstacleslowdown.The following properties are relatively easy to prove for this simple type ofdiscontinuous F in 2D problems: F =
1. Level sets of u and representative optimal trajectories computed bythe original FMM (Left) and by Algorithm 2 (Right). The latter avoids obvious numerical artifacts nearrarefying corners.(A) maze: original (B) maze: localized factoringFig. 9: Navigating a maze with F ( x , y ) = + . ( π x ) sin ( π y ) . Level sets of u and representativeoptimal trajectories computed by the original FMM (Left) and by Algorithm 2 (Right). The latter avoidsobvious numerical artifacts near rarefying corners.
1. Rarefaction fans can only arise at point sources or at rarefying corners of slowlypermeable obstacles. (E.g., there are no fans arising on non-corner parts of obsta-cle boundaries.)2. When a rarefaction fan arises at an obstacle corner ˜ x , it does not propagate intothat obstacle.3. Such rarefaction fans are always confined to a sector between the characteristicdirection ( − a ) and another vector ( − b ) found from Snell’s law.4. Suppose a makes an angle α ∈ ( , π / ) with a normal to one side of a rectangularobstacle at ˜ x and ( − b ) makes an angle β ∈ ( , π / ] with a normal to the other ynamic Factoring 17 log(1/h) -5.5-5-4.5-4-3.5-3 l og ( e rr o r) unit slopeoriginallocal cone+plane (A) maze with constant F log(1/h) -5.5-5-4.5-4-3.5-3 l og ( e rr o r) unit slopeoriginallocal cone+plane (B) maze with variable F Fig. 10: L ∞ error for maze navigation examples. side of that obstacle; see Figure 11. Then these angles must satisfysin β = min (cid:16)(cid:112) ϒ − sin α , (cid:17) . (10)and the rarefaction fan takes place in a sector of angle δ = ( α + β − π ) . For the sake of brevity, we sketch the proof of the last of these only. (A) (B)Fig. 11: A rarefaction fan at the corner of a slowly permeable obstacle. (Left) α is the incidence angleof a ray from the point source to the rarefying corner ˜ x and β is the “refracted” angle from the normalon the other side of the obstacle. The rarefaction fan appears in a sector corresponding to the angle δ be-tween ( − b ) shown in red and the yellow dash-dotted line corresponding to ( − a ) . (Right) A ray refractionhappening at a point z close to ˜ x . As z → ˜ x , θ and θ will tend to α and β respectively, with the purplesegment disappearing, and yielding the ( − a , − b ) path through ˜ x in the limit. Proof.
We can reinterpret the characteristics as light rays traveling from a pointsource and refracted at the boundary of slowly permeable obstacles. We then usethe Snell’s Law to determine their changing directions.
Consider a point z close to the corner ˜ x , whose characteristic has an incidenceangle θ and refraction angle θ ; see Figure 11(B). The incidence angle of the secondrefraction is π − θ and the second refraction angle is θ . If θ < π , by Snell’s Lawthese three angles must satisfysin θ F f ree = sin θ F ob , cos θ F ob = sin θ F f ree . (11)Eliminating θ , we obtainsin θ = (cid:115)(cid:18) F f ree F ob (cid:19) − sin θ = (cid:113) ϒ − sin θ We note that this equality only makes sense if (cid:112) ϒ − sin θ ≤ . Otherwise, we willobserve the “total internal reflection” with θ = π and Snell’s Law not holding forthe θ - θ transition. So, the more accurate version of this relationship is sin θ = min (cid:16)(cid:112) ϒ − sin θ , (cid:17) . As z → ˜ x , we have θ → α , θ → β , and we recover (10)in the limit. Remark 3.
It is easy to provide a sufficient condition for the rarefaction fan fillingthe whole region between ( − a ) and the obstacle boundary (exactly as we saw in thenon-permeable case). Whenever ϒ ≥ √ , we have sin β = and hence β = π ; so,optimal trajectories from all starting positions in Ω \ Ω ob reach the exit set Q withoutpassing through Ω ob . On the other hand, in the continuous case ( ϒ = ), formula (10) implies that β = π − α , δ = and no rarefaction fan is present. Using a and b defined at a rarefying corner ˜ x in the above properties, it is naturalto split B r ( ˜ x ) into three regions: S = { x ∈ Ω | ( x − ˜ x ) is between ( − b ) and ( − a ) } ; S = { x ∈ Ω | ( x − ˜ x ) is between ( − b ) and a } ; S = Ω \ ( S ∪ S ) . We can now build a suitable (localized) factoring function T as follows: T ( x ) = | x − ˜ x | F ( ˜ x ) , x ∈ S − b · ( x − ˜ x ) F ( ˜ x ) , x ∈ S − a · ( x − ˜ x ) F ( ˜ x ) , x ∈ S . (12)Based on the shape of the graph, we refer to this function as a “cone+2 planes”;see Figure 12(C). This formulation makes T discontinuous along a ray parallel to a = − ∇ u ( ˜ x ) / | ∇ u ( ˜ x ) | , but for a sufficiently small r all gridpoints close to this ray in B r ( ˜ x ) will be already Accepted by the time we start this factoring. Both T and ∇ T are continuous along ( − a ) and ( − b ) . ynamic Factoring 19(A) level sets of u (B) domain splitting (C) level sets of T Fig. 12: A simple example with one “permeable obstacle”. (Left) level curves of the “partially refracted”distance to a point source. (Center) a dynamic domain splitting based on the rarefaction fan. (Right) thelevel sets of a “cone + 2 planes” function T capturing the correct rarefaction behavior. The function has asmall discontinuous jump along the ray parallel to a through ˜ x . F ob = √ ≈ .
894 inside the ob-stacle Ω ob = ( , . ) × ( . , ) and F f ree = Ω \ Ω ob . Based on the propertiesdiscussed above, this will result in a rarefaction fan spreading in a sector of angle δ = π between the two white dashed lines in Figure 13(A). We test the convergenceof several methods described in section 3.1 and report the results in Figure 13(B). Thenumerical experiments are conducted using gridsizes h = − k , where k = , . . . , h = / F ob (indicated in Figure 14(A)) and F f ree = ( − b ) and use two white linesegments to indicate the rarefying region. The “ground truth” is computed using h = / h = − k , k = , . . . , . Figure 14(B) demonstrates that our method reduces the errors and recovers the first-order convergence.
We have introduced a new just-in-time factoring algorithm for Eikonal equations toreduce the numerical errors due to rarefaction fans. Prior (global and localized) fac-toring algorithms were meant to deal with rarefactions arising at point sources andwe have carefully compared their accuracy in that setting. However, our main focushas been on rarefactions arising in 2D due to nonsmothness of ∂ Ω (e.g., cornersof non-permeable obstacles) or discontinuities in the speed function (e.g., cornersof “slowly-permeable” obstacles). The locations and the geometry of such rarefac-tion fans are a priori unknown. Our algorithm uncovers them dynamically and adop-tively applies the localized factoring. This dynamic aspect makes our approach nat- u ( x ) log(1/h) -7-6.5-6-5.5-5-4.5-4 l og ( e rr o r) unit slopeoriginalswitching coneslocal cone onlylocal cone+2planes (B) L ∞ errorFig. 13: Optimal trajectories and convergence for a “single permeable obstacle” example introduced inFigure 12. (Left) The level sets of u , with purple dashed lines showing the obstacle boundaries and whitedashed lines showing the rarefaction fan boundaries. Eight representative optimal trajectories shown inblack: (1) outside any region influenced by the obstacle, taking a straight line to the point source; (2,3)starting within the rarefaction fan, coinciding after reaching the rarefying corner; (4,5,6) experiencing adouble refraction; (7) starting inside the obstacle and experiencing the “total internal reflection” describedin the proof, with two different segments inside the obstacle; (8) starting inside the obstacle and expe-riencing a single refraction. [Note: the “light rays” (5-8) enter the obstacle with small incidence angles,resulting in barely changed refracted angles, so the first refraction is difficult to identify visually.](A) level sets of u ( x ) log(1/h) -5.5-5-4.5-4-3.5-3 l og ( e rr o r) unit slopeoriginallocal cone+2planes (B) L ∞ errorFig. 14: A maze with several slowly permeable obstacles. (Left) The level sets of u with dashed linesshowing the obstacle boundaries and white line segments showing the rarefaction fan boundaries. Thespeed F ob is shown inside each obstacle. ural in the Fast Marching framework. (With Fast Sweeping, one could in principlesolve the original Eikonal on the entire domain, then identify all rarefaction fans inpost-processing and re-solve the correctly factored equation on Ω . ) Numerical testsconfirm that our method restores the full linear convergence and prevents numericalartifacts in approximating optimal trajectories once the value function is already com-puted. While we have only implemented and tested the “additive” dynamic factoring,we expect that in the “multiplicative” case the results would be qualitatively simi- ynamic Factoring 21 lar. All presented examples were in the context of time-optimal path planning, butother optimization criteria (e.g., a cumulative exposure to an enemy observer) wouldalso lead to a similarly factored Eikonal equation as long as the running cost remainsisotropic. Even though our focus so far has been on applications in robotic naviga-tion and computational geometry, we hope that the same general approach might alsobe useful in seismic imaging problems, which motivated much of the prior work onEikonal factoring.Several straightforward generalizations will make our method more useful inpractice.1. We can easily treat general polygonal obstacles by adding dynamic factoring toprior Fast Marching techniques on (obstacle-fitted) triangulated meshes [12, 24].The definition of our “additive factor” T will stay exactly the same; see alsoRemark 1 in section 3.2. The examples presented in section 4 are based on rectangular “slowly permeableobstacles” with a piecewise constant speed function. However, the same approachis also applicable for the general discontinuous speed functions as long as thediscontinuity lines are polygonal and aligned with the discretization mesh. Therarefaction fans can be determined based on a local information only (i.e., thedirectional limits of the speed function at a rarefying corner of the discontinuityline), and the definition of T in dynamic factoring will remain the same evenwhen F is not piecewise constant.3. If the speed of motion is anisotropic (i.e., dependent on the direction of motionrather than just the current location), the value function satisfies a more generalHamilton-Jacobi-Bellman PDE. Point-source-based factoring for the latter has al-ready been developed (e.g., by Fast Sweeping in [16]). Marching-type techniquesfor anisotropic problems (e.g., [26] or [18]) can be similarly modified to handlethe corner-induced rarefactions.4. Another easy extension is to treat rarefaction fans due to more general boundaryconditions (e.g., fast-varying u = g specified on ∂ Ω can result in rarefactions evenif the boundary is smooth).It will be more difficult to move to factoring suitable for higher-order accurate dis-cretizations. For point-source-induced fans, this has been addressed in [15] and [17].Similar ideas might work in our context, but higher derivatives will need to be esti-mated at rarefying corners and one would need to construct a smoother T than theversion used in this paper.Finally, the obvious limitation of our current approach is that Ω ⊂ R . We expectthat Eikonal problems in higher dimension will be much harder to factor dynamically.Even with F = rarefying edges rather than corners. Acknowledgements:
The authors are grateful to anonymous reviewers for their sug-gestions on improving this paper.
References
1. Alton, K., Mitchell, I.M.: An ordered upwind method with precomputed stenciland monotone node acceptance for solving static convex Hamilton-Jacobi equa-tions. Journal of Scientific Computing (2), 313–348 (2012)2. Bardi, M., Italo, C.D.: Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations. Springer Science+Business Media, LLC (1997)3. Benamou, J.D., Luo, S., Zhao, H.: A compact upwind second order scheme forthe eikonal equation. Journal of Computational Mathematics pp. 489–516 (2010)4. Bertsekas, D.P.: Network Optimization: Continuous and Discrete Models.Athena Scientific (1998)5. Chacon, A., Vladimirsky, A.: Fast two-scale methods for eikonal equations.SIAM Journal on Scientific Computing (2), A547–A578 (2012)6. Chacon, A., Vladimirsky, A.: A parallel two-scale method for eikonal equations.SIAM Journal on Scientific Computing (1), A156–A180 (2015)7. Crandall, M.G., Lions, P.L.: Viscosity solutions of Hamilton-Jacobi equations.Transactions of the American Mathematical Society (1), 1–42 (1983)8. Dahiya, D., Cameron, M.: Ordered line integral methods for computing thequasi-potential. J. Scientific Computing (2017). In revision; https://arxiv.org/abs/1706.07509
9. Dijkstra, E.: A note on two problems in connexion with graphs. NumerischeMathematik , 269–271 (1959)10. Fomel, S., Luo, S., Zhao, H.: Fast sweeping method for the factored eikonalequation. J. Comput. Phys. (17), 6440–6455 (2009)11. Kao, C.Y., Osher, S., Qian, J.: Lax-Friedrichs sweeping scheme for staticHamilton–Jacobi equations. Journal of Computational Physics (1), 367–391(2004)12. Kimmel, R., Sethian, J.A.: Fast marching methods on triangulated domains. Pro-ceedings of the National Academy of Science , 8341–8435 (1998)13. Kumar, A., Vladimirsky, A.: An efficient method for multiobjective optimal con-trol and optimal control subject to integral constraints. Journal of ComputationalMathematics (4), 517–551 (2010)14. Luo, S.: Numerical Methods for Static Hamilton-Jacobi Equations. Ph.D. thesis,University of California at Irvine (2009)15. Luo, S., Qian, J.: Factored singularities and high-order Lax–Friedrichs sweepingschemes for point-source traveltimes and amplitudes. Journal of ComputationalPhysics (12), 4742–4755 (2011)16. Luo, S., Qian, J.: Fast sweeping methods for factored anisotropic eikonal equa-tions: multiplicative and additive factors. Journal of Scientific Computing (2),360–382 (2012)17. Luo, S., Qian, J., Burridge, R.: High-order factorization based high-order hybridfast sweeping methods for point-source eikonal equations. SIAM Journal onNumerical Analysis (1), 23–44 (2014)18. Mirebeau, J.M.: Efficient fast marching with Finsler metrics. Numerische math-ematik (3), 515–557 (2014) ynamic Factoring 23
19. Mitchell, I.M., Sastry, S.: Continuous path planning with multiple constraints.In: 42nd IEEE International Conference on Decision and Control (IEEE Cat.No.03CH37475), vol. 5, pp. 5502–5507 (2003)20. Noble, M., Gesret, A., Belayouni, N.: Accurate 3-d finite difference computationof traveltimes in strongly heterogeneous media. Geophysical Journal Interna-tional (3), 1572–1585 (2014)21. Sethian, J.A.: A fast marching level set method for monotonically advancingfronts. Proceedings of the National Academy of Sciences (4), 1591–1595(1996)22. Sethian, J.A.: Fast marching methods. SIAM review (2), 199–235 (1999)23. Sethian, J.A., Adalsteinsson, D.: An overview of level set methods for etching,deposition, and lithography development. IEEE Transactions on SemiconductorManufacturing (1), 167–184 (1997)24. Sethian, J.A., Vladimirsky, A.: Fast methods for the eikonal and relatedHamilton–Jacobi equations on unstructured meshes. Proceedings of the NationalAcademy of Sciences (11), 5699–5703 (2000)25. Sethian, J.A., Vladimirsky, A.: Ordered upwind methods for static Hamilton-Jacobi equations. Proceedings of the National Academy of Sciences (20),11,069–11,074 (2001)26. Sethian, J.A., Vladimirsky, A.: Ordered upwind methods for static Hamilton-Jacobi equations: Theory & algorithms. SIAM Journal on Numerical Analysis (1), 325–363 (2003)27. Treister, E., Haber, E.: A fast marching algorithm for the factored eikonal equa-tion. Journal of Computational Physics , 210–225 (2016)28. Tsai, Y.H.R., Cheng, L.T., Osher, S., Zhao, H.K.: Fast sweeping algorithms for aclass of Hamilton–Jacobi equations. SIAM journal on numerical analysis (2),673–694 (2003)29. Tsitsiklis, J.N.: Efficient algorithms for globally optimal trajectories. IEEETransactions on Automatic Control (9), 1528–1538 (1995)30. Zhang, Y.T., Chen, S., Li, F., Zhao, H., Shu, C.W.: Uniformly accurate discon-tinuous galerkin fast sweeping methods for eikonal equations. SIAM Journal onScientific Computing (4), 1873–1896 (2011)31. Zhang, Y.T., Zhao, H.K., Qian, J.: High order fast sweeping methods for staticHamilton–Jacobi equations. Journal of Scientific Computing (1), 25–56(2006)32. Zhao, H.: A fast sweeping method for eikonal equations. Mathematics of com-putation74