QuickHull, in depth: farthest-point recursion, exact orientation, measured cost from linear to quadratic and the path to Qhull

The idea and the one primitive

Everything rests on one primitive, the orientation test. For points o, a and b, cross(o, a, b) = (a.x - o.x)(b.y - o.y) - (a.y - o.y)(b.x - o.x). It is positive when b is to the left of the directed line o to a, negative when to the right, and zero when the three are collinear. Its absolute value is twice the area of triangle o, a, b, which is also the distance from b to the line times the length of o to a. For a fixed chord, the length is constant, so the point with the largest absolute cross value is the farthest point. No square roots, no division.

The algorithm, in four steps:

1. Take a, the lexicographically smallest point (least x, then least y), and b, the largest. Both are hull vertices.
2. Split the rest into points strictly right of a to b (the lower side) and strictly right of b to a (the upper side). Points exactly on the chord are dropped.
3. For a side with chord p to q and candidates S, find the point f in S farthest from the chord. f is a hull vertex.
4. Points inside triangle p, f, q are discarded. Recurse on the candidates strictly right of p to f, then on those strictly right of f to q, emitting f between them.

Why is f on the hull? Draw the line through f parallel to the chord. No candidate lies beyond it, because f is the farthest, and the chord side has no candidates by construction. So that line supports the whole set at f. If several points tie for farthest, they lie on that supporting line; the two outermost are hull vertices and any between them are collinear, and the strict tests in the next level drop those middle points automatically.

A correct implementation

The implementation below is the version used for every measurement in this article. It deduplicates points first, uses integer arithmetic so every orientation sign is exact, and only ever keeps points strictly outside a chord, which is what removes collinear points and handles ties. If all points are collinear it returns the two endpoints.

def cross(o, a, b):
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def quickhull(points):
    """Counter-clockwise hull, no collinear vertices. Exact for integer coordinates."""
    pts = sorted(set(points))
    if len(pts) < 3:
        return pts
    a, b = pts[0], pts[-1]
    out = [a]
    _hull(a, b, [p for p in pts if cross(a, b, p) < 0], out)   # lower side
    out.append(b)
    _hull(b, a, [p for p in pts if cross(a, b, p) > 0], out)   # upper side
    return out

def _hull(p, q, cand, out):
    """cand are strictly right of p->q; append hull vertices between p and q in order."""
    if not cand:
        return
    far = max(cand, key=lambda x: -cross(p, q, x))
    _hull(p, far, [x for x in cand if cross(p, far, x) < 0], out)
    out.append(far)
    _hull(far, q, [x for x in cand if cross(far, q, x) < 0], out)

Test it the way you would test any geometry code: against a second, independent algorithm on many random inputs. Five thousand random sets on small grids (coordinates from 0 to 2, 3 or 6, where duplicates, collinear runs and ties for farthest are common) plus wider ranges all produced exactly the same hull as Andrew monotone chain after rotating both to start at the smallest point. Small grids matter most; random points in a large square almost never produce the degenerate cases that break hull code.

Worked example

Take the eleven points in the figure. The extremes are a=(0,0) and b=(5,2). The orientation value cross(a, b, p) = 5p.y - 2p.x puts (4,0) at -8 and (3,1) at -1, so both are below; the other seven are above. Above the chord, (2,5) scores 21, more than (0,3) at 15 or (4,4) at 12, so it is the first hull vertex found, and the four points inside triangle a, b, (2,5) are gone.

Recursing on b to (2,5) leaves only (4,4) strictly outside; on (2,5) to a, only (0,3). Below the chord, (4,0) is farthest and (3,1) falls inside triangle a, (4,0), b. The output order is a, the lower chain, b, the upper chain, which is counter-clockwise.

What it costs, measured

Each recursion node scans its candidates a constant number of times and produces one hull vertex, so total work is the sum of candidate set sizes over nodes. That gives a clean worst-case bound of O(n h), where h is the number of hull vertices, and O(n squared) when every point is on the hull. When the farthest point splits candidates roughly in half, the depth is logarithmic and the cost is O(n log n). When most points fall inside the first few triangles, it is close to linear.

Measured candidate visits (candidates scanned per recursive call, not counting the first split) for the implementation above:

