Minimum Enclosing Circle

Library Checker — Minimum Enclosing Circle: given n ≤ 10⁵ points, decide for each point whether it lies on the circumference of the minimum enclosing circle (MEC) of the whole set.

The MEC is the unique circle of smallest radius that contains every point (interior or boundary). There is a randomized algorithm that looks cubic but runs in expected O(n) time — Welzl’s algorithm (1991).

This lesson follows CP-Algorithms — Minimum enclosing circle.

Related: Basic Geometry.

Uniqueness

The MEC is unique. Intuition: if two distinct minimum-radius circles both enclose P, their intersection is a smaller-radius region that still covers P (lens argument / convex combination of centers), a contradiction. So the optimum is a single circle.

Welzl’s algorithm

Write mec(S) for the MEC of set S, and Pᵢ = {p₁,…,pᵢ}.

Apply a random permutation to the points.
Start with C = mec(p₁, p₂) (diameter circle).
For i = 3 … n, if pᵢ ∈ C, then C is already mec(Pᵢ).
Otherwise pᵢ must lie on the new boundary. Set C = mec(pᵢ, p₁) and scan j = 2 … i:
- If pⱼ ∈ C, keep going (MEC among circles through pᵢ).
- Else set C = mec(pᵢ, pⱼ) and scan k = 1 … j:
  - If pₖ ∉ C, set C = mec(pᵢ, pⱼ, pₖ) (unique circle through three non-collinear points).

Each nesting maintains an invariant: C is the MEC among circles forced to pass through 0, 1, or 2 extra boundary points. When an inner loop finishes, that invariant matches the parent’s current iteration — correctness follows.

mec enclosing_circle(vector<point> &p) {
    int n = p.size();
    shuffle(p.begin(), p.end(), gen);
    mec C{p[0], p[1]};
    for (int i = 0; i < n; i++) {
        if (!inside(C, p[i])) {
            C = mec{p[i], p[0]};
            for (int j = 0; j < i; j++) {
                if (!inside(C, p[j])) {
                    C = mec{p[i], p[j]};
                    for (int k = 0; k < j; k++) {
                        if (!inside(C, p[k])) {
                            C = mec{p[i], p[j], p[k]};
                        }
                    }
                }
            }
        }
    }
    return C;
}

(inside is O(1) — next section.)

Why expected `O(n)`?

The k-loop costs O(j) when it runs. The j-loop only enters the k-loop when pⱼ is one of the (at most two) points of Pⱼ that lie on the boundary of the MEC of Pⱼ through pᵢ and are critical (removing them shrinks the circle). After a random shuffle, that happens with probability ≤ 2/j.

Σⱼ (2/j) · O(j) = O(i)

The same backward-analysis argument lifts one level: the outer loop’s expected cost is O(n).

Point-in-MEC via `std::complex`

Represent points as complex<ftype> (z = x + y i). Multiplying adds polar angles; conj(z) negates the angle.

Two points (diameter)

mec(a,b) has diameter ab. Point z is inside iff ∠azb is not acute iff

I₀ = (b − z) · conj(a − z)

has non-positive real part (Re(I₀) ≤ 0).

Three points

Let

I₀ = (b − z) conj(a − z)
I₂ = (a − c) conj(b − c)
I₁ = I₀ · I₂

Signed angle sum ∠azb + ∠bca is the polar angle of I₁. Comparing signs of Im(I₂) and Im(I₁) (with the four side/orientation cases) tells inside vs outside. With integer input, the products stay integers (coefficients up to O(A⁴)).

Represent the MEC as the 2 or 3 defining points (not center/radius) so these predicates stay exact.

using ftype = int64_t;
using point = complex<ftype>;

struct mec { vector<point> pts; }; // size 2 or 3

ftype indicator(mec const& C, point z) {
    point a = C.pts[0], b = C.pts[1];
    point I0 = (b - z) * conj(a - z);
    if (C.pts.size() == 2) return real(I0);
    point c = C.pts[2];
    point I2 = (a - c) * conj(b - c);
    point I1 = I0 * I2;
    return imag(I2) < 0 ? -imag(I1) : imag(I1);
}

bool inside(mec const& C, point p) {
    return indicator(C, p) <= 0; // <0 inside, =0 on boundary, >0 outside
}

Library Checker asks which points have indicator == 0 on the final MEC.

🧪 Welzl on a tiny point set

Check: Welzl / MEC

Quiz