Smallest Enclosing Circle (Welzl's Algorithm)

How Welzl's algorithm works

You're handed a cloud of points and asked for the smallest circle that contains all of them. Try the obvious things and they all fail: the centroid minimizes average distance, not the farthest point; the bounding-box circle can be 41% too big; and "circle through the two most distant points" ignores that a third point off the diameter might poke out. The correct object is the 1-center — the centre that minimizes the maximum distance to any point — and its circle is unique.

The structural insight that makes the whole problem tractable: the smallest enclosing circle is pinned by at most three points sitting exactly on its boundary. Either two points form a diameter, or three points form a triangle whose circumcircle is the answer. Every other point sits strictly inside. So the "answer" is really a tiny set of 2 or 3 support points; once you know them, the circle is determined by a closed-form formula.

Emo Welzl's 1991 algorithm turns that observation into a recursion. Process points one at a time. Keep the current smallest circle D for the points seen so far. When a new point p arrives:

1. If p is already inside D, do nothing — D is still optimal.
2. If p is outside D, then p must lie on the boundary of the new answer. Recurse, recomputing the circle of all earlier points but now with p forced onto the boundary.

The forced-boundary set R never exceeds 3 points, so the recursion bottoms out fast: with 0, 1, 2, or 3 boundary points the circle is computed directly (the empty circle, a point, a diameter, or a circumcircle). The magic ingredient is a random shuffle of the input first. Shuffling makes step 2 rare — for the i-th point the chance it lands outside is at most 3/i — and that single fact collapses the expected running time from quadratic to linear.

Why it's expected linear

Here is the backwards-analysis argument, because it's short and it's the reason anyone uses Welzl over the older approaches. Fix the first i points (in random order) and ask: what is the probability that the i-th point was outside the circle of the first i−1, forcing a recomputation? Equivalently, remove a uniformly random one of the i points and ask how often the circle changes.

The circle only changes if you removed one of its (at most 3) support points. Since each of the i points is equally likely to be the one removed, that probability is at most 3/i. A recomputation costs O(i) work, so the expected cost of step i is at most i · (3/i) = 3, a constant. Summing over all n insertions gives expected O(n) total — and crucially, the n outer "is it inside?" tests are each O(1), so they sum to O(n) as well.

E[total work]  ≤  Σ_{i=1..n}  O(1)  +  (3/i)·O(i)
               =  Σ_{i=1..n}  O(1)
               =  O(n)

The expectation is over the random shuffle, not over the inputs. That's the whole trick: there is no worst-case input, because you choose the order. And the constant is small: the i-th insertion triggers a recompute with probability ≤ 3/i, so the expected number of recomputes over the whole run is Σ 3/i ≈ 3 ln n — only about 40 even for a million points — and each one is cheap on average, which is what keeps the total work linear.

When to use it

• Bounding-volume hierarchies and collision culling — a minimum circle (or sphere) is the tightest rotation-invariant bound, so it rejects more non-collisions than a bounding box.
• Facility location / the 1-center problem — place one transmitter, warehouse, or emergency depot to minimize the worst-case distance to any client.
• Tolerance and metrology — fitting the minimum circumscribed circle to a machined part to check it fits a hole (a standard GD&T measurement).
• Map and chart labelling — sizing a marker that must cover a cluster of geographic points.
• Machine learning — the support-vector "minimum enclosing ball" underlies one-class SVMs and core-set approximations for clustering.

Reach for a different tool when the inputs are dynamic (points inserted and deleted online — Welzl recomputes from scratch), when you need a min-enclosing rectangle or convex region, or when n is in the billions and a (1+ε)-approximation via core-sets is good enough and far cheaper.

Welzl vs other minimum-circle methods

	Welzl (randomized)	Megiddo (prune & search)	Skyum / convex-hull	Elzinga–Hearn	Core-set (1+ε)
Time complexity	O(n) expected	O(n) worst case	O(n log n)	O(n²) worst case	O(n / ε)
Determinism	Las Vegas (random order)	Fully deterministic	Deterministic	Deterministic	Deterministic
Constant factor	Tiny (few recomputes)	Large, branchy	Moderate (hull first)	Small per step	Tiny, ε-tunable
Lines of code	~40	Hundreds	~80 (needs a hull)	~50	~30
Exact answer	Yes	Yes	Yes	Yes	No — within ε
Extends to d dimensions	Yes (basis d+1)	Hard past 2D	2D only	Slow	Yes — best for huge d/n
Best for	The default in practice	Theoretical interest	When hull is already computed	Tiny inputs / teaching	Streaming / massive n

Megiddo's 1983 prune-and-search achieves the same linear bound deterministically, but the constant and the code complexity are so much larger that almost every production library (CGAL, Boost.Geometry's variant, miniball) ships Welzl or a Welzl-derived "move-to-front" heuristic instead. The headline is that randomization buys you a 40-line routine that beats a hundreds-of-lines deterministic one on real inputs.

