poj 3528 Ultimate Weapon (3D Convex Hull)

  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 = 555;
 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(a - o, b - o)) / 2.0;}
 33 // directed volume
 34 inline double volume(Point o, Point a, Point b, Point c) { return dot(cross(a - o, b - o), c - o);}
 35 
 36 struct _3DCH {
 37     Point pt[N];
 38     int pcnt;
 39     struct Face { // follow right-hand
 40         int a, b, c;
 41         bool ok;
 42         Face() {}
 43         Face(int a, int b, int c) : a(a), b(b), c(c) { ok = true;}
 44     } ;
 45     Face f[N << 3];
 46     int fcnt, fid[N][N];
 47     bool outside(int p, int a, int b, int c) {
 48         int dir = sgn(volume(pt[a], pt[b], pt[c], pt[p]));
 49         return dir < 0;
 50     }
 51     bool outside(int p, int x) { return outside(p, f[x].a, f[x].b, f[x].c);}
 52     void addface(int p, int a, int b, int c) {
 53         int dir = sgn(volume(pt[a], pt[b], pt[c], pt[p]));
 54         if (dir == 0) return ;
 55         if (dir > 0) {
 56             f[fcnt] = Face(a, b, c);
 57             fid[a][b] = fid[b][c] = fid[c][a] = fcnt++;
 58         } else {
 59             f[fcnt] = Face(c, b, a);
 60             fid[c][b] = fid[b][a] = fid[a][c] = fcnt++;
 61         }
 62     }
 63     bool dfs(int p, int x) {
 64         if (!f[x].ok) return true;
 65         if (outside(p, x)) {
 66             f[x].ok = false;
 67             if (!dfs(p, fid[f[x].b][f[x].a])) addface(f[x].c, f[x].b, f[x].a, p);
 68             if (!dfs(p, fid[f[x].c][f[x].b])) addface(f[x].a, f[x].c, f[x].b, p);
 69             if (!dfs(p, fid[f[x].a][f[x].c])) addface(f[x].b, f[x].a, f[x].c, p);
 70             return true;
 71         }
 72         return false;
 73     }
 74     bool build() {
 75         bool ok;
 76         fcnt = 0;
 77         sort(pt, pt + pcnt);
 78         pcnt = unique(pt, pt + pcnt) - pt;
 79         ok = false;
 80         for (int i = 2; i < pcnt; i++) {
 81             if (sgn(veclen(cross(pt[i], pt[0], pt[1])))) {
 82                 swap(pt[2], pt[i]);
 83                 ok = true;
 84                 break;
 85             }
 86         }
 87         if (!ok) return false;
 88         ok = false;
 89         for (int i = 3; i < pcnt; i++) {
 90             if (sgn(volume(pt[i], pt[0], pt[1], pt[2]))) {
 91                 swap(pt[3], pt[i]);
 92                 ok = true;
 93                 break;
 94             }
 95         }
 96         if (!ok) return false;
 97         addface(0, 1, 2, 3);
 98         addface(1, 2, 3, 0);
 99         addface(2, 3, 0, 1);
100         addface(3, 0, 1, 2);
101         for (int i = 4; i < pcnt; i++) {
102             for (int j = fcnt - 1; j >= 0; j--) {
103                 if (f[j].ok && outside(i, j)) {
104                     dfs(i, j);
105                     break;
106                 }
107             }
108         }
109         int sz = fcnt;
110         fcnt = 0;
111         for (int i = 0; i < sz; i++) if (f[i].ok) f[fcnt++] = f[i];
112         return true;
113     }
114     double vol() {
115         Point O = Point(0.0, 0.0, 0.0);
116         double ret = 0.0;
117         for (int i = 0; i < fcnt; i++) ret += volume(O, pt[f[i].a], pt[f[i].b], pt[f[i].c]);
118         return fabs(ret) / 6.0;
119     }
120     double facearea() {
121         double ret = 0.0;
122         for (int i = 0; i < fcnt; i++) ret += area(pt[f[i].a], pt[f[i].b], pt[f[i].c]);
123         return ret;
124     }
125 } ch;
126 
127 int main() {
128 //    freopen("in", "r", stdin);
129     while (cin >> ch.pcnt) {
130         for (int i = 0; i < ch.pcnt; i++) scanf("%lf%lf%lf", &ch.pt[i].x, &ch.pt[i].y, &ch.pt[i].z);
131         ch.build();
132         printf("%.3f
", ch.facearea());
133 =    }
134     return 0;
135 }