1. 程式人生 > >hdu3662 3D Convex Hull(三維凸包【三維計算幾何基本操作)

hdu3662 3D Convex Hull(三維凸包【三維計算幾何基本操作)

題目連線

分析:
三維凸包模板

瞧好了基本操作
三維和二維簡直不是一個級別的。。。orz

#include<bits/stdc++.h>

using namespace std;

const double eps=1e-8;
const int N=503;
struct node{
    double x,y,z;
    node (double xx=0,double yy=0,double zz=0)
    {
        x=xx;y=yy;z=zz;
    }
};

node operator +(const node &A,const node &B) {return
node(A.x+B.x,A.y+B.y,A.z+B.z);} node operator -(const node &A,const node &B) {return node(A.x-B.x,A.y-B.y,A.z-B.z);} node operator *(const node &A,const double &B) {return node(A.x*B,A.y*B,A.z*B);} node operator /(const node &A,const double &B) {return node(A.x/B,A.y/B,A.z/B);} double
Dot(node A,node B) {return A.x*B.x+A.y*B.y+A.z*B.z;} node Cross(node A,node B) {return node(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);} struct CH3D{ struct face{ int a,b,c; //凸包一個面上的三個點編號 bool ff; //該面是否屬於最終凸包上的面 }; int n; //初始頂點數 node P[N]; int
num; //凸包表面的三角形數 face F[8*N]; int vis[N][N]; //向量長度 double lenth(node a) { return sqrt(Dot(a,a)); } //三角形面積*2 double S(node a,node b,node c) { return lenth(Cross(b-a,c-a)); //注意要取lenth } //四面體有向體積*6 混合積 double V(node a,node b,node c,node d) { return Dot(d-a,Cross(b-a,c-a)); } //面是否可見 int cansee(node a,face f) { double x=Dot(a-P[f.a],Cross(P[f.b]-P[f.a],P[f.c]-P[f.a])); if (x>eps) return 1; else return 0; } //搜尋與該邊相鄰的另一個平面 void deal(int p,int a,int b) { int f=vis[a][b]; face add; if (F[f].ff) { if (cansee(P[p],F[f])) dfs(p,f); else { add.a=b; add.b=a; add.c=p; add.ff=1; vis[p][b]=vis[a][p]=vis[b][a]=num; F[num++]=add; } } } //遞迴搜尋所有應該從凸包內刪除的面 void dfs(int p,int now) { F[now].ff=0; deal(p,F[now].b,F[now].a); deal(p,F[now].c,F[now].b); deal(p,F[now].a,F[now].c); } //構建三維凸包 void create() { face add; num=0; if (n<4) return; //保證前四個點不共面 bool flag=1; for (int i=1;i<n;i++) { if (lenth(P[0]-P[i])>eps) { swap(P[1],P[i]); flag=0; break; } } if (flag) return; flag=1; for (int i=2;i<n;i++) { if (S(P[0],P[1],P[i])>eps) { swap(P[2],P[i]); flag=0; break; } } if (flag) return; flag=1; for (int i=3;i<n;i++) { if (fabs(V(P[0],P[1],P[2],P[i]))>eps) { swap(P[3],P[i]); flag=0; break; } } if (flag) return; for (int i=0;i<4;i++) { add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ff=1; if (cansee(P[i],add)) swap(add.b,add.c); vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=num; //加面就要維護vis F[num++]=add; } for (int i=4;i<n;i++) for (int j=0;j<num;j++) if (F[j].ff&&cansee(P[i],F[j])){ dfs(i,j); break; } int tmp=num; for (int i=num=0;i<tmp;i++) if (F[i].ff) F[num++]=F[i]; } //表面積 double Ssum() { double ans=0; if (n==3) { ans=S(P[0],P[1],P[2]); return ans/2.0; } for (int i=0;i<num;i++) ans=ans+S(P[F[i].a],P[F[i].b],P[F[i].c]); return ans/2.0; } //體積 double Vsum() { double ans=0; node tmp(0,0,0); for (int i=0;i<num;i++) ans+=V(tmp,P[F[i].a],P[F[i].b],P[F[i].c]); return fabs(ans/6.0); } bool same(int s,int t) { node &a=P[F[s].a]; node &b=P[F[s].b]; node &c=P[F[s].c]; return fabs(V(a,b,c,P[F[t].a]))<eps && fabs(V(a,b,c,P[F[t].b]))<eps && fabs(V(a,b,c,P[F[t].c]))<eps; } //表面多邊形個數 int facecnt() { int flag,ans=0; for (int i=0;i<num;i++) { flag=1; for (int j=0;j<i;j++) if (same(i,j)) { flag=0; break; } ans+=flag; } return ans; } //三維凸包重心 node centre() { node ans(0,0,0),o(0,0,0); double all=0; for (int i=0;i<num;i++) { double vol=V(o,P[F[i].a],P[F[i].b],P[F[i].c]); ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol; all+=vol; } ans=ans/all; return ans; } //點到面的距離 double dis_face(node p,int i) { return fabs(V(P[F[i].a],P[F[i].b],P[F[i].c],p)/lenth(Cross(P[F[i].b]-P[F[i].a],P[F[i].c]-P[F[i].a]))); } }; CH3D hull; int main() { while (scanf("%d",&hull.n)!=EOF) { for (int i=0;i<hull.n;i++) scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z); hull.create(); printf("%d\n",hull.facecnt()); } return 0; }