求点集的
最小包围球可以使用 Welzl
算法（也称为 Randomized Incremental Algorithm），该
算法的时间复杂度为 O(n)。
C++ 代码如下：
```c++
#include <iostream>
#include <vector>
#include <random>
#include <chrono>
#include <cmath>
using namespace std;
struct Point {
double x, y, z;
};
double dist(const Point& p1, const Point& p2) {
return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2));
}
double sphere_radius(const Point& p1, const Point& p2, const Point& p3) {
double a = dist(p1, p2);
double b = dist(p2, p3);
double c = dist(p3, p1);
double s = (a + b + c) / 2.0;
double r = a * b * c / (4.0 * sqrt(s * (s - a) * (s - b) * (s - c)));
return r;
}
bool point_in_sphere(const Point& p, const Point& center, double radius) {
return dist(p, center) <= radius;
}
Point sphere_center(const Point& p1, const Point& p2, const Point& p3) {
double A1 = p2.x - p1.x;
double B1 = p2.y - p1.y;
double C1 = p2.z - p1.z;
double A2 = p3.x - p1.x;
double B2 = p3.y - p1.y;
double C2 = p3.z - p1.z;
double D1 = A1 * (p1.x + p2.x) + B1 * (p1.y + p2.y) + C1 * (p1.z + p2.z);
double D2 = A2 * (p1.x + p3.x) + B2 * (p1.y + p3.y) + C2 * (p1.z + p3.z);
double det = 2.0 * (A1 * B2 * C2 + B1 * C2 * A2 + C1 * A2 * B2 - A2 * B1 * C1 - B2 * C1 * A1 - C2 * A1 * B1);
double x = (D1 * B2 * C2 + B1 * C2 * D2 + C1 * D2 * B2 - D2 * B1 * C1 - C2 * B1 * D1 - C1 * D1 * B2) / det;
double y = (A1 * D2 * C2 + D1 * C2 * A2 + C1 * A2 * D2 - A2 * C1 * D1 - D2 * C1 * A1 - A2 * D1 * C2) / det;
double z = (A1 * B2 * D2 + B1 * D1 * A2 + D1 * B2 * C1 - A2 * B1 * D1 - B2 * D2 * C1 - C2 * A1 * D2) / det;
return {x, y, z};
}
bool min_sphere_helper(const vector<Point>& P, const vector<Point>& R, Point& center, double& radius) {
if (R.size() == 0 || P.size() == 0) {
if (R.size() == 2) {
center = {(R[0].x + R[1].x) / 2.0, (R[0].y + R[1].y) / 2.0, (R[0].z + R[1].z) / 2.0};
radius = dist(R[0], R[1]) / 2.0;
return true;
}
else if (R.size() == 3) {
center = sphere_center(R[0], R[1], R[2]);
radius = sphere_radius(R[0], R[1], R[2]);
return true;
}
else {
return false;
}
}
int idx = rand() % P.size();
Point p = P[idx];
P.erase(P.begin() + idx);
if (min_sphere_helper(P, R, center, radius)) {
if (point_in_sphere(p, center, radius)) {
return true;
}
R.push_back(p);
return min_sphere_helper(P, R, center, radius);
}
R.push_back(p);
return min_sphere_helper(P, R, center, radius);
}
pair<Point, double> min_sphere(vector<Point> P) {
random_device rd;
srand(rd());
random_shuffle(P.begin(), P.end());
Point center;
double radius;
min_sphere_helper(P, {}, center, radius);
return {center, radius};
}
int main() {
vector<Point> points = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {-1, 0, 0}, {0, -1, 0}, {0, 0, -1}};
auto [center, radius] = min_sphere(points);
cout << "Center: (" << center.x << ", " << center.y << ", " << center.z << "), Radius: " << radius << endl;
return 0;
}
```
该代码使用了 C++11 的新特性，如结构体解构（`auto [center, radius] = min_sphere(points);`）和 if 语句中的布尔表达式（`if (min_sphere_helper(P, R, center, radius)) {...}`）。
最低0.47元/天 解锁文章
687
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
