# Geometric Algorithms Geometric Algorithms Suman Sourav Paramasiven Appavoo Anuja Meetoo Appavoo Li Jing Lu Bingxin Suhendry Effendy Dumitrel Loghin Introduction & Motivation Suman Sourav

Introduction Computational geometry is a branch of computer science devoted to the study of algorithms which can be stated in terms of geometry. Deep connections with classical mathematics and theoretical computer science on the one hand, and many ties with applications on the other. Focus is on discrete nature of geometric problems as opposed to continuous issues. Introduction

Any problem that is to be solved using a digital computer has to be discrete in form. For CG techniques to be applied to areas that involves continuous issues, discrete approximations to continuous curves or surfaces are needed. People deal more with straight or flat objects or simple

curved objects, than with high degree algebraic curves. Introduction Programming in CG is a little difficult. Libraries like LEDA and CGAL implement various data structures and geometric algorithms. CG algorithms suffer from the curse of degeneracies. So, certain simplifying assumptions at times like no three points are collinear, no four points are cocircular,

etc are needed to be made. Motivation & Application There are many areas that give rise to geometric problems. Computer graphics

Computer vision and image processing Robotics Computer-aided designing (CAD) Geographic information systems (GIS) Telecommunication VLSI design (chip layout) Example To find the nearest phone booth on campus Example The motion planning problem Example

To build a model of the terrain surface, we can start with a number of sample points where we know the height. How do we interpolate the height at other points? Nearest neighbor interpolation

Piecewise linear interpolation by a triangulation Moving windows interpolation Natural neighbor interpolation For interpolation, it is good if triangles are not long and skinny. Example Each time you click inside the editing area of a diagram or inside a web page, the application must know inside which object

the mouse pointer is. Given the high frequency of such a query, a fast algorithm must be used. Facility location Given a set of houses and farms in an isolated area. Can we place a helicopter ambulance post so that each house and farm can be reached within 15 minutes? Where should we place an antenna so that a number of locations have maximum

Reception? Art gallery problem How many cameras are needed to guard a given art gallery so that every point is seen? A Real Life Application Voronoi Voronoi diagram - Properties - Construction methods Fortunes Algorithm - How does it work?

- Implementation aspects - Data structures - Pseudocode - Complexity (time & space) Paramasiven Appavoo & Anuja Meetoo Appavoo Voronoi diagram Properties Voronoi edges: Each point on an edge of the Voronoi diagram is

equidistant from its two nearest neighbors pi and pj. Voronoi vertices: It is the center of the circle passing through these sites, and this circle contains no other sites in its interior. Degree: If we assume that no four points are cocircular, then the vertices of the Voronoi diagram all have degree three. Voronoi diagram Properties

Convex hull: A cell of the Voronoi diagram is unbounded if and only if the corresponding site lies on the convex hull. Size: If n denotes the number of sites, then the Voronoi diagram is a planar graph (if we imagine all the unbounded edges as going to a common vertex infinity) with exactly n faces.

Voronoi Construction methods Perpendicular bisector Divide and conquer Fortunes algo

Voronoi Fortunes algo Time complexity of perpendicular bisector Extra computation at the merge of divide and conquer Steven Fortune developed a plane sweep algorithm called the fortunes algorithm that uses: o a sweep line, o a beach line, o site event, and

o circle event to generate the diagrams Steven Fortune, Bell Laboratories Fortunes Algo Fortunes algo - Concepts Sweep line

The idea behind algorithms of this type is to imagine that a line is swept or moved across the plane, stopping at some points. A sweep line splits the problem domain into two regions, an explored region and an unexplored region Operations are restricted to objects that are either intersecting or are in the immediate vicinity of the sweep line whenever it stops The complete solution is available once the line has passed over all objects. It also applies an ordering to the problem

Fortunes algo - Concepts Site events A site event happens when the sweep line hits a site Property of the line perpendicular to the sweep line at the site event o All points are equidistant to: The site The closest point on the sweep line (the site site event itself, at this stage)

