3D Convex Hull, in depth: the orientation determinant, horizons, randomized incremental construction and robust Qhull use

How big a 3D hull can be

Start with what the output looks like. If every face is split into triangles, the hull surface is a planar graph, and Euler's formula V - E + F = 2 applies. Every triangle has three edges and every edge borders exactly two triangles, so 3F = 2E. Substituting gives F = 2h - 4 and E = 3h - 6, where h is the number of hull vertices. The output is therefore linear in the number of points: never more than 2n - 4 triangles. Running SciPy on 1,000 Gaussian points produced h = 38 hull vertices and exactly 72 = 2 x 38 - 4 triangles, as the formula predicts.

That linear size bound is what makes O(n log n) possible, and the 2D sorting lower bound shows it cannot be beaten in the worst case. Optimal algorithms include Preparata and Hong's divide and conquer, randomized incremental construction (expected time), and Chan's output-sensitive O(n log h).

The orientation determinant

Every 3D hull algorithm rests on one predicate. Given an oriented triangle (a, b, c) and a point d, compute the signed volume of the tetrahedron they form:

orient(a, b, c, d) = det | b-a |
                         | c-a |  = ((b - a) x (c - a)) . (d - a)
                         | d-a |

The cross product (b - a) x (c - a) is the triangle's normal; its direction follows the right-hand rule, so listing the vertices counter-clockwise as seen from outside makes the normal point outward. The dot product then measures how far d lies along that normal. Positive means d is in front of the face (outside, so the face is visible from d), negative means behind, and zero means coplanar. The value is six times the tetrahedron's volume, which the implementation below uses to compute the hull volume for free.

With integer coordinates the determinant is exact. Its magnitude grows with the cube of the coordinate range, so in C++ check that range against 64-bit limits; Python integers never overflow. Doubles are a different story, covered under robustness. The 2D version of this test is covered in the robust orientation predicate article, and everything said there about error bounds carries over.

Incremental construction and the horizon

The incremental algorithm keeps a valid hull of the points seen so far and adds one point at a time. For a new point p there are two cases. If p sees no face, it is inside the current hull and nothing changes. Otherwise the set of visible faces is a connected patch of the surface, like the part of the Moon lit by a lamp. Its boundary is a single closed cycle of edges called the horizon. Delete the visible faces and connect p to every horizon edge with a new triangle. The result is the hull of the old points plus p.

Finding the horizon needs no geometry. Store each face as an oriented triple (i, j, k), meaning directed edges (i, j), (j, k), (k, i). With consistent orientation, the face across edge (i, j) contains (j, i). A directed edge of a visible face is on the horizon exactly when its reverse is not an edge of a visible face. The new triangle (u, v, p) reuses the edge's direction, so outward orientation is preserved automatically.

A tested implementation

The implementation below is deliberately small: integer coordinates, a set of oriented faces, and a linear scan for visible faces. It was checked against scipy.spatial.ConvexHull on 200 random integer point sets of up to 200 points: volumes matched, every directed edge had its twin (a closed surface), and no input point lay in front of any face.

import random

def orient(a, b, c, d):
    ux, uy, uz = b[0]-a[0], b[1]-a[1], b[2]-a[2]
    vx, vy, vz = c[0]-a[0], c[1]-a[1], c[2]-a[2]
    wx, wy, wz = d[0]-a[0], d[1]-a[1], d[2]-a[2]
    return ux*(vy*wz - vz*wy) - uy*(vx*wz - vz*wx) + uz*(vx*wy - vy*wx)

def initial_simplex(P):
    # Raises StopIteration if all points are collinear or coplanar: handle upstream.
    a = 0
    b = next(i for i in range(len(P)) if P[i] != P[a])
    def off_line(i):
        u = [P[b][k] - P[a][k] for k in range(3)]
        v = [P[i][k] - P[a][k] for k in range(3)]
        return any([u[1]*v[2] - u[2]*v[1], u[2]*v[0] - u[0]*v[2], u[0]*v[1] - u[1]*v[0]])
    c = next(i for i in range(len(P)) if off_line(i))
    d = next(i for i in range(len(P)) if orient(P[a], P[b], P[c], P[i]) != 0)
    if orient(P[a], P[b], P[c], P[d]) > 0:   # make d lie behind face (a, b, c)
        b, c = c, b
    return a, b, c, d

def hull3d(P, seed=0):
    # P: list of distinct integer (x, y, z). Returns outward-oriented triangles.
    a, b, c, d = initial_simplex(P)
    faces = {(a, b, c), (a, d, b), (b, d, c), (c, d, a)}
    rest = [i for i in range(len(P)) if i not in (a, b, c, d)]
    random.Random(seed).shuffle(rest)          # random order: the expected-time argument
    for p in rest:
        visible = [f for f in faces if orient(P[f[0]], P[f[1]], P[f[2]], P[p]) > 0]
        if not visible:
            continue                           # p is inside (or on) the current hull
        edges = set()
        for (i, j, k) in visible:
            edges |= {(i, j), (j, k), (k, i)}
        horizon = [(u, v) for (u, v) in edges if (v, u) not in edges]
        faces.difference_update(visible)
        faces.update((u, v, p) for (u, v) in horizon)
    return faces

def volume(P, faces):
    r = P[next(iter(faces))[0]]                # any point on the hull works as apex
    return -sum(orient(P[i], P[j], P[k], r) for (i, j, k) in faces) / 6

