用给出的点求出凸包的重心，并求出重心到多边形表面的最近距离。
代码如下：
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath> using namespace std; const double EPS = 1e-;
const int N = ;
inline int sgn(double x) { return (x > EPS) - (x < -EPS);}
struct Point {
double x, y, z;
Point() {}
Point(double x, double y, double z) : x(x), y(y), z(z) {}
bool operator < (Point a) const {
if (sgn(x - a.x)) return x < a.x;
if (sgn(y - a.y)) return y < a.y;
return z < a.z;
}
bool operator == (Point a) const { return !sgn(x - a.x) && !sgn(y - a.y) && !sgn(z - a.z);}
Point operator + (Point a) { return Point(x + a.x, y + a.y, z + a.z);}
Point operator - (Point a) { return Point(x - a.x, y - a.y, z - a.z);}
Point operator * (double p) { return Point(x * p, y * p, z * p);}
Point operator / (double p) { return Point(x / p, y / p, z / p);}
} ;
inline double dot(Point a, Point b) { return a.x * b.x + a.y * b.y + a.z * b.z;}
inline double dot(Point o, Point a, Point b) { return dot(a - o, b - o);}
inline Point cross(Point a, Point b) { return Point(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);}
inline Point cross(Point o, Point a, Point b) { return cross(a - o, b - o);}
inline double veclen(Point x) { return sqrt(dot(x, x));}
inline double area(Point o, Point a, Point b) { return veclen(cross(o, a, b)) / 2.0;}
inline double volume(Point o, Point a, Point b, Point c) { return dot(cross(o, a, b), c - o) / 6.0;} struct _3DCH{
Point pt[N];
int pcnt;
struct Face {
int a, b, c;
bool ok;
Face() {}
Face(int a, int b, int c) : a(a), b(b), c(c) { ok = true;}
} fc[N << ];
int fcnt;
int fid[N][N]; bool outside(int p, int a, int b, int c) { return sgn(volume(pt[a], pt[b], pt[c], pt[p])) < ;}
bool outside(int p, int f) { return outside(p, fc[f].a, fc[f].b, fc[f].c);}
void addface(int a, int b, int c, int d) {
if (outside(d, a, b, c)) fc[fcnt] = Face(c, b, a), fid[c][b] = fid[b][a] = fid[a][c] = fcnt++;
else fc[fcnt] = Face(a, b, c), fid[a][b] = fid[b][c] = fid[c][a] = fcnt++;
} bool dfs(int p, int f) {
if (!fc[f].ok) return true;
int a = fc[f].a, b = fc[f].b, c = fc[f].c;
if (outside(p, a, b, c)) {
fc[f].ok = false;
if (!dfs(p, fid[b][a])) addface(p, a, b, c);
if (!dfs(p, fid[c][b])) addface(p, b, c, a);
if (!dfs(p, fid[a][c])) addface(p, c, a, b);
return true;
}
return false;
}
void build() {
bool ok;
fcnt = ;
sort(pt, pt + pcnt);
pcnt = unique(pt, pt + pcnt) - pt;
ok = false;
for (int i = ; i < pcnt; i++) {
if (sgn(area(pt[], pt[], pt[i]))) {
ok = true;
swap(pt[i], pt[]);
break;
}
}
if (!ok) return ;
ok = false;
for (int i = ; i < pcnt; i++) {
if (sgn(volume(pt[], pt[], pt[], pt[i]))) {
ok = true;
swap(pt[i], pt[]);
break;
}
}
if (!ok) return ;
addface(, , , );
addface(, , , );
addface(, , , );
addface(, , , );
for (int i = ; i < pcnt; i++) {
for (int j = fcnt - ; j >= ; j--) {
if (outside(i, j)) {
dfs(i, j);
break;
}
}
}
int sz = fcnt;
fcnt = ;
for (int i = ; i < sz; i++) if (fc[i].ok) fc[fcnt++] = fc[i];
}
double facearea() {
double ret = 0.0;
for (int i = ; i < fcnt; i++) ret += area(pt[fc[i].a], pt[fc[i].b], pt[fc[i].c]);
return ret;
}
bool same(int i, int j) {
int a = fc[i].a, b = fc[i].b, c = fc[i].c;
if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].a]))) return false;
if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].b]))) return false;
if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].c]))) return false;
return true;
}
int facecnt() {
int cnt = ;
for (int i = ; i < fcnt; i++) {
bool ok = true;
for (int j = ; j < i; j++) {
if (same(i, j)) {
ok = false;
break;
}
}
if (ok) cnt++;
}
return cnt;
}
Point gravity() {
Point ret = Point(0.0, 0.0, 0.0);
Point O = Point(0.0, 0.0, 0.0);
double vol = 0.0;
for (int i = ; i < fcnt; i++) {
double tmp = volume(pt[fc[i].a], pt[fc[i].b], pt[fc[i].c], O);
ret = ret + (pt[fc[i].a] + pt[fc[i].b] + pt[fc[i].c]) / 4.0 * tmp;
vol += tmp;
}
// cout << ret.x << ' ' << ret.y << ' ' << ret.z << endl;
return ret / vol;
}
} ch; inline double pt2plane(Point o, Point a, Point b, Point c) { return 3.0 * fabs(volume(o, a, b, c) / area(a, b, c));} double getMin() {
double ret = 1e100;
Point g = ch.gravity();
for (int i = ; i < ch.fcnt; i++) {
ret = min(ret, pt2plane(g, ch.pt[ch.fc[i].a], ch.pt[ch.fc[i].b], ch.pt[ch.fc[i].c]));
}
// cout << ret << endl;
return ret;
} int main() {
// freopen("in", "r", stdin);
double x, y, z;
while (true) {
double ans = 0.0;
for (int i = ; i < ; i++) {
if (scanf("%d", &ch.pcnt) == -) return ;
for (int j = ; j < ch.pcnt; j++) {
scanf("%lf%lf%lf", &x, &y, &z);
ch.pt[j] = Point(x, y, z);
}
ch.build();
ans += getMin();
}
printf("%.8f\n", ans);
}
return ;
}
——written by Lyon