Convex Polygon Intersection, in depth: the half-plane view, a linear-time exact algorithm and the degenerate cases

The contract

Before writing code, fix the contract, because most bugs here are disagreements about it. Inputs are simple convex polygons given counter-clockwise with no repeated vertices. Three collinear consecutive vertices are tolerated but are best removed. The output is the closed intersection as a counter-clockwise vertex list. Decide what to return when the overlap has zero area: polygons that share only an edge or a vertex touch, but their overlap is a segment or a point. For area, IoU and most clipping uses, zero area should return an empty list, and that is the convention used below. If you need contact information, as a collision system does, compute it separately.

If you only need a yes or no answer, do not compute the polygon. The separating axis theorem says two convex polygons are disjoint exactly when some edge normal of either polygon separates their projections. That test is simple and fast for the small polygons typical in games. Point-in-polygon tests cover the containment half of the question.

Convex polygons as half-planes

A convex polygon with n edges is exactly the set of points on the left of every one of its directed edges, an intersection of n half-planes. The intersection of P and Q is therefore the intersection of all n + m half-planes. General half-plane intersection sorts the half-planes by direction angle and sweeps them with a deque, so it costs O(k log k) for k half-planes. Here the sort is free: walking a counter-clockwise convex polygon visits its edge directions in increasing angle, so each polygon's edge list is already sorted up to a rotation. Rotate each list to start at its smallest angle, merge the two, and the sweep runs in O(n + m).

The sweep keeps a deque of half-planes whose boundaries form the current partial polygon. When a new half-plane arrives, any half-plane at the back whose corner with its neighbour now lies outside the new one is redundant, so it is popped. The same check runs at the front, because the polygon closes on itself. Two parallel half-planes facing the same way keep only the tighter one. If, after the pops, the new half-plane points opposite to the one at the back of the deque, nothing survives between them and the answer is empty.

A tested implementation

Coordinates are integers or Python Fractions, so every orientation test is exact. Intersection points become rationals. I compared the area of every result against an exact Sutherland-Hodgman reference on 30,994 random pairs of convex hulls of integer points, on grids from 3 x 3, where degeneracies dominate, up to 1,000 x 1,000. There were no mismatches, and every non-empty output was strictly convex and counter-clockwise.

from collections import deque
from fractions import Fraction

def edges(poly):                      # (point, direction); keep the left side
    n = len(poly)
    return [(poly[i], (poly[(i+1) % n][0] - poly[i][0], poly[(i+1) % n][1] - poly[i][1]))
            for i in range(n)]

def half(d):                          # 0 for angles in [0, pi), 1 for [pi, 2 pi)
    return 0 if d[1] > 0 or (d[1] == 0 and d[0] > 0) else 1

def angle_less(d1, d2):
    h1, h2 = half(d1), half(d2)
    return h1 < h2 or (h1 == h2 and d1[0]*d2[1] - d1[1]*d2[0] > 0)

def by_angle(poly):
    E = edges(poly)
    k = min(range(len(E)), key=lambda i: half(E[i][1]))
    while angle_less(E[(k - 1) % len(E)][1], E[k][1]):
        k = (k - 1) % len(E)
    return E[k:] + E[:k]

def merge(A, B):
    out, i, j = [], 0, 0
    while i < len(A) or j < len(B):
        if j == len(B) or (i < len(A) and not angle_less(B[j][1], A[i][1])):
            out.append(A[i]); i += 1
        else:
            out.append(B[j]); j += 1
    return out

def outside(h, q):
    (p, d) = h
    return d[0]*(q[1]-p[1]) - d[1]*(q[0]-p[0]) < 0

def meet(h1, h2):
    (p, d), (q, e) = h1, h2
    t = Fraction((q[0]-p[0])*e[1] - (q[1]-p[1])*e[0]) / (d[0]*e[1] - d[1]*e[0])
    return (p[0] + t*d[0], p[1] + t*d[1])

def area2(poly):
    return sum(poly[i-1][0]*poly[i][1] - poly[i][0]*poly[i-1][1] for i in range(len(poly)))

def convex_intersection(P, Q):
    dq = deque()
    for h in merge(by_angle(P), by_angle(Q)):
        while len(dq) > 1 and outside(h, meet(dq[-1], dq[-2])):
            dq.pop()
        while len(dq) > 1 and outside(h, meet(dq[0], dq[1])):
            dq.popleft()
        if dq:
            d, e = dq[-1][1], h[1]
            if d[0]*e[1] - d[1]*e[0] == 0:            # parallel boundaries
                if d[0]*e[0] + d[1]*e[1] < 0:
                    return []                         # opposite after pops: empty
                if outside(h, dq[-1][0]):
                    dq.pop()                          # h is tighter
                else:
                    continue
        dq.append(h)
    while len(dq) > 2 and outside(dq[0], meet(dq[-1], dq[-2])):
        dq.pop()
    while len(dq) > 2 and outside(dq[-1], meet(dq[0], dq[1])):
        dq.popleft()
    if len(dq) < 3:
        return []
    pts = [meet(dq[i], dq[(i + 1) % len(dq)]) for i in range(len(dq))]
    res = [q for k, q in enumerate(pts) if q != pts[k - 1]]
    return res if area2(res) > 0 else []      # zero area counts as empty

