Randomized Bounding ball

randomized_bounding_ball.cpp

#include <iostream>
#include <cmath>
#include <complex>
#include <algorithm>
using namespace std;
typedef complex<double> point;
#define x real()
#define y imag()
#define EPS 1e-7
#define PRECISION 12
#define cpoint const point&

/**
 * The class that represents a circle.
 */
struct circle {
    point center;
    double radius;
    circle() {}
    circle(cpoint center, double radius):
        center(center),
        radius(radius) {}
};

/**
 * Converts from barycentric coordinates to cartesian coordinates.
 *
 * @param (A, B, C)  the cartesian coordinates of the base triangle
 * @param (a, b, c)  the barycentric coordintes of the wanted point
 * @return           the cartesian coordinates of the wanted point
 */
point bary(cpoint A, cpoint B, cpoint C, double a, double b, double c) {
    return (A*a + B*b + C*c) / (a + b + c);
}

/**
 * Gets the circumcenter of a triangle (A, B, C).
 *
 * @param (A, B, C)  the triangle
 * @return           the circumcenter of the triangle
 */
point circumcenter(cpoint A, cpoint B, cpoint C) {
    double a = norm(B-C), b = norm(C-A), c=norm(A-B);
    return bary(A, B, C, a*(b + c - a), b*(c + a - b), c*(a + b - c));
}

/**
 * Computes for the minimum bounding ball for the "shaky-cam" radius for
 * error analysis.
 *
 * @param points  the array of points
 * @param n       the size of the array
 * @return        the circle representing the bounding ball
 */
circle bounding_ball(point points[], int n) {
    point *p = new point[n];
    copy(points, points + n, p);
    // randomize
    random_shuffle(p, p + n);
    circle ball;
    for (int i = 0; i < n; ++i) {
        if (abs(ball.center - p[i]) > ball.radius + EPS) {
            ball = circle(p[i], 0);
            for (int j = 0; j < i; ++j) {
                if (abs(ball.center - p[j]) > ball.radius + EPS) {
                    // get the midpoint
                    ball.center = (p[i] + p[j]) / 2.0;
                    ball.radius = abs(ball.center - p[i]);
                    for (int k = 0; k < j; ++k) {
                        if (abs(ball.center - p[k]) > ball.radius + EPS) {
                            ball.center = circumcenter(p[i], p[j], p[k]);
                            ball.radius = abs(ball.center - p[i]);
                        }
                    }
                }
            }
        }
    }
    delete[] p;
    return ball;
}

/**
 * Input format:
 *  First line contains N, the number of sample points.
 *  The next N lines contain two real numbers X, Y separated by a single space.
 * Output format:
 *  Outputs the minimum enclosing circle's X, Y, and radius separated by single
 *  spaces.
 */
int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    int N;
    cin >> N;
    point *p = new point[N];
    for (int i = 0; i < N; ++i) {
        double X, Y;
        cin >> X >> Y;
        p[i] = point(X, Y);
    }
    circle ball = bounding_ball(p, N);
    cout.precision(PRECISION);
    cout << fixed
         << ball.center.x
         << ' '
         << ball.center.y
         << ' '
         << ball.radius
         << endl;
    delete[] p;
    return 0;
}