三维凸包
来源:互联网 发布:百度的人工智能 编辑:程序博客网 时间:2024/05/12 06:20
三维凸包模板 hdu4266
#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <cmath>#include <algorithm>using namespace std;const int MAXN = 1100;const double eps = 1e-9;struct Point{ double x,y,z; Point(){} Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){} Point operator -(const Point &P) { return Point(x-P.x, y-P.y, z-P.z); } Point operator ^ (const Point &P)//叉积 { return Point(y*P.z - z * P.y, z * P.x - x * P.z, x*P.y - y * P.x); } double operator * (const Point &P)//点积 { return (x * P.x + y * P.y + z * P.z); } };struct CH3D{ struct face { int a,b,c; bool ok ;// 判断是否是凸包的面 face(){} face(int _a,int _b,int _c, bool _ok):a(_a),b(_b),c(_c),ok(_ok){} void init(int _a,int _b,int _c, bool _ok){ a = _a; b = _b; c =_c; ok = _ok; } }; int n,num; Point p[MAXN]; int g[MAXN][MAXN];//凸包表面的三角形 face F[MAXN * 8]; double vlen(const Point &a) //向量长度 { return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); } 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)//四面体体积 { return ((b-a)^(c-a)) * (d-a); } double dblcmp(Point &P,face &f) { Point m=p[f.b]-p[f.a]; Point n=p[f.c]-p[f.a]; Point t=P-p[f.a]; return (m ^ n) * t; } void deal(int P,int a,int b) { int f=g[a][b]; if(F[f].ok) { if(dblcmp(p[P],F[f])>eps) dfs(P,f); else { face add(b,a,P,1); g[P][b]=g[a][P]=g[b][a]=num; F[num++]=add; } } } void dfs(int p,int now) { F[now].ok=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); } bool same(int s,int t) { Point &a=p[F[s].a]; Point &b=p[F[s].b]; Point &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; } void solve()//构建三维凸包 { face add; bool flag=true; num=0; if(n<4) return; for(int i=1;i<n;i++)//此段是为了保证前四个点不共面,若以保证,则可去掉 { if(vlen(p[0]-p[i])>eps){ swap(p[1],p[i]); flag=false; break; } } if(flag) return; flag=true; for(int i=2;i<n;i++)//使前三点不共线 { if(vlen((p[0]-p[1])^(p[1]-p[i]))>eps) { swap(p[2],p[i]); flag=false; break; } } if(flag) return; flag=true; for(int i=3;i<n;i++)//使前四点不共面 { if(fabs(((p[0]-p[1])^(p[1]-p[2]))*(p[0]-p[i]))>eps){ swap(p[3],p[i]); flag=false; break; } } if(flag) return; for(int i=0;i<4;i++) { add.init((i+1)%4,(i+2)%4,(i+3)%4,true); if(dblcmp(p[i],add)>0) swap(add.b,add.c); g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; F[num++]=add; } for(int i=4;i<n;i++){ for(int j=0;j<num;j++) { if(F[j].ok && dblcmp(p[i],F[j])>eps){ dfs(i,j); break; } } } int tmp=num; num = 0; for(int i=0;i<tmp;i++) if(F[i].ok) F[num++]=F[i]; } double area()//表面积 { double res=0.0; if(n==3){ //Point p=cross(P[0],P[1],P[2]); res=vlen((p[1]-p[0])^(p[2]-p[0]))/2.0; return res; } for(int i=0;i<num;i++) res+=vlen((p[F[i].b]-p[F[i].a]) ^ (p[F[i].c]-p[F[i].a])); return res/2.0; } double volume() //体积 { double res=0.0; Point temp(0,0,0); for(int i=0;i<num;i++) res+=((p[F[i].b]-p[F[i].a])^(p[F[i].c]-p[F[i].a])) * (temp-p[F[i].a]); return fabs(res/6.0); } int triangle()//表面三角形个数 { return num; } int polygon()//表面多边形个数 { int res = 0; for(int i=0,flag = 1; i<num;i++) { for(int j=0;j<i;j++) if(same(i,j)){ flag=0;break; } res+=flag; } return res; } double getMinDis(Point Q) { double ans = 0x7fffffff; for (int i = 0; i < num;++i){ ans = min(ans, volume(Q,p[F[i].a],p[F[i].b],p[F[i].c]) /area(p[F[i].a],p[F[i].b],p[F[i].c])); } return ans; } }hull; int main(){ while (scanf("%d",&hull.n)!=EOF) { if (hull.n ==0) break; 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.solve(); int q; cin>>q; for (int i = 0; i < q;++i) { Point Q; scanf("%lf%lf%lf",&Q.x,&Q.y,&Q.z); printf("%.4lf\n", hull.getMinDis(Q)); } } return 0;}
2013 ACM/ICPC 长沙邀请赛- I 题 Throw the stones
#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <cmath>#include <algorithm>#include <map>using namespace std;const int MAXN = 10002;const double eps = 1e-9;int child;double last,D;struct Point{ double x,y,z; Point(){} Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){} Point operator -(const Point &P) { return Point(x-P.x, y-P.y, z-P.z); } Point operator ^ (const Point &P)//叉积 { return Point(y*P.z - z * P.y, z * P.x - x * P.z, x*P.y - y * P.x); } double operator * (const Point &P)//点积 { return (x * P.x + y * P.y + z * P.z); } };struct CH3D{ struct face { int a,b,c; bool ok ;// 判断是否是凸包的面 face(){} face(int _a,int _b,int _c, bool _ok):a(_a),b(_b),c(_c),ok(_ok){} void init(int _a,int _b,int _c, bool _ok){ a = _a; b = _b; c =_c; ok = _ok; } }; int n,num; Point p[MAXN]; map<int,int> g[MAXN];//凸包表面的三角形 face F[MAXN * 10]; double vlen(const Point &a) //向量长度 { return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); } 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)//四面体体积 { return ((b-a)^(c-a)) * (d-a); } double dblcmp(Point &P,face &f) { Point m=p[f.b]-p[f.a]; Point n=p[f.c]-p[f.a]; Point t=P-p[f.a]; return (m ^ n) * t; } void deal(int P,int a,int b) { int f=g[a][b]; if(F[f].ok) { if(dblcmp(p[P],F[f])>eps) dfs(P,f); else { D += (p[a]^p[P]) * p[b]; face add(b,a,P,1); g[P][b]=g[a][P]=g[b][a]=num; F[num++]=add; } } } void dfs(int P,int now) { D-=(p[F[now].b] ^ p[F[now].c] )* p[F[now].a]; F[now].ok=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); } bool same(int s,int t) { Point &a=p[F[s].a]; Point &b=p[F[s].b]; Point &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; } void init() { face add; for(int i=0;i<4;i++) { add.init((i+1)%4,(i+2)%4,(i+3)%4,true); if(dblcmp(p[i],add)>0) swap(add.b,add.c); g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; F[num++]=add; } } void update(int i) { for(int j=0;j<num;j++) if(F[j].ok && dblcmp(p[i],F[j])>eps){ dfs(i,j); break; } } void solve()//构建三维凸包 { face add; bool flag=true; num=0; if(n<4) return; for(int i=1;i<n;i++)//此段是为了保证前四个点不共面,若以保证,则可去掉 { if(vlen(p[0]-p[i])>eps){ swap(p[1],p[i]); flag=false; break; } } if(flag) return; flag=true; for(int i=2;i<n;i++)//使前三点不共线 { if(vlen((p[0]-p[1])^(p[1]-p[i]))>eps) { swap(p[2],p[i]); flag=false; break; } } if(flag) return; flag=true; for(int i=3;i<n;i++)//使前四点不共面 { if(fabs(((p[0]-p[1])^(p[1]-p[2]))*(p[0]-p[i]))>eps){ swap(p[3],p[i]); child = i; init(); for (int j = 4; j<=child;++j) update(j); last = volume(); flag=false; break; } } if(flag) return; for(int i=child+1;i<n;i++) { D = 0;update(i);D/=6.0; if (D>last+eps){ child = i; last = D; } } int tmp=num; num = 0; for(int i=0;i<tmp;i++) if(F[i].ok) F[num++]=F[i]; } double area()//表面积 { double res=0.0; if(n==3){ //Point p=cross(P[0],P[1],P[2]); res=vlen((p[1]-p[0])^(p[2]-p[0]))/2.0; return res; } for(int i=0;i<num;i++) res+=vlen((p[F[i].b]-p[F[i].a]) ^ (p[F[i].c]-p[F[i].a])); return res/2.0; } double volume() //体积 { double res=0.0; Point temp(0,0,0); for(int i=0;i<num;i++) if (F[i].ok) res+=((p[F[i].b]-p[F[i].a])^(p[F[i].c]-p[F[i].a])) * (temp-p[F[i].a]); //res +=(p[F[i].b] ^ p[F[i].c]) *p[F[i].a]; return fabs(res/6.0); } int triangle()//表面三角形个数 { return num; } int polygon()//表面多边形个数 { int res = 0; for(int i=0,flag = 1; i<num;i++) { for(int j=0;j<i;j++) if(same(i,j)){ flag=0;break; } res+=flag; } return res; } double getMinDis(Point Q) { double ans = 0x7fffffff; for (int i = 0; i < num;++i){ ans = min(ans, volume(Q,p[F[i].a],p[F[i].b],p[F[i].c]) /area(p[F[i].a],p[F[i].b],p[F[i].c])); } return ans; } }hull; int main(){ int test =0; while (scanf("%d",&hull.n)!=EOF) { test ++; child = 0;last = 0; 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.solve(); printf("Case #%d:\n%d %.2lf\n",test, child+1,last); } return 0;}
- 三维凸包
- 三维凸包模板
- 三维凸包模板
- 三维凸包
- 三维凸包模板
- 三维凸包模板
- 三维凸包模板
- 三维凸包大全
- uvalive5090(三维凸包)
- 三维凸包模板
- 【增量算法】三维凸包
- poj3528 三维凸包面积
- Asteroids (三维凸包+重心)
- POJ3528 (三维凸包模板)
- hdu 4273 Rescue(三维凸包 三维中心)
- 求取模型的三维凸包问题
- hdu4273Rescue三维凸包的重心
- 解题报告-HDU 3662 (三维凸包)
- 13:java并发编程总结
- UcHome二次开发调试技巧
- linux shell编程控制结构:expr、let、for、while、until、shift、if、case、break、continue、函数、select 学习笔记
- Cocos2d-x Tiled Map Editor(一)
- UCHome二次开发 模板基础语法
- 三维凸包
- Cocos2d-x Tiled Map Editor(二)
- Ext中左侧tree与右侧grid,grid分页问题
- 创业的野性
- UCHome二次开发 规范
- 软件更新提示(检查app store 软件是否有新的版本)---需要app id
- 计算机的因果机制
- 免费 python_books
- OpenGL学习笔记(二):视图