POJ 3528--Ultimate Weapon(三維凸包)
阿新 • • 發佈：2018-08-23
cep 一個 view swa acc esc 維護 3.1 lec 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
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 }
View Code
POJ 3528--Ultimate Weapon(三維凸包)