hdu3662 3D Convex Hull(三維凸包【三維計算幾何基本操作)
阿新 • • 發佈:2019-01-09
分析:
三維凸包模板
瞧好了基本操作
三維和二維簡直不是一個級別的。。。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;
}