CS 101.4 Algorithms in Geometry and Topology

1. CS 101.4Algorithms in Geometry and Topology http://www.ist.caltech.edu/~sthite/cs101/ https://courses.caltech.edu/course/view.php?id=162 Shripad Thite shripad@caltech.edu

2. About This Course Lecture plus seminar style course on algorithms in geometry processing “Computational geometry” and “Computational topology” Quick and comprehensive introduction to topics, followed by presentation and discussion of open research problems My goals: to excite you about research in the topics in this course, to give you enough background to follow-up on topics by reading papers, to encourage collaborative thinking about problems of active research interest, to learn from you about applications in your research and your other courses Emphasis on collaboration. Actively seek other sources of knowledge and inspiration. Follow the Caltech Honor Code, especially in written work. Identify and share due credit with collaborators and colleagues.

3. Topics We will study different stages of the geometry processing pipeline Part 1. Mesh generation (2 weeks) Voronoi diagrams, Delaunay triangulations in 2D and 3D Part 2. Surface reconstruction (3 weeks) Building a triangulation of a surface in R3 from point samples (laser scans) Part 3. Smoothing and mesh simplification (2 weeks) Removing unwanted noise and unnecessary detail for efficiency Part 4. Topological analysis and simplification (3 weeks) Discovering tunnels and loops, removing unwanted topological detail

4. Prerequisites No previous background in computational geometry or topology is assumed This is an algorithms course. So we will constantly talk about the complexity of algorithms, which means the running time and memory usage “Running time” usually means the number of elementary operations (comparisons, arithmetic operations, evaluating a simple function, etc.) performed by an algorithm when presented with an input instance of size n You should be aware of O() (“Big-Oh”) notation which is used to denote the asymptotic complexity of an algorithm Thus, O(n) is faster than O(n log n) is faster than O(n log2n) is faster than O(n2) is faster than O(n3) is faster than O(2n) …

5. Organization We meet TueThu 2:30-3:55pm in 287 Jorgensen TA: Keenan Crane (keenan@cs.caltech.edu) I encourage you to meet with me and/or Keenan often My office is Moore 335, but I am more likely to be at the Red Door café My likely office hours: TueThu after class 4-6pm Feel free to call me on my cellphone: (626) 487-5786

6. Homework and Grades Units: 3-0-6 Weekly homework: reading 2-3 research papers and writing a 1-page critical review, typically assigned Tuesday, due next Tuesday before class. Occasional exercises for extra credit, due end of term In-class presentation: give one lecture in class on the current topic, presenting research paper(s) and original ideas, with help from me/Keenan Term project: focusing on a theoretical research problem, 10-20 page report surveying previous work and detailing a solution or attempts at a solution, including conjectures and promising directions for future research. Project chosen and conducted in consultation with me, preferably in a team. Each of the above is worth 1/3rd of final grade {A,B,C,F}; A or B to pass Choose project topic and volunteer for presentations early!

7. Questions?

8. Voronoi Diagram a.k.a. Thiessen polygons, Dirichlettesselation, etc. Partition of space into cells, points in a cell are closer to one site than others The canonical geometric concept and data structure used in applications!

9. Voronoi Diagram: Applications Facility location: place new facility to maximize the area of its Voronoi cell to capture largest possible market share Clustering: partition space into groups or categories, each group is the set of all items that are closer to one representative item than the others Nearest-neighbor query data structure: given a query point q, find its nearest site; build a point location data structure on top of the Voronoi diagram … and many more!

10. General Position We frequently assume that the input is in general position or non-degenerate E.g., No three points are collinear General position assumption is different for each problem, and must always be clearly stated Purpose: Degenerate configurations are combinatorially complicated, e.g., a Voronoi vertex may have arbitrarily large degree. Algorithms that correctly handle degenerate inputs are much more complicated and a lot less efficient. It is simpler to explain an algorithm working on non-degenerate input Degenerate configurations are “unlikely” (measure zero) Standard perturbation techniques exist to handle degenerate inputs

11. Voronoi Diagram General position: No three points collinear and no four points co-circular Every Voronoi vertex is equidistant from exactly three sites

12. Delaunay Triangulation Empty circumcircle property: Triangle pqr if and only if circumscribing disk of pqr is empty of other points Delaunay triangulation of P partitions the convex hull of P http://en.wikipedia.org/wiki/Delaunay_triangulation

13. Duality Delaunay triangulation of P is the dual of the Voronoi diagram of P Edge pqiffVor(p),Vor(q) intersect Triangle pqriffVor(p),Vor(q),Vor(r) intersect Theorem:Δpqr has empty circumcircle iffVor(p),Vor(q),Vor(r) intersect

