15-211 Fundamental Data Structures and Algorithms

1. 15-211Fundamental Data Structures and Algorithms Klaus Sutner April 22, 2004 Computational Geometry

2. Announcements • HW7 the clock is ticking ... • Quiz 3: Today, after class • Final Exam on Tuesday May 4, 5:30 pm • review session April 29

3. The Basics

4. Computational Geometry Lots of applications. - computer graphics - CAD/CAM - motion planning Could be done in dimension D, but we will stick to D=2. D=3 is already much harder.

5. Airplane wing triangulation

6. The Basics What are the basic objects in plane geometry, how do we represent them as data structures? - point two floats - line two points - ray two points, constraints - line segment two points, constraints and triangles, rectangles, polygons, ... Never mind curves such as circles and ellipses.

7. Disclaimer There are two annoying issues one has to confront in CompGeo algorithms. - Degenerate Cases Points are collinear, lines parallel, … Dealing with degeneracy is a huge pain, not to mention boring. - Floating point numbers are not reals. Floating point arithmetic introduces errors and hence leads to accuracy problems.

8. High Precision Arithmetic If 32 or 64 bits are not enough, use 128, or 512, or 1024, or … Problems: We have to be able to analyze how many bits are required. Slows down computation significantly. See GNU gmp.

9. Interval Arithmetic Replace “real number” x by an interval xL x  xU Problems: Arithmetic becomes quite a bit more complicated. Slows down computation significantly.

10. Symbolic Computation Use (arbitrary precision) rationals and whatever algebraic operations that are necessary symbolically: compute with the formal expressions, not their numerical values. Sqrt[2] Problems: Manipulating these expressions is complicated. Slows down computation significantly.

11. Warm-Up Let's figure out how to check if a point P lies on a line L. P = (p1,p2) L = (A,B) two-points representation Need P = A + l(B - A) where l is a real parameter. P B P A

12. Linear Equations The vector equation P = A + l (B – A) = (1- l) A + l B is just short-hand for two linear equations: p1 = a1 + l (b1 – a1) p2 = a2 + l (b2 – a2) Will have a common solution if P lies on L.

13. Subroutines Solving systems of equations such as p1 = a1 + l (b1 – a1) p2 = a2 + l (b2 – a2) is best left to specialized software: send all linear equation subproblems to an equation solver - that will presumably handle all the tricky cases properly and will produce reasonable accuracy.

14. Rays and Line Segments What if L =[A,B) is the ray from A through B? No problem. Now we need P = A + l (B - A) where l is a non-negative real. Likewise, if L = [A,B] is the line segment from A to B we need 0 l 1. Both l and (1 – l) must be non-negative.

15. Intersections How do we check if two lines intersect? A + l (B - A) = C + k (D – C) Two linear equations, two unknowns l , k. Add constraints for intersection of rays, line segments. B C A D

16. General Placement Here is a slightly harder problem: Which side of L is the point P on? P P left on B P right A

17. left on B right A Left/Right Turns March from A to B, then perform a left/right turn.

18. Calculus: Cross Product The cross product of two (3-D) vectors is another vector perpendicular to the given ones Z = X  Y Length and direction of Z express the signed area of the parallelogram spanned by the vectors. antisymmetric: X  Y = - Y  X

19. Cross Product Simple expression for cross product: (x1,x2,x3)  (y1,y2,y3) = (-x3y2 + x2 y3, x3 y1 - x1y3, -x2y1 + x1y2 ) Compute w = (b1-a1,b2-a2,0)  (p1-a1,p2-a2) = (0,0,-a2 b1+a1 b2+a2 p1-b2 p1-a1 p2+b1 p2 ) P is to the right of ray A to B iff w[3] > 0.

20. Determinants Alternatively one can compute the determinant of the matrix a1 a2 1 b1 b2 1 p1 p2 1 This may be a better solution since presumably the determinant subroutine presumably deals better with numerical problems.

21. Some Applications

22. Membership How do we check if a point X belongs to some region R in the plane (triangle, rectangle, polygon, ... ). R X R X X X

23. Membership Clearly the difficulty of a membership query depends a lot on the complexity of the region. For some regions it is not even clear that there is an inside and an outside.

24. Simple Polgyons A polygon is simple if its lines do not intersect except at the vertices of the polygon (and then only two). By the Jordan curve theorem any simple polygon partions the plane into inside and outside.

