Numerical geometry of nonrigid shapes Shortest Path Problems

Numerical geometry of non-rigid shapes Shortest Path Problems 1 As regards obstacles, the shortest distance between two points can be a curve. B. Brecht Shortest Path Problems Alexander Bronstein, Michael Bronstein © 2008 All rights reserved. Web: tosca. cs. technion. ac. il

2 Numerical geometry of non-rigid shapes Shortest Path Problems How to compute the intrinsic metric? n So far, we represented itself. n Our model of rigid shapes as metric spaces involves the intrinsic metric n Sampling procedure requires as well. n We need a tool to compute geodesic distances on .

Numerical geometry of non-rigid shapes Shortest Path Problems Shortest path problem Brussels 183 Paris 346 943 566 Prague 504 1542 194 407 902 285 Bern 271 Munich 1146 Vienna 3

Numerical geometry of non-rigid shapes Shortest Path Problems Shapes as graphs n Sample the shape at vertices n Represent shape as an undirected graph n set of edges representing adjacent vertices. n Define length function measuring local distances as Euclidean ones, . 4

Numerical geometry of non-rigid shapes Shortest Path Problems Shapes as graphs n Path between where is an ordered set of connected edges and . n Path length = sum of edge lengths 5

Numerical geometry of non-rigid shapes Shortest Path Problems Geodesic distance n Shortest path between n Length metric in graph n Approximates the geodesic distance n Shortest path problem: compute between any and . n Alternatively: given a source point distance map on the shape. , compute the. 6

Numerical geometry of non-rigid shapes Shortest Path Problems 7 Bellman’s principle of optimality n Let be shortest path between and n Then, a point on the path. and are shortest sub-paths between , and . Richard Bellman (1920 -1984) n Suppose there exists a shorter path n Contradiction to . being shortest path.

8 Numerical geometry of non-rigid shapes Shortest Path Problems Dynamic programming n How to compute the shortest path between source n Bellman principle: there exists n has to minimize path length n Recursive dynamic programming equation. such that and on ?

Numerical geometry of non-rigid shapes Shortest Path Problems Edsger Wybe Dijkstra (1930– 2002) [‘ɛtsxər ‘wibə ‘dɛɪkstra] 9

10 Numerical geometry of non-rigid shapes Shortest Path Problems Dijkstra’s algorithm n Initialize and for the rest of the graph; Initialize queue of unprocessed vertices . n While n Find vertex with smallest value of , n For each unprocessed adjacent vertex n Remove from n Return distance map , . .

11 Numerical geometry of non-rigid shapes Shortest Path Problems Dijkstra’s algorithm Brussels 183 566 183 Paris Prague 749 0 407 504 346 194 271 346 Bern 617 679 Munich 285 904 Vienna

12 Numerical geometry of non-rigid shapes Shortest Path Problems Dijkstra’s algorithm – complexity n While there are still unprocessed vertices n Find and remove minimum n For each unprocessed adjacent vertex n Perform update n Every vertex is processed exactly once: outer iterations. n Minimum extraction straightforward complexity: n Can be reduced to using binary or Fibonacci heap. n Updating adjacent vertices is in general . n In our case, graph is sparsely connected, update in . n Total complexity: .

Numerical geometry of non-rigid shapes Shortest Path Problems Troubles with the metric n Grid with 4 -neighbor connectivity. n True Euclidean distance n Shortest path in graph (not unique) n Increasing sampling density does not help. 13

Numerical geometry of non-rigid shapes Shortest Path Problems 14 Metrication error 4 -neighbor topology 8 -neighbor topology Manhattan distance Continuous Euclidean distance n Graph representation induces an inconsistent metric. n Increasing sampling size does not make it consistent. n Neither does increasing connectivity.

Numerical geometry of non-rigid shapes Shortest Path Problems Metrication error n How to approximate the metric consistently? Solution 1 n Stick to graph representation. n Change connectivity and sampling. n Under certain conditions consistency is guaranteed. Solution 2 n Stick to given sampling (and connectivity). n Compute distance map on a surface in some representation (e. g. , mesh). n Requires a new algorithm. 15

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching algorithms Imagine a forest fire… 16

Numerical geometry of non-rigid shapes Shortest Path Problems Forest fire n Fire starts at a source at . n Propagates with constant velocity n Arrives at time to a point n Fermat’s (least action) principle: The fire chooses the quickest path to travel. n Governs refraction laws in optics (Snell’s law) and acoustics. n Fire arrival time distance map = from source. . 17

Numerical geometry of non-rigid shapes Shortest Path Problems Distance maps on surfaces n Distance map on surface n Mapped locally to the tangent space n A small step in the direction n changes the distance by is directional derivative in the direction . 18

Numerical geometry of non-rigid shapes Shortest Path Problems Intrinsic gradient n For some direction , n The perpendicular direction is the direction of steepest change of the distance map. n is referred to as the intrinsic gradient. n Formally, the intrinsic gradient of function at a point is a map satisfying for any 19