The final filter matters. When polygons share only an edge or a corner, the sweep can end with three or more half-planes whose corners coincide or are collinear. Removing duplicate corners and rejecting zero area turns every touching case into the empty result the contract promises.

Worked example: a square and a triangle

Let P be the square with corners (0, 0), (4, 0), (4, 4), (0, 4) and Q the triangle (2, -1), (5, 3), (-1, 3). Q's lower-right edge has direction (3, 4). It cuts the bottom of the square at x = 11/4 and the right side at y = 5/3. Its lower-left edge, from (-1, 3) to (2, -1), cuts the bottom at x = 5/4 and the left side at y = 5/3. Its top edge, y = 3, cuts both sides. The square's top edge, y = 4, is redundant, because Q lies entirely below it, so the sweep pops it.

The function returns the six corners (11/4, 0), (4, 5/3), (4, 3), (0, 3), (0, 5/3), (5/4, 0), with area 119/12, about 9.917. As a sanity check, the square minus the band above y = 3 leaves 12, and the two corner triangles below the slanted edges have legs 5/4 and 5/3, so each has area 25/24. 12 minus 25/12 is 119/12.

The classic edge chase and its degenerate cases

The O'Rourke edge chase walks one edge of each polygon at a time. At every step it tests whether the two current edges cross, emits the crossing, and advances whichever edge is aiming at the other, emitting vertices of whichever polygon is currently inside. It is elegant and uses O(1) extra memory. Its weak point is that every advance rule assumes general position.

To make that concrete, I transcribed the textbook rules in floating point and compared the result with a Sutherland-Hodgman reference. On 3,000 random pairs with real coordinates it matched every time. On 19,243 pairs of hulls of points on a 7 x 7 integer grid, where shared edges, collinear edges and vertices lying on edges are common, 424 areas were wrong. My first transcription also looped forever on two polygons touching at one corner, because it reset its step counters on every crossing it found. The polygons shared no interior, so no inside flag was ever set. With that fixed, the same case returned the whole triangle as the answer, because the containment fallback counted a boundary point as inside. Published implementations handle these cases with additional rules. If you use one, run it against a simple reference on degenerate inputs first.

Measured cost

These are pure Python timings with rational arithmetic, on two regular polygons with coordinates snapped to a 10^6 grid. Sutherland-Hodgman, described in the polygon algorithms overview, clips the subject by each clip edge in turn, so it is O(nm). Ten times more vertices cost the merge ten times more time and Sutherland-Hodgman about a hundred times more. For triangles and boxes the quadratic method is simpler and fast enough. Above a few dozen vertices, use the linear method.

Operational guidance

• Snap to integers. Scale coordinates to a fixed grid, such as millimetres or 1e-6 degrees, and store integers. Orientation tests on integers are exact. In a language with 64-bit integers, check that cross products of coordinate differences fit, or use 128-bit intermediates.
• Keep exact results or round once. Intersection points are rational. Either keep them as rationals for further exact steps, or round to the grid once at the end, then re-validate convexity, because rounding can create tiny reflex corners.
• Validate inputs at the boundary. Check orientation by signed area, reverse clockwise inputs, remove duplicate and collinear vertices, and reject non-convex input. Build convex inputs with Andrew's monotone chain, which returns counter-clockwise order with collinear points removed.
• Use a library for general shapes. CGAL, Boost.Geometry and Clipper2 handle non-convex polygons and holes. Use the convex routine when convexity is guaranteed and speed matters.

Failure modes

• Float orientation errors. A cross product near zero can take the wrong sign, which flips an inside test and yields self-intersecting output. Use exact arithmetic, or an adaptive-precision predicate.
• Clockwise input. Every half-plane faces the wrong way and the answer is empty or nonsense. Normalise orientation first.
• Containment without crossings. Methods that only collect crossings must detect nesting separately. The half-plane sweep handles it naturally.
• Touching treated as overlap. A zero-area result returned as a two-vertex polygon breaks downstream area and IoU code. Apply the contract's zero-area rule.
• Repeated vertices. Zero-length edges have no direction. Remove them before building half-planes.

Trade-offs

What to do next

1. Write the contract down: orientation, zero-area policy, coordinate grid and output type.
2. Implement the code above and reproduce the worked example's area of 119/12.
3. Build a randomised test against exact Sutherland-Hodgman on small integer grids, where degeneracies are frequent, and run it in CI.
4. Snap production coordinates to integers before calling the routine and validate convexity after any rounding.
5. Read Graham scan and the art gallery theorem for the neighbouring convexity and visibility problems.