25. Point Membership The proof of the JCT provides a membership test. Let P be a simple polygon and X a point. Draw a line through X and count the intersections of that line with the edges of P.

26. Odd One can show that X lies inside of P iff the number of intersections (on either side) is odd. Note that there are bothersome degenerate cases: X lies on an edge, an edge lies on the line through X. Maybe wiggle the line a little.

27. Algorithm Theorem: One can test in linear time whether a point lies inside a simple polygon. Linear here refers to the number of edges of the polygon. And, of course, we assume that these edges are given as a list of points, say, in counterclockwise order.

28. Convexity Here is one type of simple region. A region R in the plane is convex if for all A, B in R: the line segment [A,B] is a subset of R. In particular convex polygons: boundary straight line segments.

29. Membership in Convex Polygon Traverse the boundary of the convex polygon in counterclockwise order, check that the point in question is always to the left of the corresponding line. This is clearly linear in the number of boundary points.

30. Intersections of LSs Suppose we are given a collection of LSs and want to compute all the intersections between them. Note that the number of intersections may be quadratic in the number n of LSs. The brute force algorithms is (n2): consider all pairs of LSs and check each pair for intersection. Bad if there are only a few intersections.

31. Sweep-Line Suppose there are s intersections. We would like an algorithm with running time O( (n+s) ??? ) where ??? is hopefully small. The key idea is to use a sweep-line: moves from left to right. At each point, maintain the intersection of the sweep line and the given LSs.

32. Variant It helps to think about the following easier problem first: check if there is at least one intersection between the given LSs. Brute force is still quadratic. Note how the SL needs to consider only immediate neighbors.

33. Again The sweep line needs to consider only immediate neighbors. t s r

34. SL Operations The SL has to keep track of all the LSs that currently intersect it. So we need operations insert(s) delete(s) where s is a given LS. The “time” when the segments are inserted/deleted are given by the end points of the input LSs. (Static events, preprocess by sorting.)

35. Dealing With Intersections When a new LS is entered, we have to check if it interacts with any of the neighboring LSs. So we need additional operations above(s) below(s) which return the immediate neighbors above/below on the SL. And we need to stop the sweep at any intersections found that way. Dynamic events, no known ahead of time.

36. Data Structures What is the right data structure for all of this? Sweep line: dictionary Perhaps threaded to get O(1) above/below operations. Event structure: priority queue. Theorem: Can compute all s intersections of n line segments in O((n + s) log n) steps.

37. Testing Simplicity How do we check if a polygon is simple? Use the line segment intersection algorithm on the edges.

38. Polgyon Intersection How do we check if two simple polygon intersect? Note: two cases!

39. Polgyon Intersection The first case can be handled with LS intersection testing. For the second, pick any vertex in P and check wether it lies in (the interior of) Q. Total running time: O(n log n).

40. Convex Hulls

41. A Hull Operation Suppose P is a set of points in the plane. The convex hull of P is the least set of points Q such that: - P is a subset of Q, - Q is convex. Written CH(P). This pattern should look very familiar by now (reachability in graphs, equivalence relations, ...)

42. Convex Combinations Abstractly it is easy to describe CH(P): is the set of all points X =  i Pi where 0 i and  i = 1. X is a convex combination of the Pi. So the convex combinations of A and B are the line segment [A,B]. And the convex combinations of three points (in general position) form a triangle. And so on.

43. So? Geometrically this is nice, but computationally this characterization of the convex hull is not too useful. We want an algorithm that takes as input a simple polygon P and returns as output the simple polygon CH(P). CH(P) P

44. Extremal Points Point X in region R is extremal iff X is not a convex combination of other points in R. For a polygon, the extremal points make up the CH.

45. An Observation So let's deal with a set of points P rather than regions. Lemma: Point X in P is not extremal iff X lies in the interior of a triangle spanned by three other points in P.

46. An Algorithm Hence we can find the convex hull of a simple polygon: for all n points we eliminate non-extremal points by trying all possible triangles using the membership test for convex polygons. In the end we sort the remaining extremal points in counterclockwise order. Unfortunately, this is O(n4). Could it be faster? When do we get n4?

47. Next Time Clearly, O(n4) is way too slow to be of any practical use. Next time we will see a number of simple algorithms that push this down to O(n log n). This turns out to be optimal: sorting can be reduced to convex hull. Note, though, that in dimension 3 and higher things get much more messy.