poj 3862 && LA 4589 Asteroids （三维凸包+多面体重心）
用给出的点求出凸包的重心，并求出重心到多边形表面的最近距离。
代码如下：
1 #include <cstdio> 2 #include <iostream> 3 #include <algorithm> 4 #include <cstring> 5 #include <cmath> 6 7 using namespace std; 8 9 const double EPS = 1e-10; 10 const int N = 333; 11 inline int sgn(double x) { return (x > EPS) - (x < -EPS);} 12 struct Point { 13 double x, y, z; 14 Point() {} 15 Point(double x, double y, double z) : x(x), y(y), z(z) {} 16 bool operator < (Point a) const { 17 if (sgn(x - a.x)) return x < a.x; 18 if (sgn(y - a.y)) return y < a.y; 19 return z < a.z; 20 } 21 bool operator == (Point a) const { return !sgn(x - a.x) && !sgn(y - a.y) && !sgn(z - a.z);} 22 Point operator + (Point a) { return Point(x + a.x, y + a.y, z + a.z);} 23 Point operator - (Point a) { return Point(x - a.x, y - a.y, z - a.z);} 24 Point operator * (double p) { return Point(x * p, y * p, z * p);} 25 Point operator / (double p) { return Point(x / p, y / p, z / p);} 26 } ; 27 inline double dot(Point a, Point b) { return a.x * b.x + a.y * b.y + a.z * b.z;} 28 inline double dot(Point o, Point a, Point b) { return dot(a - o, b - o);} 29 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);} 30 inline Point cross(Point o, Point a, Point b) { return cross(a - o, b - o);} 31 inline double veclen(Point x) { return sqrt(dot(x, x));} 32 inline double area(Point o, Point a, Point b) { return veclen(cross(o, a, b)) / 2.0;} 33 inline double volume(Point o, Point a, Point b, Point c) { return dot(cross(o, a, b), c - o) / 6.0;} 34 35 struct _3DCH{ 36 Point pt[N]; 37 int pcnt; 38 struct Face { 39 int a, b, c; 40 bool ok; 41 Face() {} 42 Face(int a, int b, int c) : a(a), b(b), c(c) { ok = true;} 43 } fc[N << 4]; 44 int fcnt; 45 int fid[N][N]; 46 47 bool outside(int p, int a, int b, int c) { return sgn(volume(pt[a], pt[b], pt[c], pt[p])) < 0;} 48 bool outside(int p, int f) { return outside(p, fc[f].a, fc[f].b, fc[f].c);} 49 void addface(int a, int b, int c, int d) { 50 if (outside(d, a, b, c)) fc[fcnt] = Face(c, b, a), fid[c][b] = fid[b][a] = fid[a][c] = fcnt++; 51 else fc[fcnt] = Face(a, b, c), fid[a][b] = fid[b][c] = fid[c][a] = fcnt++; 52 } 53 54 bool dfs(int p, int f) { 55 if (!fc[f].ok) return true; 56 int a = fc[f].a, b = fc[f].b, c = fc[f].c; 57 if (outside(p, a, b, c)) { 58 fc[f].ok = false; 59 if (!dfs(p, fid[b][a])) addface(p, a, b, c); 60 if (!dfs(p, fid[c][b])) addface(p, b, c, a); 61 if (!dfs(p, fid[a][c])) addface(p, c, a, b); 62 return true; 63 } 64 return false; 65 } 66 void build() { 67 bool ok; 68 fcnt = 0; 69 sort(pt, pt + pcnt); 70 pcnt = unique(pt, pt + pcnt) - pt; 71 ok = false; 72 for (int i = 2; i < pcnt; i++) { 73 if (sgn(area(pt[0], pt[1], pt[i]))) { 74 ok = true; 75 swap(pt[i], pt[2]); 76 break; 77 } 78 } 79 if (!ok) return ; 80 ok = false; 81 for (int i = 3; i < pcnt; i++) { 82 if (sgn(volume(pt[0], pt[1], pt[2], pt[i]))) { 83 ok = true; 84 swap(pt[i], pt[3]); 85 break; 86 } 87 } 88 if (!ok) return ; 89 addface(0, 1, 2, 3); 90 addface(1, 2, 3, 0); 91 addface(2, 3, 0, 1); 92 addface(3, 0, 1, 2); 93 for (int i = 4; i < pcnt; i++) { 94 for (int j = fcnt - 1; j >= 0; j--) { 95 if (outside(i, j)) { 96 dfs(i, j); 97 break; 98 } 99 } 100 } 101 int sz = fcnt; 102 fcnt = 0; 103 for (int i = 0; i < sz; i++) if (fc[i].ok) fc[fcnt++] = fc[i]; 104 } 105 double facearea() { 106 double ret = 0.0; 107 for (int i = 0; i < fcnt; i++) ret += area(pt[fc[i].a], pt[fc[i].b], pt[fc[i].c]); 108 return ret; 109 } 110 bool same(int i, int j) { 111 int a = fc[i].a, b = fc[i].b, c = fc[i].c; 112 if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].a]))) return false; 113 if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].b]))) return false; 114 if (sgn(volume(pt[a], pt[b], pt[c], pt[fc[j].c]))) return false; 115 return true; 116 } 117 int facecnt() { 118 int cnt = 0; 119 for (int i = 0; i < fcnt; i++) { 120 bool ok = true; 121 for (int j = 0; j < i; j++) { 122 if (same(i, j)) { 123 ok = false; 124 break; 125 } 126 } 127 if (ok) cnt++; 128 } 129 return cnt; 130 } 131 Point gravity() { 132 Point ret = Point(0.0, 0.0, 0.0); 133 Point O = Point(0.0, 0.0, 0.0); 134 double vol = 0.0; 135 for (int i = 0; i < fcnt; i++) { 136 double tmp = volume(pt[fc[i].a], pt[fc[i].b], pt[fc[i].c], O); 137 ret = ret + (pt[fc[i].a] + pt[fc[i].b] + pt[fc[i].c]) / 4.0 * tmp; 138 vol += tmp; 139 } 140 // cout << ret.x << ' ' << ret.y << ' ' << ret.z << endl; 141 return ret / vol; 142 } 143 } ch; 144 145 inline double pt2plane(Point o, Point a, Point b, Point c) { return 3.0 * fabs(volume(o, a, b, c) / area(a, b, c));} 146 147 double getMin() { 148 double ret = 1e100; 149 Point g = ch.gravity(); 150 for (int i = 0; i < ch.fcnt; i++) { 151 ret = min(ret, pt2plane(g, ch.pt[ch.fc[i].a], ch.pt[ch.fc[i].b], ch.pt[ch.fc[i].c])); 152 } 153 // cout << ret << endl; 154 return ret; 155 } 156 157 int main() { 158 // freopen("in", "r", stdin); 159 double x, y, z; 160 while (true) { 161 double ans = 0.0; 162 for (int i = 0; i < 2; i++) { 163 if (scanf("%d", &ch.pcnt) == -1) return 0; 164 for (int j = 0; j < ch.pcnt; j++) { 165 scanf("%lf%lf%lf", &x, &y, &z); 166 ch.pt[j] = Point(x, y, z); 167 } 168 ch.build(); 169 ans += getMin(); 170 } 171 printf("%.8f\n", ans); 172 } 173 return 0; 174 }
——written by Lyon