14. ComputingVoronoi diagrams andDelaunay triangulations

15. Geometric Predicates: Orientation Test Given three points (x1,y1), (x2,y2), (x3,y3), are they in counterclockwise orientation? > 0, then counterclockwise = 0, then collinear < 0, then clockwise (x3,y3) (x2,y2) The above determinant is 2 x signed area of the triangle (It is the magnitude of the cross-product of two edges) (x1,y1)

16. Geometric Predicates: In-circle Test Given four points (x1,y1), (x2,y2), (x3,y3), (x4,y4), is the last point outside the disc through the first three? > 0, then outside = 0, then co-circular < 0, then inside The above determinant is 6 x signed area of the tetrahedron with vertices (x1,y1,z1=x12+y12), (x2,y2,z2=x22+y22), (x3,y3,z3=x32+y32), (x4,y4,z4=x42+y42)

17. Geometric Predicates: Robustness Computation depends on computing sign of determinants Roundoff error may lead to an incorrect result when the true determinant is near zero; error might be unpredictable and inconsistent Exact arithmetic using arbitrary precision floating-point arithmetic is slow Solution: Use adaptive precision—compue a sequence of successively more accurate approximations to the determinant, stopping as soon as the accuracy of the sign is guaranteed Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates. Jonathan Richard Shewchuk. Discrete & Computational Geometry 18:305-363, 1997 http://www.cs.cmu.edu/~quake/robust.html

18. Convex Hull The convex hull of a set P of points, CH(P), is the intersection of all convex sets that contain P Let P be a set of n points in R3 General position: No four points are coplanar The boundary of CH(P) is a triangulation Triangles = “facets” “1-skeleton” = Graph of vertices P and edges, a planar graph The lower hull or lower convex hull of P is the convex hull of P together with the point (0,∞)

19. Delaunay Triangulation = Convex Hull Delaunay triangulation of P is (the projection of) the lower convex hull of P lifted to the standard paraboloid z=x2+y2 http://cse.iitkgp.ac.in/geomalgos/slides/SCNandy.ppt Voronoi diagram of P is (the projection of) the upper envelope of planes tangent to the paraboloid at the lifted versions of P

20. 2D Convex Hull: Melkman’s algorithm 1. Sort the n points from left to right, in O(n log n) time 2. Interpolate a simple polygon through the points by doubling each edge between consecutive points 3. Repeatedly convexify the polygon at each vertex pi by testing whether p0pipi+1 is counterclockwise Data structure: Push and pop vertices from a stack Analysis: Since each vertex is pushed and popped exactly once, convexifying a simple polygon takes O(n) time and O(n) space Convex hulls in O(n log h) time and O(n) space: Optimal output-sensitive convex hull algorithms in two and three dimensions. Timothy Chan. Discrete & Computational Geometry, 16:361-368, 1996. http://www.cs.uwaterloo.ca/~tmchan/pub.html

21. 2D Convex Hull: Melkman’s algorithm Convexifying a simple polygon : O(n) time

22. 3D Convex Hull: Gift Wrapping Imagine wrapping paper tightly around a 3D point set P Find an edge ab of the 3D convex hull from an edge of the convex hull of the 2D projection of P : O(n log n) time Find the vertex c such that all other points are on the same side of Δabc; Sort all points by angle made with edge ab and choose c making the smallest angle : O(n log n) time Starting from Δabc, repeatedly find the next adjacent triangle by looking for its third vertex Running time = O(nh) where h is the number of triangles on CH(P) h = “output complexity”

23. “Locally Delaunay” An edge e of a triangulation is “locally Delaunay” if: e Theorem: A triangulation is Delaunay if and only if every edge is locally Delaunay

24. Edge Flip Replace one diagonal of a convex quadrilateral by the other diagonal e’ e Theorem: The diagonal e is locally Delaunay if and only if the diagonal e’ is not locally Delaunay.

25. Randomized Incremental Algorithm for 2D Delaunay Triangulation(and for 3D Convex Hull)

26. Randomized Incremental Algorithm Randomly permute the n input points to get {p1,p2,p3,…,pi-1,pi,…,pn} Insert each point in random order into previous triangulation / convex hull Repeatedly flip edges to restore local Delaunay property everywhere Equivalently, replace facets of old convex hull by those of new hull

