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
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