Smallest Enclosing Circle Algorithm: Interview Story, Theory, and C Implementation

The author recounts an interview where a candidate was asked to compute the smallest enclosing circle for a set of planar points, discussing the importance of understanding problem characteristics and using appropriate algorithms.

The problem is defined: given a set P of n (n>2) distinct points, find the unique smallest circle that contains all points. Special cases for n=3 are illustrated, showing that the optimal circle is either the circumcircle of the three points or the circle with the longest side as diameter, depending on whether the third point lies inside that diameter circle.

For n>3, the solution is reduced to the convex hull of P, because only the outermost points can affect the minimal circle. The convex hull is described intuitively as the shape formed by stretching a rubber band around the set of points.

A randomized incremental approach is introduced: points are processed in a random order P₁,…,Pₙ. Assuming the minimal circle D_{i‑1} for the first i‑1 points is known, if the new point P_i lies inside D_{i‑1} the circle remains unchanged; otherwise a new minimal circle D_i must be computed, which will have P_i on its boundary. The algorithm proceeds by examining at most three boundary points to construct the new circle.

The detailed step‑by‑step algorithm is presented, including how to handle the cases where a point lies outside the current circle, how to form circles using two points as a diameter, and how to compute the circumcircle of three points when necessary. The process repeats until all points are processed, yielding D_n, the minimal enclosing circle of the whole set.

A correctness argument is given: the minimal enclosing circle is unique, and any point not contained in the minimal circle of the remaining set must lie on the boundary of the overall minimal circle. This ensures the incremental construction is valid.

The algorithm runs in expected O(n) time and O(n) space because each insertion requires only constant‑time work after the convex hull preprocessing, and the overall loop iterates n times.

The full C source code implementing the algorithm is provided below. The code follows the described steps, defining point and circle structures, utility functions for distance, midpoint, orientation, and the core functions min_disc_2point and min_disc_points that maintain the minimal circle incrementally.

#include
#include
#define MAX_N       ( 1001 )
#define EPS_V       ( 1e-6 )

typedef struct { double x; double y; } point;

typedef struct { point o; double r; } circle;

point p_list[MAX_N];    // input point set
circle c_res;           // output minimal enclosing circle

int dcmp(double x) {
    if (fabs(x) < EPS_V) return 0;
    return x < 0 ? -1 : 1;
}

point midpoint(point p1, point p2) {
    point p; p.x = (p1.x + p2.x) / 2; p.y = (p1.y + p2.y) / 2; return p;
}

double multiply(point sp, point ep, point op) {
    return (sp.x - op.x)*(ep.y - op.y) - (ep.x - op.x)*(sp.y - op.y);
}

double distance(point p1, point p2) {
    return sqrt((p2.x - p1.x)*(p2.x - p1.x) + (p2.y - p1.y)*(p2.y - p1.y));
}

void min_disc_2point(point p, point q, int n) {
    c_res.o = midpoint(p, q);
    c_res.r = distance(p, q) / 2;
    for (int i = 0; i < n; ++i) {
        if (distance(c_res.o, p_list[i]) <= c_res.r) continue;
        if (dcmp(multiply(p, q, p_list[i]))) {
            double c1 = ((p.x*p.x + p.y*p.y) - (q.x*q.x + q.y*q.y)) / 2;
            double c2 = ((p.x*p.x + p.y*p.y) - (p_list[i].x*p_list[i].x + p_list[i].y*p_list[i].y)) / 2;
            c_res.o.x = (c1*(p.y - p_list[i].y) - c2*(p.y - q.y)) / ((p.x - q.x)*(p.y - p_list[i].y) - (p.x - p_list[i].x)*(p.y - q.y));
            c_res.o.y = (c1*(p.x - p_list[i].x) - c2*(p.x - q.x)) / ((p.y - q.y)*(p.x - p_list[i].x) - (p.y - p_list[i].y)*(p.x - q.x));
            c_res.r = distance(c_res.o, p_list[i]);
        } else {
            double t1 = distance(p, q);
            double t2 = distance(q, p_list[i]);
            double t3 = distance(p, p_list[i]);
            if (t1 >= t2 && t1 >= t3) { c_res.r = distance(p, q); c_res.o = midpoint(p, q); }
            else if (t2 >= t1 && t2 >= t3) { c_res.r = distance(p_list[i], q); c_res.o = midpoint(p_list[i], q); }
            else { c_res.r = distance(p_list[i], p); c_res.o = midpoint(p_list[i], p); }
        }
    }
}

void min_disc_points(point p, int n) {
    c_res.o = midpoint(p, p_list[0]);
    c_res.r = distance(p, p_list[0]) / 2;
    for (int i = 1; i < n; ++i) {
        if (distance(c_res.o, p_list[i]) <= c_res.r) continue;
        min_disc_2point(p, p_list[i], i);
    }
}

int main() {
    int n;
    while (scanf("%d", &n) != EOF) {
        if (n < 1) { printf("Invalid data, please enter n > 0.\n"); continue; }
        for (int i = 0; i < n; ++i) scanf("%lf%lf", &p_list[i].x, &p_list[i].y);
        if (n == 1) { printf("%.2lf %.2lf %.2lf\n", p_list[0].x, p_list[0].y, 0); continue; }
        c_res.o = midpoint(p_list[0], p_list[1]);
        c_res.r = distance(p_list[0], p_list[1]) / 2;
        for (int i = 2; i < n; ++i) {
            if (distance(c_res.o, p_list[i]) <= c_res.r) continue;
            min_disc_points(p_list[i], i);
        }
        printf("%.2lf %.2lf %.2lf\n", c_res.o.x, c_res.o.y, c_res.r);
    }
    return 0;
}

The article concludes that the algorithm is correct, runs in linear expected time, and can be further optimized by better point selection strategies; it also reflects on the joy of programming and the beauty of algorithms.