What happens when the sweep line moves further away? Fortunes algo - Concepts What happens The set of points that are equidistant to both the site and the closest point on the sweep line degenerate from

a line to form a parabola... The parabola defines the set of points... Set of points that belong to this sites voronoi cell Why are we doing this? Important property implying that the cross section of such parabolas, from a number of sites will define the edges of a voronoi cell.. Fortunes algo - Concepts

2 sites - Crossing parabolas P1 P2 q P1 p1 and p2 are two sites

As the sweep line goes down, the intersection of the parabolas, breakpoints, are the endpoints for a partial edge... At the site event of p2, point q is equidistant from p1 and p2 P1

P2 site event part of voronoi diagrams edge breakpoints Fortunes algo - Concepts More sites & beach line

This dividing line is termed the beach line The algorithm for computing the Voronoi diagram of a set of points depends entirely on how this beach line changes as the sweep moves through the space. The beach lines topology changes when a new arc is added or another is absorbed. parts of the edges of the voronoi diagrams... breakpoints What is the stopping condition for an edges growth?

Fortunes algo - Concepts Circle event... One endpoint is formed with a circle event o a vertex in the voronoi diagram. Each relating site is equidistant from the vertex. Vertex is the centre of a circle on which the three points lie vertex

circle event Either by cross sections of parabolas or simply by using the intersection of each or the perpendicular bisector of the inner triangle Fortunes algo - Concepts From the vertex... vertex q a b c

Now from the vertex... If a circumcircle is empty in its interior then, in a Voronoi diagram: a, b, c would be Voronoi sites q would be a Voronoi vertex The perpendicular bisectors of abc would be Voronoi edges. Fortunes algo Implementation Aspects...

Parabolas defining beach line not computed as sweep line moves through the problem space o Computationally expensive and unnecessary o Calculations required only when the beach line changes topology Site events Circle events Sweep line makes discrete steps, rather than a continuous sweep Fortunes algo

Data Structures Priority queue o State of sweep line o Indicates when topology of beach line can change Site events Circle events Doubly connected edge list (DCEL) o State of Voronoi diagram o Stores

Vertex records Edge records Face records Balanced binary search tree ( AVL or red-black tree) o State of beach line Fortunes algo Priority Queue of Events Priority based on their y-coordinates Site events are represented by their

coordinates Circle events o Computed on the fly o Represented by coordinates of lowest point of empty circle touching 3 sites o Anticipated and may be false s1 (x1, y1) s2 (x2, y2) s3 (x3, y3) s5 (x5, y5) s4 (x4, y4) c1 (x, y) circle event added to queue Fortunes algo DCEL

Vertex records Vertex Coordinates Incident edge v1 (x1, y1) e12

v2 (x2, y2) e23 ... ... ... v8

(x8, y8) e87 Edge records Half edge Source vertex Twin

e23 v2 .. ... Face records Incident face Prev

edge Next edge e32 s3 e52 e37 ...

... Face Incident edge s1 e32 s2 e56

... ... Fortunes algo Balanced Binary Search Tree Internal nodes o Represent breakpoints between two arcs (tuple) o Contains pointer to record of edge being traced Leaf nodes o Stores site defining arc on beach line

o Contains a pointer to a potential circle event Fortunes algo Pseudocode 1. 2. 3. 4. 5. Initialise data structures

Event queue, Q all site events Binary search tree, T DCEL, D while(queue not empty) event = pop first event from queue process(event, T, D) finish all edges in binary search tree How to process site event? - Add the event and breakpoints to BST - Add new edge in DCEL

- Delete false alarms for circle events - Add potential circle event(s) to queue How to process circle event? - Update breakpoints and leaf in BST - Record new vertex - Delete false alarms - Add possible circle event to queue Fortunes algo Processing Site Event

Locate existing arc (if any) that is above the new site o x-coordinate of new site is used for binary search o x-coordinate of each breakpoint along root to leaf path is computed on the fly Q s5 s4 Fortunes algo Processing Site Event

