A Tutorial on Computational Geometry

A Tutorial on Computational Geometry Pham Minh Tri Ph.D. Candidate and Project Officer School of Computer Engineering 1 Mar 2008 presented by

What will be covered? Mainly popular static problems: Point inside polygon Convex hull Line segment intersection Delaunay triangulation Closest pair of points

Point Inside Polygon Test Given a point, determine if it lies inside a polygon or not

Ray Test Fire ray from point Count intersections Odd = inside polygon Even = outside polygon

Problems With Rays Fire ray from point Count intersections Odd = inside polygon Even = outside polygon Problems Ray through vertex

Problems With Rays Fire ray from point Count intersections Odd = inside polygon Even = outside polygon Problems Ray through vertex Ray parallel to edge

A Better Way One winding = inside

A Better Way zero winding = outside

Requirements Oriented edges Edges can be processed in any order

Computing Winding Number Given unit normal n =0 For each edge ( p 1 , p 2 ) for some integer k k = winding number If k == 1 , then inside If k == 0 , then outside

Advantages Extends to 3D! Numerically stable Even works on models with holes: Odd k : inside Even k : outside No ray casting

Convex Hulls Subset of S of the plane is convex, if for all pairs p,q in S the line segment pq is completely contained in S. The C onvex Hull CH(S) is the smallest convex set, which contains S. p q pq p q pq

Convex hull of a set of points in the plane The convex hull of a set P of points is the unique convex polygon whose vertices are points of P and which contains all points from P. Rubber band experiment

Convex Hull Input: Set S = {s 1 , …, s n } of n points Output: Find its convex hull Many algorithms: Naïve – O ( n 3 ) Insertion – O ( n log n ) Divide and Conquer – O ( n log n ) Gift Wrapping – O ( nh ), h = no of points on the hull Graham Scan – O ( n log n )

Graham Scan Polar sort the points around a point inside the hull Scan points in counter-clockwise (CCW) order Discard any point that causes a clockwise (CW) turn If CCW, advance If !CCW, discard current point and back up

