Computational Geometry

Chapter 1: Cross Products & Segment Intersection

Before we can detect intersections, build convex hulls, or find closest pairs, we need a single primitive operation: given three points p₁, p₂, p₃, is the path p₁ → p₂ → p₃ a left turn, a right turn, or collinear?

This is the orientation test, and it is built on the cross product of two 2D vectors.

The Cross Product in 2D

Given two vectors a = (a_x, a_y) and b = (b_x, b_y), the 2D cross product (really the z-component of the 3D cross product) is:

This single number tells you everything about the spatial relationship between the two vectors:

The magnitude |a × b| equals the area of the parallelogram formed by a and b. Half of that is the area of the triangle. So the cross product simultaneously tells you orientation AND area.

Why Does This Formula Work?

The full 3D cross product of a = (a_x, a_y, 0) and b = (b_x, b_y, 0) is the vector (0, 0, a_xb_y − a_yb_x). The z-component is the only non-zero part. Its sign tells you whether b is counterclockwise from a (positive z = counterclockwise when viewed from above) or clockwise (negative z).

Another way to see it: the cross product is the determinant of a 2×2 matrix:

Determinants measure "signed area." When the two vectors span a parallelogram that opens counterclockwise, the area is positive. When clockwise, negative. When the vectors are parallel (zero area), the determinant is zero. This is the geometric content of the cross product.

The Area Connection

The triangle formed by points p₁, p₂, p₃ has signed area:

Positive area = counterclockwise (left turn). Negative area = clockwise (right turn). Zero area = collinear (degenerate triangle). This connects orientation testing to area computation: they are the same operation.

The Orientation Test

To determine the turn direction at p₁ → p₂ → p₃, we compute the cross product of the vectors (p₂ − p₁) and (p₃ − p₁):

If d > 0, the triple (p₁, p₂, p₃) makes a left turn (counterclockwise). If d < 0, it is a right turn (clockwise). If d = 0, the three points are collinear.

Worked Example

Let p₁ = (0, 0), p₂ = (4, 0), p₃ = (4, 3). What is the orientation?

Makes sense: from (0,0) we go right to (4,0), then up to (4,3). That is a left (counterclockwise) turn.

Now try p₃ = (4, −3):

From (0,0) right to (4,0), then down to (4,−3). A right turn. Correct.

Segment Intersection Test

Now the payoff. Two line segments s₁ = (p₁, p₂) and s₂ = (p₃, p₄) intersect if and only if one of these conditions holds:

General case: The endpoints of each segment straddle the line containing the other segment. Formally:

"Opposite signs" means the two endpoints of s₁ are on different sides of the line through s₂, and vice versa. If both conditions hold, the segments must cross.

Implementation

python
def cross(o, a, b):
    """Cross product of vectors (a - o) and (b - o)."""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def on_segment(p, q, r):
    """Check if point q lies on segment pr, given that p, q, r are collinear."""
    return (min(p[0], r[0]) <= q[0] <= max(p[0], r[0]) and
            min(p[1], r[1]) <= q[1] <= max(p[1], r[1]))

def segments_intersect(p1, p2, p3, p4):
    """Do segments (p1,p2) and (p3,p4) intersect?"""
    d1 = cross(p3, p4, p1)
    d2 = cross(p3, p4, p2)
    d3 = cross(p1, p2, p3)
    d4 = cross(p1, p2, p4)

    # General case: endpoints straddle each other's line
    if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
       ((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0)):
        return True

    # Collinear boundary cases
    if d1 == 0 and on_segment(p3, p1, p4): return True
    if d2 == 0 and on_segment(p3, p2, p4): return True
    if d3 == 0 and on_segment(p1, p3, p2): return True
    if d4 == 0 and on_segment(p1, p4, p2): return True

    return False

This is O(1) per test. Four cross products, at most four boundary checks. Every single computational geometry algorithm in this chapter calls this function (or the orientation test it depends on) as a subroutine.

Worked Example: Segment Intersection

Do segments (1,1)-(5,5) and (1,5)-(5,1) intersect? These form an "X" shape.

Now try parallel segments: (0,0)-(4,0) and (0,2)-(4,2):

Both endpoints of segment 1 are on the same side of segment 2's line. They can never cross.

The On-Segment Test in Detail

When a cross product is zero, we know three points are collinear. But that does not mean the point is on the segment — it might be on the line extension beyond the endpoints. The on-segment test checks whether the point's coordinates fall within the bounding box of the segment:

This handles the case where segments touch at endpoints (T-junctions), overlap along a line, or one segment's endpoint lies exactly on the other segment. Without this test, collinear touching segments are missed.

Interactive: Drag Segments

