三维凸包之增量算法

三维凸包之增量算法

算法工作过程：

选取点集中的四个点构成一个初始的四面体：

选两个点p1, p2，再选一个点p3，p3和p1, p2不在同一条直线上，这样形成了一个面。再选一个点，不在同一个面上，形成一个四面体。

然后继续研究剩余的点：

点在四面体内部，略过。

点在四面体外部，删除此点能“看见”的面，扩大凸包体积。

POJ 3528 Ultimate Weapon

​​#include #include using namespace std;const double eps=1e-7;const int N=5e2+10;struct point{ double x,y,z; point(){ x=0; y=0; z=0; } point(double _x,double _y,double _z):x(_x),y(_y),z(_z) {} point operator -(const point other){ return point(x-other.x,y-other.y,z-other.z); } point operator *(const point other){ return point(y*other.z-z*other.y,-x*other.z+z*other.x,x*other.y-y*other.x); } double operator ^(const point other){ //点乘 return x*other.x+y*other.y+z*other.z; }};struct ch3D{ struct fac{ int a,b,c; // 面上三点的编号 bool ok; // 是否属于凸包 }; int n; //初始点数 point P[N]; //初始点 int cnt; //凸包上的三角形 fac F[N*6]; // 被创建的面不超过6N int vis[N][N]; double vlen(point v){ return sqrt(v.x*v.x+v.y*v.y+v.z*v.z); } double area(point p1,point p2,point p3){ return vlen((p2-p1)*(p3-p1)); // 三角形面积*2 } double volumn(point p1,point p2,point p3,point p4){ return fabs((p2-p1)*(p3-p1)^(p4-p1)); //V的6倍 } double pToPlane(point p,fac ff){ // 点在面的同向: + point m=P[ff.b]-P[ff.a]; point n=P[ff.c]-P[ff.a]; point v=p-P[ff.a]; return m*n^v; } void deal(int p,int a,int b){ int f=vis[a][b]; fac add; if(F[f].ok){ if(pToPlane(P[p],F[f])>eps) dfs(p,f); //在同侧，能看见 else { add.a=b, add.b=a, add.c=p, add.ok=1; vis[p][b]=vis[a][p]=vis[b][a]=cnt; F[cnt++]=add; } } } void dfs(int p,int num){ // 维护凸包 F[num].ok=0; //删除能看见的面 deal(p,F[num].b,F[num].a); deal(p,F[num].c,F[num].b); deal(p,F[num].a,F[num].c); } bool same(int s,int e){ //判断两个面是否是同一个面 point a=P[F[s].a],b=P[F[s].b],c=P[F[s].c]; return fabs(volumn(a,b,c,P[F[e].a]))eps){ jud=0; swap(P[1],P[i]); break; } } if(jud) return ; jud=1; // 三点在面 for(int i=2;ieps){ swap(P[2],P[i]); jud=0; break; } } if(jud) return ; jud=1; // 四点成体 for(int i=3;ieps){ swap(P[3],P[i]); jud=0; break; } } if(jud) return ; /*************************************/ fac add; for(int i=0;i<4;i++){ add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=1; if(pToPlane(P[i],add)>0) swap(add.b,add.c); vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=cnt; F[cnt++]=add; } for(int i=4;ieps){ dfs(i,j); break; } } } int tmp=cnt; cnt=0; for(int i=0;i

3D Convex Hull

​​#include #include using namespace std;const double eps=1e-7;const int N=5e2+10;struct point{ double x,y,z; point(){ x=0; y=0; z=0; } point(double _x,double _y,double _z):x(_x),y(_y),z(_z) {} point operator -(const point other){ return point(x-other.x,y-other.y,z-other.z); } point operator *(const point other){ return point(y*other.z-z*other.y,-x*other.z+z*other.x,x*other.y-y*other.x); } double operator ^(const point other){ //点乘 return x*other.x+y*other.y+z*other.z; }};struct ch3D{ struct fac{ int a,b,c; // 面上三点的编号 bool ok; // 是否属于凸包 }; int n; //初始点数 point P[N]; //初始点 int cnt; //凸包上的三角形 fac F[N*6]; // 被创建的面不超过6N int vis[N][N]; double vlen(point v){ return sqrt(v.x*v.x+v.y*v.y+v.z*v.z); } double area(point p1,point p2,point p3){ return vlen((p2-p1)*(p3-p1)); // 三角形面积*2 } double volumn(point p1,point p2,point p3,point p4){ return fabs((p2-p1)*(p3-p1)^(p4-p1)); //V的6倍 } double pToPlane(point p,fac ff){ // 点在面的同向: + point m=P[ff.b]-P[ff.a]; point n=P[ff.c]-P[ff.a]; point v=p-P[ff.a]; return m*n^v; } void deal(int p,int a,int b){ int f=vis[a][b]; fac add; if(F[f].ok){ if(pToPlane(P[p],F[f])>eps) dfs(p,f); //在同侧，能看见 else { add.a=b, add.b=a, add.c=p, add.ok=1; vis[p][b]=vis[a][p]=vis[b][a]=cnt; F[cnt++]=add; } } } void dfs(int p,int num){ // 维护凸包 F[num].ok=0; //删除能看见的面 deal(p,F[num].b,F[num].a); deal(p,F[num].c,F[num].b); deal(p,F[num].a,F[num].c); } bool same(int s,int e){ //判断两个面是否是同一个面 point a=P[F[s].a],b=P[F[s].b],c=P[F[s].c]; return fabs(volumn(a,b,c,P[F[e].a]))eps){ jud=0; swap(P[1],P[i]); break; } } if(jud) return ; jud=1; // 三点在面 for(int i=2;ieps){ swap(P[2],P[i]); jud=0; break; } } if(jud) return ; jud=1; // 四点成体 for(int i=3;ieps){ swap(P[3],P[i]); jud=0; break; } } if(jud) return ; /*************************************/ fac add; for(int i=0;i<4;i++){ add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=1; if(pToPlane(P[i],add)>0) swap(add.b,add.c); vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=cnt; F[cnt++]=add; } for(int i=4;ieps){ dfs(i,j); break; } } } int tmp=cnt; cnt=0; for(int i=0;i

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。