Graham-Scan : (1/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12

Graham-Scan :(1/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 1.Calculate polar angle 2.Sorted by polar angle

Graham-Scan : (2/11) p 2 p 1 p 0 Stack S : p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12

Graham-Scan : (3/11) Stack S : p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 3 p 1 p 0

Graham-Scan : (4/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 4 p 3 p 1 p 0 Stack S :

Graham-Scan (5/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 5 p 3 p 1 p 0 Stack S :

Graham-Scan (6/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 8 p 7 p 6 p 5 p 3 p 1 p 0 Stack S :

Graham-Scan (7/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 9 p 6 p 5 p 3 p 1 p 0 Stack S :

Graham-Scan (8/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 10 p 3 p 1 p 0 Stack S :

Graham-Scan (9/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 11 p 10 p 3 p 1 p 0 Stack S :

Graham-Scan (10/11) p 1 p 0 p 3 p 4 p 2 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 12 p 10 p 3 p 1 p 0 Stack S :

Time complexity Analysis Graham-Scan Sorting in step 2 needs O ( n log n ) . Time complexity of stack operation is O(2 n) The overall time complexity in Graham-Scan is O ( n log n ) . Demo: http://www.cse.unsw.edu.au/~lambert/java/3d/hull.html

Line segment intersection Input: Set S = {s 1 , …, s n } of n line segments, s i = (x i , y i ) Output: k = All intersection points among the segments in S

Line segment intersection Worst case: k = n ( n –1)/2 = O ( n 2 ) intersections Sweep line algorithm (near optimal algorithm): O ( n log n + k ) time and O ( n ) space O( n ) space

Sweep Line Algorithm Avoid testing pairs of segments that are far apart . Idea : imagine a vertical sweep line passes through the given set of line segments, from left to right. Sweep line

Assumption on Non-degeneracy No segment is vertical. // the sweep line always hits a segment at // a point. If a segment is vertical, imagine we rotate it clockwise by a tiny angle. This means: For each vertical segment, we will consider its lower endpoint before upper point.

Sweep Line Status The set of segments intersecting the sweep line. It changes as the sweep line moves, but not continuously . Updates of status happen only at event points . left endpoints right endpoints intersections event points A G C T

Ordering Segments A total order over the segments that intersect the current position of the sweep line: Based on which parts of the segments we are currently interested in A B C D E B > C > D (A and E not in the ordering) C > D (B drops out of the ordering) At an event point, the sequence of segments changes:  Update the status.  Detect the intersections. D > C (C and D swap their positions)

Status Update (1) A new segment L intersecting the sweep line L M K Check if L intersects with the segment above ( K ) and the segment below ( M ). new event point Intersection(s) are new event points. Event point is the left endpoint of a segment. N K, M, N K, L, M, N O

Status Update (2) The two intersecting segments ( L and M ) change order. L M K Check intersection with new neighbors ( M with O and L with N ). Intersection(s) are new event points. Event point is an intersection. N O O, L, M, N O, M, L, N

Status Update (3) The two neighbors ( O and L ) become adjacent. L M K Check if they ( O and L ) intersect. Intersection is new event point. Event point is a lower endpoint of a segment. N O, M, L, N O, L, N O

Data Structure for Event Queue Data structure: balanced binary search tree (e.g., red - black tree). Ordering of event points: by x -coordinates by y -coordinates in case of a tie in x -coordinates. Supports the following operations on a segment s .  inserting an event  fetching the next event Every event point p is stored with all segments starting at p . // O (log m ) // O (log m ) m = # event points in the queue

Data Structure for Sweep-line Status Describes the relationships among the segments intersected by the sweep line. Use a balanced binary search tree T to support the following operations on a segment s . Insert( T , s ) Delete( T , s ) Above( T , s ) // segment immediately above s Below( T , s ) // segment immediately below s e.g , Red - black trees, splay trees ( key comparisons replaced by cross-product comparisons ) . O (log n ) for each operation.

An Example L K M N O K L O N M K L N O  The bottom-up order of the segments correspond to the left-to-right order of the leaves in the tree T .  Each internal node stores the segment from the rightmost leaf in its left subtree.

The Algorithm FindIntersections( S ) Input : a set S of line segments Ouput : all intersection points and for each intersection the segment containing it. 1 . Q  // initialize an empty event queue 2. Insert the segment endpoints into Q // store with every left endpoint // the corresponding segments 3. T   // initialize an empty status structure while Q   do extract the next event point p 6. Q  Q – { p } 7. HandleEventPoint( p )

Handling Event Points Status updates (1) – (3) presented earlier. Degeneracy : several segments are involved in one event point (tricky). T: (a) Delete D, E, A, C (b) Insert B , A, C p A B C C A D G E H l D D A G E H C A C G E G C H A B A G C B

Motivation Transformation of a topographic map into a perspective view

Terrains Given: A number of sample points p 1 ..., p n Requir ed: A triangulation T of the points resulting in a “realistic” terrain. "Flipping" of an edge: 900 930 20 50 900 930 20 50 Goa l: Maximise the minimum angle in the triangulation

Triangulation of Planar Point Sets Given: Set P of n points in the plane (not all collinear ). A triangulation T(P) of P is a planar subdivision of the convex hull of P into triangles with vertices from P . T(P) is a maximal planar subdivision. For a given point set there are only fin itel y many different triangulations.

Size of Triangulations Theorem : Let P be a set of n points in the plane, not all collinear and let k denote the number of points in P that lie on the boundary of the convex hull for P . Then any trianglation of P has 2n-2-k triangles and 3n-3-k edges. Proof : Let T be triangulation of P , and let m denote the # of triangles of T . Each triangle has 3 edges, and the unbounded face has k edges .  n f = # of faces of t riangulation = m + 1 every edge is incident to exactly 2 faces . Hence, # of edges n e = (3m +k)/2 . Euler ‘s formula : n - n e + n f = 2 . Substituting values of n e and n f , we obtain: m = 2n – 2 – k and n e = 3n – 3 – k .

Angle Vector Let T(P) be a triangulation of P ( set of n points ) . Suppose T(P) has m triangles. Consider the 3m angles of triangles of T(P) , sorted by increasing value . A(T) = { a 1 ..., a 3m } is called angle- vector of T . Triangulations can be sorted in lexicographical order according to A(T). A triangulation T(P) is called angle-optimal if A(T(P))  A(T´(P)) for all triangulations T´ of P .

Illegal Edge a 1 ‘ a 2 ‘ a 4 ‘ a 3 ‘ a 5 ‘ a 6 ‘ P k P j P i Edge flip  The edge p i p j is illegal if  ‘ i Note: Let T be a triangulation with an illegal edge e . Let T´ be the triangulation obtained from T by flipping e . Then, A(T´)  A(T) .  i  Consider a quadrilateral: a 1 a 2 a 6 a 5 a 3 a 4 P j P i

Legal Triangulation Definition : A triangulation T(P) is called a legal triangulation , if T(P) does not contain any illegal edges. Test for illegality Lemma : Let edge p i p j be incident to triangles p i p j p k and p i p j p l , and let C be the circle thru p i , p j and p k . The edge p i p j is illegal iff the point p l lies in the interior of C . Furthermore, if the points p i , p j , p k , p l form a convex quadrilateral and do not lie on a common circle, then exactly one of p i p j or p k p l is an illegal edge. P i P k P j P l

Test of Illegality Observation: p l lies inside the circle through p i , p j and p k iff p k lies inside the circle through p i , p j , p l . When all four points lie on circle, both p i p j and p k p l are legal. P i P k P j P l P i P k P j P l

Thales Theorem asb   aqb =  apb   arb illegal Lemma: Let C be the circle through the triangle p i , p j , p k and let the point p l be the fourth point of a quadrilateral. The edge p i p j is illegal iff p l lies in the interior of C. a b p q r s p i p j p k p l

Thales Theorem Consider the quadrilateral with p l in the interior of the circle that goes through p i , p j , p k . Claim: The minimum angle does not occur at p k ! (likewise: Minimum angle does not occur at p l ) G oal: Show that p i p j is illegal p i p j p l p k p i p j p l p k

W.l.o.g. a 4 minimal Thales Theorem: Proof a 1 ‘ > a 4 a 2 ‘ > a 2 p i p j p k p l a 1 a 2 a 3 a 4 p i p j p k p l a 1 a 2 a 3 a 4 a 1 ‘ a 2 ‘ a 3 ‘ a 4 ‘ p i p j p k p l a 1 ‘ a 2 ‘ a 3 ‘ a 4 ‘

Circle criterion violated  illegal edge Thales Theorem: Proof a 4 ‘ > a 1 a 3 ‘ > a 3 Hence, min{a i ‘} > min{a i } p i p j p k p l a 1 a 2 a 3 a 4 a 1 ‘ a 2 ‘ a 3 ‘ a 4 ‘

Circle Criterion Definition: A triangulation fulfills the circle criterion if and only if the circumcircle of each triangle of the triangulation does not contain any other point in its interior .

Theorems Theorem: A triangulation T(P) of a set P of points does not contain an illegal edge if and only if nowhere the circle criterion is violated . Theorem: Every triangulation T(P) of a set P of points can be finally transformed into an angle-optimal triangulation in a finite number of steps.

Definition of the Delaunay Triangulation A triangulation T(P) is a Delaunay Triangulation of P, denoted as DT(P), if and only if the circumcircle of any triangle of T does not contain any other point of P in its interior (i.e. T fulfills the circle criterion).

Equivalent Characterisations of the Delaunay Triangulation DT(P) is the straight - line - dual of Voronoi Diagram VD(P) . DT(P) is a triangulation of P such that all edges are legal ( local angle - optimal ) . DT(P) is a triangulation of P such that for each triangle the circle criterion is fulfilled . DT(P) is global angle - optimal triangulation . DT(P) is a triangulation of P such that for each edge p i p j ther e is a circle , on wh ich p i and p j lie and which does not contain any other point from P .

Algorithm DT(P) (randomized, incremental) Given: Point set P = {p 1 ..., p n } Initially : Compute triangle (x, y, z) , which includes the points p 1 ..., p n . m z (0,3m) y (3m,0) x (-3m,-3m)

Algorithm DT(P) m = max { |x i |,|y i | } T = ((3m, 0), (3m , 3m), (0, 3m)) initialize DT(P) as T . permutate the points in P random ly . for r = 1 to n do find the triangle in DT(P), which contains p r ; insert new edges in DT(P) to p r ; legalize new edges. 4. remove all edges, which are connected with x, y or z .

Inserting a Point p i p j p r 2 cases : p r is inside a triangle p r is on an edge p i p j p k p r p l p i p j p k p r Legalize (p r ,p i p j ,T): if p i p j is illegal then let p i p j p k be the triangle adjacent to p r p i p j along p i p j . flip p i p j , i.e. replace p i p j by p r p k Legalize (p r , p i p k , T) Legalize (p r , p k p j , T)

Algorithm Delaunay Triangulation Input: A set of p oint s P = {p 1 ..., p n } in general position Output : The Delaunay triangulation DT(P) of P DT(P) = T = (x, y, z) for r = 1 to n do find a triangle p i p j p k  T , th at contains p r . if p r lies in the interior of the triangle p i p j p k then split p i p j p k Legalize( p r, p i p j ), Legalize( p r, p i p k ), Legalize( p r, p j p k ) if p r lies on an edge of p i p j p k (say p i p j ) then split p i p j p k and p i p j p l Legalize (p r , p i p l ) , Legalize (p r , p i p k ), Legalize (p r , p j p l ) , Legalize (p r , p j p k ) Delete (x, y, z) with all incident edges to P

Correctness Lemma : Every new edge created in the algorithm for constructing DT during the intersection of p r is an edge of the Delaunay graph of   { p 1 ,...,p r } . pq is a Delaunay edge iff there is a (empty) circle, which contains only p and q on the circumference . Proof idea : Shrink a circle which was empty before addition of p r !

Correctness Observation: After insertion of p r , every new edge produced by edge-flips is incident to p r !  Correctness of the algorithm : Consider newly produced edges: p r p r

Edge Flips Edge-flips produce only legal edges. Before inserting p r , circle that goes through p i , p j , p k was empty! Edge- flips produce edges that are always in c ident to p r ! p r p i p j p k

Data Structure for Point Location t 1 t 2 t 4 t 5 t 3 t 6 t 7  Split t 1 p i p j p i p k t 1 t 2 t 3 t 2 t 3 t 4 t 5 t 4 t 6 t 7 t 1 t 2 t 3 t 1 t 2 t 3 t 1 t 2 t 4 t 5 t 3  flip p i p j  flip p i p k

Analysis of the Algorithm Lemma : The expected number of triangles created by the incremental algorithm for constructing DT(P) is atmost 9n + 1 .

Analysis of the Runtime Theorem : The Delaunay triangulation of a set of P of n points in the plane can be computed in O(n log n) expected time, using O(n) expected storage. Proof : Running time without Point Location : Proportional to the number of created triangles = O(n). Point Location : The time to locate the point p r in the current triangulation is linear in the number of nodes of D that we visit.

Relation to Euclidean MST Problem: Find a spanning tree on a set P of N points in a plane such that the total length of all the edges of the tree is minimized (Euclidean MST) Relation of Delaunay Triangulation A Euclidean MST is a sub-graph of a Delaunay Triangulation Blindly applying MST finding methods like Kruskal’s or Prim’s on this problem gives O ( m ) where m = n(n-1)/2 . If one builds the DT(P) in O ( n log n ) and then does Kruskal’s or Prim’s, one get an O ( n log n ) solution.

Geometry and Air Traffic Control » Avoid collisions at all costs. » Start by f inding the two aeroplanes closest to each other, at specific elevations. » During peak times over 5000 planes, per hour, around Frankfurt airspace.

Geometry and Astronomy » Find the two nearest objects in space.

Search Quickly… and Why ? » Find the two nearest objects in space.

Problem Description Input: A set S of n points in the plane. Output: The pair of points with the shortest Euclidean distance between them. Other applications in fundamental geometric primitive: Graphics, computer vision, geographic information systems, molecular modeling Special case of nearest neighbour, Euclidean MST, Voronoi

Closest Pair of Points Naïve method: Brute force. Check all pair of points p and q with Θ ( n 2 ) comparisons. 1D version Simple if all points are on a single line. Run time O ( n log n ) Key to reducing quadratic runtime of the na ï ve method: Simply don´t visit each point n times.

Algorithm: DAC – Divide Sort S in ascending order of x -coordinates. Find the median line L , so that we have ½ n points in S L and ½ n points S R . S L S R

Algorithm: DAC – Conquer Recursively, find the closest pair in S L and S R . S L S R

Algorithm: DAC – Merge (Strip-Margins) Assuming that distance <  S = min(  SL ,  SR ). Observation: Only need to consider the points within  S of line L . S L S R  S  S = min(12,21)

Handling Merge Sort points in 2  S -strip in ascending y -coordinate order. ( First observation ) Only check those points within 11 positions in sorted list. S L S R  S  S = min(12,21)

Comparing Efficiently Definition: Let s i be the point in the 2  S -strip with the i th smallest coordinate. Claim: If | i - j |  12, then the distance between s i and s j is at least  S . Proof: The points s i and s j cannot lie in same ½  S -by- ½  S box. s i and s j are at least 2 rows apart, and have distance  2(½  S ). Fact: The above still holds if we replace 12 with 7.  S  S

Closest Pair: DAC Algorithm Algorithm: Closest-Pair ( S = { p 1 , …, p n } ) Sort all n points in ascending x -coordinate order Compute the median line L , and divide S into S L and S R Set  SL  Closest-Pair ( S L ) Set  SR  Closest-Pair ( S R ) Set  S  min(  SL ,  SR ) Delete all points further than  S  from separation line L Sort remaining points in ascending y -coordinate order Scan points in y -order and compare distance between each point and next 7 neighbours. Update  S if new shorter distance is found. Return  S

Closest Pair: DAC Analysis Run time: T ( N )  2 T ( N /2) + O ( N log N )  T ( N ) = O ( N log 2 N ) Can we achieve a better run time? Yes. Don’t sort points in strip from scratch each time. Each recursive step returns two lists: all points sorted by y -coordinate, and all points sorted by x -coordinate. Sort by merging two pre-sorted list. T ( N )  2 T ( N /2) + O ( N )  T ( N ) = O ( N log N ) Space: O ( n )

Algorithm: Scan-Line Sort S in ascending order of x -coordinates. Place sorted points p 1 , …, p n in the event-queue. Initialise   distance( p 1 , p 2 ). Start scanning from p 3 . » Event-points? » Status structure?

The Status-Structure Maintain the status-structure (  -slab): Insert the new event-point p i . Remove all points that are greater than  distance to the left of the current event-point p i . All points in the status structure are maintained in ascending y -coordinate order.

The Status-Structure Maintain the status-structure (  -slab): Insert the new event-point p i . Remove all points that are greater than  distance to the left of the current event-point p i . All points in the status structure are maintained in ascending y -coordinate order. Perform the ‘ Efficient Comparison ’, as described earlier, BUT only between the new event point p i and at most 7 other nearest points to p i . Update  if a shorter distance is found.

Closest Pair: Scan-Line Analysis Invariant: Whenever an event-point p i has been handled, then  is the minimum distance between any pair of points encountered from p 1 up to p i . Run time: O ( n log n ) Space: O ( n )

Related Problems Nearest neighbour query: Given a set S of n points, and a query point q , determine which s  S is closest to q . Brute force: O ( n ) Goal: O ( n log n ) pre-processing, O (log n ) per query. Voronoi region, VR ( S ): A set of points closest to a given point. Voronoi diagram, VD ( S ): Planar subdivision deliniating Voronoi regions. Delaunay Triangulation (the dual of VD ( S )).

Acknowledgment Slides prepared mainly from the materials of Prof. Dr. Thomas Ottman

A Tutorial on Computational Geometry

More Related Content

What's hot

Viewers also liked

Similar to A Tutorial on Computational Geometry

Recently uploaded

In this document

A Tutorial on Computational Geometry