The canvas below shows two segments. Drag the endpoints to move them. The display shows the four cross-product values (d₁ through d₄) and whether the segments intersect. Watch how the signs flip as endpoints cross the other segment's line.

Chapter 2: Sweep Line

You have n line segments in the plane. You want to know: do any two segments intersect? The brute-force approach checks all &binom;(n,2) = O(n²) pairs. For n = 10,000, that is 50 million intersection tests. Can we do better?

Yes. The sweep line technique processes the plane from left to right, maintaining a vertical line that sweeps across all segments. At any moment, the sweep line intersects a subset of the segments. The key insight: two segments can only intersect if they are adjacent in the vertical ordering at the sweep line's current position.

How the Sweep Works

Think of a vertical ruler sliding from left to right across the plane. As it moves, it "enters" segments (at their left endpoint) and "exits" segments (at their right endpoint). At any x-coordinate, the ruler intersects some subset of segments, and these segments appear in some vertical order on the ruler.

Why O(n log n)?

There are 2n events (one per endpoint). Each event involves one BST operation (insert or delete) which takes O(log n). Each event also involves at most two intersection tests (checking neighbors), each O(1). Total: O(n log n) for sorting + O(n log n) for BST operations = O(n log n).

The Key Invariant

If two segments intersect, there must exist some sweep-line position where those two segments are vertically adjacent in the BST. Why? Before the intersection point, one segment is above the other. After the intersection, their order reverses. At the exact intersection, they swap — and to swap, they must be adjacent. The algorithm catches this swap either when one of them is inserted (and becomes a neighbor of the other) or when an intervening segment is deleted (making them newly adjacent).

Worked Example: Sweep Line Trace

Three segments: s₀ = (1,4)-(6,1), s₁ = (2,1)-(5,5), s₂ = (3,3)-(7,4).

This example shows the sweep line's early termination advantage. In the best case, it finds an intersection after just a few events. Brute force always checks all O(n²) pairs regardless.

Animated Sweep Line

The canvas below shows a set of line segments. Press "Sweep" to watch the vertical sweep line process events from left to right. When it detects an intersection, it highlights the intersecting pair.

Pseudocode

Implementation: Sweep Line

python
from sortedcontainers import SortedList

def any_segments_intersect(segments):
    """Sweep line: detect if any two segments in the set intersect.
    Each segment is ((x1,y1), (x2,y2)) with x1 <= x2."""

    # Build events: (x, type, segment_index)
    # type: 0 = left endpoint (insert), 1 = right endpoint (delete)
    events = []
    for i, ((x1,y1),(x2,y2)) in enumerate(segments):
        if x1 > x2:
            (x1,y1),(x2,y2) = (x2,y2),(x1,y1)
        events.append((x1, 0, i))  # left endpoint
        events.append((x2, 1, i))  # right endpoint
    events.sort()  # O(n log n)

    # Active segments ordered by y-coordinate at current sweep x
    # In practice, use a balanced BST with y-at-sweep-x as key
    active = SortedList(key=lambda idx: y_at_x(segments[idx], sweep_x))

    for x, typ, seg_idx in events:
        sweep_x = x  # update global sweep position

        if typ == 0:  # LEFT endpoint: insert
            active.add(seg_idx)
            pos = active.index(seg_idx)
            # Check neighbors
            if pos > 0 and intersects(segments[active[pos-1]], segments[seg_idx]):
                return True
            if pos < len(active)-1 and intersects(segments[active[pos+1]], segments[seg_idx]):
                return True

        else:  # RIGHT endpoint: delete
            pos = active.index(seg_idx)
            # Check if removing this makes its neighbors adjacent
            if pos > 0 and pos < len(active)-1:
                if intersects(segments[active[pos-1]], segments[active[pos+1]]):
                    return True
            active.remove(seg_idx)

    return False

Bentley-Ottmann: Finding ALL Intersections

The sweep line above only detects whether any intersection exists. The Bentley-Ottmann algorithm (1979) finds all k intersections in O((n + k) log n) time. The key modification: when two segments swap their vertical ordering at an intersection point, we add the intersection as a new event. We then swap the two segments in the BST and check their new neighbors for additional intersections.

For n segments with k intersections, brute force is O(n²). Bentley-Ottmann is O((n + k) log n). When k is small (few intersections), this is near-linear. When k = Θ(n²) (many intersections), it is no better than brute force — but at least it does not waste time when intersections are sparse.

Why Sweep Line Matters Beyond Segments

The sweep line is not just for segment intersection. It is a general paradigm that transforms 2D problems into dynamic 1D problems. Here are other problems solved by sweep line:

