Prinston___GeometricPrimitives-2x2

Algorithms in Java, 4th Edition · Robert Sedgewick and Kevin Wayne · Copyright © 2009 · January 26, 2010 8:51:53 AM 6.1 Geometric Primitives ‣ primitive operations ‣ convex hull ‣ closest pair ‣ voronoi diagram 2 Geometric algorithms Applications. • Data mining. • VLSI design. • Computer vision. • Mathematical models. • Astronomical simulation. • Geographic information systems. • Computer graphics (movies, games, virtual reality). • Models of physical world (maps, architecture, medical imaging). History. • Ancient mathematical foundations. • Most geometric algorithms less than 25 years old. http://www.ics.uci.edu/~eppstein/geom.html airflow around an aircraft wing 3 ‣ primitive operations ‣ convex hull ‣ closest pair ‣ voronoi diagram 4 Geometric primitives Point: two numbers (x, y). Line: two numbers a and b. [ax + by = 1] Line segment: two points. Polygon: sequence of points. Primitive operations. • Is a polygon simple? • Is a point inside a polygon? • Do two line segments intersect? • What is Euclidean distance between two points? • Given three points p1, p2, p3, is p1!p2!p3 a counterclockwise turn? Other geometric shapes. • Triangle, rectangle, circle, sphere, cone, … • 3D and higher dimensions sometimes more complicated. any line not through origin 5 Geometric intuition Warning: intuition may be misleading. • Humans have spatial intuition in 2D and 3D. • Computers do not. • Neither has good intuition in higher dimensions! Q. Is a given polygon simple? we think of this algorithm sees this no crossings x y 1 6 5 8 7 2 7 8 6 4 2 1 x y 1 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 2 18 4 18 4 19 4 19 4 20 3 20 2 20 x y 1 10 3 7 2 8 8 3 4 6 5 15 1 11 3 14 2 16 Jordan curve theorem. [Jordan 1887, Veblen 1905] Any continuous simple closed curve cuts the plane in exactly two pieces: the inside and the outside. Q. Is a point inside a simple polygon? Application. Draw a filled polygon on the screen. 6 Polygon inside, outside Puzzle. Are A and B inside or outside the maze? 7 Fishy maze http://britton.disted.camosun.bc.ca/fishmaze.pdf 8 Polygon inside, outside Jordan curve theorem. [Jordan 1887, Veblen 1905] Any continuous simple closed curve cuts the plane in exactly two pieces: the inside and the outside. Q. Is a point inside a simple polygon? Application. Draw a filled polygon on the screen. http://www.ics.uci.edu/~eppstein/geom.html Q. Does line segment intersect ray? 9 public boolean contains(double x0, double y0) { int crossings = 0; for (int i = 0; i < N; i++) { double slope = (y[i+1] - y[i]) / (x[i+1] - x[i]); boolean cond1 = (x[i] <= x0) && (x0 < x[i+1]); boolean cond2 = (x[i+1] <= x0) && (x0 < x[i]); boolean above = (y0 < slope * (x0 - x[i]) + y[i]); if ((cond1 || cond2) && above) crossings++; } return crossings % 2 != 0; } Polygon inside, outside: crossing number y0 = yi+1 - yi xi+1 - xi (x0 - xi) + yi xi ! x0 ! xi+1 (xi, yi) (xi+1, yi+1) (x0, y0) 10 CCW. Given three point a, b, and c, is a-b-c a counterclockwise turn? • Analog of compares in sorting. • Idea: compare slopes. Lesson. Geometric primitives are tricky to implement. • Dealing with degenerate cases. • Coping with floating-point precision. Implementing ccw a b yes a c no a b Yes (!-slope) a b ??? (collinear) b a ??? (collinear) a c ??? (collinear) cc b c c b CCW. Given three point a, b, and c, is a!b!c a counterclockwise turn? • Determinant gives twice signed area of triangle. • If area > 0 then a!b!c is counterclockwise. • If area < 0, then a!b!c is clockwise. • If area = 0, then a!b!c are collinear. < 0> 0 11 Implementing ccw ! 2 " Area(a, b, c) = ax ay 1 bx by 1 cx cy 1 = (bx # ax )(cy # ay ) # (by # ay )(cx # ax ) (ax, ay) (bx, by) (cx, cy) (ax, ay) (bx, by) (cx, cy) 12 Immutable point data type public class Point { private final int x; private final int y; public Point(int x, int y) { this.x = x; this.y = y; } public double distanceTo(Point that) { double dx = this.x - that.x; double dy = this.y - that.y; return Math.sqrt(dx*dx + dy*dy); } public static int ccw(Point a, Point b, Point c) { int area2 = (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x); if (area2 < 0) return -1; else if (area2 > 0) return +1; else return 0; } public static boolean collinear(Point a, Point b, Point c) { return ccw(a, b, c) == 0; } } cast to long to avoid overflowing an int l1.p1l2.p1 13 Intersect. Given two line segments, do they intersect? • Idea 1: find intersection point using algebra and check. • Idea 2: check if the endpoints of one line segment are on different "sides" of the other line segment (4 calls to ccw). Sample ccw client: line intersection not handled l1.p2 l2.p2 public static boolean intersect(LineSegment l1, LineSegment l2) { int test1 = Point.ccw(l1.p1, l1.p2, l2.p1) * Point.ccw(l1.p1, l1.p2, l2.p2); int test2 = Point.ccw(l2.p1, l2.p2, l1.p1) * Point.ccw(l2.p1, l2.p2, l1.p2); return (test1 <= 0) && (test2 <= 0); } 14 ‣ primitive operations ‣ convex hull ‣ closest pair ‣ voronoi diagram 15 Convex hull A set of points is convex if for any two points p and q in the set, the line segment pq is completely in the set. Convex hull. Smallest convex set containing all the points. Properties. • "Simplest" shape that approximates set of points. • Shortest perimeter fence surrounding the points. • Smallest area convex polygon enclosing the points. convex not convex convex hull p q p q 16 Mechanical solution Mechanical convex hull algorithm. Hammer nails perpendicular to plane; stretch elastic rubber band around points. http://www.dfanning.com/math_tips/convexhull_1.gif 17 An application: farthest pair Farthest pair problem. Given N points in the plane, find a pair of points with the largest Euclidean distance between them. Fact. Farthest pair of points are on convex hull. 18 Brute-force algorithm Observation 1. Edges of convex hull of P connect pairs of points in P. Observation 2. p-q is on convex hull if all other points are counterclockwise of pq. O(N3) algorithm. For all pairs of points p and q: • Compute ccw(p, q, x) for all other points x. • p-q is on hull if all values are positive. p q 19 Package wrap (Jarvis march) Package wrap. • Start with point with smallest (or largest) y-coordinate. • Rotate sweep line around current point in ccw direction. • First point hit is on the hull. • Repeat. 20 Package wrap (Jarvis march) Implementation. • Compute angle between current point and all remaining points. • Pick smallest angle larger than current angle. • "(N) per iteration. 21 Jarvis march: demo http://www.cs.princeton.edu/courses/archive/fall08/cos226/demo/ah/JarvisMarch.html 22 Jarvis march: demo http://www.cs.princeton.edu/courses/archive/fall08/cos226/demo/ah/JarvisMarch.html 23 Jarvis march: demo http://www.cs.princeton.edu/courses/archive/fall08/cos226/demo/ah/JarvisMarch.html 24 How many points on the hull? Parameters. • N = number of points. • h = number of points on the hull. Package wrap running time. "(N h). How many points on hull? • Worst case: h = N. • Average case: difficult problems in stochastic geometry. - uniformly at random in a disc: h = N1/3 - uniformly at random in a convex polygon with O(1) edges: h = log N 25 Graham scan Graham scan. • Choose point p with smallest (or largest) y-coordinate. • Sort points by polar angle with p to get simple polygon. • Consider points in order, and discard those that would create a clockwise turn. p 26 Graham scan: demo http://www.cs.princeton.edu/courses/archive/fall08/cos226/demo/ah/GrahamScan.html 27 Graham scan: demo http://www.cs.princeton.edu/courses/archive/fall08/cos226/demo/ah/GrahamScan.html 28 Graham scan: implementation Implementation. • Input: p[1], p[2], …, p[N] are distinct points. • Output: M and rearrangement so that p[1], p[2], …, p[M] is convex hull. Running time. O(N log N) for sort and O(N) for rest. // preprocess so that p[1] has smallest y-coordinate // sort by polar angle with respect to p[1] p[0] = p[N]; // sentinel int M = 2; for (int i = 3; i <= N; i++) { while (Point.ccw(p[M-1], p[M], p[i]) <= 0) M--; M++; swap(p, M, i); } why? discard points that would create clockwise turnadd i to putative hull 29 Quick elimination Quick elimination. • Choose a quadrilateral Q or rectangle R with 4 points as corners. • Any point inside cannot be on hull. - 4 ccw tests for quadrilateral - 4 compares for rectangle Three-phase algorithm. • Pass through all points to compute R. • Eliminate points inside R. • Find convex hull of remaining points. In practice. Eliminates almost all points in linear time. Q these points eliminated R Asymptotic cost to find h-point hull in N-point set. 30 Convex hull algorithms costs summary t assumes "reasonable" point distribution output sensitive algorithm running time package wrap N h Graham scan N log N quickhull N log N mergehull N log N sweep line N log N quick elimination N t marriage-before-conquest N log h output sensitive 31 Convex hull: lower bound Models of computation. • Compare-based: compare coordinates. (impossible to compute convex hull in this model of computation) • Quadratic decision tree model: compute any quadratic function of the coordinates and compare against 0. Proposition. [Andy Yao, 1981] In quadratic decision tree model, any convex hull algorithm requires #(N log N) ops. higher constant-degree polynomial tests don't help either [Ben-Or, 1983] even if hull points are not required to be output in counterclockwise order (a.x < b.x) || ((a.x == b.x) && (a.y < b.y))) (a.x*b.y - a.y*b.x + a.y*c.x - a.x*c.y + b.x*c.y - c.x*b.y) < 0 32 ‣ primitive operations ‣ convex hull ‣ closest pair ‣ voronoi diagram 33 Closest pair Closest pair problem. Given N points in the plane, find a pair of points with the smallest Euclidean distance between them. Fundamental geometric primitive. • Graphics, computer vision, geographic information systems, molecular modeling, air traffic control. • Special case of nearest neighbor, Euclidean MST, Voronoi. fast closest pair inspired fast algorithms for these problems 34 Closest pair Closest pair problem. Given N points in the plane, find a pair of points with the smallest Euclidean distance between them. Brute force. Check all pairs with N2 distance calculations. 1-D version. Easy N log N algorithm if points are on a line. Degeneracies complicate solutions. [assumption for lecture: no two points have same x-coordinate] • Divide: draw vertical line L so that ~ !N points on each side. 35 Divide-and-conquer algorithm L 36 Divide-and-conquer algorithm • Divide: draw vertical line L so that ~ !N points on each side. • Conquer: find closest pair in each side recursively. L 12 21 37 Divide-and-conquer algorithm • Divide: draw vertical line L so that ~ !N points on each side. • Conquer: find closest pair in each side recursively. • Combine: find closest pair with one point in each side. • Return best of 3 solutions. seems like !(N2) L 12 21 8 Find closest pair with one point in each side, assuming that distance < $. L 12 21 38 How to find closest pair with one point in each side? $ = min(12, 21) Find closest pair with one point in each side, assuming that distance < $. • Observation: only need to consider points within $ of line L. 39 How to find closest pair with one point in each side? L 12 21 $ = min(12, 21) $ Find closest pair with one point in each side, assuming that distance < $. • Observation: only need to consider points within $ of line L. • Sort points in 2$-strip by their y coordinate. L 12 21 $ = min(12, 21) $ 40 How to find closest pair with one point in each side? 1 2 3 4 5 6 7 41 How to find closest pair with one point in each side? Find closest pair with one point in each side, assuming that distance < $. • Observation: only need to consider points within $ of line L. • Sort points in 2$-strip by their y coordinate. • Only check distances of those within 11 positions in sorted list! L 12 21 $ = min(12, 21) $ 1 2 3 4 5 6 7 why 11? 42 How to find closest pair with one point in each side? Def. Let si be the point in the 2$-strip, with the ith smallest y-coordinate. Claim. If |i – j| % 12, then the distance between si and sj is at least $. Pf. • No two points lie in same !$-by-!$ box. • Two points at least 2 rows apart have distance % 2(!$). ▪ Fact. Claim remains true if we replace 12 with 7. $ 27 29 30 31 28 26 25 $ !$ 2 rows !$ !$ 39 i j 43 Divide-and-conquer algorithm O(N log N) 2T(N / 2) O(N) O(N log N) O(N) Closest-Pair(p1, …, pn) { Compute separation line L such that half the points are on one side and half on the other side. $1 = Closest-Pair(left half) $2 = Closest-Pair(right half) $ = min($1, $2) Delete all points further than $ from separation line L Sort remaining points by y-coordinate. Scan points in y-order and compare distance between each point and next 11 neighbors. If any of these distances is less than $, update $. return $. } 44 Divide-and-conquer algorithm: analysis Running time recurrence. T(N) ! 2T(N/2) + O(N log N). Solution. T(N) = O(N (log N)2). Remark. Can be improved to O(N log N). Lower bound. In quadratic decision tree model, any algorithm for closest pair requires #(N log N) steps. sort by x- and y-coordinates once (reuse later to avoid re-sorting) (x1 - x2) 2 + (y1 - y2) 2 45 ‣ primitive operations ‣ convex hull ‣ closest pair ‣ voronoi diagram Life-or-death question. Given a new cholera patient p, which water pump is closest to p’s home? 46 1854 cholera outbreak, Golden Square, London http://content.answers.com/main/content/wp/en/c/c7/Snow-cholera-map.jpg 47 Voronoi diagram Voronoi region. Set of all points closest to a given point. Voronoi diagram. Planar subdivision delineating Voronoi regions. Fact. Voronoi edges are perpendicular bisector segments. Voronoi of 2 points (perpendicular bisector) Voronoi of 3 points (passes through circumcenter) 48 Voronoi diagram Voronoi region. Set of all points closest to a given point. Voronoi diagram. Planar subdivision delineating Voronoi regions. 49 Voronoi diagram: more applications Anthropology. Identify influence of clans and chiefdoms on geographic regions. Astronomy. Identify clusters of stars and clusters of galaxies. Biology, Ecology, Forestry. Model and analyze plant competition. Cartography. Piece together satellite photographs into large "mosaic" maps. Crystallography. Study Wigner-Setiz regions of metallic sodium. Data visualization. Nearest neighbor interpolation of 2D data. Finite elements. Generating finite element meshes which avoid small angles. Fluid dynamics. Vortex methods for inviscid incompressible 2D fluid flow. Geology. Estimation of ore reserves in a deposit using info from bore holes. Geo-scientific modeling. Reconstruct 3D geometric figures from points. Marketing. Model market of US metro area at individual retail store level. Metallurgy. Modeling "grain growth" in metal films. Physiology. Analysis of capillary distribution in cross-sections of muscle tissue. Robotics. Path planning for robot to minimize risk of collision. Typography. Character recognition, beveled and carved lettering. Zoology. Model and analyze the territories of animals. http://voronoi.com http://www.ics.uci.edu/~eppstein/geom.html 50 Scientific rediscoveries Reference: Kenneth E. Hoff III year discoverer discipline name 1644 Descartes astronomy "Heavens" 1850 Dirichlet math Dirichlet tesselation 1908 Voronoi math Voronoi diagram 1909 Boldyrev geology area of influence polygons 1911 Thiessen meteorology Thiessen polygons 1927 Niggli crystallography domains of action 1933 Wigner-Seitz physics Wigner-Seitz regions 1958 Frank-Casper physics atom domains 1965 Brown ecology area of potentially available 1966 Mead ecology plant polygons 1985 Hoofd et al. anatomy capillary domains 51 Fortune's algorithm Industrial-strength Voronoi implementation. • Sweep-line algorithm. • O(N log N) time. • Properly handles degeneracies. • Properly handles floating-point computations. Try it yourself! http://www.diku.dk/hjemmesider/studerende/duff/Fortune/ Remark. Beyond scope of this course. algorithm preprocess query brute 1 N Fortune N log N log N 52 Fortune's algorithm in practice Def. Triangulation of N points such that no point is inside circumcircle of any other triangle. 53 Delaunay triangulation circumcircle of 3 points Proposition 1. It exists and is unique (assuming no degeneracy). Proposition 2. Dual of Voronoi (connect adjacent points in Voronoi diagram). Proposition 3. No edges cross & O(N) edges. Proposition 4. Maximizes the minimum angle for all triangular elements. Proposition 5. Boundary of Delaunay triangulation is convex hull. Proposition 6. Shortest Delaunay edge connects closest pair of points. 54 Delaunay triangulation properties Delaunay Voronoi 55 Delaunay triangulation application: Euclidean MST Euclidean MST. Given N points in the plane, find MST connecting them. [distances between point pairs are Euclidean distances] Brute force. Compute N2 / 2 distances and run Prim's algorithm. Ingenuity. • MST is subgraph of Delaunay triangulation. • Delaunay has O(N) edges. • Compute Delaunay, then use Prim (or Kruskal) to get MST in O(N log N) ! Ingenious algorithms enable solution of large instances for numerous fundamental geometric problems. Note. 3D and higher dimensions test limits of our ingenuity. 56 asymptotic time to solve a 2D problem with N points Geometric algorithms summary problem brute clever convex hull N2 N log N farthest pair N2 N log N closest pair N2 N log N Delaunay/Voronoi N4 N log N Euclidean MST N2 N log N