先决条件：给定圆上的三个点时的圆方程，最小封闭圆。
给定一个数组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) 。