In every case, the pattern is identical: sort events, maintain a dynamic structure, update at each event, check only local neighbors. The sweep line converts an unstructured 2D search into a sequence of structured 1D updates.

Chapter 3: Convex Hull

The convex hull of a set of points Q is the smallest convex polygon that contains all points in Q. "Convex" means that for any two points inside the polygon, the line segment between them is also entirely inside. Equivalently, at every vertex the boundary turns in the same direction (all left turns, or all right turns).

Why is this useful? Convex hulls are the foundation for collision detection (if two convex hulls do not intersect, the objects do not collide), shape analysis (the hull gives you the "footprint" of a point cloud), and as a preprocessing step for more complex algorithms (halfplane intersection, Voronoi diagrams).

Graham Scan — O(n log n)

The most elegant convex hull algorithm. Invented by Ronald Graham in 1972. It works in two phases:

Phase 1: Sort by polar angle. Pick the point with the lowest y-coordinate (break ties by x). Call it p₀. Sort all other points by the polar angle they make with p₀. This means if you stand at p₀ and sweep your gaze counterclockwise from the positive x-axis, the sorted points appear in the order you would see them.

Phase 2: Build the hull with a stack. Initialize a stack with p₀, p₁, p₂ (the first three sorted points). For each remaining point p_i:

When all points are processed, the stack contains the convex hull vertices in counterclockwise order.

Why It Works

The polar angle sort guarantees we process points in a counterclockwise sweep around p₀. The stack invariant is: the points on the stack always form a convex chain (all left turns). When we encounter a right turn, the previous point is "inside" the hull being built — it creates a dent — so we remove it. This is greedy and correct because the polar angle ordering ensures we never need to reconsider a point once we have swept past its angle.

Animated Graham Scan

Watch the algorithm in action. The points are first sorted by polar angle from the bottom-most point (shown with dashed sweep lines). Then the stack builds the hull, popping points that cause right turns.

Jarvis March (Gift Wrapping) — O(nh)

An alternative algorithm that is simpler but potentially slower. Start at the leftmost point. At each step, find the point that makes the smallest counterclockwise angle with the current direction. This "wraps" the hull like wrapping a gift with string.

Each wrapping step scans all n points to find the next hull vertex — using the cross product to determine which candidate makes the most counterclockwise turn. Since there are h hull vertices, the total is O(nh). When h is small (e.g., h = O(log n)), Jarvis March beats Graham Scan. When h = O(n) (all points on the hull), it degrades to O(n²).

Worked Example: Graham Scan by Hand

Points: A(0,0), B(4,0), C(4,4), D(2,5), E(0,4), F(2,2), G(1,1), H(3,1). The bottom-most point is A(0,0).

Notice how H was pushed and then immediately popped when C arrived. The cross product caught the right turn — H was "inside" the triangle ABC. This is exactly the bouncer analogy: H tried to join the rope line, but C revealed that H did not belong.

Implementation: Graham Scan

python
def convex_hull(points):
    """Graham Scan. Returns hull vertices in CCW order."""
    # Find bottom-most point (lowest y, then leftmost x)
    p0 = min(points, key=lambda p: (p[1], p[0]))

    # Sort by polar angle from p0
    import math
    def polar_angle(p):
        return math.atan2(p[1] - p0[1], p[0] - p0[0])

    sorted_pts = sorted(points, key=polar_angle)

    # Build hull with stack
    stack = [sorted_pts[0], sorted_pts[1]]

    for i in range(2, len(sorted_pts)):
        # Pop while the last turn is clockwise (right turn)
        while len(stack) > 1 and cross(stack[-2], stack[-1], sorted_pts[i]) <= 0:
            stack.pop()
        stack.append(sorted_pts[i])

    return stack

# Example
pts = [(0,0), (1,1), (2,0), (1,2), (0,2), (2,2), (1,0.5)]
hull = convex_hull(pts)
# hull = [(0,0), (2,0), (2,2), (0,2)] — the four corners
# Interior points (1,1), (1,2), (1,0.5) were popped

Quickhull: The Practical Favorite

Quickhull is to convex hull what Quicksort is to sorting: O(n log n) expected time, O(n²) worst case, but typically the fastest in practice. It is the default algorithm in most geometry libraries (CGAL, scipy.spatial, MATLAB).

The algorithm works recursively:

The "farthest point" at each step is found using — you guessed it — the cross product. The cross product of the line direction and the vector to each candidate point gives the perpendicular distance (proportional to it). Take the maximum.

Why is Quickhull fast in practice? Because most random point sets have very few hull vertices (O(log n) expected for uniform distributions). At each recursion level, the triangle discards a large fraction of remaining points. Only the points "outside" the current hull approximation survive to the next level. This aggressive pruning means the average-case depth is O(log n), giving expected O(n log n).