What the numbers actually say

• Expected circle recomputations ≈ 3 ln n — that's Θ(log n), not Θ(n). The i-th insertion forces a recompute with probability ≤ 3/i, so the expected count of recomputes sums to Σ 3/i ≈ 3 ln n (about 40 even for a million points). Total work stays O(n) because the i-th recompute costs O(i) but happens with probability 3/i — the i's cancel, leaving O(1) expected work per insertion.
• The bounding-box circle can be 41% too large. Take the four edge-midpoints of a unit square — they lie on a circle of radius 0.5, so that is the smallest enclosing circle. But the axis-aligned bounding box is the whole unit square, and the circle circumscribing that box has radius half the diagonal = (√2)/2 ≈ 0.707. The box-circle's radius is √2 ≈ 1.414× the optimum here, a 41% overshoot in radius and ~100% in area. (Square corners, by contrast, are concyclic with the box-circle, so they show no overshoot — the gap needs points that don't reach the corners.)
• One million 2D points: a few milliseconds. A clean Welzl implementation does ~1M inside-tests (each a squared-distance compare, no sqrt) plus a handful of 3-point circumcircle solves. On a modern CPU that's single-digit milliseconds; the shuffle dominates only because it touches memory randomly.
• No square roots on the hot path. Compare squared distances against the squared radius. You only take a square root once, at the very end, if the caller actually needs the radius rather than radius².

JavaScript implementation

The full recursive Welzl, plus the three base-case circle builders. Points are {x, y}; a circle is {x, y, r}.

const EPS = 1e-9;

const dist2 = (a, b) => (a.x - b.x) ** 2 + (a.y - b.y) ** 2;
const inCircle = (c, p) => dist2(c, p) <= c.r * c.r + EPS;

function circleFrom2(a, b) {
  const cx = (a.x + b.x) / 2, cy = (a.y + b.y) / 2;
  return { x: cx, y: cy, r: Math.sqrt(dist2(a, b)) / 2 };
}

// Circumcircle of three points (assumes they're not collinear).
function circleFrom3(a, b, c) {
  const ax = a.x, ay = a.y, bx = b.x, by = b.y, cx = c.x, cy = c.y;
  const d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by));
  if (Math.abs(d) < EPS) return null;              // collinear
  const a2 = ax * ax + ay * ay, b2 = bx * bx + by * by, c2 = cx * cx + cy * cy;
  const ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d;
  const uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d;
  const center = { x: ux, y: uy };
  return { x: ux, y: uy, r: Math.sqrt(dist2(center, a)) };
}

// Smallest circle with the ≤3 boundary points in R already fixed.
function trivial(R) {
  if (R.length === 0) return { x: 0, y: 0, r: 0 };
  if (R.length === 1) return { x: R[0].x, y: R[0].y, r: 0 };
  if (R.length === 2) return circleFrom2(R[0], R[1]);
  // 3 points: try the circle on each pair; if it covers the third, use it
  // (handles the obtuse-triangle case). Otherwise use the circumcircle.
  for (let i = 0; i < 3; i++)
    for (let j = i + 1; j < 3; j++) {
      const c = circleFrom2(R[i], R[j]);
      if (R.every(p => inCircle(c, p))) return c;
    }
  return circleFrom3(R[0], R[1], R[2]);
}

// Iterative move-to-front Welzl — the practical variant.
function smallestEnclosingCircle(points) {
  const pts = points.slice();
  for (let i = pts.length - 1; i > 0; i--) {       // Fisher–Yates shuffle
    const j = (Math.random() * (i + 1)) | 0;
    [pts[i], pts[j]] = [pts[j], pts[i]];
  }
  let c = trivial([]);
  for (let i = 0; i < pts.length; i++) {
    if (inCircle(c, pts[i])) continue;
    c = { x: pts[i].x, y: pts[i].y, r: 0 };         // pts[i] on boundary
    for (let j = 0; j < i; j++) {
      if (inCircle(c, pts[j])) continue;
      c = circleFrom2(pts[i], pts[j]);              // two boundary points
      for (let k = 0; k < j; k++) {
        if (inCircle(c, pts[k])) continue;
        c = circleFrom3(pts[i], pts[j], pts[k]);    // three boundary points
      }
    }
  }
  return c;
}