Break the arc o Replace the leaf node with a sub tree representing the new arc and its breakpoints Q s4 Different arcs can be identified by the same site!!! Fortunes algo Processing Site Event

Add new edge record in DCEL Create new face record with pointer to new edge New half-edge record Q s4 Fortunes algo Processing Site Event Check for potential circle events o Scan for triple of consecutive arcs and determine if breakpoints converge

Add potential circle event to to priority queue!!! Potential circle event Store a pointer to circle event in leaf node for s1!!! Q c1 s4 Fortunes algo Processing Site Event

Converging breakpoints may not always yield a circle event o Appearance of a new site before the circle event makes the potential circle non-empty Original circle event becomes a false alarm! Fortunes algo Processing Circle Event

Create new vertex record Add vertex to corresponding edge record Delete disappearing arc Q s4 Fortunes algo Processing Circle Event Delete disappearing arc

Q s4 Fortunes algo Processing Circle Event Q s4 Fortunes algo Processing Circle Event Create new edge record New half-edge record New edge being traced by new breakpoint

Check new triple edges for potential circle events Fortunes algo Algorithm Termination When Q is empty, beach line and its breakpoints continue to trace edges o Terminate half-infinite edges via a bounding box Fortunes algo Complexity Steps in handling site events

Running Time Locate leaf representing existing arc above new site o Delete potential circle event from queue O(log n) Break arc by replacing leaf node with a sub

tree representing new arc and break points O(1) Add new edge and face records to DCEL O(1) Check for potential circle event(s) o Add event to queue if they exist o Store pointer to event in proper leaf

O(log n) Fortunes algo Complexity Steps in handling circle events Running Time Delete disappearing BST leaf node and its associated circle events from event queue

Add vertex record in DCEL O(1) Create new edge record in DCEL O(1)

Check for potential circle event(s) O(1) O(log n) Fortunes algo Complexity Time complexity for each event: O(log n) How many events are there? o Number of site events + Number of circle events

o How many site events? n o How many circle events? Each circle event corresponds to a vertex 2n - 5 (Eulers formula) False alarms deleted before they are processed Total number of events = 3n - 5 Overall running time: O(n log n)

Storage complexity: O(n) Delaunay Triangulation Delaunay and Voronoi Delaunay properties Randomized Algorithm - Idea - Implementation aspects - Pseudocode - Data structure - Complexity (time & space) Li Jing & Lu Bingxin Delaunay and Voronoi

An intuitive conception (a) Voronoi diagram (b) Delaunay triangulation General position assumption: no 4 points are co-circular Delaunay and Voronoi Delaunay and Voronoi complexes are dual to each other Dual correspondence Voronoi complexes Delaunay triangulation

cells (regions) vertices edges edges vertices faces Delaunay Triangulation (DT)

Whats the difference? Properties Circumcircle property The circumcircle of any triangle in DT{S) is empty. (It contains no points of S in its interior.) Proof By general position assumption, the degree of all Voronoi vertex is 3, edge vu exists

Consider S4, S1S4 is perpendicular to l and divided half by l . S4 is outside the circle In fact, one definition of Delaunay Triangulation is: Delaunay Triangulation is a triangulation that circumcircle of each triangle is empty. Properties Empty circle property Two points are connected by an edge in the Delanuay triangulation.

There is an empty circle passing through these two points. Proof Trivial from the circumcircle property. There are a series of circumcircles that pass through Si and Sj Note down the centres of these circumcircles until they pass through another point

The track is the bound for both cell(Si) and cell(Sj) Si Sj is an edge in the DT(S) Locally Delaunay edge An edge ab is locally Delaunay if it belongs to only one triangle, or it belongs to two triangles, abc and abd, and d lies outside the circumcircle of abc.

If an edge is locally Delaunay,we also call it legal, otherwise illegal. Delaunay Lemma If every edge in Ts is locally Delaunay, then Ts is the Delaunay triangulation of S. Let x be an arbitrary point in Proof abc Let abc=0, A1, A2, Ak be the sequence of triangles that intersect xp Let di(p) = |p ai|2 ri2

