source: http://acm.hdu.edu.cn/showproblem.php?pid=3662
title : 3D Convex Hull
eclipse C++ 编译通过
#include <stdio.h>
#include <algorithm>
#include <cmath>
#include <list>
using namespace std;
const double eps = 1e-8;
typedef list<int>::iterator liit;
inline int sign(double d){
if(d < -eps) return -1;
return (d > eps) ? 1 : 0;
}
struct point{
double x, y, z;
point(double _x=0, double _y=0, double _z=0):x(_x), y(_y), z(_z) {}
void read(){ scanf("%lf%lf%lf", &x, &y, &z); }
bool operator==(const point& tp){
return (!sign(x-tp.x)) && (!sign(y-tp.y)) && (!sign(z-tp.z));
}
bool isOrg(){
return !sign(x) && !sign(y) && !sign(z);
}
point operator-(point tp){
return point(x-tp.x, y-tp.y, z-tp.z);
}
};
struct face{
list<int> pos;
point v; //法向量
double d; //常量
};
inline bool cmp(point p1, point p2){
return (p1.x < p2.x) || (p1.x == p2.x && p1.y < p2.y)
|| (p1.x == p2.x && p1.y == p2.y && p1.z < p2.z);
}
//向量st-->ed1和st-->ed2的叉积
inline point xmul3d(point st, point ed1, point ed2){
ed1 = ed1 - st;
ed2 = ed2 - st;
return point(ed1.y*ed2.z-ed2.y*ed1.z, ed1.z*ed2.x-ed2.z*ed1.x,
ed1.x*ed2.y-ed2.x*ed1.y);
}
//向量(0,0)-->p1和(0,0)-->ed2的叉积
inline double dmul3d(point p1, point p2){
return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z;
}
inline double dist3d(point p1, point p2){
return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y)
+ (p1.z-p2.z)*(p1.z-p2.z));
}
//三维凸包类
struct convex3d{
static const int N = 305; //点的数目的最大值
point ps[N]; //求解凸包的点集
list<face> convex; //convex存储三维的凸包的各个面,这些面是三维空间的凸多边形
int pn; //点的个数
int fsign[N][N];
//添加由ps里的第a,b,c个点组成的面
void addFace(int a, int b, int c){
face tf;
tf.pos.push_back(a);
tf.pos.push_back(b);
tf.pos.push_back(c);
tf.v = xmul3d(ps[a], ps[b], ps[c]);
tf.d = -dmul3d(ps[a], tf.v);
convex.push_back(tf);
}
//插入第i个点是,对面f进行处理
int handleFace(face f, int i){
int s = sign(dmul3d(f.v, ps[i]) + f.d);
liit now, nxt;
now = f.pos.begin();
nxt = f.pos.begin();
nxt++;
for(; nxt != f.pos.end(); nxt++, now++){
fsign[*now][*nxt] = s;
}
fsign[*now][*f.pos.begin()] = s;
return s;
}
//判断第i个点是否在f里面
bool inFace(face f, int i){
liit now, nxt;
int pn, nn, s;
now = f.pos.begin();
nxt = f.pos.begin();
nxt++;
for(pn = nn = 0; nxt != f.pos.end(); nxt++, now++){
s = sign(dmul3d(xmul3d(ps[*now], ps[*nxt], ps[i]), f.v));
if(s == 1) pn++;
else if(s == -1) nn++;
}
s = sign(dmul3d(xmul3d(ps[*now], ps[*f.pos.begin()], ps[i]), f.v));
if(s == 1) pn++;
else if(s == -1) nn++;
if(pn >= 1 && nn >= 1) return false;
return true;
}
//扩展面f,返回true表示需要删除当前的面
bool extFace(face& f, int i){
liit now, nxt;
bool flag = false;
now = f.pos.begin();
nxt = f.pos.begin();
nxt++;
if(fsign[*now][*nxt] == 0){
list<int> tpos;
while(true){
if(sign(dmul3d(xmul3d(ps[i], ps[*now], ps[*nxt]), f.v)) >= 1){
break;
}
now++;
if(now == f.pos.end()) now = f.pos.begin();
nxt++;
if(nxt == f.pos.end()) nxt = f.pos.begin();
}
tpos.push_back(*now);
int st = *now;
while(*nxt != st){
if(sign(dmul3d(xmul3d(ps[*now], ps[i], ps[*nxt]), f.v)) >= 0){
break;
}
tpos.push_back(*nxt);
now++;
if(now == f.pos.end()) now = f.pos.begin();
nxt++;
if(nxt == f.pos.end()) nxt = f.pos.begin();
}
tpos.push_back(i);
while(true){
now++;
if(now == f.pos.end()) now = f.pos.begin();
nxt++;
if(nxt == f.pos.end()) nxt = f.pos.begin();
if(*now == st) break;
if(sign(dmul3d(xmul3d(ps[i], ps[*now], ps[*nxt]), f.v)) >= 1){
tpos.push_back(*now);
}
}
f.pos = tpos;
}else if(fsign[*now][*nxt] > 0){
for(; nxt != f.pos.end(); now++, nxt++){
if(fsign[*nxt][*now] < 0){
addFace(*now, *nxt, i);
}
}
if(fsign[*f.pos.begin()][*now] < 0){
addFace(*now, *f.pos.begin(), i);
}
flag = true;
}
return flag;
}
//对ps里的pn个点求最小包围多面体,结果放在convex中
void initConvex(){
sort(ps, ps+pn, cmp);
pn = unique(ps, ps+pn) - ps;
convex.clear();
if(pn <= 2){
return;
}
int a, b, c;
double ab, bc, ac;
a = 0;
b = 1;
for(c = 2; c < pn; c++){
ab = dist3d(ps[a], ps[b]);
bc = dist3d(ps[b], ps[c]);
ac = dist3d(ps[a], ps[c]);
if(sign(ab+bc-ac) == 0){
b = c;
}else if(sign(ab+ac-bc) == 0){
a = c;
}else if(sign(ac+bc-ab) != 0){
break;
}
}
if(c == pn){
return;
}
int i, size, j;
list<face>::iterator it;
addFace(a, b, c);
addFace(a, c, b);
for(i = c+1; i < pn; i++){
size = convex.size();
for(it = convex.begin(), j = 0; j < size; j++, it++){
if(handleFace(*it, i) == 0 && inFace(*it, i)){
break;
}
}
if(j < size) continue;
for(it = convex.begin(), j = 0; j < size; j++){
if(extFace(*it, i)){
it = convex.erase(it);
}else{
it++;
}
}
}
}
};
convex3d c;
int main(){
#ifndef ONLINE_JUDGE
// freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
#endif
while(~scanf("%d", &c.pn)){
int i;
for(i = 0; i < c.pn; i++) c.ps[i].read();
c.initConvex();
printf("%d\n", c.convex.size());
}
return 0;
}