hdu4273Rescue(三维凸包重心)

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<stdlib.h>
#include<vector>
#include<cmath>
#include<queue>
#include<set>
using namespace std;
#define N 510
#define INF 1e20
#define max(a,b) (a>b?a:b)
#define min(a,b) (a<b?a:b)
#define eps 1e-8
#define MAXV 505
const double pi = acos(-1.0);
const double inf = ~0u>>;
//三维点
struct point3
{
double x, y,z;
point3() {}
point3(double _x, double _y, double _z): x(_x), y(_y), z(_z) {}
point3 operator +(const point3 p1)
{
return point3(x+p1.x,y+p1.y,z+p1.z);
}
point3 operator - (const point3 p1)
{
return point3(x - p1.x, y - p1.y, z - p1.z);
}
point3 operator * (point3 p)
{
return point3(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x); //叉乘
}
point3 operator *(double d)
{
return point3(x*d,y*d,z*d);
}
point3 operator /(double d)
{
return point3(x/d,y/d,z/d);
}
double operator ^ (point3 p)
{
return x*p.x+y*p.y+z*p.z; //点乘
} } pp[N],rp[N];
struct point
{
double x,y;
point(double x=,double y=):x(x),y(y) {}
point operator -(point b)
{
return point(x-b.x,y-b.y);
}
} p[N],ch[N];
struct _3DCH
{
struct fac
{
int a, b, c; //表示凸包一个面上三个点的编号
bool ok; //表示该面是否属于最终凸包中的面
}; int n; //初始点数
point3 P[MAXV]; //初始点 int cnt; //凸包表面的三角形数
fac F[MAXV*]; //凸包表面的三角形 int to[MAXV][MAXV];
double vlen(point3 a)
{
return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
} //向量长度
double area(point3 a, point3 b, point3 c)
{
return vlen((b-a)*(c-a));
} //三角形面积*2
double volume(point3 a, point3 b, point3 c, point3 d)
{
return (b-a)*(c-a)^(d-a); //四面体有向体积*6
}
//正：点在面同向
double point3of(point3 &p, fac &f)
{
point3 m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
return (m * n) ^ t;
}
void deal(int p, int a, int b)
{
int f = to[a][b];
fac add;
if (F[f].ok)
{
if (point3of(P[p], F[f]) > eps)
dfs(p, f);
else
{
add.a = b, add.b = a, add.c = p, add.ok = ;
to[p][b] = to[a][p] = to[b][a] = cnt;
F[cnt++] = add;
}
}
}
void dfs(int p, int cur)
{
F[cur].ok = ;
deal(p, F[cur].b, F[cur].a);
deal(p, F[cur].c, F[cur].b);
deal(p, F[cur].a, F[cur].c);
}
bool same(int s, int t)
{
point3 &a = P[F[s].a], &b = P[F[s].b], &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 construct()
{
cnt = ;
if (n < )
return;
bool sb = ;
//使前两点不公点
for (int i = ; i < n; i++)
{
if (vlen(P[] - P[i]) > eps)
{
swap(P[], P[i]);
sb = ;
break;
}
}
if (sb)return;
sb = ;
//使前三点不公线
for (int i = ; i < n; i++)
{
if (vlen((P[] - P[]) * (P[] - P[i])) > eps)
{
swap(P[], P[i]);
sb = ;
break;
}
}
if (sb)return;
sb = ;
//使前四点不共面
for (int i = ; i < n; i++)
{
if (fabs((P[] - P[]) * (P[] - P[]) ^ (P[] - P[i])) > eps)
{
swap(P[], P[i]);
sb = ;
break;
}
}
if (sb)return;
fac add;
for (int i = ; i < ; i++)
{
add.a = (i+)%, add.b = (i+)%, add.c = (i+)%, add.ok = ;
if (point3of(P[i], add) > )
swap(add.b, add.c);
to[add.a][add.b] = to[add.b][add.c] = to[add.c][add.a] = cnt;
F[cnt++] = add;
}
for (int i = ; i < n; i++)
{
for (int j = ; j < cnt; j++)
{
if (F[j].ok && point3of(P[i], F[j]) > eps)
{
dfs(i, j);
break;
}
}
}
int tmp = cnt;
cnt = ;
for (int i = ; i < tmp; i++)
{
if (F[i].ok)
{
F[cnt++] = F[i];
}
}
}
//表面积
double area()
{
double ret = 0.0;
for (int i = ; i < cnt; i++)
{
ret += area(P[F[i].a], P[F[i].b], P[F[i].c]);
}
return ret / 2.0;
}
double ptoface(point3 p,int i)
{
return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
}
//体积
double volume()
{
point3 O(, , );
double ret = 0.0;
for (int i = ; i < cnt; i++)
{
ret += volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
}
return fabs(ret / 6.0);
}
//表面三角形数
int facetCnt_tri()
{
return cnt;
} //表面多边形数
int facetCnt()
{
int ans = ;
for (int i = ; i < cnt; i++)
{
bool nb = ;
for (int j = ; j < i; j++)
{
if(same(i, j))
{
nb = ;
break;
}
}
ans += nb;
}
return ans;
}
//三维凸包重心
point3 barycenter()
{
point3 ans(,,),o(,,);
double all=;
for(int i=;i<cnt;i++)
{
double vol=volume(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;
} }hull; void solve()
{
double ans = INF;
int i;
int cnt = hull.cnt;
point3 pp = hull.barycenter();
for(i = ; i < cnt ; i++)
{
ans = min(ans,hull.ptoface(pp,i));
}
printf("%.3f\n",ans);
}
int main()
{
int n,i;
while(scanf("%d",&n)!=EOF)
{
hull.n = n;
for(i = ; i < n; i++)
{
scanf("%lf%lf%lf",&pp[i].x,&pp[i].y,&pp[i].z);
hull.P[i] = pp[i];
}
hull.construct();
solve();
}
}