先决条件：给定圆上三个点时的圆方程，最小封闭圆。
给定一个数组arr[][]在具有整数坐标的二维平面中包含N个点。任务是找到最小包围圆（MEC）的中心和半径。最小封闭圆是所有点都位于圆内或其边界上的圆。
例子：
Input: arr[][] = {{0, 0}, {0, 1}, {1, 0}}
Output: Center = {0.5, 0.5}, Radius = 0.7071
Explanation:
On plotting the above circle with radius 0.707 and center (0.5, 0.5), it can be observed clearly that all the mentioned points lie either inside or on the circle.
Input: arr[][] = {{5, -2}, {-3, -2}, {-2, 5}, {1, 6}, {0, 2}}
Output: Center = {1.0, 1.0}, Radius = 5.000
方法：在上一篇文章中，讨论了一种朴素的方法和一种优化的方法，即首先使点集充满凸面，然后再执行朴素的方法。尽管优化的解决方案对于某些输入非常有效，但优化后的最坏情况时间复杂度仍然是O(N 4 ) 。在本文中，已经讨论了一种优化方法。
这个想法是使用 Welzl 的递归算法。使用该算法，可以在O(N) 中找到 MEC。算法的工作取决于从上一篇文章中得出的观察和结论。该算法的思想是从给定的输入集中随机去除一个点，形成一个圆方程。一旦方程成立，检查被删除的点是否被方程包围。如果不是，则该点必须位于 MEC 的边界上。因此，该点被视为边界点并递归调用该函数。该算法的具体工作如下：
该算法将一组点P和一组R初始为空并用于表示 MEC 边界上的点作为输入。
该算法的基本情况是当P变为空或集合R的大小等于3 时：
- 如果 P 为空，则所有点都已处理。
- 如果|R| = 3 ，那么已经找到 3 个点位于圆的边界上，并且由于仅使用 3 个点就可以唯一确定一个圆，因此可以停止递归。
当算法达到上述基本情况时，它返回 R 的平凡解，即：
- 如果|R| = 1 ，我们返回一个以 R[0] 为中心、半径 = 0 的圆
- 如果|R| = 2 ，我们返回 R[0] 和 R[2] 的 MEC
- 如果|R| = 3 ，我们通过尝试 3 对 (R[0], R[1]), (R[0], R[2]), (R[1], R[2]) 来返回 MEC
- 如果这些对都无效，我们返回由 R 中的 3 个点定义的圆
如果尚未达到基本情况，我们执行以下操作：
- 从 P 中选取一个随机点 p 并将其从 P 中删除
- 在 P 和 R 上调用算法得到圆 d
- 如果 p 被 d 包围，那么我们返回 d
- 否则，p 必须位于 MEC 的边界上
- 将 p 添加到 R
- 返回算法在 P 和 R 上的输出
下面是上述方法的实现：
// C++ program to find the minimum enclosing
// circle for N integer points in a 2-D plane
#include
#include
#include
#include
#include
using namespace std;
// Defining infinity
const double INF = 1e18;
// Structure to represent a 2D point
struct Point {
double X, Y;
};
// Structure to represent a 2D circle
struct Circle {
Point C;
double R;
};
// Function to return the euclidean distance
// between two points
double dist(const Point& a, const Point& b)
{
return sqrt(pow(a.X - b.X, 2)
+ pow(a.Y - b.Y, 2));
}
// Function to check whether a point lies inside
// or on the boundaries of the circle
bool is_inside(const Circle& c, const Point& p)
{
return dist(c.C, p) <= c.R;
}
// The following two functions are used
// To find the equation of the circle when
// three points are given.
// Helper method to get a circle defined by 3 points
Point get_circle_center(double bx, double by,
double cx, double cy)
{
double B = bx * bx + by * by;
double C = cx * cx + cy * cy;
double D = bx * cy - by * cx;
return { (cy * B - by * C) / (2 * D),
(bx * C - cx * B) / (2 * D) };
}
// Function to return a unique circle that
// intersects three points
Circle circle_from(const Point& A, const Point& B,
const Point& C)
{
Point I = get_circle_center(B.X - A.X, B.Y - A.Y,
C.X - A.X, C.Y - A.Y);
I.X += A.X;
I.Y += A.Y;
return { I, dist(I, A) };
}
// Function to return the smallest circle
// that intersects 2 points
Circle circle_from(const Point& A, const Point& B)
{
// Set the center to be the midpoint of A and B
Point C = { (A.X + B.X) / 2.0, (A.Y + B.Y) / 2.0 };
// Set the radius to be half the distance AB
return { C, dist(A, B) / 2.0 };
}
// Function to check whether a circle
// encloses the given points
bool is_valid_circle(const Circle& c,
const vector& P)
{
// Iterating through all the points
// to check whether the points
// lie inside the circle or not
for (const Point& p : P)
if (!is_inside(c, p))
return false;
return true;
}
// Function to return the minimum enclosing
// circle for N <= 3
Circle min_circle_trivial(vector& P)
{
assert(P.size() <= 3);
if (P.empty()) {
return { { 0, 0 }, 0 };
}
else if (P.size() == 1) {
return { P[0], 0 };
}
else if (P.size() == 2) {
return circle_from(P[0], P[1]);
}
// To check if MEC can be determined
// by 2 points only
for (int i = 0; i < 3; i++) {
for (int j = i + 1; j < 3; j++) {
Circle c = circle_from(P[i], P[j]);
if (is_valid_circle(c, P))
return c;
}
}
return circle_from(P[0], P[1], P[2]);
}
// Returns the MEC using Welzl's algorithm
// Takes a set of input points P and a set R
// points on the circle boundary.
// n represents the number of points in P
// that are not yet processed.
Circle welzl_helper(vector& P,
vector R, int n)
{
// Base case when all points processed or |R| = 3
if (n == 0 || R.size() == 3) {
return min_circle_trivial(R);
}
// Pick a random point randomly
int idx = rand() % n;
Point p = P[idx];
// Put the picked point at the end of P
// since it's more efficient than
// deleting from the middle of the vector
swap(P[idx], P[n - 1]);
// Get the MEC circle d from the
// set of points P - {p}
Circle d = welzl_helper(P, R, n - 1);
// If d contains p, return d
if (is_inside(d, p)) {
return d;
}
// Otherwise, must be on the boundary of the MEC
R.push_back(p);
// Return the MEC for P - {p} and R U {p}
return welzl_helper(P, R, n - 1);
}
Circle welzl(const vector& P)
{
vector P_copy = P;
random_shuffle(P_copy.begin(), P_copy.end());
return welzl_helper(P_copy, {}, P_copy.size());
}
// Driver code
int main()
{
Circle mec = welzl({ { 0, 0 },
{ 0, 1 },
{ 1, 0 } });
cout << "Center = { " << mec.C.X << ", " << mec.C.Y
<< " } Radius = " << mec.R << endl;
Circle mec2 = welzl({ { 5, -2 },
{ -3, -2 },
{ -2, 5 },
{ 1, 6 },
{ 0, 2 } });
cout << "Center = { " << mec2.C.X << ", " << mec2.C.Y
<< " } Radius = " << mec2.R << endl;
return 0;
}
输出：
Center = { 0.5, 0.5 } Radius = 0.707107
Center = { 1, 1 } Radius = 5
时间复杂度：该算法的预期时间和空间复杂度为O(N) ，其中 N 是点数。空间是由于正在使用递归的事实。要理解为什么时间复杂度是线性的，我们可以观察不同状态的数量，以了解递归函数可以发生多少次调用。每次调用时， P的大小都会减少1 。此外，R 的大小可以保持不变或可以增加 1。由于|R|不能超过3 ，那么不同状态的数量将是3N 。因此，这使得时间复杂度为O(N) 。