POJ 3528 三维立体空间凸包模板
POJ 3528 三维凸包模板
#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) { //判点是否在面的同侧，> 同侧
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;
}