【HDU4273】Rescue（三维凸包）传送门

89 阅读 0 评论 59 点赞

我是靠谱客的博主聪明羽毛，最近开发中收集的这篇文章主要介绍【HDU4273】Rescue（三维凸包）传送门，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

传送门

全部都是常规操作，求一个三维凸包之后算一下重心然后求到每个面的距离就行了。

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define db double
#define cs const

using std::cerr;
using std::cout;

cs db eps=1e-8;

int sign(db x){return x<-eps?-1:(x>eps?1:0);}

struct Pnt{
	db x,y,z;Pnt(){}Pnt(db _x,db _y,db _z):x(_x),y(_y),z(_z){}
	friend Pnt operator+(cs Pnt &a,cs Pnt &b){return Pnt(a.x+b.x,a.y+b.y,a.z+b.z);}
	friend Pnt operator-(cs Pnt &a,cs Pnt &b){return Pnt(a.x-b.x,a.y-b.y,a.z-b.z);}
	friend Pnt operator*(cs Pnt &a,cs Pnt &b){return Pnt(a.y*b.z-b.y*a.z,a.z*b.x-b.z*a.x,a.x*b.y-b.x*a.y);}
	friend Pnt operator*(cs Pnt &a,cs db &b){return Pnt(a.x*b,a.y*b,a.z*b);}
	friend Pnt operator/(cs Pnt &a,cs db &b){return Pnt(a.x/b,a.y/b,a.z/b);}
	friend db dot(cs Pnt &a,cs Pnt &b){return a.x*b.x+a.y*b.y+a.z*b.z;}
	friend bool operator==(cs Pnt &a,cs Pnt &b){return !sign(a.x-b.x)&&!sign(a.y-b.y)&&!sign(a.z-b.z);}
	friend bool operator!=(cs Pnt &a,cs Pnt &b){return !(a==b);}
	inline db len()cs{return sqrt(x*x+y*y+z*z);}
};

cs Pnt O=Pnt(0,0,0);

db area(cs Pnt &a,cs Pnt &b,cs Pnt &c){
	return ((b-a)*(c-a)).len();
}

db vol(cs Pnt &a,cs Pnt &b,cs Pnt &c,cs Pnt &d){
	return dot((b-a)*(c-a),d-a);
}

cs int N=1e2+7;

int n,ct,id[N][N];
bool ok[N<<1|1];

Pnt p[N];

struct Face{
	int a,b,c;Face(){}Face(int _a,int _b,int _c):a(_a),b(_b),c(_c){}
	Pnt normal()cs{return (p[b]-p[a])*(p[c]-p[a]);}
	db area()cs{return normal().len();}
	db vol(cs Pnt &d)cs{return dot(normal(),d-p[a]);}
	Pnt sm()cs{return p[a]+p[b]+p[c];}
}f[N<<1|1];

db dist(cs Face &f,cs Pnt &p){
	return fabs(f.vol(p))/f.area();
}

void dfs(int,int);
void deal(int,int,int);

void deal(int x,int a,int b){
	if(ok[id[a][b]]){
		if(f[id[a][b]].vol(p[x])>eps)dfs(x,id[a][b]);
		else {
			id[b][a]=id[a][x]=id[x][b]=++ct;
			f[ct]=Face(b,a,x);ok[ct]=true;
		}
	}
}

void dfs(int p,int x){
	ok[x]=false;
	deal(p,f[x].b,f[x].a);
	deal(p,f[x].a,f[x].c);
	deal(p,f[x].c,f[x].b);
}

void Convex_3D(){
	ct=0;bool flag=true;
	for(int re i=2;i<=n;++i)if(p[1]!=p[i])
	{std::swap(p[2],p[i]);flag=false;break;}
	if(flag)return ;flag=true;
	for(int re i=3;i<=n;++i)if(sign(area(p[1],p[2],p[i])))
	{std::swap(p[3],p[i]);flag=false;break;}
	if(flag)return ;flag=true;
	for(int re i=4;i<=n;++i)if(sign(vol(p[1],p[2],p[3],p[i])))
	{std::swap(p[4],p[i]);flag=false;break;}
	if(flag)return ;flag=true;
	
	for(int re i=1;i<=4;++i){
		int a=(i&3)+1,b=(a&3)+1,c=(b&3)+1;
		if(vol(p[a],p[b],p[c],p[i])>eps)std::swap(b,c);
		id[a][b]=id[b][c]=id[c][a]=++ct;
		f[ct]=Face(a,b,c);ok[ct]=true;
	}
	for(int re i=5;i<=n;++i){
		for(int re j=1;j<=ct;++j)
			if(f[j].vol(p[i])>eps)
				{dfs(i,j);flag=false;break;}
		if(flag)continue;flag=true;int cur=0;
		for(int re j=1;j<=ct;++j)if(ok[j]){
			if((++cur)==j)continue;f[cur]=f[j],ok[cur]=true;
			id[f[cur].a][f[cur].b]=id[f[cur].b][f[cur].c]=id[f[cur].c][f[cur].a]=cur;
		}ct=cur;
	}
}

Pnt calc_G(){
	Pnt ans=O;db all=0;
	for(int re i=1;i<=ct;++i){
		db vol=f[i].vol(O);
		all+=vol;ans=ans+f[i].sm()*vol;
	}return ans/(all*4);
}

void Main(){
	while(~scanf("%d",&n)){
		for(int re i=1;i<=n;++i)
			scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
		Convex_3D();Pnt G=calc_G();db ans=1e20;
		for(int re i=1;i<=ct;++i)ans=std::min(ans,dist(f[i],G));
		printf("%.3fn",ans);
	}
}

void file(){
#ifdef zxyoi
	freopen("rescue.in","r",stdin);
#endif
}
signed main(){file();Main();return 0;}