python
def quickhull(points):
    """Quickhull algorithm. Expected O(n log n)."""
    if len(points) < 2: return points
    left = min(points); right = max(points)
    hull = []
    upper = [p for p in points if cross(left, right, p) > 0]
    lower = [p for p in points if cross(left, right, p) < 0]
    _qh(left, right, upper, hull)
    _qh(right, left, lower, hull)
    return hull

def _qh(a, b, pts, hull):
    if not pts:
        hull.append(a); return
    # Farthest point from line a-b
    far = max(pts, key=lambda p: abs(cross(a, b, p)))
    # Split: outside triangle(a, far) and outside triangle(far, b)
    left_of_af = [p for p in pts if cross(a, far, p) > 0]
    left_of_fb = [p for p in pts if cross(far, b, p) > 0]
    _qh(a, far, left_of_af, hull)
    _qh(far, b, left_of_fb, hull)

Chapter 4: Closest Pair

Given n points in the plane, find the two points that are closest to each other. This is one of the most elegant applications of divide and conquer in all of computer science.

Brute force: check all &binom;(n,2) pairs, take the minimum distance. That is O(n²). For a million points, that is 500 billion distance computations. We can do O(n log n).

The Divide-and-Conquer Strategy

The genius is in the combine step. It seems like checking the strip could be O(n²) — the strip might contain many points. But a remarkable geometric argument shows that each point needs to be compared with at most 7 other points in the strip.

The Strip Argument (Why Only 7 Points)

Consider a point p in the strip. We only need to check points within distance δ of p. These points lie in a rectangle of width 2δ (the strip) and height δ (above p). Divide this rectangle into 8 cells, each δ/2 × δ/2.

Each cell can contain at most one point. Why? If two points were in the same cell, their distance would be at most δ/2 · √2 = δ · √2/2 ≈ 0.707δ, which is less than δ. But δ is the minimum distance within each half — contradiction. So each cell has at most one point, the rectangle has 8 cells, and p itself occupies one. That leaves at most 7 points to check.

Worked Example: Closest Pair by Hand

Points sorted by x: A(1,2), B(2,8), C(3,1), D(5,7), E(6,3), F(7,6), G(8,1), H(9,5).

Notice that the strip check did not find a closer cross-boundary pair. This happens often — the recursive halves usually contain the closest pair. But not always. When the closest pair straddles the dividing line, the strip check catches it, and the 7-point bound keeps the cost linear.

Implementation Detail: Sorting by Y

To efficiently find which points are within δ of each other in the strip, we sort the strip points by y-coordinate. Then for each point, we scan forward in the sorted list until the y-distance exceeds δ. This scanning never goes more than 7 steps, so it is O(1) per point amortized.

The total sort takes O(n log n). But we only sort once at the beginning if we maintain the y-sorted order through the recursion (merge the y-sorted halves back together, like merge sort).

Animated Closest Pair

Watch the divide-and-conquer algorithm split the point set, recurse on each half, and then check the strip. The closest pair is highlighted at the end.

Implementation

python
import math

def dist(p, q):
    return math.sqrt((p[0]-q[0])**2 + (p[1]-q[1])**2)

def brute_force(pts):
    """O(n^2) base case for small inputs."""
    best = float('inf')
    pair = None
    for i in range(len(pts)):
        for j in range(i+1, len(pts)):
            d = dist(pts[i], pts[j])
            if d < best:
                best, pair = d, (pts[i], pts[j])
    return best, pair

def closest_pair(pts):
    """O(n log n) divide and conquer."""
    pts.sort()  # sort by x
    return _closest(pts)

def _closest(pts):
    n = len(pts)
    if n <= 3:
        return brute_force(pts)

    mid = n // 2
    mid_x = pts[mid][0]

    # Recurse on each half
    d_left, pair_left = _closest(pts[:mid])
    d_right, pair_right = _closest(pts[mid:])

    # delta = min distance from each half
    if d_left < d_right:
        delta, best_pair = d_left, pair_left
    else:
        delta, best_pair = d_right, pair_right

    # Build strip: points within delta of the dividing line
    strip = [p for p in pts if abs(p[0] - mid_x) < delta]
    strip.sort(key=lambda p: p[1])  # sort by y

    # Check strip pairs — at most 7 comparisons per point
    for i in range(len(strip)):
        j = i + 1
        while j < len(strip) and strip[j][1] - strip[i][1] < delta:
            d = dist(strip[i], strip[j])
            if d < delta:
                delta, best_pair = d, (strip[i], strip[j])
            j += 1

    return delta, best_pair

Chapter 5: Geometry Showcase