Because the edges along xp are locally delaunay, d0(p)> d1(p)> > dk(p) d (p) = 0, so d (p)>0 Edge Flipping Flip all edges in a triangulation until they are all locally Delaunay edges Randomized Incremental Algorithm Randomized the points as p1, p2, pn Find a sufficiently large triangle that contains P Insert p1, then p2... and finally pn suppose we have computed DT(Pi-1)

insert pi which splits a triangle into three perform edge flips until no illegal edge remains we have just computed DT(Pi) Repeat the process until i = n First step Find a sufficiently large triangle that contains P First step First step Example

insert p Example split abc into abp, bcp, and acp Example check edges ab, bc, and ac Example edge ab is illegal

flip it Example edge ab is flipped into pd edge ad and bd are to be checked edge ad is legal, keep it Example edge bd is illegal flip it Example

edge bd is flipped into pe edge ed and be are legal, keep them Example edge bc is illegal flip it Example edge bc is flipped into fp edge bf and cf are legal, keep them

Example edge ac is illegal flip it Example edge ac is flipped into pg no more edge to flip: we are done Pseudocode Algorithm DelaunayTriangulation(P) Input: a suitably shuffled (permuted uniformly at random) set of points P = (p1, p2, p3,., pn)

Output: DT(P) (* use a global DCEL to store DT(P) *) 1. Find a sufficiently large triangle T(p-3p-2p-1) containing P 2. for i = 1 to n do 3. Insert(pi) 4. Endfor Algorithm Insert(p) 5. Discard the triangle T(p-3p-2p-1) Input: a point p, a set of point P and T = DT(P) Algorithm SwapTest(ab) Output: DT(P u {p}) 1. if ab is an edge of the exterior face of DT(P) 1. Find the triangle T(abc) of DT(P) 2. do return

containing p 3. d <- the vertex (other than a,b) of the triangle (* use conflict lists *) adjacent to triangle T(pab) along edge ab 2. Insert edges pa,pb and pc 4. if inCircle(p, a, b, d) < 0 (* update conflict lists *) 5. do Flip edge ab for pd 3. SwapTest(ab) (* update conflict lists *) 4. SwapTest(bc) 6. SwapTest(ad) 5. SwapTest(ca) 7. SwapTest(db)

Conflict list Conflict --- a non-inserted point is inside a triangle in the bi-directional pointer current triangulation Current triangle Non-inserted points T(p1p2p3) p7, p8 ...

... T(p4p5p6) p9 ... ... Triangles in the current Delaunay Triangulation non-inserted points

Non-inserted point Current triangle p7 T(p1p2p3) p8 T(p1p2p3)

p9 T(p4p5p6) ... ... Each triangle of the current triangulation --- Bucket Time complexity

Major steps in the algorithm Find a sufficiently large triangle Find the triangle containing a non-inserted point Update the triangulation Update conflict lists Find a sufficiently large triangle M: maximum absolute value of either x or y coordinate of all the points in P Time cost: O(1)

Find the triangle containing a non-inserted point Query the conflict list Non-inserted point Current triangle p7 T(p1p2p3) p8

T(p1p2p3) p9 T(p4p5p6) ... ... Time cost for one iteration: O(1)

Time cost for all the n iterations: O(n) Backward analysis Imagine that the algorithm is run backwards starting from the delaunay triangulation we have at the end In analyzing the ith step, imagine that we are deleting one of the i points in the current triangulation and then update the triangulation The work done in this case is the same as running the algorithm forward Assume that each of the i points is equally likely to be deleted at the ith step, since the points were

added randomly in the original algorithm Time to update triangulation Consider the ith step of the algorithm Pi: the set of the first i points (p1, p2, p3,., pi) in the whole point set P, i >3, or the set of points in DT(Pi) Run the ith step backward (deleting a random point p in Pi ) Note that each new edge added in DT(Pi) due to the insertion of p is incident to p The total number of edge changes made in the triangulation for the insertion of p is proportional to