Note this is the iterative three-loop form, which is what people actually ship — it avoids recursion-depth limits and is identical in expectation to Welzl's recursion. The trivial helper above is shown for the recursive variant; the iterative loop inlines the same logic by trying pairs before the circumcircle.

Python implementation

import math, random

EPS = 1e-9

def dist2(a, b):
    return (a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2

def in_circle(c, p):
    cx, cy, r = c
    return (p[0] - cx) ** 2 + (p[1] - cy) ** 2 <= r * r + EPS

def circle_from2(a, b):
    cx, cy = (a[0] + b[0]) / 2, (a[1] + b[1]) / 2
    return (cx, cy, math.sqrt(dist2(a, b)) / 2)

def circle_from3(a, b, c):
    ax, ay, bx, by, cx, cy = a[0], a[1], b[0], b[1], c[0], c[1]
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < EPS:                       # collinear
        return None
    a2, b2, c2 = ax*ax + ay*ay, bx*bx + by*by, cx*cx + cy*cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return (ux, uy, math.sqrt((ux - ax) ** 2 + (uy - ay) ** 2))

def smallest_enclosing_circle(points):
    pts = points[:]
    random.shuffle(pts)                    # the line that makes it O(n)
    c = (0.0, 0.0, 0.0)
    for i, pi in enumerate(pts):
        if in_circle(c, pi):
            continue
        c = (pi[0], pi[1], 0.0)            # pi must be on the boundary
        for j in range(i):
            pj = pts[j]
            if in_circle(c, pj):
                continue
            c = circle_from2(pi, pj)
            for k in range(j):
                pk = pts[k]
                if in_circle(c, pk):
                    continue
                c = circle_from3(pi, pj, pk)
    return c                               # (center_x, center_y, radius)

The triple-nested loop looks like O(n³), and in the worst-case ordering it is O(n²); the random shuffle is what guarantees each inner loop runs to completion only a constant number of times in expectation. Delete the random.shuffle line and feed sorted-by-distance points and you'll measure the quadratic blow-up directly.

Variants worth knowing

Minimum enclosing ball in d dimensions. The same recursion solves the problem in any fixed dimension; the boundary basis grows from 3 (a circle) to d+1 (a ball needs up to 4 points in 3D). Gärtner's miniball library is the canonical fast implementation with numerically robust base-case solving up to ~30 dimensions.

Move-to-front and pivot heuristics. Bernd Gärtner's practical version keeps the support points near the front of the list, so future insertions hit them first and rarely trigger a recompute. It's still expected linear but with a markedly smaller constant on clustered data.

Megiddo's prune-and-search. A deterministic O(n) method that, at each round, discards a constant fraction of points that provably can't be support points. Theoretically beautiful, practically slower than Welzl because of its large constant and branchy control flow.

Core-set (1+ε)-approximation. For streaming or enormous n, maintain a small "core set" of O(1/ε) points whose minimum circle is within (1+ε) of the true answer. Time is O(n/ε), independent of dimension's exponential blowup — the method of choice for high-dimensional minimum-enclosing-ball in ML.

Smallest enclosing circle with k outliers. Allow up to k points to fall outside (robust fitting). Much harder — the unique-3-point structure breaks — and typically solved with parametric search or LP-type methods, not vanilla Welzl.

Common bugs and edge cases

• Forgetting the shuffle. The code is correct without it but degrades to O(n²) on adversarial or already-sorted input. The shuffle is not optional for the complexity guarantee.
• The obtuse-triangle trap. Three points don't always give a circumcircle answer. If the triangle is obtuse, the smallest enclosing circle is built on the longest edge (a diameter), and the opposite vertex sits inside. The pair-first loop in the code handles this; a naive "always circumcircle of the 3 support points" is wrong.
• Collinear or duplicate points. circleFrom3 divides by a determinant that is zero for collinear points — guard with an epsilon and fall back to the pair circle. Duplicate points are harmless but waste an inside-test.
• Strict vs non-strict inside test. Use ≤ r² + ε, not < r². Boundary points are inside; a strict comparison makes the algorithm loop or recompute on the very points that define the circle, thanks to floating-point round-off.
• Comparing distances with a square root. Calling Math.sqrt in the inside-test is both slower and a source of round-off. Keep everything squared until the final radius.
• Empty or single-point input. Define the base cases explicitly: zero points → radius 0 at the origin (or raise), one point → radius 0 at that point. Skipping these crashes the inside-test on an undefined circle.