Half-Plane Intersection, in depth: sort by angle and sweep a deque in O(n log n), bounding boxes, parallel and empty cases, tested code and kernels

Half-planes and their intersection

Represent each half-plane by a point P on its boundary and a direction D along it. The feasible side is on the left: a point X is inside when cross(D, X - P) >= 0. This form needs no division and no normalisation, and the sign test is the same orientation predicate used everywhere else in computational geometry. A constraint written as ax + by <= c converts to it directly. Take D = (-b, a), which keeps the feasible side on the left, and take P as any point on the line.

The intersection is convex because each half-plane is convex, and an intersection of convex sets is convex. That has a useful consequence. Walk the boundary of the result counter-clockwise and the boundary lines appear in increasing order of angle, each at most once. The algorithm exploits exactly this. Most of the n inputs usually turn out to be redundant, meaning they never touch the boundary.

The algorithm: sort by angle, sweep with a deque

The algorithm has five steps.

1. Add four half-planes forming a large bounding box. This makes the result bounded, so every pair of neighbouring boundary lines meets at a vertex.
2. Sort all half-planes by the angle of D, using atan2(D.y, D.x). Among half-planes with the same angle only the tightest one can matter. The tightest one is the one whose boundary lies on the inside of the others.
3. Sweep in sorted order and keep a deque of the half-planes that are currently on the boundary. Before appending a new half-plane h, pop from the back while the vertex formed by the last two entries lies outside h. Then pop from the front while the vertex formed by the first two entries lies outside h.
4. When the sweep ends, the front and back must be reconciled, because the boundary wraps around. Pop the back while its last vertex lies outside the front half-plane, and pop the front while its first vertex lies outside the back half-plane.
5. If fewer than three half-planes remain, the region is empty or degenerate. Otherwise the vertices are the intersections of neighbouring deque entries, in counter-clockwise order.

Why pop from both ends? A new half-plane has the largest angle seen so far, so it can only cut off vertices formed by the most recent lines, at the back, or vertices near the start of the angular order, at the front, because the boundary is a closed loop. Each half-plane is appended once and popped at most once, so the sweep is O(n) and the sort dominates. The total is O(n log n). If the input arrives already sorted by angle, for example the edges of a convex polygon, the whole computation is linear.

A tested implementation

import math
from collections import deque

EPS = 1e-9

class HalfPlane:
    """Points X with cross(d, X - p) >= 0: the region LEFT of p -> p + d."""
    __slots__ = ("px", "py", "dx", "dy", "ang")
    def __init__(self, px, py, dx, dy):
        self.px, self.py, self.dx, self.dy = px, py, dx, dy
        self.ang = math.atan2(dy, dx)

    def out(self, x, y):
        return self.dx * (y - self.py) - self.dy * (x - self.px) < -EPS

def cross_dir(a, b):
    return a.dx * b.dy - a.dy * b.dx

def meet(a, b):                       # caller guarantees a and b are not parallel
    t = ((b.px - a.px) * b.dy - (b.py - a.py) * b.dx) / cross_dir(a, b)
    return a.px + t * a.dx, a.py + t * a.dy

def bounding_box(lo, hi):
    return [HalfPlane(lo, lo, 1, 0), HalfPlane(hi, lo, 0, 1),
            HalfPlane(hi, hi, -1, 0), HalfPlane(lo, hi, 0, -1)]

def intersect(planes, box=1e9):
    """CCW vertex list of the intersection, or [] if empty or degenerate."""
    hs = sorted(planes + bounding_box(-box, box), key=lambda h: h.ang)
    dq = deque()
    for h in hs:
        while len(dq) > 1 and h.out(*meet(dq[-1], dq[-2])):
            dq.pop()
        while len(dq) > 1 and h.out(*meet(dq[0], dq[1])):
            dq.popleft()
        if dq and abs(cross_dir(h, dq[-1])) < EPS:
            if h.dx * dq[-1].dx + h.dy * dq[-1].dy < 0:
                return []                    # antiparallel neighbours: empty
            if h.out(dq[-1].px, dq[-1].py):  # same direction: keep the tighter
                dq[-1] = h
            continue
        dq.append(h)
    while len(dq) > 2 and dq[0].out(*meet(dq[-1], dq[-2])):
        dq.pop()
    while len(dq) > 2 and dq[-1].out(*meet(dq[0], dq[1])):
        dq.popleft()
    if len(dq) < 3:
        return []
    return [meet(dq[i], dq[(i + 1) % len(dq)]) for i in range(len(dq))]

def through(ax, ay, bx, by):          # half-plane left of the directed edge a -> b
    return HalfPlane(ax, ay, bx - ax, by - ay)

Same-direction duplicates are handled inside the sweep rather than in a separate pass. When the new half-plane is parallel to the back of the deque and points the same way, the code keeps whichever is tighter. When the two are antiparallel and still neighbours after the pops, nothing lies between them and the region is empty. The bounding box guarantees that a box edge with an angle between them would otherwise separate them.

How it was tested. The oracle intersects every pair of boundary lines, keeps the points that satisfy every constraint, and takes their convex hull. That is O(n^3) but obviously correct. Over 3,000 random instances of 1 to 12 half-planes in a box of half-width 100, with 15% of the inputs deliberately parallel, antiparallel or duplicated, the two areas never disagreed. 1,799 of the instances were empty, which is typical for random constraints and is why empty-region handling deserves tests of its own.