Time to update triangulation The total degree of the vertices in Pi is less than 6i There are at most 3i edges in DT(Pi) (planar graph) The sum of vertex degrees is twice the number of edges in a graph The expected degree of a random point in Pi is at most 6 The expected number of edge changes for inserting a random point p is O(1) Summing over all the n steps, the time for updating the triangulation is O(n) Update conflict lists Each non-inserted point is assigned to a triangle/bucket For a non-inserted point, if the triangle which it lies

within is destroyed, we have to find another new triangle containing this non-inserted point The expected time to update conflict lists is the expected time to rebucket non-inserted points is proportional to the expected number of noninserted points required to be rebucketed Rebucket points Consider the ith step of the algorithm Pi: the set of the first i points (p1, p2, p3, , pi) in the whole point set P, i >3, or the set of points in DT(Pi) P\Pi : the set of non-inserted points (pi+1, pi+2, pi+3, , pn) Run the ith step backward (deleting a random point in Pi) some triangles in DT(Pi) are destroyed

some points in P\Pi are required to be rebucketed Time to rebucket points Assume a random point p in Pi is deleted For a random point q in P\Pi, suppose that q is bucketed in the triangle T(abc) of DT(Pi) If p is one of the three vertices of T(abc), T(abc) will be destroyed and q will be rebucketed Pr (p is deleted) = 1/i Pr (T(abc) is destroyed) = 3/i Pr (q needs to be rebucketed) = 3/i The expected number of points in P\Pi required to be rebucketed is 3(n-i)/i

Time to rebucket points Summing over all the n steps The total expected number of non-inserted points required to be rebucketed is O(nlogn) The expected time to update conflict lists is O(nlogn) Time complexity Major steps in the algorithm Running time

o Find a sufficiently large triangle O(1) o Find the triangle containing a point O(n) o Update the triangulation o Update conflict lists O(n) O(nlogn)

Complexity Expected time complexity: O(nlogn) Expected space complexity: O(n) DCEL, Conflict list Trapezoidal Decomposition Concepts Randomized Algorithm Motivation - Point Location Complexity Analysis Dumitrel Loghin & Suhendry Effendy

Defining the problem Given a set S of n segments in the plane, with no two distinct end-points having the same x coordinate (general position) These segments intersect in k points There is a bounding rectangle R containing all n segments If we draw vertical lines through each end-point or intersection point, till they intersect a segment or R, we obtain a set of trapezoids T(S) which is the trapezoidal decomposition of S Example R

q2 p3 k1 k2 q1 p1 p2 k3

q3 Details Construction Example Randomized Algorithm Start with a random permutation S = { s1, , sn }, with a bounding rectangle R and with T(S)= FOR i = 1 TO n DO # add segment si and update T(S) Get the trapezoid which contains left end-point of si and draw a vertical line. Move to the right on the planar graph and get all the

intersected trapezoids: If the segment intersect the bottom or top segment of the trapezoid, the point is one of the k intersection points, hence, draw a vertical line. Draw a vertical line at the right end-point. Rescan the intersected vertical lines and trim them based on their starting point; merge the trapezoids which share the deleted line. Point Location Given a query point q, find in which trapezoid it lies Data structure DAG internal nodes are end-points, intersection points or segments

leaves are trapezoids and they may be shared (in-degree of a leaf may be more than one) Can be constructed simultaneously with trapezoidal decomposition Algorithm key idea: At step i, some of the leaves (trapezoids) become internal nodes (end-points of the segment si or intersection points) and new trapezoids are added. Point Location Example 1

2 7 k1 s1 q2 8 5 s2

p2 3 6 4 q1 q p1 10

9 1 2 7 s1 p1 5

k1 8 10 q2 4 q1 q Point Search in DAG

s2 p2 3 6 9 p1 q1 1 s1

4 p2 q2 k1 3 s2 6 k1 9

5 10 s2 3 8 7 Complexity Analysis Assumption: points (end of segment) are in general position

= no two points lies in a same vertical line. = all x-coordinate are different. Later we will see how to lift this assumption. Complexity Analysis What is the complexity of adding one segment s? find the left end point of segment s. trace (and update) all trapezoids intersected by s.