Two caveats. The visibility scan makes this version O(n F), quadratic in the worst case; the conflict graph below removes it. And the strict > 0 test lets a point lying in a face's plane become a vertex when it is outside a neighbouring face. Across the 200 test sets, 18 such boundary points were kept that SciPy omits, and on a 4 x 4 x 4 grid the code kept 25 vertices where the cube has 8 corners (the volume was still exactly 27). If you need only extreme points, merge coplanar neighbouring faces or filter vertices afterwards.

Worked example: six points

Take six points: the tetrahedron O = (0,0,0), X = (4,0,0), Y = (0,4,0), Z = (0,0,4), an interior point (1,1,1) and an exterior point Q = (3,3,3). The initial simplex uses O, X, Y and Z. Orienting gives faces (O,Y,X) on the floor, (O,Z,Y), (O,X,Z) and the slanted face (X,Y,Z) on the plane x + y + z = 4, each listed counter-clockwise from outside.

Insert (1,1,1). For the slanted face, orient is proportional to 1 + 1 + 1 - 4 = -1, and the axis faces are negative too, so the point is discarded. Insert Q. Only the slanted face sees it (3 + 3 + 3 - 4 = 5), and its three edges form the horizon. Replacing it with (X,Y,Q), (Y,Z,Q) and (Z,X,Q) gives h = 5 and 2 x 5 - 4 = 6 faces. The code reports volume 24; SciPy agrees on the volume and the five vertices.

Randomization and the conflict graph

To reach expected O(n log n), stop scanning. Maintain a conflict graph: a bipartite graph linking every point not yet inserted to every current face it can see, stored as a list on both sides. When p is inserted, its conflict list already is the visible set, so finding it costs nothing extra. Only the new faces need conflict lists. The key lemma: if a later point q sees a new face (u, v, p), then q already saw at least one of the two old faces that shared the horizon edge (u, v), the deleted visible one and the surviving one beyond it. So the candidates for a new face are just the union of those two lists, each tested once with orient.

Clarkson and Shor's backwards analysis shows that, for a uniformly random insertion order, the expected total number of faces created is O(n) and the expected total size of all those candidate lists is O(n log n). The shuffle is the whole guarantee: insert points in sorted or adversarial order and the bound is gone. Qhull's Quickhull uses a close relative of this structure: each outside point is assigned to one facet's outside set, and the farthest point of a facet is processed next instead of a random one.

build tetrahedron; for each remaining point q, for each face f: link(q, f) if q sees f
for p in random order:
    V = conflict_faces(p)                             # exactly the faces p sees
    if V is empty: continue                           # inside, discard
    H = directed edges of V whose twin face is not in V
    for (u, v) in H:
        f_in, f_out = face of V on (u, v), face beyond (u, v)
        g = new face (u, v, p); link neighbours along the cone
        for q in conflict_points(f_in) | conflict_points(f_out), q != p:
            if orient(g, q) > 0: link(q, g)
    unlink and delete every face in V

Robustness and Qhull

In 3D, floating-point error corrupts topology, not just coordinates. One wrong orient sign can disconnect the visible set, split the horizon into several cycles, and leave edges bordering three faces, which every later step then trusts. Three standard defences:

• Exact arithmetic. Snap inputs to an integer grid, or use adaptive predicates that fall back to exact arithmetic only near zero.
• Facet merging. Qhull by default merges facets that are coplanar within a computed tolerance.
• Joggling. Qhull's QJ option slightly perturbs the input so every facet is simplicial, trading exactness for simple output.

In Python, scipy.spatial.ConvexHull wraps Qhull. Option Qt (triangulated output) is always enabled. The object gives simplices, neighbors, vertices, volume, area (surface area in 3D) and equations: rows of [normal, offset] with outward unit normals, so a point x is inside when equations[:, :3] @ x + equations[:, 3] <= tol for every row. That one line is the fastest correct containment test once the hull exists.

The lifting map: hulls make Delaunay triangulations

One reason 3D hulls matter outside graphics is the lifting map. Send each planar point (x, y) to (x, y, x^2 + y^2) on a paraboloid. The faces of the lower hull of the lifted points, projected back down, are exactly the triangles of the Delaunay triangulation, because a point lies inside a triangle's circumcircle exactly when its lifted image lies below the plane through the lifted triangle. So a 3D hull routine is also a Delaunay routine; this is how Qhull computes it. See the Delaunay triangulation article for the in-circle test that the lifting map turns into a 3D orientation test.

Failure modes

• Degenerate input. All points collinear or coplanar means no initial tetrahedron exists. Detect it and fall back to a 2D hull in the containing plane.
• Duplicate points. These make zero-length edges. Deduplicate before building, as the tests above do.
• Float orientation. A wrong sign gives a non-manifold surface. Validate by checking that every directed edge has exactly one twin and that no input point is strictly outside any face.
• Sorted insertion. Inserting along a scan line removes the randomness the expected bound depends on. Always shuffle.
• Mixed winding. If downstream code assumes clockwise faces, every normal points inward and containment tests silently invert.

Trade-offs

What to do next

1. Run the code on the six-point example and confirm a volume of 24 and six faces.
2. Add the validation checks (edge twins, no point outside) and run them on 1,000 random integer sets against SciPy's volume.
3. Replace the visibility scan with conflict lists and time both versions as n doubles.
4. Lift 2D points onto the paraboloid, take the lower hull, and compare the result with scipy.spatial.Delaunay.
5. For production meshes, call Qhull through SciPy, decide between merged and joggled output, and record that choice next to the call.
6. Read the 2D convex hull article and the QuickHull article to see which ideas carry across dimensions.