Ordered Line Integral Methods for Solving the Eikonal Equation Samuel F. Potter Maria K. Cameron July 17, 2019
The eikonal equation Solve �∇ u ( x ) � = s ( x ) for u , where s is the slowness function ( c = 1 /s is speed of sound/light ) Where does it come from? Reduced wave (Helmholtz) equation: ∆ + ω 2 c 2 w = 0 Ansatz (WKB): w ( x ) ∼ A ( x ) e iωu ( x ) → yields eikonal (and transport equation) 1/25
Ordered line integral methods (OLIMs) Higher accuracy through use of midpoint quadrature rules and increased directional coverage Semi-Lagrangian = ⇒ parametrize local characteristics Use local characteristic information to skip updates (important in 3D) More tricks that work for the eikonal equation but not the quasipotential Related work OLIMs for computing the quasipotential [DC ‘18] concurrent : quasipotential OLIM in 3D [YPC ‘19] 2/25
History of related algorithms eikonal solvers label-setting algorithms iterative solvers Dijkstra-like Dial-like Bou´ e & Dupuis [‘99] Tsitsiklis [‘95] fast sweeping OLIMs for Sethian [‘96] method [Zhao] eikonal [PC ‘19] 3/25
Semi-Lagrangian updates x ∂ Ω ϕ � � � u ( y ) = min x,ϕ u ( x ) + ϕ s y Ω Approximate and minimize Approximate using simple quadrature rules over straight line segments Minimize using linear algebra + geometry or numerical optimization 4/25
Stencils, neighborhoods, other methods Tsitsiklis (2D) FMM (2D) Tsitsiklis (3D) FMM (3D) 5/25
A generic Dijkstra-like eikonal solver Grid of N d nodes ( d = 2 , 3 ) Compute numerical eikonal U Boundary data: U ≡ 0 on a subset of grid points States: far , trial , valid Use a priority queue to sort trial nodes by U value trial Until done: → Remove minimum node from priority queue → Set node state to valid → Set far node neighbors’ states to trial → Update neighbors our focus 6/25
Contributions Approximate with simple quadrature rules: mp0 , mp1 , rhr Update using efficient algorithms: top-down and bottom-up Skip redundant updates Exact semi-Lagrangian solution for mp0 / rhr in all dimensions Theorem relating mp0 and mp1 Numerical results + code profiling 7/25
OLIM Classification eikonal OLIMs update alg. bottom-up top-down quad. rule mp0 mp1 rhr mp0 mp1 rhr num. neighb. 26 26 26 6 18 26 6 18 26 6 18 26 o o o OLIM name o o o o o o o o o l l l l l l l l l l l l i i i i i i i i i i i i m m m m m m m m m m m m 3 3 3 6 1 2 6 1 2 6 1 2 d d d 8 6 8 6 8 6 m m r m m r p m m p m m h r r p p h 0 p p 1 p p r h h 0 1 r 0 0 1 1 r r 8/25
OLIM Classification eikonal OLIMs update alg. bottom-up top-down quad. rule mp0 mp1 rhr mp0 mp1 rhr num. neighb. 26 26 26 6 18 26 6 18 26 6 18 26 o o o OLIM name o o o o o o o o o l l l l l l l l l l l l i i i i i i i i i i i i m m m m m m m m m m m m 3 3 3 6 1 2 6 1 2 6 1 2 d d d 8 6 8 6 8 6 m m r m m r p m m p m m h r r p p h 0 p p 1 p p r h h 0 1 r 0 0 1 1 r r FMM 3D Tsitsiklis 3D 8/25
OLIM Classification eikonal OLIMs update alg. bottom-up top-down quad. rule mp0 mp1 rhr mp0 mp1 rhr num. neighb. 26 26 26 6 18 26 6 18 26 6 18 26 o o o OLIM name o o o o o o o o o best l l l l l l l l l l l l i i i i i i i i i i i i m m m m m m m m m m m m 3 3 3 6 1 2 6 1 2 6 1 2 d d d 8 6 8 6 8 6 m m r m m r p m m p m m h r r p p h 0 p p 1 p p r h h 0 1 r 0 0 1 1 r r 8/25
Quadrature Rules: Setup approximate U linearly �� L � �� U (ˆ x ) = min 0 ≤ λ ≤ 1 (1 − λ ) U ( x 0 ) + λU ( x 1 ) + Q 0 s ( ϕ ) dl ˆ U U λ s, ˆ x, ˆ ˆ U Approximate integral with quadrature rule Q : rhr s � ˆ ˆ x − x λ � ϕ � ˆ s + s λ � � ˆ x − x λ � mp1 x 1 , s 1 , U 1 2 � � s + s 0+ s 1 ˆ mp0 � ˆ x − x λ � 2 x λ , s λ , U λ 2 x 0 , s 0 , U 0 9/25
Quadrature Rules: Cost Functions With ∆ d = { λ ∈ R d : λ i ≥ 0 for i = 1 , . . . , d and � d i =1 λ i ≤ 1 } : �� L � �� U (ˆ x ) = min λ ∈ ∆ d U λ + Q 0 s ( ϕ ) dl =: F ( λ ) For each quadrature rule ( rhr , mp0 , mp1 ): F rhr ( λ ) = u λ + ˆ s � ˆ x − x λ � � � � d 1 s + ˆ i =0 s i F mp0 ( λ ) = u λ + d +1 � ˆ x − x λ � 2 � ˆ s + s λ � F mp1 ( λ ) = u λ + � ˆ x − x λ � 2 ˆ U = min λ ∈ ∆ d F ( λ ) 10/25
Quadrature Rules: θ -rules Cost functions F rhr , F mp0 , F mp1 generalized by: � � � d θ F 0 ( λ ) = u λ + (1 − θ )ˆ s + i =0 s i � ˆ x − x λ � d +1 F 1 ( λ ) = u λ + [(1 − θ )ˆ s + θs λ ] � ˆ x − x λ � θ = 0 Specifically: s ˆ s θ θ = 1 / 2 F rhr = F 0 , θ = 0 = ⇒ F mp0 = F 0 , θ = 1 / 2 s 1 s 0 + s 1 2 F mp1 = F 1 , θ = 1 / 2 s λ s 0 λ 11/25
E.g., setup for tetrahedron update p ˆ s = s (ˆ ˆ p ) s θ s θ λ δp 2 p 0 p 2 s 0 s 2 δp 1 p λ = p 0 + δPλ s λ = s ( p λ ) s 1 p 1 s i = s ( p i ) � p 1 − p 0 � � δp 1 � p 2 − p 0 δp 2 δP = = 12/25
Exact solution for F 0 ( F rhr and F mp0 ) � p 1 − p 0 � Theorem : Let δP = · · · p d − p 0 , and let QR = qr ( δP, 0 ) and λ ∗ = argmin λ ∈ R d F 0 ( λ ) . Then: � � p ⊤ � 0 ( I − QQ ⊤ ) p 0 � � p λ ∗ � = � 2 � R −⊤ δU � � 1 − s θ h � � Q ⊤ p 0 + � p λ ∗ � R −⊤ δU λ ∗ = − R − 1 s θ h ˆ U = F ( λ ∗ ) Works in any dimension for any simplex, parametrizes local characteristic in addition to computing ˆ U . Proof: geometry & linear algebra. 13/25
Problem with continuity for F mp0 p ˆ p 0 p 2 p 2 p 1 14/25
Problem with continuity for F mp0 p ˆ p 0 p 2 p 2 p 1 Set ˆ U = F mp1 ( λ ∗ mp0 ) 14/25
Validating mp0 Theorem : Let h be small enough that F 1 is strictly convex, let λ ∗ 0 and λ ∗ 1 be (unconstrained) minimizers of θ -rules F 0 and 0 ) | = O ( h 3 ) . F 1 . Then | F 1 ( λ ∗ 1 ) − F 1 ( λ ∗ Proof : Kantorovich’s theorem + 3 lemmas for verification of conditions. Compute λ ∗ 0 using QR decomposition, then set ˆ U ← F mp1 ( λ ∗ 0 ) (i.e., use F mp1 to evaluate). ∆ 2 15/25
Solution for mp1 Nothing fancy: − → For tri updates, use hybrid method − → For tetra updates, use sequential quadratic programming Our implementation is carefully tested and reasonably optimized ∆ 2 16/25
A few other things... Some relevant things in the paper which I won’t go into here Convexity of θ -rules F 0 and F 1 Causality (of individual updates) for F 0 and F 1 General formula for the update gap of a simplex in d dimensions (for Dial-like algorithm) Simple idea : ∇ 2 F 1 = ∇ 2 F 0 + “small indefinite perturbation” Hence, � perturbation � = O ( h 2 ) lets proofs for F 0 apply to F 1 for h small enough 17/25
Local Factoring At point sources and rarefaction fans, order of convergence can degrade from O ( h ) to O ( h log 1 h ) — see also Qi’s talk Very easy to modify our cost functions to work with the additively factored eikonal equation rarefaction fan obstacle point source Qi, Dongping, and Alexander Vladimirsky. ”Corner cases, singularities, and dynamic factoring.” Journal of Scientific Computing 79.3 (2019): 1456-1476. 18/25
Enumerating Updates (top-down) Manually enumerate update tetrahedra and triangles = olim6 neighbors + = olim18 neighbors + + = olim26 neighbors Consider only one octant at a time (by symmetry) Can then enumerate small number of p ˆ tetrahedra (hence triangle) updates 19/25
Enumerating Updates (top-down) p ˆ IVa IVb p ˆ I II III V VIa VIb 19/25
Enumerating Updates (top-down) olim6 p ˆ IVa IVb p ˆ I II III V VIa VIb 19/25
Enumerating Updates (top-down) olim18 p ˆ IVa IVb p ˆ I II III V VIa VIb 19/25
Enumerating Updates (top-down) p ˆ IVa IVb p ˆ I II III V VIa VIb olim26 19/25
Skipping Updates (top-down) Case 1: F 0 (i.e. F rhr and mp0 ) = ⇒ exact solution via QR ∈ ∆ d , check its If λ ∗ 0 / 1 λ 2 location For d = 2 , easy to divide p 2 R 2 \ ∆ 2 into disjoint “skip 2 2 λ ∗ zones” 0 ∆ 2 Skip zone tells us which λ 1 lower-dimensional p 1 p 0 1 updates to skip 1 2 skip Case 2: F 1 (i.e. F mp1 ) = ⇒ constrained optimization Skip all incident lower-dimensional updates (i.e., after a tetra update, skip all incident tri and line updates (6 total)) 20/25
Enumerating Updates (bottom-up) ( olim3d ) Start with p 0 = p new , do line update Check gradient to see if line update “minimizer” is minimizer of adjacent triangle updates p ˆ Check if tri update minimizer λ ∗ is minimizer of for adjacent tetrahedron updates (using KKT conditions) λ ∈ ∆ d is a simple linear matrix inequality p 2 Check if tri update argmin is argmin of tetra update by computing only one Lagrange multiplier p 0 = p new p 1 21/25
Recommend
More recommend