Adding 1 Segment Adding 1 Segment How to find this trapezoid? Finding Trapezoid How to find trapezoid which contain the left end-point of segment s? bi-directional couple of wayspointer Trapezoid Segment

Segment Trapezoid s5 T1 s6 T3 T1 s5, s7, ...

s7 T1 T2 s8, ... s8 T2 T3

s6, s9, s10, ... s9 T3 ... ... s10 T3

Bi-directional pointer bi-directional pointer query find in O(1) need to update the pointer for each changed trapezoid (only update trapezoids which are intersected by new segment s). Adding 1 Segment Adding 1 Segment Adding 1 Segment Complexity Analysis For each segment s, we need to:

1. find the left end-point of s. 2. trace intersected trapezoids. 3. update the trapezoids. 4. update the bi-directional pointers. For each trapezoid: Overall: t(f) : the complexity of trapezoid f. p(f) : update the bi-directional pointer for trapezoid f. Backward Analysis Imagine the algorithm run backwards, deleting segment one at a time.

When we delete a segment s from Hi, only trapezoids in Hi which adjacent to s will be affected. Complexity Analysis Since we insert (or delete in backward analysis) segment s in random order, then every remaining segment is equally likely to be chosen. E(i): the expected complexity of deleting the ith segment. (in backward analysis) Observation Each trapezoid is adjacent to at most 4 segments. (using general position assumption)

Observation Trapezoid can have an arbitrary number of adjacent segments in non general position. We will deal with this case later. Complexity Analysis Each trapezoid is adjacent to at most 4 segments. = Each trapezoid appears in at most 4 segments adjacency list

Complexity Analysis = all the left endpoints of segment that have yet to be added (or have been deleted in backward analysis). Complexity Analysis = proportional to decomposition the k = the number of intersection

number of vertices in all Complexity Analysis Complexity Analysis

Randomized vs. Deterministic Non General Position How to handle non general position? Rotation Transformation or

Shear Transformation x = x + y Summary Voronoi Diagrams -- Use Fortune Algorithm Delaunay Triangulation -- Randomized Incremental Construction Dual of Voronoi Trapezoidal Decomposition -- Randomized Incremental Construction Conclusion

Widely used in various other areas. We use it sometimes without even realising it. Lot of potential of further development. Numerous interesting open problems. http://compgeom.cs.uiuc.edu/~jeffe/open/

Thank You References Voronoi Derek Johns, "An Optimal Algorithm for Computing 2D Voronoi Diagrams Fortune's Sweep Algorithm", Available at http://cgm.cs.mcgill.ca/~mcleish/644/Projects/DerekJohns/Sweep.htm [Accessed February 2014]

Dheeraj Kumar Singh, "Lecture 20: Voronoi Diagrams and Fortunes Algorithm", Available at http://intinno.iitkgp.ernet.in/courses/91/wfiles/37906 [Accessed February 2014] Voronoi Diagrams, Available at http://ima.udg.edu/~sellares/ComGeo/Vor2D_1.ppt [Accessed February 2014] Delaunay Triangulation http://www.cs.umd.edu/~mount/754/Lects/754lects.pdf http://www.cs.uu.nl/geobook/interpolation.pdf http://www.comp.nus.edu.sg/~hcheng/academic/courses/cs5237/notes/04.pdf http://www.comp.nus.edu.sg/~hcheng/academic/courses/cs5237/notes/05.pdf http://www.comp.nus.edu.sg/~tantc/ioi_training/CG/l9cs4235.pdf http://www.comp.nus.edu.sg/~tantc/ioi_training/CG/l10cs4235.pdf http://groups.csail.mit.edu/graphics/classes/6.838/F01/lectures/Delaunay/Delaunay2D.ppt

Trapezoidal Decomposition Rajeev Motwani, Prabhakar Raghavan, Randomized Algorithms, 1995 Subhash Suri, Point Location, Available at http://www.cs.ucsb.edu/~suri/cs235/Location.pdf