Worked example

The constraints are y >= 0, x <= 4, y <= 3, x >= 0 and y >= x + 1. The code produces the triangle (0, 1), (2, 3), (0, 3), with area 2. Figure 1 lists the deque after each step. Two moments matter. Adding y >= x + 1 evicts y >= 0, because the corner (0, 0) violates the new constraint. Adding y <= 3 evicts x <= 4, because the corner (4, 5) lies above y = 3. Neither evicted constraint touches the final region, so both were redundant. The bounding-box planes never enter the deque, because each shares an angle with a user constraint and loses the tightness tie. In the code these are the same-direction branches.

What it is used for

Polygon kernel. The kernel of a simple polygon is the intersection of the half-planes to the left of its counter-clockwise edges. A guard placed anywhere in the kernel sees the whole polygon, and a polygon is star-shaped exactly when its kernel is non-empty. For the L-shape (0,0), (4,0), (4,1), (1,1), (1,4), (0,4), the code returns the unit square [0, 1] x [0, 1]. The general method costs O(n log n). Lee and Preparata's specialised algorithm finds a kernel in linear time.

Largest inscribed circle. Shift every half-plane inward by r and ask whether the intersection is still non-empty. Then bisect on r. For the 6-8-10 right triangle, 60 bisection steps returned r = 2.0000000000, which matches area / semiperimeter = 24 / 12. The centre is any point of the final, tiny region.

Voronoi cells. The cell of site s is the intersection of the half-planes closer to s than to each other site, bounded by perpendicular bisectors. Computing one cell this way costs O(n log n), which is fine for a few cells. To get all of them, use a proper Voronoi algorithm.

Linear programming and duality. Maximising a linear objective over the region means evaluating it at the vertices. If you only need the optimum, Seidel's randomised incremental algorithm solves two-variable LPs in expected O(n) time, and Megiddo's prune-and-search does so deterministically, without building the polygon. The LP duality view also explains a classic equivalence. Under point-line duality, intersecting upper half-planes is the same problem as computing the lower convex hull of the dual points.

Measured costs

Measured with the code above in CPython 3.13. The inputs were n tangent lines to the unit circle at random angles, so nearly every half-plane is on the boundary and the result approximates the circle.

At 100,000 inputs two vertices disappeared. Two tangents had angles so close that the |cross| < EPS test treated them as parallel and kept only the tighter one. The area is unaffected, but if you count vertices the tolerance is part of the specification. A redundancy-heavy input of 100,000 tangents at radii between 1 and 2 took 1,027 ms and left 127 vertices. The sweep pops more there, but its work is still linear. In this Python code the per-element sweep costs far more than the sort.

Operational guidance

Choose the bounding box from your data, not from the largest float. Coordinates near 1e9 carry an absolute rounding error of about 1e-7 in double precision, which swamps an EPS of 1e-9. A box about ten times your data extent is a good default. Then check whether any box plane is still in the final deque. If one is, the true region is unbounded, and you should report that rather than a polygon with arbitrary corners.

Scale EPS with the input. The test compares a cross product, whose magnitude grows with the lengths of D and with distance, so either normalise directions or set EPS relative to the coordinate scale. With integer inputs you can go exact. Sort by which half of the circle each direction points into (upper or lower), then by cross product, instead of using atan2, and evaluate out on rational vertices using Python's fractions.Fraction or 128-bit integers. That removes EPS entirely at a constant-factor cost.

Keep the brute-force oracle in your test suite. It is twenty lines long, catches every ordering bug, and its random generator should deliberately produce parallel, antiparallel, duplicate and empty cases. For small inputs, roughly a few dozen, clipping a convex polygon by each half-plane in turn is O(n^2) but simpler and very hard to get wrong. The polygon algorithms page covers the clipping step.

Failure modes

• No bounding box. Neighbouring antiparallel lines then have no meeting point, meet divides by zero, and unbounded results come out as garbage.
• A box that is too small silently clips a real region. A box that is too large destroys precision.
• Equal angles not deduplicated. Two parallel lines become neighbours, their meeting point is undefined, and the sweep pops valid entries.
• The wrong side. Mixing "left of D" with "right of D" conventions, or feeding clockwise polygon edges, makes kernels come back empty. Assert the orientation with a signed-area check.
• Degenerate results. Constraints that meet in a single point or a segment return [] with this code. If those cases are meaningful to you, for example a tight feasibility check, decide what to return and test it explicitly.
• Skipping the final clean-up loops. This leaves a wrap-around vertex outside the first half-plane, so the polygon gains a spurious spike.

Trade-offs

What to do next

1. Write each constraint as a point plus a direction with the feasible side on the left, and assert the convention on a known example.
2. Pick the bounding box and EPS from your coordinate scale. If inputs are integers, consider exact arithmetic.
3. Copy the code above together with the brute-force oracle, and run at least a few thousand random cases that include parallel, antiparallel and empty inputs.
4. Decide the policy for unbounded, point and segment results, and test each one.
5. If you need only an optimum, switch to a 2D LP routine. If you need many Voronoi cells, use a Voronoi algorithm instead of one intersection per cell.
6. Benchmark at your real n and check whether the sort or the sweep dominates before you optimise either one.