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);
    }
}

　　

转载于:https://www.cnblogs.com/jie-dcai/p/3902543.html