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

Posted 2020-12-31 flyerbird

Ultimate Weapon

Time Limit: 2000MS		Memory Limit: 131072K
Total Submissions: 2430		Accepted: 1173

Description

In year 2008 of the Cosmic Calendar, the Aliens send a huge armada towards the Earth seeking after conquest. The humans now depend on their ultimate weapon to retain their last hope of survival. The weapon, while capable of creating a continuous, closed and convex lethal region in the space and annihilating everything enclosed within, unfortunately exhausts upon each launch a tremendous amount of energy which is proportional to the surface area of the lethal region.

Given the positions of all battleships in the Aliens‘ armada, your task is to calculate the minimum amount of energy required to destroy the armada with a single launch of the ultimate weapon. You need to report the surface area of the lethal region only.

Input

The first line contains one number N -- the number of battleships.(1 ≤ N ≤ 500) 
Following N lines each contains three integers presenting the position of one battleship.

Output

The minimal area rounded to three decimal places.

Sample Input

4
0 0 0
4 0 0
2 3 0
1 1 2

Sample Output

19.137

Hint

There are no four coplaner battleships.

Source

POJ Founder Monthly Contest – 2008.03.16, Jiang Liyang

 

三维凸包裸题，求这些点组成的凸包的表面积，数据最大为500，用增量算法，且由于只求表面积，可以不用合并同一平面的三角形。

 

  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
", hull.area());
173     }
174     return 0;
175 }

View Code