首页 > 代码库 > HDU 4273
HDU 4273
计算凸包重心到各面的最短距离。
若知道重心,按四面体用体积法即可求出高。
关键在于,多面体重心的求法。这必须把多面体分割成多个四面体来求。下面从多边形的重心说起。
一般来用,对于一个多边形(p0,p1,p2....pn-1),其重心一般为pc.x=(p0.x+p1.x+....)/n对于y也一样。
但这其实是不正确的。反例以梯形为例。上面的式子当各点的权值均匀时是正确的。(三角形是一个特例)
但在多边形上,由于面的密度一样,所以,应当是与面积有关的。于是,把多边形分割成多个三角形,求出其重心。这样重心组成一个新的多边形与原多边形重心相同。于是,就把质量都集中在了重心上。而质量与面积相关。
于是,可由代码求重心:
以P0为顶点划分三角形,求得是有向面积,因为可以正负抵消for(多边形上的点){ //逆时针 与p0组成三角形。 有向面积V=(p1-p0)*(p2-p0)/2; Vtot+=V; sum_x+=(p1.x+p2.x+p0.x)*V; sum_y+=(p1.y+p2.y+p0.y)*V; C.x=sum_x/3/Vtot; C.y=sum_y/3/Vtot}
对于求多面体,只需划分成四面体来求即可。增加一个z坐标,同时
sum_x+=(p1.x+p2.x+p0.x+p3.x)*V;C.x=sum_x/4/Vtot;
/*增量法求凸包。选取一个四面体,同时把它各面的方向向量向外,增加一个点时,若该点与凸包上的某些面的方向向量在同一侧,则去掉那些面,并使某些边与新增点一起连成新的凸包上的面。 */ #include <iostream>#include <cstring>#include <cstdio>#include <algorithm>#include <cmath> using namespace std;const int MAXN=110;const double eps=1e-8;const double inf=1e10;struct point { double x,y,z;};struct face { int a,b,c; bool ok;};int n; //初始点数 point p[MAXN]; //空间点int trianglecnt; //凸包上三角形数face tri[6*MAXN]; //凸包上被创建的三角形int vis[MAXN][MAXN]; //点i到点j是属于哪一个三角形。此处是有方向 point operator -(const point &x, const point &y){ point ret; ret.x=x.x-y.x; ret.y=x.y-y.y; ret.z=x.z-y.z; return ret;} point operator * (const point &u,const point &v){ //叉积 point ret; ret.x=u.y*v.z-u.z*v.y; ret.y=u.z*v.x-u.x*v.z; ret.z=u.x*v.y-u.y*v.x; return ret;} double operator ^(const point &u,const point &v){ //点积 return (u.x*v.x+u.y*v.y+u.z*v.z);} double dist(point t){ return sqrt(t.x*t.x+t.y*t.y+t.z*t.z);} double ptoplane(point &tmp,face &f){ //若结果大于0,证明点面的同向,即法向量方向 point m=p[f.b]-p[f.a]; point n=p[f.c]-p[f.a]; point t=tmp-p[f.a]; return (m*n)^t;} double farea(point a,point b,point c ){ point t1=a-c; point t2=b-c; return fabs(dist(t1*t2));}void dfs(int pt, int ct);void deal(int pt,int a,int b){ int f=vis[a][b]; //所属三角形,即原来的ab。 face add; if(tri[f].ok){ if((ptoplane(p[pt],tri[f]))>eps) dfs(pt,f); //若点同样在该f三角形方向一侧,继续调整 else { add.a=b; add.b=a; add.c=pt; add.ok=1; vis[pt][b]=vis[a][pt]=vis[b][a]=trianglecnt; tri[trianglecnt++]=add; } }} void dfs(int pt, int ct){ tri[ct].ok=0; //去掉该面 deal(pt,tri[ct].b,tri[ct].a); //因为有向边ab所属三角形去掉,则反方向边必定属于另一个三角形. deal(pt,tri[ct].c,tri[ct].b); deal(pt,tri[ct].a,tri[ct].c);} void construct (){ int i,j; trianglecnt=0; if(n<4) return ; //不可能构成一个多面体 bool tmp=true; for(i=1;i<n;i++){ //不共点两点 if(dist(p[0]-p[i])>eps){ swap(p[1],p[i]); tmp=false; break; } } if(tmp) return ; tmp=true; for(i=2;i<n;i++){ //不共线 if(dist((p[0]-p[1])*(p[1]-p[i]))>eps){ swap(p[2],p[i]); tmp=false; break; } } if(tmp) return ; tmp=true; for(i=3;i<n;i++){ //四点不共面K if(fabs((p[0]-p[1])*(p[1]-p[2])^(p[0]-p[i]))>eps){ swap(p[3],p[i]); tmp=false; break; } } if(tmp) return ; face add; for(i=0;i<4;i++){ //使各三角形的方向向量向外,同时记录下三角形的序号 add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=1; //等于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]=trianglecnt; tri[trianglecnt++]=add; } for(i=4;i<n;i++){ //构建凸包 for(j=0;j<trianglecnt;j++){ if(tri[j].ok&&(ptoplane(p[i],tri[j]))>eps){ //增加点可见该平,即在面方向一侧 dfs(i,j); break; } } } int cnt=trianglecnt; trianglecnt=0; for(i=0;i<cnt;i++){ //只有ok为1的才属于凸包上的三角形 if(tri[i].ok){ tri[trianglecnt++]=tri[i]; } }}double cdis(point p0){ double ans=inf; point p1,p2,p3; for(int i=0;i<trianglecnt;i++){ p1=p[tri[i].a]; p2=p[tri[i].b]; p3=p[tri[i].c]; double V=fabs(((p0-p1)^((p2-p1)*(p3-p1)))/6); // printf("%lf\n",V); ans=min(ans,V*3*2/dist((p2-p1)*(p3-p1))); } return ans; } point Cenconstruct(){ point p0=p[0]; point p1,p2,p3; double sum_area=0,sum_x=0,sum_y=0,sum_z=0; for(int i=0;i<trianglecnt;i++){ p1=p[tri[i].a]; p2=p[tri[i].b]; p3=p[tri[i].c]; double V=((p0-p1)^((p2-p1)*(p3-p1)))/6; sum_area+=V; sum_x+=(p0.x+p1.x+p2.x+p3.x)*V; sum_y+=(p0.y+p1.y+p2.y+p3.y)*V; sum_z+=(p0.z+p1.z+p2.z+p3.z)*V; } point ret; ret.x=sum_x/4/sum_area; ret.y=sum_y/4/sum_area; ret.z=sum_z/4/sum_area; return ret;}int main(){ while(scanf("%d",&n)!=EOF){ memset(vis,0,sizeof(vis)); for(int i=0;i<n;i++) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z); construct(); point centroid=Cenconstruct(); double ans; ans = cdis(centroid); printf("%.3lf\n",ans); }}
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。