POJ3528 (三维凸包模板)

                              Ultimate Weapon

 POJ - 3528 

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

40 0 04 0 02 3 01 1 2

Sample Output

19.137

Hint

There are no four coplaner battleships.

这是一道裸的三维凸包模板题，直接套用模板就行了，下面给出三维凸包的超强模板。

#include <cstdio>#include <cstring>#include <cmath>#include <algorithm>using namespace std;const double eps = 1e-8;const int maxn = 505;struct Point {double x, y, z;Point(double x = 0, double y = 0, double z = 0) :x(x), y(y), z(z) {}Point operator+(const Point &t) const {return Point(x + t.x, y + t.y, z + t.z);}Point operator-(const Point &t) const {return Point(x - t.x, y - t.y, z - t.z);}Point operator*(const Point &t) const {return Point(y*t.z - z*t.y, z*t.x - x*t.z, x*t.y - y*t.x);}double operator^(const Point &t) const {return  x*t.x + y*t.y + z*t.z;}Point operator*(const double &t) const {return Point(x*t, y*t, z*t);}Point operator/(const double &t) const {return Point(x / t, y / t, z / t);}void in() {scanf("%lf%lf%lf", &x, &y, &z);}};double vlen(Point t) {return sqrt(t.x*t.x + t.y*t.y + t.z*t.z);}struct CH3D {struct face {int a, b, c;bool ok;face(int a, int b, int c) :a(a), b(b), c(c), ok(1) {}face() {}}f[8 * maxn];int cnt, n, to[maxn][maxn];Point p[maxn];bool fix() { //为了保证前四个点不公面, 若有保证就可以不写fix函数int i;bool ok = true;for (i = 1; i < n; i++) //使前两点不公点if (vlen(p[i] - p[0]) > eps) {swap(p[i], p[1]);ok = false;break;}if (ok) return 0; ok = true;for (i = 2; i < n; i++)  //使前三点不公线if (vlen((p[1] - p[0])*(p[i] - p[0])) > eps) {swap(p[i], p[2]);ok = false;break;}if (ok) return 0; ok = true;for (i = 3; i < n; i++)  //使前四点不共面if (fabs((p[1] - p[0])*(p[2] - p[0]) ^ (p[i] - p[0])) > eps) {swap(p[i], p[3]);ok = false;break;}if (ok) return 0;return 1;}double ptoface(Point &t, face &f) { //判点是否在面的同侧，>0 同侧return volume(p[f.a], p[f.b], p[f.c], t);}void dfs(int i, int j) {f[j].ok = 0;deal(i, f[j].b, f[j].a);deal(i, f[j].c, f[j].b);deal(i, f[j].a, f[j].c);}void deal(int i, int a, int b) {int j = to[a][b];if (f[j].ok) {if (ptoface(p[i], f[j]) > eps)dfs(i, j);else {face add(b, a, i);to[b][a] = to[a][i] = to[i][b] = cnt;f[cnt++] = add;}}}void creat() { //构造凸包if (n < 4 || !fix()) return;int i, j;cnt = 0;for (i = 0; i < 4; i++) {face add((i + 1) % 4, (i + 2) % 4, (i + 3) % 4);if (ptoface(p[i], add) > eps)swap(add.b, add.c);to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;f[cnt++] = add;}for (i = 4; i < n; i++) {for (j = 0; j < cnt; j++)if (f[j].ok && ptoface(p[i], f[j]) > eps) {dfs(i, j);break;}}int t = cnt; cnt = 0;for (i = 0; i < t; i++)if (f[i].ok) f[cnt++] = f[i];}double area() { //凸包表面积double ret = 0.0;for (int i = 0; i < cnt; i++)ret += area(p[f[i].a], p[f[i].b], p[f[i].c]);return ret*0.5;}double volume() { //凸包体积Point o;  //o是原点double ret = 0.0;for (int i = 0; i < cnt; i++)ret += volume(o, p[f[i].a], p[f[i].b], p[f[i].c]);return fabs(ret / 6.0);}bool same(int s, int t) { // 判两个平面是否为同一平面Point &a = p[f[s].a], &b = p[f[s].b], &c = p[f[s].c];return fabs(volume(a, b, c, p[f[t].a])) < eps &&fabs(volume(a, b, c, p[f[t].b])) < eps &&fabs(volume(a, b, c, p[f[t].c])) < eps;}int faceCnt() {  //凸包的面数（除去相同的平面）int i, j;int ans = 0;for (i = 0; i < cnt; i++) {bool ok = 1;for (j = 0; j < i; j++) {if (same(i, j)) {ok = 0;break;}}ans += ok;}return ans;}double area(Point &a, Point &b, Point &c) { // 三角形面积*2;return vlen((b - a)*(c - a));}double volume(Point &a, Point &b, Point &c, Point &d) { // 四面体的有向面积*6;return (b - a)*(c - a) ^ (d - a);}}hull;int main() {int i;while (~scanf("%d", &hull.n)) {for (i = 0; i < hull.n; i++) hull.p[i].in();hull.creat();printf("%.3f\n", hull.area());}return 0;}