Fast Marching on Triangulated Domains

1. Computer Science Department • Technion-Israel Institute of Technology Fast Marching on Triangulated Domains Ron Kimmel www.cs.technion.ac.il/~ron • Geometric Image Processing Lab

2. Brief Historical Review • Upwind schemes: Godunov 59 • Level sets: Osher & Sethian 88 • Viscosity SFS: Rouy & Tourin 92, (Osher & Rudin) • Level sets SFS: Kimmel & Bruckstein 92 • Continuous morphology: Brockett & Maragos 92,Sapiro et al. 93 • Minimal geodesics: Kimmel, Amir & Bruckstein 93 • Fast marching method: Sethian 95 • Fast optimal path: Tsitsiklis 95 • Level sets on triangulated domains:Barth & Sethian 98 • Fast marching on triangulated domains: Kimmel & Sethian 98 • Applications based on joint works with: Elad, Kiryati, Zigelman

3. 1D Distance: Example 1 • Find distance T(x), given T(x0)=0. • Solution: T(x)=|x-x0|. • except at x0. T(x) x x0

4. Find the distance T(x), given T(x0)=T(x1)=0 Solution: T(x) = min{|x-x0|,|x-x1|}, Again, , except x0,x1and(x0+x1)/2. 1D Distance: Example 2 T(x) x x0 x1

5. 1D Eikonal Equation with boundary conditions T(x0)=T(x1)=0. • Goal: Compute T that satisfies the equation `the best'. T(x) x x0 x1

6. Numerical Approximation • Restrict , where h= grid spacing. • Possible solutions for are T(x) x x0 x1

7. Approximation II • Updated i has always • `upwind' from where the `wind blows' T(x) x x0 x1

8. Update Procedure Set , and T(x0)=T(x1)=0. REPEAT UNTIL convergence, • FOR each i T i T i+1 T i-1 i i-1 i+1 h

9. Update Order What is the optimal order of updates? Solution I: Scan the line successively left to right. N scans, i.e. O(N ) Solution II: Left to right followed by right to left. Two scans are sufficient. (Danielson`s distance map 1980) Solution III: Start from x0, update its neighboring points, accept updated values, and update their neighbors, etc. 2 1 2 3 1 2

10. Weighted Domains • Local weight ,Arclength • Goal: distance function characterized by: • By the chain rule: • The Eikonal equation is T(x) F(x) x x0 x1

11. 2D Rectangular Grids Isotropic inhomogeneous domains Weighted arclength: the weight is Goal: Compute the distance T(x,y) from p0 where

12. Upwind Approximation in 2D T ij T i+1,j T i,j-1 i+1,j T i,j-1 i+1,j T i-1,j ij i,j+1 i-1,j

13. 2D Approximation Initialization: • given initial value or Update: Fitting a tilted plane with gradient , and two values anchored at the relevant neighboring grid points. T 1 i+1,j T 2 i,j-1 ij i,j+1 i-1,j

14. Computational Complexity • T is systematically constructed from smaller to larger T values. • Update of a heap element is O(log N). • Thus, upper bound of the total is O(N log N). www.math.berkeley.edu/~sethian

15. Shortest Path on Flat Domains Why do graph search based algorithms (like Dijkstra's) fail?

16. Edge Integration Cohen-Kimmel, IJCV, 1997. Solve the 2D Eikonal equation given T(p)=0 Minimal geodesic w.r.t.

17. Shape from Shading Rouy-Tourin SIAM-NU 1992, Kimmel-Bruckstein CVIU 1994, Kimmel-Sethian JMIV 2001. Solve the 2D Eikonal equation where Minimal geodesic w.r.t.

18. Path Planning 3 DOF Solve the Eikonal Eq. in 3D {x,y,j}-CS given T(x0,y0,j0)=0, Minimal geodesic w.r.t.

19. Path Planning 3 DOF

20. Path Planning 4 DOF Solve the Eikonal Eq. in 4D Minimal geodesic w.r.t.

21. Update Acute Angle Given ABC, update C. Consistency and monotonicity: Update only `from within the triangle' hinABC Find t=EC that satisfies the gradient approximation (t-u)/h= F. c c

22. Update Procedure We end up with: t must satisfy u<t, and hinABC. The update procedure is IF (u<t) AND (a cos q < b(t-u)/t < a/cosq) THEN T(C) = min {T(C),t+T(A)}; ELSE T(C)=min {T(C),bF+T(A),aF+T(B)}. u C B A

23. Obtuse Problems This front first meets B, next A, and only then C. A is `supported’ by a single point. The supported section of incoming fronts is a limited section.

24. Solution by splitting Extend this section and link the vertex to one within the extended section. Recursive unfolding: Unfold until a new vertex Q is found. Initialization step!

25. Recursive Unfolding: Complexity e = length of longest edge The extended section maximal area is bounded bya<= e /(2a ). The minimal area of any unfolded triangle is bounded below a >= (h a ) q /2, The number of unfolded triangles before Q is found is bounded by m<= a /a = e /(q h a ). max 2 max min 2 min min min min 2 2 3 max max min min min min

26. 1st Order Accuracy The accuracy for acute triangles isO(e ) Accuracy for the obtuse case O(e /(p-q )) max max max

27. Minimal Geodesics

28. Minimal Geodesics

29. Linear Interpolation ODE ‘back tracking’

30. Quadratic Interpolation

31. Voronoi Diagrams and Offsets • Given n points, { p D, j 0,..,n-1} • Voronoi region: • G = {p D| d(p,p ) < d(p,p ), V j = i}. j i i j

32. Geodesic Voronoi Diagrams and Geodesic Offsets

33. Geodesic Voronoi Diagrams and Geodesic Offsets

34. Marching Triangles • The intersection set of two functions is linearly interpolated via `marching triangle'

35. Voronoi Diagrams and Offsets on Weighted Curved Domains

36. Voronoi Diagrams and Offsets on Weighted Curved Domains

37. A frame grabber built at the Technion by Y Grinberg Lego Mindstorms rotates the laser (E. Gordon) Cheap and Fast 3D Scanner • PC + video frame grabber. • Video camera. • Laser line pointer. Joint with G. Zigelman motivated by simple shape from structure light methods, like Bouguet-Perona 99, Klette et al. 98

38. Cheap and Fast 3D Scanner

39. Detection and Reconstruction

40. Examples of Decimation Decimation - 3% of vertices Sub-grid sampling

41. Results

42. Results

43. Texture Mapping • Environment mapping: Blinn, Newell (76). • Environment mapping: Greene, Bier and Sloan (86). • Free-form surfaces: Arad and Elber (97). • Polyhedral surfaces: Floater (96, 98), Levy and Mallet (98). • Multi-dimensional scaling: Schwartz, Shaw and Wolfson (89).

44. Difficulties • Need for user intervention. • Local and global distortions. • Restrictive boundary conditions. • High computational complexity.

45. Flattening via MDS • Compute geodesic distances between pairs of points. • Construct a square distance matrix of geodesic distances^2. • Find the coordinates in the plane via multi-dimensional scaling. The simplest is `classical scaling’. • Use the flattened coordinates for texturing the surface, while preserving the texture features. Zigelman, Kimmel, Kiryati, IEEE T. on Visualization and Computer Graphics (in press). 0 0 0

46. Flattening

47. Flattening

48. Distances - comparison

49. Texture Mapping

50. Texture Mapping