Computational Geometry
COMPSCI 164
- Baseline
- Convex Hull
- Line Segment Intersections
- Doubly Connected Edge List (DCEL)
- Polygon Triangulation (Art Gallery Guarding Problem)
- Lower Envelope
- Linear Programming
- Orthogonal Range Searching
- Point Location
- Voronoi Diagrams
- Point Line Duality
- Line Arrangements
- Delaunay Triangulation
- 3D Convex Hull
- Windowing
Baseline
-
Real Random Access Machine (RAM) Model is the basis for algorithmic analysis
- Stored program that has a memory and CPU
- Each memory cell can store a real number; supports all data structures
- Performs basic arithmetic operations (add, subtract, multiply, divide) in O(1) time
- Depending on definition, can also perform more complicated ones like powers or rounding
Convex Hull
- Convexity: A shape S is convex if for any two points $p, q \in S$ the line $\bar{pq}$ is covered wholly by the shape
- Convex hull: Any set of points that are convex and contains a set of constraining points
- Problem statement: Given n points, find a convex hull that covers the points
Line Sweep Algorithm
- Stop at event points to form the upper and lower convex hulls separately
- Add points to hull linearly, and if the last three points form an “incorrect turn” (i.e. a left turn for the upper hull, or a right turn for the lower hull), remove offending points
- O(n log n) time; must sort points by x-coordinate first (O(n log n)) and then form the convex hulls linearly (O(n))
- Start with point with the smallest x-coordinate, end with point with largest x-coordinate
- These two are guaranteed to be on the hull
Gift Wrapping (Jarvis March)
- Start at lowest point, draw a horizontal line of infinite length
- Look at all points to determine which one will be reached first when rotating the line
- O(hn) runtime, where h is the size of the hull
- Each point on hull requires n comparisons to the other points
Divide and Conquer
- Separate the set of points S into two sets, S1 and S2
- Calculate the upper hull for S1 and S2, and then merge them by adding a bridge edge
- No points will lie above this bridge edge, thus keeping the property of the upper hull containing all points
- Can use binary search for a O(log n) merge step
Quick Hull
- Select the leftmost and rightmost points in S, p and q
- Select the point furthest from $\bar{pq}$, r
- Prune all points in the triangle formed by pqr and call this algorithm on $\bar{pr}$ and $\bar{rq}$
- O(nh) runtime
Randomized Quickhull
- Select the leftmost and rightmost points in S, p and r
- Select a random point and shoot a ray perpendicular to $\bar{pr}$
- The bridge edge (AKA edge furthest away from $\bar{pr}$) defined by points s and t will be in the convex hull
- Add s and t to the hull
- Call randomized quickhull on $\bar{ps}$ and $\bar{tq}$
- Expected runtime is O(n log h)
Line Segment Intersections
- Problem: Given a set S of closed line segments, return all intersections
- Calculating if two lines intersect is O(1); can use matrix operations
- Brute force: O(n2); compare all lines with each other
- Run time can be no better than O(n log n + k)
- Element uniqueness problem can be reduced into line segment intersections, and that problem is provably O(n log n)
Plane Sweep Algorithm
- Cleanliness property: all points to the left of the sweep line have been processed
- Sweep Line Status: store segments that currently intersect with the sweep line in order of their y-coordinate intersections
- Events
- Points in time when sweep line status changes
- e.g. status goes from A, B to B, A
- Endpoints of segments
- Intersection points
- Points in time when sweep line status changes
- Event handling
- Left endpoint: Add segment to sweep line status, check for intersection with adjacent segments, add new intersection point(s) to queue
- Intersection event: swap order of intersecting segments, check intersections of swapped segments with new adjacencies
- Right endpoint: delete segment from status, check for intersections in newly adjacent segments
- Works by intersection lemma: all intersecting lines must become adjacent at some point
- Runs in O((n + k)log n), where k is the size of the output
Doubly Connected Edge List (DCEL)
- Representation of a map
- Records each edge, face, and vertex
- Stores geometric information, topological information, attribute information, etc
- Used for traversing all faces of a planar intersection, visiting all edges around a given vertex or face
- Main design paradigm: utilizing half-edges
- Each edge borders two faces, so use two half-edges to represent one edge with both half-edges going in different directions
- Each face’s half-edges are oriented counter-clockwise, excluding the outer face
- Object information
- Vertex stores coordinates, pointer to an incident half-edge that considers the vertex its origin
- Face stores a pointer to one of its half-edges to begin traversal
- Half-edge stores origin vertex, twin half-edge (half-edge that combines with this one to create a single edge), next half-edge on the boundary of the incident face, the previous edge, the face to its left (AKA incident face)
- All of this information can be represented with functions (e.g. IncidentFace(v))
- Available operations:
- Walk around the boundary of a given face
- Access a face from an adjacent one
- Visit all edges around a vertex
- With n vertices, there are <= 3n-6 edges and <= 2n-4 faces
Computing Map Overlay
- Input: Two DCELs; Output: A DCEL that is a combination of both
- Can use the plane sweep line intersection algorithm to construct this DCEL by finding new intersections (vertices)
- Can use some algorithms for boolean operations: union, intersection, difference of DCELs
- Runs in O((n + k)log n), where n is the total size of the two DCELS and k is the number of intersections
Polygon Triangulation (Art Gallery Guarding Problem)
- Given a simple polygon P with n vertices, place a number of guards/cameras on vertices such that every point in P can be seen by some guard/camera
- NP-hard to find minimum number of guards (NP-complete for decision problem)
- Simple polygon: A polygon that does not intersect itself and has no holes
- Guard using triagnulations
- Triangulate: Get diagonals between all vertices
- Triangulation: Maximal set of non-crossing diagonals
- Creates many triangles inside the polygon
- Simple solution: Place one guard per triangle, either in the middle of it or on a vertex
-
Theorem: Every simple polygon has a triangulation with n-2 triangles
- First bound: n-2 guards necessary
- 3-coloring: assign each vertex a color such that adjacent vertices have different colors
-
Lemma: 3-coloring is possible for vertices in a triangulation
- Create a dual graph: vertex for each triangle, edge between traingles that share an edge
- Max degree of 3 (3 sides), acylcic (tree)
-
Lemma: 3-coloring is possible for vertices in a triangulation
-
Theorem: You need the floor of n/3 guards to guard the polygon
- Derived from 3-coloring: since each triangle has one of each color, you can choose 1 color and place guards on those vertices
Triangulating a Polygon
- O(n log n) algorithm
- Split polygon into monotone polygons, O(n log n)
- Triangulate each monotone polygon, O(n)
- A polygon is monotone with respect to a line l iff for every line l’ perpendicular to l is connected with P
- x-monotone, y-monotone: monotone with respect to the x or y axis
- Can check if triangle is monotone in O(n): check upper and lower half, see if x-coordinates are increasing in values
Plane Sweep Algorithm
- Triangulate in O(n) time, with everything to the left of the sweep line being triangulated
- Use left/right turns to determine when to draw triangles
Creating Monotone Polygons
- From each vertex, shoot a ray up/down until it hits a boundary
- AKA trapezoidal decomposition
- Takes O(n log n) time to sweep n events, O(log n) time per event bc of binary tree
- Identify split/merge vertices as ones with two edges going the same direction
- Add diagonals to each by joining it to the polygon vertex of the trapezoid in the opposite direction of the edges
- Both edges going left = split vertex; both going right = merge vertex
Lower Envelope
- Problem: Given n real-valued functions, find the function f(x) that gives the minimum of all n functions; known as the lower envelope
- Other parameters: s is the maximum number of intersection points between any two functions
- Assumptions
- All functions are x-monotone
- Function evaluation and intersection finding takes O(1)
- Functions take constant space
- Lower envelope represented by list-like object (array, linked list); record the current function being used and the range of x-values for which it’s being used
D&C Algorithm
- Divide the set S of n functions in half into S1 and S2
- Compute the lower envelopes L1 and L2 for S1 and S2, respectively
- Use a sweep-line algorithm to merge the lower envelopes
Sweep Line Merge
- Event points: All vertices of L1 and L2, all intersection points
- Queue
- Next (right) endpoint of L1
- Next (right) endpoint of L2
- Intersection point (if it exists)
- Sweep status: y-order of L1 and L2
- Runtime is O(λs(n) log n)
- λs has an upper bound of n log* n
- log* n is defined as the number of times you have to log a number for it to be less than 1
- Inverse Ackermann function: slowest growing function possible
- λ3(n) has an average runtime of n * α(n), where α(n) is the Inverse Ackermann function
Linear Programming
- Problem: Optimize an objective function subject to some constraints
- Linear program: Defined by d variables and n constraints
- Each constraint can be considered a half-space in Rd; e.g. the 2D case with lines
- The set of points satisfying all constraints is $\cap_{i=1}^{n} h_i$
- Want to maximize $f_{\vec{c}} = \vec{x} \cdot \vec{c}$, i.e. finding a point $\vec{x}$ that is extreme in the direction of $\vec{c}$
Half-Plane Intersection
- Given a set of n half-planes in 2D, return the set of points in the interesction of the set
- The output will be a convex polygon with at most n edges
- Can use the same method employed in the lower envelope algorithm; divide and conquer
- O(n log n) runtime
Incremental Linear Programming
- Assumes that the linear program is bounded and that there’s a unique solution
- Add half-planes one by one and keep track of the optimal solution
- Definitions: $H_i = [h_1, h_2, …, h_3], C_n = \cap_{i=1}^{n} h_i, v_i = \text{ Optimal vertex for } C_i$
- Lemma
- If $v_{i-1} \in h_i$, then $v_i = v_{i-1}$
- If $v_{i-1} \notin h_i$, then $C_i = \emptyset$ or $v_i \in l_i$, the line bounding $h_i$
- Takes O(n2) time in worst case because of the case where $v_{i-1} \notin h_i$
- Takes O(n) time on average when randomizing the input
Seidel’s Algorithm
- Randomly permute constraints
- Choose an infinity for an optimal solution
- For each constraint:
- Check if solution point obeys new constraint
- If not, solve the new constraint’s (d-1)-dimensional linear program recursively and update optimal solution
- Exepected runtime is O(d!n) which is O(n) in lower dimensions
Orthogonal Range Searching
- Given n points in d dimensions, query an axis algined box
- Can report points in box, the number inside the box, etc.
- Goal: preprocess this list into a data structure that supports fast queries
1D Case
- Naive: store points in an array and query by binary searching on endpoints
- O(k + log n) runtime to list all points in query
- Better: use a balanced binary search tree
- Store points in the leaves of tree
- Internal nodes store copies of leaves for binary search
- key[x] is the maximum key of any leaf in the left subtree
- Faster than binary searching because you can “share” the two binary search queries up until a split node
- Can make counting queries (i.e. how many points lie in this range) faster by storing number of nodes in subtree for each node
- BST: Space is O(n), preprocessing is O(n log n)
Naive 2D Case
- Create a 1D binary tree for x-coordinates, and for every node, store another 1D tree for y-coordinates
- Query both trees to find all points
- Query time: $O(log^2n)$ or $O(k + log^2n)$ for reporting list of points
- Preprocessing: $O(n log n)$
- Space: O(n log n)
Better 2D: Fractional Cascading
- Wikipedia link
- Results in O(log n) for counting queries or O(k + log n) for listing queries
d-Dimension Case
- For each node in the d-dimension tree, add a tree for the (d-1)-dimension
- In the 2D case, use fractional cascading
- Query time is $O(log^{d-1}n)$ or $O(log^{d-1}n + k)$ depending on query type
- log factor is saved by fractional cascading; without it would be $O(log^d n)$
- Preprocessing: $O(n log^{d-1} n)$
- Space: $O(log^{d-1}n)$
KD tree
- Split odd layers by x-coordinate and even layers by y-coordinate in the 2D case
- Higher dimensions will use modulo to determine which coordinate to discriminate on
- Similar to a decision tree
- O(n log n) preprocessing time
- $O(n^{1 - \frac{1}{d}})$ or $O(n^{1 - \frac{1}{d}} + k)$ query time
Point Location
- An undirected graph is considered planar if it can be embedded in the plane without edge crossings
- A planar embedding induces a planar subdivision; represented by DCEL
- Edges, faces, and vertices are all linear compared to each other; O(n)
- Problem: Given a planar subdivision, find which face a point p lies in
- Naive: Use DCEL to check each face in O(n) time
Slab Method
- Draw a vertical line through each point, dividing the subdivision into “slabs”
- Binary search for the slab p lies in and binary search the slab to determine the face
- $O(n^2 log n)$ preprocessing, $O(n^2)$ space, $O(log n)$ query
Plane Sweep
- Events: Intersections and endpoints
- Status: Ordering of lines, lines are removed/added/reordered based on events
- Runtime is O(n log n), can be used to preprocess the DCEL
-
Persistent data structures allow for queries on past versions of themselves
- BSTs do this via path copying; never overwrite old nodes, simply create new copies
- Path copying for persistence maintains the O(n log n) runtime but extends space to O(n log n)
- Fat nodes, as defined in the Sarnak and Tarjan paper, reduce space to O(n)
- Point location is done in O(log n)
Kirkpatrick’s Point Location Algorithm
- Takes traingulation as input
- Can convert a planar subdivision to a triangulation in linear time
- Triangulate each face, maintaining label, in O(n) time
- Turn outer face into triangle by drawing a triangle around the convex hull
- Size is O(n), conversion is O(n)
Kirkpatrick’s Hierarchy
- Compute a sequence $T_0, …, T_k$ of increasingly coarse triangulations
- $T_0$ is the input triangulation, $T_k$ is the outer triangle
- k is O(log n) and each triangle in $T_{i+1}$ intersects a constant number of triangles in $T_i$
- Build this sequence by deleting vertices and re-triangulating
- Delete a constant fraction of vertices O(log n) times
- Only delete vertices with a constant degree (≤8)
- Choose an independent set of vertices to delete; don’t delete adjacent ones
- O(n) run time
- Removes n/18 vertices each time; $log_{18/17} n$ steps
- Query: Store hierarchy as DAG (or tree) and path find to the corresponding triangle in $T_0$
- Complete algorithm has O(n) preprocessing, O(n) space, and O(log n) query but with large constant factors
- O(216 log n), O(18n)
Voronoi Diagrams
- Given a set of point sites $P = {p_1, …, p_n} \in \R^2$
- Voronoi cell: $V(p_i) = {q \in \R^2 \vert d(p_i, q) < d(p_j, q) \text{ for all } j \neq i}$
- Voronoi diagram: The planar division derived from removing all V cells
- Voronoi edges are bisectors of points; equidistant from two points p and q
-
Voronoi vertex: the intersection of three edges
- Guaranteed to be the center of an empty disk defined by three sites
- Voronoi cells are convex, as they are the intersection of half planes
- Cells have at most n-1 cells
- Number of vertices ≤ 2n - 5, Number of edges ≤ 3n - 6
- A Voronoi cell $V(p_i)$ will be unbounded if and only if $p_i$ is on the convex hull of P
- Can use non-Euclidean distance formulas
Fortune’s Sweep
- Status: “beach line”
- Can be thought of as the lower envelope of parabolas above the sweep line
- Parabolas are equidistant from their origin site and line
- Intersections of parabolas will be on a Voronoi edge
- Can be thought of as the lower envelope of parabolas above the sweep line
- Events
- New point (Site Event)
- A new parabola must be drawn for the new point
- Create a new Voronoi edge by tracing the intersection of parabolas
-
Circle Event
- An arc shrinks to a point and disappears; “eaten up” by two other parabolas
- Defines a Voronoi vertex
- New point (Site Event)
- Voronoi diagram will be stored in a DCEL
- Status will be stores in a BST
- The tree’s leaves store sitres that define the arcs above the beach line
- Internal nodes refer to break points; intersections of arcs
- Event queue: Priority queue ordered by y-coordinate
- New sites are stored as site events
- Circle events
- Store lowest point of circle
- Store pointer to the leaf/arc in the tree that will disappear
- Size: O(n) because there are n sites and O(n) circles
- Detecting circle events
- Store the circles created by any three consecutive points
- Calculate point on circle
- Handling site events
- Locate parabola directly above new point
- If a circle event related to the parabola exists, delete it
- Replace the prior leaf node corresponding to the above parabola with a three-subtree containing the new point’s leaf node in the middle and leaf nodes for the old point on the outside
- Add new half-edges separated by the new breakpoints to the DCEL
- Create new circle events with the newly consecutive points and add them if they are below the sweep line
- Handling circle events
- Delete the point corresponding to the disappearing arc
- Add intersection point to Voronoi diagram and add half-edges for the breakpoints
- Check newly consecutive parabolas for circle events
- Total runtime: O(n log n)
- Insert and delete is O(log n), O(n) events
Point Line Duality
- Let $P = {p_1, …, p_n} \in \R ^2$ be a set of points
- $P* = {p_1^, …, p_n^}$ is defined as the dual set of lines
- Primal: Point p: (px, py), Line l: y = mx + b
- Dual: Line p* = y = pxx - py, Point l* = (m, -b)
- Properties
- (p) = p
- $p \in l \iff l^* \in p^*$ (Incidence preservation)
- $\text{p lies above l} \iff l^* \text{ lies above } p^*$ (Order preservation)
- From the properties, it follows that the intersection of lines in the dual plane defines the line between two points in the primal plane
- Used for mapping algorithms between planes
- e.g. the lower convex hull in primal is the same as the upper envelope in dual
- Conversion takes O(n) time
Line Arrangements
- Let L be a set of lines; then, A(L) is the arrangement of L
- Defined as the planar subdivision induced by all lines in L; represented by DCEL
- Vertices <= n choose 2 (O(n2)), edges <= n2, faces <= 1 + n(n+1)/2
- Using a sweep line algorithm would take O(n2 log n) time
- Incremental construction takes this down to O(n2)
- Add lines one at a time
- Each line is guaranteed to intersect/split O(i) faces by the Zone Theorem
- O(i) < O(n)
Delaunay Triangulation
- The triangulation of a set of point is guaranteed to have an outer face bounded by the convex hull
- The dual graph G* has a node for each face of G and an edge connecting nodes if and only if the faces share an edge
- G would be a planar subdivision
- The Delaunay triangulation is the dual graph of a Voronoi diagram for a set of points
- Making each edge straight in the Delaunay triangulation turns it into a triangulation
- Only occurs if P is in general position
- Can be stored as an abstract graph; less space than a DCEL
- Straight line edges don’t intersect due to the fact that an edge exists in the Delaunay graph if and only if an empty disk exists that has both points on its boundary
- For any three points in P:
- Their triangle is in the Delaunay triangulation
- The center of the circle is a vertex in the Voronoi diagram
- Equivalent statements
- Can be used to compute the convex hull of a 3D set of points
3D Convex Hull
Gift Wrapping
- Find a face guaranteed to be on the convex hull
- Choose an edge on this face whose other face has not been found
- Find the point that creates the minimum angle with the current face and add it to the convex hull
- Note that the new point creates a triangle when joined with the original edge - Runtime is O(nF) where F is the number of faces
Divide and Conquer
- Same paradigm as 2D; divide points into two sets and merge the calculated convex hulls
- Merge step
- Option 1: Find each new face by projecting the set of points into 2D and finding the new hull using the same 2D binary search method
- Option 2: Find the bottom plane that merges the two hulls, and gift wrap it around
- Since you can ignore points if they don’t belong on the hull when processed, this takes O(n) time
- Runtime: O(n log n)
Windowing
- Input: A set S of line segments
- Query: Report all segments that intersect a query window
- Subproblem (1D): Report all intervals that contain a query point
Interval Trees
- Input: A set of I intervals on a line
- Uses D&C paradigm: Partition I at each step into Ileft, Imid, and Iright
- Imid contains intervals containing the median of the 2n endpoints
- Imid stores intervals in two lists; one sorted by left endpoint, and one stored by right endpoint
- O(n) space (each interval stored in two lists), O(n log n) preprocessing (Master’s theorem), O(log n + k) query
- Each node takes O(1 + kv) time, as we must report the kv intervals hit at each node
Priority Search Trees
- Allows for 3-sided range queries; find all points in a range where one of the sides is unbounded
- e.g. two y constraints and one x constraint
- Acts like a heap on x and a BST on y
- Each vertex stores the minimum x-coordinate of its set of points and the median y-coordinate
- For children nodes: split the set of points by median y-coordinate, and then intialize them with the point with the minimum x-coordinate
- O(log n + k) 3-sided query, O(n) space, built in O(n log n)
Segment Tree
- Partition line segments by all endpoints
- $\text{Line segments} = {(-\infty, p_1), [p_1, p_1], (p_1, p_2), …, (p_n , \infty)}$
- Store intervals in leaves of tree
- Parent nodes are unions of their children’s intervals
- O(n) space; each interval can be stored in a maximum of two nodes
- O(log n + k) query, O(n log n) preprocess
- Built bottom up from the leaf nodes
Solving the Windowing Problem
- Create slabs using the segment tree and search for the slabs contained by the query
- In each slab, use a 1D range query (taken from point location) to detect segments that intersect
- O(log2 n + k) query; query the segment tree in O(log n) time, and for each of those queries, query the 1D range tree in O(log n) time
- O(n log n) preprocessing
- O(n log n) space because each node contains a 1D range tree