20 Numerical geometry of non-rigid shapes Shortest Path Problems Extrinsic gradient n Consider the distance map as a function n The extrinsic gradient of at a point . is a map satisfying for any direction n In the standard Euclidean basis n Usually called “the gradient” of . n What is the connection between intrinsic and extrinsic gradients?

Numerical geometry of non-rigid shapes Shortest Path Problems Intrinsic and extrinsic gradients n Intrinsic gradient = projection of extrinsic gradient on tangent plane n In coordinates of a parametrization , n is the Jacobian matrix whose columns span . 21

22 Numerical geometry of non-rigid shapes Shortest Path Problems Eikonal equation n Let be a minimal geodesic between and . n The derivative is the fire front propagation direction. n In arclength parametrization . n Fermat’s principle: n Propagation direction = direction of steepest increase of n Geodesic is perpendicular to the level sets of on . .

Numerical geometry of non-rigid shapes Shortest Path Problems Eikonal equation n Eikonal equation (from Greek εικων) n Hyperbolic PDE with boundary condition n Minimal geodesics are characteristics. n Describes propagation of waves in medium. 23

24 Numerical geometry of non-rigid shapes Shortest Path Problems Eikonal equation n Let be a minimal geodesic between and . n The derivative is the fire front propagation direction. n In arclength parametrization . n Fermat’s principle: n Propagation direction = direction of steepest increase of n Geodesic is perpendicular to the level sets of on . .

Numerical geometry of non-rigid shapes Shortest Path Problems Uniqueness of solution n In classic PDE theory, a solution is a continuous differentiable function satisfying n PDE theory guarantees existence and uniqueness of solution. n Distance map is not everywhere differentiable. n Solution is not unique! 1 D example 25

Numerical geometry of non-rigid shapes Shortest Path Problems We need assistance of a super-creature… 26

Numerical geometry of non-rigid shapes Shortest Path Problems Sub- and super-derivatives (1 D case) n Superderivative: the set of all slopes above the graph n Subderivative: the set of all slopes below the graph n where is differentiable. 27

Numerical geometry of non-rigid shapes Shortest Path Problems 28 Viscosity solution n is a viscosity solution of the 1 D eikonal equation if n Monotonicity: viscosity solution Not a viscosity solution does not have local maxima. n The largest among all n Existence and uniqueness guaranteed. Viscosity solution

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching methods (FMM) n A family of numerical methods for n solving eikonal equation. n Finds the viscosity solution = distance map. n Simulates wavefront propagation from a source set. n A continuous variant of Dijkstra’s algorithm. n Consistently approximate the intrinsic metric on the surface. 29

30 Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching algorithm n Initialize and mark it as black. n Initialize for other vertices and mark them as green. n Initialize queue of red vertices . n Repeat n Mark green neighbors of black vertices as red (add to ) n For each red vertex n For each triangle sharing the vertex Update n Mark from the triangle. with minimum value of as black (remove from n Until there are no more green vertices. n Return distance map . )

31 Numerical geometry of non-rigid shapes Shortest Path Problems Update step Dijkstra’s update n Vertex updated from adjacent vertex n Distance computed from n Path restricted to graph edges Fast marching update n Vertex n Distance from updated from triang computed and n Path can pass on mesh faces

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching update step n Update from triangle n Compute from and n Model wave front propagating from planar source n unit propagation direction n source offset n Front hits n Hits at time n When does the front arrive to Planar source ? 32

33 Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching update step n Assume w. l. o. g. n and . is given by the point-to-plane distance n Solve for parameters and using the point-to-plane distance n In vector notation where , n In a non-degenerate triangle matrix , and is full-rank .

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching update step n Apparently, we have two equations with three variables. n However, is a unit vector, hence where . n Substitute . and obtain a quadratic equation 34

Numerical geometry of non-rigid shapes Shortest Path Problems Causality condition n Quadratic equation is satisfied by both and . n Two solutions for n Causality: front can propagate only forward in time. n Causality condition 35

Numerical geometry of non-rigid shapes Shortest Path Problems Causality condition n Causality condition In other words n has to form obtuse angles with both triangle edges . n Causality is required to obtain consistent approximation of the distance map. n Smallest solution for is inconsistent and is discarded. n If largest solution is consistent, live the largest solution! 36

37 Numerical geometry of non-rigid shapes Shortest Path Problems Monotonicity condition n Viscosity solution has to be a monotonically increasing function. n Monotonicity condition: increase when In other words: n Differentiate w. r. t obtaining or increase.

38 Numerical geometry of non-rigid shapes Shortest Path Problems Monotonicity condition n Substitute n Monotonicity satisfied when both coordinates of have the same sign. n is positive definite n Causality condition: n Monotonicity condition: At least one coordinate of is negative

Numerical geometry of non-rigid shapes Shortest Path Problems Monotonicity condition n Since n Rows of we have are orthogonal to triangle edges n Monotonicity condition: Geometric interpretation: n must form obtuse angles with normals to triangle edges. Said differently: n must come from within the triangle. 39

Numerical geometry of non-rigid shapes Shortest Path Problems One-sided update n Monotonicity condition: update direction must come from within the triangle. n If it does not, project n inside the triangle. will coincide with one of the edges. n Update will reduce to Dijkstra’s update or 40

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching update n Solve for the quadratic equation n Compute propagation direction n If monotonicity condition n Set is violated, 41

