POJ 3528--Ultimate Weapon(三维凸包)

  1 #include<iostream>
  2 #include<algorithm>
  3 #include<cmath>
  4 #include<cstdio>
  5 using namespace std;
  6 const int N = 510;
  7 const double eps = 1e-8;
  8 typedef struct point3 {
  9     double x, y, z;
 10     point3() {    }
 11     point3(double a, double b, double c) :x(a), y(b), z(c) {  }
 12     point3 operator -(const point3 &b)const {        //返回减去后的新点
 13         return point3(x - b.x, y - b.y, z - b.z);
 14     }
 15     point3 operator +(const point3 &b)const {        //返回加上后的新点
 16         return point3(x + b.x, y + b.y, z + b.z);
 17     }
 18     //数乘计算
 19     point3 operator *(const double &k)const {    //返回相乘后的新点
 20         return point3(x * k, y * k, z*k);
 21     }
 22     point3 operator /(const double &k)const {    //返回相除后的新点
 23         return point3(x / k, y / k, z / k);
 24     }
 25     double operator *(const point3 &b)const {    //点乘
 26         return (x*b.x + y*b.y + z*b.z);
 27     }
 28     point3 operator ^(const point3 &p) const {    //叉积
 29         return point3(y*p.z - p.y*z, z*p.x - x*p.z, x*p.y - y*p.x);
 30     }
 31     double vlen()const {    //向量的模
 32         return sqrt(x*x + y*y + z*z);
 33     }
 34 }point3;
 35 struct fac {
 36     int a, b, c;//凸包一个面上的三个点的编号
 37     bool ok;    //该面是否是最终凸包中的面
 38 };
 39 struct T3dhull {
 40     int n;                    //初始点数
 41     point3 ply[N];            //初始点
 42     int trianglecnt;        //凸包上三角形数
 43     fac tri[N];                //凸包三角形可证明被创建的面不会超过6N
 44     int vis[N][N];            //点i到点j是属于哪个面
 45     double dist(point3 a) { return sqrt(a.x*a.x + a.y*a.y + a.z*a.z); }            //两点长度
 46     double area(point3 a, point3 b, point3 c){return dist((b - a) ^ (c - a));}  //三角形面积*2
 47 
 48     //返回四面体有向体积*6
 49     //在储存面时，保证面的法线方向朝向凸包外部，如果在某一平面和点p所组成的四面体的有向体积为正，则p点在凸包外部，并且此点可以被p点看见。
 50     double volume(point3 a, point3 b, point3 c, point3 d){
 51         return ((b - a) ^ (c - a))* (d - a);
 52     }
 53     double ptoplane(point3 &p, fac &f){     //点到平面距离，体积法
 54         point3 m = ply[f.b] - ply[f.a], n = ply[f.c] - ply[f.a], t = p - ply[f.a];
 55         return (m^n) * t;
 56     }
 57     void deal(int p, int a, int b) {
 58         int f = vis[a][b];
 59         fac add;
 60         if (tri[f].ok)
 61         {
 62             if ((ptoplane(ply[p], tri[f])) > eps)
 63                 dfs(p, f);
 64             else
 65             {
 66                 add.a = b, add.b = a, add.c = p, add.ok = 1;
 67                 vis[p][b] = vis[a][p] = vis[b][a] = trianglecnt;
 68                 tri[trianglecnt++] = add;
 69             }
 70         }
 71     }
 72     void dfs(int p, int cnt) {//维护凸包，如果点p在凸包外侧则更新凸包
 73         tri[cnt].ok = 0;
 74         deal(p, tri[cnt].b, tri[cnt].a);
 75         deal(p, tri[cnt].c, tri[cnt].b);
 76         deal(p, tri[cnt].a, tri[cnt].c);
 77     }
 78     bool same(int s, int e) {
 79         point3 a = ply[tri[s].a], b = ply[tri[s].b], c = ply[tri[s].c];
 80         return fabs(volume(a, b, c, ply[tri[e].a])) < eps
 81             && fabs(volume(a, b, c, ply[tri[e].b])) < eps
 82             && fabs(volume(a, b, c, ply[tri[e].c])) < eps;
 83     }
 84     void construct()//构造凸包
 85     {
 86         int i, j;
 87         trianglecnt = 0;
 88         if (n<4) return;
 89         bool tmp = true;
 90         for (i = 1; i < n; i++)    //前两点不共点
 91         {
 92             if ((dist(ply[0] - ply[i])) > eps)
 93             {
 94                 swap(ply[1], ply[i]);
 95                 tmp = false;
 96                 break;
 97             }
 98         }
 99         if (tmp)return;
100         tmp = true;
101         for (i = 2; i < n; i++)        //前三点不共线
102         {
103             if ((dist((ply[0] - ply[1]) ^ (ply[1] - ply[i]))) > eps)
104             {
105                 swap(ply[2], ply[i]);
106                 tmp = false;
107                 break;
108             }
109         }
110         if (tmp) return;
111         tmp = true;
112         for (i = 3; i < n; i++)    //前四点不共面
113         {
114             if (fabs(((ply[0] - ply[1]) ^ (ply[1] - ply[2]))* (ply[0] - ply[i]))>eps)
115             {
116                 swap(ply[3], ply[i]);
117                 tmp = false;
118                 break;
119             }
120         }
121         if (tmp)return;
122         fac add;
123         for (i = 0; i < 4; i++)    //构建初始四面体
124         {
125             add.a = (i + 1) % 4, add.b = (i + 2) % 4, add.c = (i + 3) % 4, add.ok = 1;
126             if ((ptoplane(ply[i], add))>0)
127                 swap(add.b, add.c);
128             vis[add.a][add.b] = vis[add.b][add.c] = vis[add.c][add.a] = trianglecnt;
129             tri[trianglecnt++] = add;
130         }
131         for (i = 4; i < n; i++)        //构建更新凸包
132         {
133             for (j = 0; j < trianglecnt; j++)
134             {
135                 if (tri[j].ok && (ptoplane(ply[i], tri[j])) > eps)
136                 {
137                     dfs(i, j); break;
138                 }
139             }
140         }
141         int cnt = trianglecnt;
142         trianglecnt = 0;
143         for (i = 0; i < cnt; i++)
144         {
145             if (tri[i].ok)
146                 tri[trianglecnt++] = tri[i];
147         }
148     }
149     double area()        //表面积
150     {
151         double ret = 0;
152         for (int i = 0; i < trianglecnt; i++)
153             ret += area(ply[tri[i].a], ply[tri[i].b], ply[tri[i].c]);
154         return ret / 2.0;
155     }
156     double volume()
157     {
158         point3 p(0, 0, 0);
159         double ret = 0;
160         for (int i = 0; i < trianglecnt; i++)
161             ret += volume(p, ply[tri[i].a], ply[tri[i].b], ply[tri[i].c]);
162         return fabs(ret / 6);
163     }
164 }hull;
165 
166 int main() {
167     while (~scanf("%d", &hull.n)) {
168         int i;
169         for (i = 0; i < hull.n; i++)
170             scanf("%lf %lf %lf", &hull.ply[i].x, &hull.ply[i].y, &hull.ply[i].z);
171         hull.construct();
172         printf("%.3lf\n", hull.area());
173     }
174     return 0;
175 }