27. Randomized Incremental Algorithm Start with the convex hull of the first 4 points = a tetrahedron At the beginning of the ith step, where i=5,6,…,n, we have computed the convex hull CHi-1 of the first i-1 points • In the ith step, we insert point pi and update the convex hull to get CHi: • If pi in the interior of CHi-1 then nothing to be done, set Chi:=CHi-1 • Otherwise, delete all facets of CHi-1 that are visible from pi and add new facets that join pi to each edge on the horizon of the facets visible from pi to get CHi

28. Data Structures Store the convex hull boundary triangulation in a half-edge data structure, which also allows us to traverse the dual graph For each point not yet inserted, maintain a pointer to a visible facet of the current convex hull (the facet directly above it, for 2D Delaunay triangulation) For each facet of the current convex hull, maintain a list of uninserted points that point to it

29. Analysis: Update Cost Cost of deleting facets visible from point pi is proportional to their number Start with the representative visible facet of pi and walk in the dual graph to identify all visible facets Each facet is deleted at most once So, we can count the number of facets created and double it to measure the cost of both insertion and deletion Number of facets created by inserting pi = deg(pi, CHi)

30. Backwards Analysis Number of facets created by inserting pi = deg(pi, CHi) Each of the first i points {p1,p2,…,pi-1,pi} is equally probable to be the ith (last) vertex inserted to obtain CHi CHi is independent of the order of insertion of its points Therefore, the expected degree of pi in CHi is: We used the fact that the 1-skeleton of CHi is a planar graph with ≤ i vertices, so it has at most 3i-6 edges

31. Analysis: Update Cost Cost of deleting facets visible from point pi is proportional to their number Each facet is deleted at most once So, we can count the number of facets created and double it to measure the cost of both insertion and deletion Expected number of facets created by inserting pi = O(1) Therefore, expected total update cost of n insertions = O(n)

32. Analysis: Visibility Maintenance Cost We need to bound the number of times that the visible facet of point pi changes in steps 1 through i-1 Suppose the visible facet changes in the jth step, where 1 ≤ j ≤ i-1 Let Δj be the new visible facet of pi Backwards analysis: Triangle Δj must have been “completed” in step j, i.e., the last of the 3 vertices of Δj must have been inserted in step j Each of the 3 vertices of Δj is equally likely to be the last of its vertices inserted Pr[Δj created in step j] = 1/j+1/j+1/j = 3/j

33. Analysis: Visibility Maintenance Cost We need to bound the number of times that the visible facet of point pi changes in steps 1 through i-1 Pr[visible facet of pi changes in jth step] = 3/j Expected update cost for pi: Therefore, expected total update cost = O(n log n)

34. Randomized Incremental Algorithm Randomly permute the n input points to get {p1,p2,p3,…,pi-1,pi,…,pn} Start with the convex hull of the first 4 points = a tetrahedron At the beginning of the ith step, where i=5,6,…,n, we have computed the convex hull CHi-1 of the first i-1 points • In the ith step, we insert point pi and update the convex hull to get CHi: • If pi in the interior of CHi-1 then nothing to be done, set Chi:=CHi-1 • Otherwise, delete all facets of CHi-1 that are visible from pi and add new facets that join pi to each edge on the horizon of the facets visible from pi to get CHi Expected running time = O(n) + O(n log n) = O(n log n)

35. Mesh Generation

36. Mesh A mesh is a partition of a domain into simple elements, such as triangles, tetrahedra etc. A Delaunay triangulation of points P is a mesh of the convex hull of P A constrained or conforming Delaunay triangulation of a domain D with piecewise linear boundary is a mesh of D that includes or contains the boundary of D in its facets http://www.cse.ohio-state.edu/~tamaldey/papers.html

37. Mesh Quality The quality of a triangle or tetrahedron is measured by its ratio of circumradius to shortest edge length The quality of a mesh is that of its poorest quality element Often conflicting goals: To minimize gradient interpolation error, large angles are bad but small angles are okay To minimize condition number of stiffness matrix (minimize largest eigenvalue), small angles are bad but large angles are okay What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measures. Jonathan Shewchuk. Unpublished preprint, 2002. http://www.cs.cmu.edu/~jrs/jrspapers.html#quality

38. Delaunay Mesh in 2D is “Best Possible” A Delaunay triangulation maximizes the smallest angle For every point set P, among all triangulations of CH(P), the Delaunay triangulation of P has the maximum smallest angle. A Delaunay triangulation minimizes the largest circumradius A Delaunay triangulation minimizes the largest smallest-enclosing-circle radius But … A Delaunay triangulation does not minimize the largest angle