Numerical geometry of non-rigid shapes Shortest Path Problems 42 Consistency and monotonicity encore cy Co Co n si st en ns Monotonicity ist en cy Acute triangle Obtuse triangle All directions in the triangle Some directions in the triangle satisfy consistency and violate consistency condition! monotonicity conditions.

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching on obtuse meshes n Inconsistent solution if the mesh contains obtuse triangles n Remeshing is costly n Solution: split obtuse triangles by adding virtual connections to non-adjacent vertices n Done as a pre-processing step in 43

Numerical geometry of non-rigid shapes Shortest Path Problems ® MATLAB intermezzo Fast marching 44

Numerical geometry of non-rigid shapes Shortest Path Problems 45

Numerical geometry of non-rigid shapes Shortest Path Problems Eikonal equation on parametric surfaces n Parametrization of n Compute distance map from source over , . n Chain rule n Extrinsic gradient in parametrization coordinates n Intrinsic gradient in parametrization coordinates . 46

Numerical geometry of non-rigid shapes Shortest Path Problems Eikonal equation on parametric surfaces n Eikonal equation in parametrization coordinates 47

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching on parametric surfaces n Solve eikonal equation in parametrization domain n March on discretized parametrization domain. n We need to express update step in parametrization coordinates. 48

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching on parametric surfaces n Cartesian sampling of with unit step. n Some connectivity (e. g. 4 - or 8 -neighbor). n Vertex updated from triangle n Assuming w. l. o. g. or in matrix form 49

Numerical geometry of non-rigid shapes Shortest Path Problems Fast marching on parametric surfaces n Inner product matrix n Describes triangle geometry. n lengths of the edges. n n Substitute cosine of the angle. into the update quadratic equation n Only first fundamental form coefficients and grid connectivity are required for update. n Can measure distances when only surface gradients are known. 50

51 Numerical geometry of non-rigid shapes Shortest Path Problems Heap-based grid update n Fast marching and Dijkstra’s algorithm use heap-based grid update. n Next vertex to be updated is decided by extracting the smallest . n Update order is unknown and data-dependent. n Inefficient use of memory system and cache. n Inherently sequential algorithm – next update depends on previous one. n Can we do better? n Regular access to memory (known in advance). n Vectorizable (parallelizable) algorithm.

Numerical geometry of non-rigid shapes Shortest Path Problems Marching even faster n Danielsson’s algorithm: update the grid in a raster scan order n In Euclidean case, parametrization is trivial. n Geodesics are straight lines in parametrization domain. n Each raster scan covers ¼ of the possible directions of the geodesics. n Euclidean distance map computed by four alternating raster scans. 52

53 Numerical geometry of non-rigid shapes Shortest Path Problems Raster scan fast marching n Generally, geodesics are curved in parametrization domain. n Raster scans have to be repeated to 1 iteration 4 iterations 2 iterations 5 iterations 3 iterations 6 iterations produce a convergent solution. n Iterative algorithm. n Number of iterations depends on geometry and parametrization. n Practically, few iterations are required.

Numerical geometry of non-rigid shapes Shortest Path Problems Raster scan fast marching n What we lost: n No more a one-pass algorithm. n Computational complexity is data-dependent. n What we found: n Coherent memory access, efficient use of cache. n No heap, each iteration is . n Raster scans can be parallelized. BBK, "Parallel algorithms for approximation of distance maps on parametric surfaces”, 2007 54

Numerical geometry of non-rigid shapes Shortest Path Problems Parallellization n Rotate scan directions by 450. n All updates performed along a row or column can be parallelized. n Constant CPU load – suitable for SIMD architecture and GPUs. 55

Numerical geometry of non-rigid shapes Shortest Path Problems Parallel marching n Rotate scan directions by 450. n All updates performed along a row or column can be parallelized. n Constant CPU load. n Suitable for SIMD architecture and GPUs. n GPU implementation computes geodesic on grid with 10, 000 vertices in less than 50 msec. n About 200 million distances per second! 56

Numerical geometry of non-rigid shapes Shortest Path Problems Minimal geodesics n We have a numerical tool to compute geodesic distance. n Sometimes, the shortest path itself is needed. n Minimal geodesics are characteristics of the eikonal equation. In other words: n Along geodesic, eikonal equation becomes an ODE with initial condition n Solve the ODE for . . 57

Numerical geometry of non-rigid shapes Shortest Path Problems 58 Minimal geodesics n To find a minimal geodesic between two points n Compute distance map from n Starting at to all other points. , follow the direction of n Steepest descent on the distance map. In the parametrization coordinates n Let be the preimage of in until is reached.

Numerical geometry of non-rigid shapes Shortest Path Problems Minimal geodesics n Substitute into characteristic equation n Steepest descent on surface = scaled steepest descent in parametrization domain. 59

60 Numerical geometry of non-rigid shapes Shortest Path Problems Uses of fast marching Geodesic distances Minimal geodesics Voronoi tessellation & sampling Offset curves