Time to put it all together. The canvas below is your computational geometry workbench. Click to add points, then run any of the three core algorithms. Watch convex hull wrap around your points, closest pair highlight the nearest neighbors, or test line segment intersections.

This is the payoff simulation. Every algorithm from the previous chapters runs live on your point set.

Combining Algorithms

In real applications, you rarely use just one algorithm. Consider this pipeline for a 2D robot planner:

Each step uses a different algorithm from this chapter. The total pipeline runs in O(n log n) where n is the total number of obstacle vertices. Without these algorithms, the planner would need O(n²) brute-force checks at every step.

Algorithm Complexity Comparison

Here is a complete comparison of all methods for each of the three core problems:

Every single improvement from brute force to optimal uses one of three paradigms: sorting, sweep line, or divide and conquer. These three paradigms, combined with the cross product primitive, are the entire toolkit of 2D computational geometry.

Numeric Precision: The Hidden Enemy

All algorithms in this chapter assume exact arithmetic. In practice, floating-point numbers are not exact. A cross product that should be exactly zero might compute to 1.2 × 10^-15. This causes subtle bugs:

Chapter 6: Voronoi Diagrams & Delaunay Triangulation

You have n cell towers scattered across a city. For any location, which tower provides the strongest signal? Assuming signal strength decreases with distance, the answer is: the nearest tower. The Voronoi diagram divides the entire plane into regions, one per tower, such that every point in a region is closer to its tower than to any other.

Voronoi Diagram

Given n points (called sites), the Voronoi diagram partitions the plane into n Voronoi cells. The cell for site s_i contains all points p such that dist(p, s_i) ≤ dist(p, s_j) for all j ≠ i. The boundary between two adjacent cells is a segment of the perpendicular bisector of the line connecting the two sites.

Key properties:

Fortune's Algorithm — O(n log n)

Fortune's sweep algorithm (1987) constructs the Voronoi diagram using a sweep line — but with a twist. Instead of a vertical line, it uses a beach line: a chain of parabolic arcs that advances across the plane.

Each site that the sweep line has already passed contributes a parabola to the beach line. The parabola for site s is the set of points equidistant from s and the sweep line. As the sweep line advances, these parabolas grow and merge. The intersections of adjacent parabolas trace out the Voronoi edges.

Two types of events drive the algorithm:

Site events: The sweep line reaches a new site. A new parabolic arc appears on the beach line, potentially splitting an existing arc. This creates a new Voronoi edge.

Circle events: Three adjacent arcs on the beach line converge. The middle arc shrinks to zero width and disappears. The point where it vanishes becomes a Voronoi vertex (equidistant from three sites). Two Voronoi edges end here, and a new one begins.

Fortune's Algorithm: Data Structures

Fortune's algorithm maintains two data structures:

Beach line (balanced BST): The beach line is stored as a balanced BST of arcs. Each internal node represents a breakpoint between two adjacent arcs — this breakpoint traces out a Voronoi edge as the sweep advances. Leaves represent the arcs themselves (each associated with a site). Insert, delete, and search are all O(log n).

Event queue (priority queue): Contains site events (sorted by y-coordinate) and circle events (predicted disappearances of arcs). Each circle event is associated with three consecutive arcs on the beach line. When the sweep line reaches the bottom of the circumscribed circle of those three arcs' sites, the middle arc disappears. Priority queue operations are O(log n).

The algorithm processes at most n site events and at most 2n − 5 circle events (since there are at most that many Voronoi vertices). Each event involves O(log n) BST and priority queue operations. Total: O(n log n).

Voronoi Diagram Visualization

Delaunay Triangulation

The Delaunay triangulation is the dual of the Voronoi diagram. Connect two sites with an edge if and only if their Voronoi cells share an edge. The result is a triangulation of the point set with a remarkable property: no point lies inside the circumcircle of any triangle.

This "empty circumcircle" property means Delaunay triangulations maximize the minimum angle among all possible triangulations. In other words, Delaunay avoids skinny, degenerate triangles — making it ideal for finite element mesh generation, terrain modeling, and computational physics.

Applications in ML and Robotics

k-Nearest Neighbors: Voronoi cells directly encode 1-NN classification boundaries. For k-NN, higher-order Voronoi diagrams partition the plane by the k closest sites. The decision boundary between two classes is exactly the perpendicular bisector of two nearby training points — which is a Voronoi edge.

Navigation meshes: Games and robots use Delaunay triangulation to create a walkable mesh. The robot plans paths through triangle edges, avoiding obstacles that lie outside the mesh. The Delaunay property (maximizing minimum angle) ensures triangles are "fat" rather than skinny, which improves path quality and numerical stability.