The square and disk cases are about 1.5 visits per point: the first triangles swallow almost everything. Points on a circle are all hull vertices, but the farthest point sits near the middle of each arc, so splits are balanced and the count grows like n log n. The adversarial input is a convex chain whose x coordinates halve at each step (x = 2 to the power i, y = x squared, plus the origin). The farthest point from each chord is always next to one end, so every level peels one vertex, and the count is exactly n(n-1)/2 for the chain. Note what that input needs: coordinates that double per point, which only exact big integers can hold. Real data does not look like that, but long smooth boundaries with uneven point density, such as a scanned outline sampled densely on one side, push QuickHull toward its O(n h) bound. If h can be large and latency matters, monotone chain is a guaranteed O(n log n).

Pre-filtering with Akl-Toussaint

Because QuickHull is a filter at heart, a cheap pre-filter helps every hull algorithm. The Akl-Toussaint heuristic finds the extreme points in a few directions, typically min and max of x, y, x+y and x-y, and discards every point strictly inside the polygon they form. One linear pass with a handful of orientation tests per point removes most of a uniform cloud before the real algorithm starts. It changes nothing for inputs where most points are on the hull, so treat it as a constant-factor tool.

Robustness with floating point

With floating-point coordinates, cross() can return the wrong sign when points are nearly collinear. In QuickHull that has two effects. A point that is actually outside a chord may be classified as inside and discarded, losing a hull vertex. Worse, the farthest-point choice and the partition tests can disagree with each other, so a point is kept by both children or neither. The result can be a polygon that is not convex or a hull missing a vertex, and neither is detected without a check.

The options, from simplest to most general: scale and round inputs to integers when the domain allows it, as in pixel or map-tile coordinates; use an adaptive exact predicate as described in the robust orientation predicate article; or accept approximate output and validate it. A cheap validator is to check that every consecutive triple of the output turns left and that no input point lies strictly outside any edge. The second check is O(n h), so run it in tests and debug builds rather than every call.

Recursion depth and an iterative version

The recursive version can reach depth h on adversarial input, which overflows the default Python recursion limit around a thousand frames and can overflow a native thread stack for very large h. An explicit stack fixes it. Push (p, q, candidates) frames, and to keep counter-clockwise order, push the right-hand subproblem before the left so the left one is processed first, emitting the farthest point between them:

def _hull_iter(p, q, cand, out):
    stack = [("solve", p, q, cand)]
    while stack:
        op, *args = stack.pop()
        if op == "emit":
            out.append(args[0])
            continue
        p, q, cand = args
        if not cand:
            continue
        far = max(cand, key=lambda x: -cross(p, q, x))
        right = [x for x in cand if cross(far, q, x) < 0]
        left = [x for x in cand if cross(p, far, x) < 0]
        stack.append(("solve", far, q, right))
        stack.append(("emit", far))
        stack.append(("solve", p, far, left))

Three dimensions and Qhull

In three dimensions and above, the same idea generalises to d-dimensional Quickhull, the algorithm implemented by Qhull (Barber, Dobkin and Huhdanpaa, 1996): maintain a hull of facets, assign each remaining point to a facet it is outside of, repeatedly take the farthest point of some facet, find all facets it can see, delete them, and stitch new facets to the horizon. SciPy exposes this as scipy.spatial.ConvexHull. Two practical notes: in 2D its volume attribute is the area and its area attribute is the perimeter, a frequent source of wrong numbers; and degenerate inputs such as coplanar points in 3D raise errors unless you pass Qhull options, for example QJ to joggle the input. For how QuickHull compares with the other common hull algorithms, see the convex hull overview.

Trade-offs against other hull algorithms

Rule of thumb: for one-off hulls of large random-looking clouds, QuickHull with an Akl-Toussaint pre-filter is hard to beat; for a library function whose callers you do not control, monotone chain gives a guarantee. Hulls often feed other algorithms, such as rotating calipers for diameter and width, or the pruning step in closest pair style geometry problems, so output order and collinear policy should be documented.

What to do next

1. Implement the integer version above and reproduce the eleven-point example exactly.
2. Add a randomized test that compares it with monotone chain on grids of size 2 to 6 and on duplicates.
3. Count candidate visits on your real data and compare them with n log n and n h.
4. Add an Akl-Toussaint pre-filter and measure how many points it removes.
5. Switch to the explicit-stack version before shipping, or cap input size.
6. If coordinates are floats, scale to integers or adopt an exact predicate, and add an output validator in tests.
7. Document whether collinear boundary points are included and in which order vertices are returned.