39. Bad Quality Tetrahedra Wedge Spade Spire Spear Spindle Spike Splinter Cap Sliver Most can be eliminated by Delaunay refinement

40. Slivers Take four evenly spaced points along the equator of a sphere Perturb one point slightly off the equator but still on the sphere Slivers are not removed by Delaunay refinement algorithms because they have good (=small) circumradius to shortest edge length ratio Sliver exudation removes slivers by building a weighted Delaunay triangulation

41. Weighted Delaunay Triangulation Point p=(x,y) with real weight W lifts to (x,y,x2+y2-W) Weighted Delaunay triangulation is the lower hull of lifted points, now displaced vertically from the standard paraboloid A sliver can be destroyed by “pumping” the weight of one of its vertices

42. Worst-case Complexity of Delaunay Triangulations Delaunay triangulation of n points in 2D has O(n) size, because its 1-skeleton is a planar graph A set of n points on each of two skew lines in 3D: {p1,p2,p3,…,pi,…,pn} {q1,q2,q3,…,qi,…,qn} The circumsphere of every tetrahedron {pi,pi+1,qj,qj+1} is empty of other points O(n2) tetrahedra Delaunay triangulation of n points in d dimensions can have simplices, e.g., npoints on the moment curve (t,t2,t3,…,td)

43. Summary • A Voronoi diagram of n sites partitions space into cells based on proximity to the sites • A Delaunay triangulation of n points partitions their convex hull into triangles • Voronoi diagram and Delaunay triangulation are dual to each other • Delaunay triangulation is computed by computing the convex hull of the points lifted up to the standard paraboloid

44. Summary • A Delaunay triangulation is a good quality mesh • Slivers present a challenge in 3D • Delaunay refinement can improve the quality • Sliver exudation, using weighted Delaunay triangulation, can remove slivers

45. A Preview of What’s to Come

46. Surface Reconstruction Given a sampling of points from an underlying smooth surface that is “sufficiently dense”, construct a triangulated surface that is guaranteed to be geometrically close to the original surface and topologically accurate Delaunay-based reconstruction: build a 3D Delaunay triangulation of the points and extract a subset of its facets as the reconstruction Moving Least Squares reconstruction: we are also given surface normals at the sample points. Build an implicit function I whose zero level set is geometrically close and topologically correct. Think of this as a generalization of linear regression or line fitting. Given certain ways to query I, build a triangulation that approximates its level set.

47. Smoothing and Simplification Laplaciansmoothing: Remove high-frequency noise from a (reconstructed) surface triangulation Reduce the number of triangles for efficient processing, trading-off accuracy and efficiency depending on the application (Think of lossy image compression) “Local” method: Quadric-based pair contraction algorithm due to Garland and Heckbert. Minimize the number of triangles subject to acceptable error. “Global” method: Variational shape approximation, due to Cohen-Steiner, Alliez, and Desbrun. Minimize a global energy functional representing the geometric error while keeping within a budget on the number of triangles.

48. Topology Introduction to homotopy and homology Computing optimal homotopy and homology generators = systems of loops of minimum total length Finding handles and tunnels Introduction to topological persistence Identifying topological features that are born and die in the process of varying a parameter (“time”). Assigning importance to such features. Removing topological noise = least important features.

49. Homework Due before class on Thursday, January 15 Read the following paper and write a critical analysis: Anisotropic Voronoi Diagrams and Guaranteed-Quality Anisotropic Mesh Generation. François Labelle and Jonathan Richard Shewchuk. Proceedings of the ACM Symposium on Computational Geometry (SoCG), 2003. http://www.cs.berkeley.edu/~flab/ • Your 1-page report should contain four parts: • Problem definition and summary of results: what is the problem being solved? Summarize the key concepts and results being proved • Pros: what are the novel contributions, key insights, and connections made? • Cons: what are the drawbacks of the solution strategy? what ideas have been overlooked? what new approaches should have been considered? • Open problem: identify one unsolved problem motivated by this paper and describe an approach to solve it, discussing why the attempt might and might not work

50. Extra Credit Give an algorithm to decide if a given convex partition of the plane into n cells is a Voronoi diagram, in O(n) time. Prove that for each point p its nearest neighbor is a Delaunay neighbor. (The nearest-neighbor graph is a subgraph of the Delaunay graph/1-skeleton.) The sorted angle vector of a triangulation is a vector of all 3m angles of a triangulation with m triangles, such that the angles are arranged in non-decreasing order. Prove that the Delaunay triangulation of a point set P lexicographically maximizes the sorted angle vector.