Sensor placement: To cover an area with minimum sensors, place sensors at the Voronoi vertices of existing coverage gaps. Each sensor covers its Voronoi cell. Iterative algorithms move sensors toward their cell centroids (Lloyd's algorithm), converging to a locally optimal placement.

Natural neighbor interpolation: Given scattered data points with values (e.g., temperature measurements), interpolate at a new location using the Voronoi diagram. The "natural neighbors" of the query point are the sites whose Voronoi cells would shrink if the query point were inserted as a new site. Weight each neighbor by how much its cell shrinks. This gives smooth, data-driven interpolation without choosing a kernel or bandwidth.

Chapter 7: Applications

Computational geometry is not abstract mathematics locked in textbooks. It runs in every robot, every game engine, every mapping application, and increasingly in machine learning pipelines. Here is where the algorithms from this chapter show up in the real world.

Robotics: Path Planning & Collision Detection

A robot arm with 6 joints has a configuration space (C-space) with 6 dimensions. Each obstacle in physical space maps to a forbidden region in C-space. The robot needs to find a path from start configuration to goal configuration that avoids all forbidden regions.

Convex hull is used to approximate obstacle boundaries. Convex shapes are easier to check for intersection than arbitrary shapes — the GJK algorithm (Gilbert-Johnson-Keerthi, 1988) exploits this to do collision detection in O(n) per convex pair using support functions.

The key insight of GJK: two convex shapes overlap if and only if their Minkowski difference contains the origin. The Minkowski difference A − B = {a − b : a ∈ A, b ∈ B} is also convex. GJK iteratively builds a simplex (point, line, triangle) inside the Minkowski difference, always picking the point farthest toward the origin (via the support function). If the simplex encloses the origin — collision. If it cannot get closer to the origin — no collision. Typically converges in 3-4 iterations.

This is why convex hull computation matters in robotics: approximate each obstacle with its convex hull, then GJK handles all collision queries in near-constant time per pair. For non-convex objects, decompose into convex parts first.

Line segment intersection determines whether a proposed trajectory (a line in C-space) passes through an obstacle. Sweep line algorithms check thousands of potential trajectories efficiently.

Voronoi diagrams in the free space give the robot maximum-clearance paths — following Voronoi edges keeps the robot as far from obstacles as possible.

Computer Graphics: Rendering & Clipping

Sutherland-Hodgman clipping: To render a polygon on screen, you clip it against the viewport boundaries. Each clip against one edge is a cross-product test: is each vertex inside or outside the edge? The algorithm outputs a new polygon that is the intersection of the original polygon with the half-plane defined by the clip edge. Four clips (one per viewport edge) render any convex polygon correctly.

Triangulation for rasterization: GPUs render triangles. Every polygon must be decomposed into triangles. Delaunay triangulation produces the most numerically stable decomposition. A convex polygon with n vertices can be triangulated in O(n) (fan from one vertex); a general polygon requires O(n log n) via ear clipping or sweep line.

GIS: Map Overlay & Spatial Queries

Geographic Information Systems layer map data: roads, buildings, rivers, administrative boundaries. Map overlay computes the intersection of two polygon layers — which buildings are in which school district? This is a sweep line problem over the edges of both polygon sets.

Point-in-polygon tests answer "which country is this GPS coordinate in?" The classic algorithm casts a ray from the point and counts crossings with the polygon boundary. Odd crossings = inside. Even = outside. Each crossing test uses the cross product.

Machine Learning: Geometry of Optimization

Support Vector Machines find the maximum-margin hyperplane separating two classes. The margin is determined by the convex hulls of the two classes — the optimal hyperplane touches both hulls. The support vectors are points on the hull boundary.

Linear programming: The feasible region is a convex polytope (higher-dimensional convex hull). The simplex algorithm walks along edges of this polytope from vertex to vertex.

Gradient descent: Projected gradient descent onto a convex set requires computing the projection of a point onto the nearest point of a convex hull — a direct geometric operation. For constrained optimization with linear constraints, each projection step reduces to finding which facet of the convex polytope is nearest — the same support function used in GJK.

Dimensionality reduction: PCA finds the directions of maximum variance, which correspond to the "widest" directions of the data's convex hull. Random projection (Johnson-Lindenstrauss lemma) preserves distances — including the closest pair distance — with high probability. The geometric structure of high-dimensional convex bodies underpins much of theoretical ML.

Game Physics: Broad & Narrow Phase

Modern physics engines use a two-phase collision detection system:

The Shoelace Formula: Area of Any Simple Polygon

Given a polygon with vertices (x₁,y₁), (x₂,y₂), ..., (x_n,y_n) listed in order, the shoelace formula computes the signed area:

where indices wrap around (x_n+1 = x₁, y_n+1 = y₁). The name comes from the criss-cross pattern of multiplication — like lacing a shoe.

Worked example: Triangle with vertices (0,0), (4,0), (2,3):

Check: base 4, height 3, area = ½·4·3 = 6. Correct.

The shoelace formula is really just the sum of cross products of consecutive edge vectors. The signed version (without absolute value) gives positive area for counterclockwise vertices and negative for clockwise — which tells you the polygon's orientation. This is used in winding-number algorithms for point-in-polygon tests.

python
def polygon_area(vertices):
    """Shoelace formula. Vertices in order (CW or CCW)."""
    n = len(vertices)
    area = 0
    for i in range(n):
        j = (i + 1) % n
        area += vertices[i][0] * vertices[j][1]
        area -= vertices[j][0] * vertices[i][1]
    return abs(area) / 2

Interactive: Point-in-Polygon Tester

Draw a polygon by clicking vertices (close it by clicking near the first vertex). Then click inside the canvas to test whether the point is inside or outside the polygon. The ray-casting algorithm runs live.

python
def point_in_polygon(point, polygon):
    """Ray casting algorithm. Cast ray to the right, count crossings."""
    x, y = point
    n = len(polygon)
    inside = False

    j = n - 1
    for i in range(n):
        xi, yi = polygon[i]
        xj, yj = polygon[j]

        # Does the edge (i,j) cross the horizontal ray from (x,y) rightward?
        if ((yi > y) != (yj > y)) and \
           (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
            inside = not inside
        j = i

    return inside

Chapter 8: Interview Arsenal

Computational geometry appears in interviews at robotics companies (Waymo, Boston Dynamics, NVIDIA), game studios (Epic, Unity), mapping companies (Google Maps, Mapbox), and any company doing spatial computing (Apple Vision Pro). Here is your cheat sheet and drill set.

Complexity Cheat Sheet

The Five Interview Patterns

Pattern 1: Orientation / Cross Product. "Given three points, determine if they are clockwise, counterclockwise, or collinear." Compute the cross product. This is the foundation — if you cannot do this in your sleep, nothing else works.

Pattern 2: Convex Hull as Preprocessing. "Find the pair of points with maximum distance." The answer must be on the convex hull (rotating calipers). "Is this set of points convex?" Check if the convex hull has n vertices. "What is the smallest enclosing circle?" Run Welzl's algorithm on the convex hull vertices, not all points.

Pattern 3: Sweep Line Reduction. "Find all intersections among n rectangles." Convert to a sweep line problem: each rectangle has a left-edge event and a right-edge event. Use an interval tree for the active set. "Find the closest pair of points." The 1D version is trivial (sort, check adjacent). The 2D version uses D&C with a strip that is essentially a 1D problem.

Pattern 4: Point-in-Region. "Is a point inside a polygon?" Ray casting (O(n)). "Is a point inside a convex polygon?" Binary search on the fan triangulation from one vertex (O(log n)). "Which of n regions contains this point?" Build a Voronoi diagram or use a point location data structure.

Pattern 5: Geometric D&C. "Find the closest pair of points." Divide by x, recurse, merge with strip check. "Find the convex hull in expected O(n)." Quickhull: pick two extreme points, recursively find the farthest point from the line. D&C reduces 2D geometry to simpler 1D subproblems.

Coding Drill: Implement from Scratch

In interviews, you will be asked to implement these from memory. Practice until you can write each in under 5 minutes:

python
# DRILL 1: Cross product / orientation test
def cross(o, a, b):
    return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])

