welzl.cpp

// Minimum enclosing circle, Welzl's algorithm
// Expected linear time.
// If there are any duplicate points in the input, be sure to remove them first.
struct point {
	double x;
	double y;
};
struct circle {
	double x;
	double y;
	double r;
	circle() {}
	circle(double x, double y, double r): x(x), y(y), r(r) {}
};
circle b_md(vector<point> R) {
	if (R.size() == 0) {
		return circle(0, 0, -1);
	} else if (R.size() == 1) {
		return circle(R[0].x, R[0].y, 0);
	} else if (R.size() == 2) {
		return circle((R[0].x+R[1].x)/2.0,
		              (R[0].y+R[1].y)/2.0,
					  hypot(R[0].x-R[1].x, R[0].y-R[1].y)/2.0);
	} else {
		double D = (R[0].x - R[2].x)*(R[1].y - R[2].y) - (R[1].x - R[2].x)*(R[0].y - R[2].y);
		double p0 = (((R[0].x - R[2].x)*(R[0].x + R[2].x) + (R[0].y - R[2].y)*(R[0].y + R[2].y)) / 2 * (R[1].y - R[2].y) - ((R[1].x - R[2].x)*(R[1].x + R[2].x) + (R[1].y - R[2].y)*(R[1].y + R[2].y)) / 2 * (R[0].y - R[2].y))/D;
		double p1 = (((R[1].x - R[2].x)*(R[1].x + R[2].x) + (R[1].y - R[2].y)*(R[1].y + R[2].y)) / 2 * (R[0].x - R[2].x) - ((R[0].x - R[2].x)*(R[0].x + R[2].x) + (R[0].y - R[2].y)*(R[0].y + R[2].y)) / 2 * (R[1].x - R[2].x))/D;
		return circle(p0, p1, hypot(R[0].x - p0, R[0].y - p1));
	}
}
circle b_minidisk(vector<point>& P, int i, vector<point> R) {
	if (i == P.size() || R.size() == 3) {
		return b_md(R);
	} else {
		circle D = b_minidisk(P, i+1, R);
		if (hypot(P[i].x-D.x, P[i].y-D.y) > D.r) {
			R.push_back(P[i]);
			D = b_minidisk(P, i+1, R);
		}
		return D;
	}
}

// Call this function.
circle minidisk(vector<point> P) {
	random_shuffle(P.begin(), P.end());
	return b_minidisk(P, 0, vector<point>());
}