POJ 2225 / ZOJ 1438 / UVA 1438 Asteroids --三维凸包,求多面体重心

题意: 两个凸多面体,可以任意摆放,最多贴着,问他们重心的最短距离。

解法: 由于给出的是凸多面体,先构出两个三维凸包,再求其重心,求重心仿照求三角形重心的方式,然后再求两个多面体的重心到每个多面体的各个面的最短距离,然后最短距离相加即为答案,因为显然贴着最优。

求三角形重心见此: http://www.cnblogs.com/whatbeg/p/4234518.html

代码:(模板借鉴网上模板)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <string>
#include <vector>
#include <set>
#define Mod 1000000007
#define eps 1e-8
#define lll __int64
#define ll long long
using namespace std;
#define N 100007
#define MAXV 505 //三维点
struct pt{
double x, y, z;
pt(){}
pt(double _x, double _y, double _z): x(_x), y(_y), z(_z){}
pt operator - (const pt p1){return pt(x - p1.x, y - p1.y, z - p1.z);}
pt operator * (pt p){return pt(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);} //叉乘
double operator ^ (pt p){return x*p.x+y*p.y+z*p.z;} //点乘
}; //pt operator - (const pt p,const pt p1){return pt(p.x - p1.x, p.y - p1.y, p.z - p1.z);}
//pt operator ** (pt p,pt p1){return pt(p.y*p1.z-p.z*p1.y, p.z*p1.x-p.x*p1.z, p.x*p1.y-p.y*p1.x);} //叉乘
//double operator ^^ (pt p1,pt p){return p1.x*p.x+p1.y*p.y+p1.z*p.z;} struct _3DCH{
struct fac{
int a, b, c; //表示凸包一个面上三个点的编号
bool ok; //表示该面是否属于最终凸包中的面
}; int n; //初始点数
pt P[MAXV]; //初始点 int cnt; //凸包表面的三角形数
fac F[MAXV*]; //凸包表面的三角形 int to[MAXV][MAXV]; double vlen(pt a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);} //向量长度
double area(pt a, pt b, pt c){return vlen((b-a)*(c-a));} //三角形面积*2
double volume(pt a, pt b, pt c, pt d){return (b-a)*(c-a)^(d-a);} //四面体有向体积*6 //正:点在面同向
double ptof(pt &p, fac &f){
pt m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
return (m * n) ^ t;
}
pt pvec(fac s) {
pt k1 = (P[s.a]-P[s.b]), k2 = (P[s.b]-P[s.c]);
return (k1*k2);
}
double ptoplane(pt p,fac s){
return fabs(pvec(s)^(p-P[s.a]))/vlen(pvec(s));
} void deal(int p, int a, int b){
int f = to[a][b];
fac add;
if (F[f].ok){
if (ptof(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){
pt &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 (ptof(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 && ptof(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 volume(){
pt 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);
} pt BaryCenter() {
pt O(, , );
double ret = 0.0,sumvolume = 0.0, sumx = 0.0, sumy = 0.0, sumz = 0.0;
for(int i=;i<cnt;i++) {
double Vol = volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
sumvolume += Vol;
sumx += (P[F[i].a].x + P[F[i].b].x + P[F[i].c].x)*Vol;
sumy += (P[F[i].a].y + P[F[i].b].y + P[F[i].c].y)*Vol;
sumz += (P[F[i].a].z + P[F[i].b].z + P[F[i].c].z)*Vol;
}
return pt(sumx/sumvolume/, sumy/sumvolume/, sumz/sumvolume/);
} //表面三角形数
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;
}
}; _3DCH hull,hull2; //内有大数组,不易放在函数内 int main()
{
while (scanf("%d", &hull.n)!=EOF){
for (int i = ; i < hull.n; i++)
scanf("%lf%lf%lf", &hull.P[i].x, &hull.P[i].y, &hull.P[i].z);
hull.construct();
pt bc1 = hull.BaryCenter();
scanf("%d",&hull2.n);
for (int i = ; i < hull2.n; i++)
scanf("%lf%lf%lf", &hull2.P[i].x, &hull2.P[i].y, &hull2.P[i].z);
hull2.construct();
pt bc2 = hull2.BaryCenter();
//printf("BARY1: %.2f %.2f %.2f\n",bc1.x,bc1.y,bc1.z);
//printf("BARY2: %.2f %.2f %.2f\n",bc2.x,bc2.y,bc2.z);
double dis1 = Mod, dis2 = Mod;
for (int i = ; i < hull.cnt; i++)
dis1 = min(dis1,fabs(hull.ptoplane(bc1,hull.F[i])));
for (int i = ; i < hull2.cnt; i++)
dis2 = min(dis2,fabs(hull2.ptoplane(bc2,hull2.F[i])));
printf("%.6f\n",dis1+dis2);
}
return ;
}
上一篇:首次发布App,In-App Purchase 无法submit for review 问题的解决方案


下一篇:【原创】kafka server源代码分析(一)