# DRILL 2: Convex hull (Graham Scan, simplified)
def convex_hull(pts):
    pts = sorted(set(pts))  # sort by (x, y), deduplicate
    if len(pts) <= 1: return pts
    # Build lower hull
    lower = []
    for p in pts:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    # Build upper hull
    upper = []
    for p in reversed(pts):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return lower[:-1] + upper[:-1]

# DRILL 3: Point in polygon (ray casting)
def inside(pt, poly):
    x, y = pt
    n, hit = len(poly), False
    j = n - 1
    for i in range(n):
        yi, yj = poly[i][1], poly[j][1]
        if (yi > y) != (yj > y):
            xi, xj = poly[i][0], poly[j][0]
            if x < (xj-xi)*(y-yi)/(yj-yi)+xi:
                hit = not hit
        j = i
    return hit

# DRILL 4: Segment intersection
def seg_intersect(a, b, c, d):
    d1, d2 = cross(c,d,a), cross(c,d,b)
    d3, d4 = cross(a,b,c), cross(a,b,d)
    if ((d1>0)!=(d2>0)) and ((d3>0)!=(d4>0)):
        return True
    # Collinear boundary cases omitted for brevity
    return False

Rotating Calipers: Diameter in O(n)

The diameter of a point set is the maximum distance between any two points. Brute force is O(n²). But the farthest pair must be on the convex hull (a point interior to the hull is always closer to every other point than some hull vertex is). So compute the hull in O(n log n), then find the diameter among the h hull vertices.

Even checking all &binom;(h,2) hull pairs is O(h²). The rotating calipers technique does it in O(h). Imagine placing two parallel lines (calipers) on opposite sides of the hull. Rotate them around the hull, maintaining contact with one vertex on each side. At each rotation step, the farthest pair is the pair of vertices touching the calipers. One full rotation visits all O(h) antipodal pairs.

python
def diameter(hull):
    """Rotating calipers on a convex hull (CCW order). O(n)."""
    n = len(hull)
    if n <= 1: return 0
    if n == 2: return dist(hull[0], hull[1])

    # Find the point farthest from edge (hull[0], hull[1])
    j = 1
    while cross(hull[0], hull[1], hull[(j+1)%n]) > cross(hull[0], hull[1], hull[j]):
        j += 1

    best = 0
    for i in range(n):
        best = max(best, dist(hull[i], hull[j]))
        # Advance j while it moves farther from edge (i, i+1)
        while cross(hull[i], hull[(i+1)%n], hull[(j+1)%n]) > \
              cross(hull[i], hull[(i+1)%n], hull[j]):
            j = (j + 1) % n
            best = max(best, dist(hull[i], hull[j]))
    return best

Drill 5: Convex Hull — Step Through It

Given points: (0,0), (1,0), (2,1), (1,2), (0,2), (0.5,1). Trace the monotone chain algorithm:

Common Interview Questions

Chapter 9: Connections

Computational geometry does not exist in isolation. Every algorithm in this chapter builds on techniques from earlier CLRS chapters and connects to advanced topics.

Where These Ideas Come From

What Comes Next

Higher dimensions: All algorithms in this chapter work in 2D. Extending to 3D is straightforward for some (convex hull via gift wrapping or incremental algorithms) but much harder for others (Voronoi diagrams in 3D have O(n²) complexity). k-d trees (not in CLRS but critical in practice) handle nearest neighbor queries in arbitrary dimensions.

Randomized algorithms: Randomized incremental construction builds Delaunay triangulations and convex hulls by inserting points in random order. Expected O(n log n) with simpler code than deterministic algorithms. Trapezoidal decomposition for point location uses randomized incremental construction.

Robust geometric computation: Floating-point arithmetic causes orientation tests to give wrong answers near degeneracies. Industrial geometry libraries (CGAL, Shapely, JTS) use exact arithmetic or epsilon-based predicates. This is one of the deepest practical challenges in computational geometry.

k-d trees and range trees: For repeated nearest-neighbor queries, build a k-d tree in O(n log n) and answer each query in O(√n) expected time (O(log n) in low dimensions with good data). Range trees answer orthogonal range queries ("find all points in this rectangle") in O(log^d n + k) time. These are the data structure counterparts to the algorithmic techniques in this chapter.

Computational topology: Beyond convex hulls and Voronoi diagrams, modern geometry studies the topology of point clouds. Persistent homology detects "holes" at multiple scales — a point cloud shaped like a donut has a 1-dimensional hole. These methods (TDA — Topological Data Analysis) build on Delaunay complexes and are increasingly used in ML for understanding data shape.

Geometry in deep learning: Attention mechanisms in transformers compute pairwise dot products — a geometric operation. Embedding spaces in LLMs have convex cones corresponding to semantic concepts. The geometry of high-dimensional spaces (concentration of measure, JL lemma) explains why random projections work for approximate nearest neighbor search (LSH). Computational geometry is not just a CLRS chapter — it is the mathematical language of modern AI.

Recommended Further Reading

The Lower Bound: Why O(n log n)?

Can we build a convex hull faster than O(n log n)? Not in the comparison model. Here is the elegant proof: reduce sorting to convex hull.

Given n numbers x₁, ..., x_n to sort, create the points (x₁, x₁²), ..., (x_n, x_n²). These points all lie on the parabola y = x², which is convex. Therefore all points are on the convex hull. The convex hull, traversed in order, gives the points sorted by x-coordinate. If we could compute the hull in o(n log n), we could sort in o(n log n), contradicting the Ω(n log n) lower bound for comparison-based sorting.

The same argument shows that closest pair is Ω(n log n): you can reduce element uniqueness (which is Ω(n log n)) to closest pair (the closest pair has distance 0 if and only if there are duplicates).

Summary of All Algorithms