poj3862 Asteroids【模板题】

题目链接：http://poj.org/problem?id=3862
题意：就是让你求两个多面体的质心到面的距离之和的最小值
解析：模板题，三维凸包质心，然后枚举每一个面算质心到面的距离取最小值

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <cstring>
#include <queue>
#include <cmath>
#include <map>
using namespace std;
const int maxn = 1e6+100;
const int inf = 0x7ffffff;
const double eps = 1e-8;
int n;
struct point
{
    double x,y,z;
    point() {}
    point(double _x,double _y,double _z)
    {
        x = _x;
        y = _y;
        z = _z;
    }
    point operator - (const point &b)const
    {
        point ret;
        ret.x = x-b.x;
        ret.y = y-b.y;
        ret.z = z-b.z;
        return ret;
    }
    point operator + (const point &b)const
    {
        point ret;
        ret.x = x+b.x;
        ret.y = y+b.y;
        ret.z = z+b.z;
        return ret;
    }
    point operator * (const double &b)const
    {
        point ret;
        ret.x = x*b;
        ret.y = y*b;
        ret.z = z*b;
        return ret;
    }
    point operator / (const double &b)const
    {
        point ret;
        ret.x = x/b;
        ret.y = y/b;
        ret.z = z/b;
        return ret;
    }
};
int dblcmp(double x)
{
    if(fabs(x)<eps)
        return 0;
    return x>0?1:-1;
}
point x_mul(point u,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 dot_mul(point u,point v)
{
    return u.x*v.x+u.y*v.y+u.z*v.z;
}
double dis(point p1,point p2)
{
    double t1 = (p1.x-p2.x)*(p1.x-p2.x);
    double t2 = (p1.y-p2.y)*(p1.y-p2.y);
    double t3 = (p1.z-p2.z)*(p1.z-p2.z);
    return sqrt(t1+t2+t3);
}
double vlen(point p)
{
    return sqrt(dot_mul(p,p));
}
double area(point a,point b,point c,point d)
{
    point t1 = d-a;
    point t2 = b-a;
    point t3 = c-a;
    point tmp = x_mul(b-a,c-a);
    double ans = dot_mul(d-a,x_mul(b-a,c-a));
    return dot_mul(d-a,x_mul(b-a,c-a));
}
double area2(point a,point b,point c)
{
    return vlen(x_mul(b-a,c-a));
}
point slove(point a,point b,point c,point d)
{
    return (a+b+c+d)/4.0;
}
struct plane
{
    int v[3];
    plane(int a, int b, int c)
    {
        v[0] = a;
        v[1] = b;
        v[2] = c;
    }
    point Normal(const vector<point>& P) const
    {
        return x_mul(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
    }
    int CanSee(const vector<point>& P, int i) const
    {
        return dot_mul(P[i]-P[v[0]], Normal(P)) > 0;
    }
};
vector<point> P;
vector<plane> faces;
int vis[100][100];
vector<plane> CH3D(const vector<point>& P)
{
    int n = P.size();
    memset(vis,0,sizeof(vis));
    vector<plane> cur;
    cur.push_back(plane(0, 1, 2));
    cur.push_back(plane(2, 1, 0));
    for(int i = 3; i < n; i++)
    {
        vector<plane> nex;
        for(int j = 0; j < cur.size(); j++)
        {
            plane& f = cur[j];
            int res = f.CanSee(P, i);
            if(!res)
                nex.push_back(f);
            for(int k = 0; k < 3; k++)
                vis[f.v[k]][f.v[(k+1)%3]] = res;
        }
        for(int j = 0; j < cur.size(); j++)
        {
            for(int k = 0; k < 3; k++)
            {
                int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
                if(vis[a][b] != vis[b][a] && vis[a][b])
                    nex.push_back(plane(a, b, i));
            }
        }
        cur = nex;
    }
    return cur;
}
double rand01()
{
    return rand() / (double)RAND_MAX;
}
double randeps()
{
    return (rand01() - 0.5) * eps;
}
point add_noise(const point& p)
{
    return point(p.x + randeps(), p.y + randeps(), p.z + randeps());
}
point centroid()
{
    point c = P[0];
    double totv = 0;
    point tot(0,0,0);
    for(int i=0;i<faces.size();i++)
    {
        point p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
        double v = -area(p1, p2, p3, c);
        totv += v;
        tot = tot + slove(p1, p2, p3, c)*v;
    }

    return tot/totv;
}
double mindis(point c)
{
    double ans = 1e30;
    for(int i=0;i<faces.size();i++)
    {
        point p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
        ans = min(ans,fabs(-area(p1,p2,p3,c))/area2(p1,p2,p3));
    }
    //printf("%.2f\n",ans);
    return ans;
}
double ss(int n)
{
    P.resize(n);
    for(int i=0;i<n;i++)
    {
        scanf("%lf %lf %lf",&P[i].x,&P[i].y,&P[i].z);
        P[i] = add_noise(P[i]);
    }
//    for(int i=0;i<n;i++)
//        printf("%.2f %.2f %.2f\n",P[i].x,P[i].y,P[i].z);
    faces = CH3D(P);
    point tt = centroid();
    //printf("%.2f %.2f %.2f\n",tt.x,tt.y,tt.z);
    double ans = mindis(tt);
    return ans;
}
int main(void)
{
    freopen("asteroids.in","r",stdin);
    freopen("asteroids.out","w",stdout);
    scanf("%d",&n);
    double ans = ss(n);
    //printf("%.2f\n",ans);
    scanf("%d",&n);
    ans += ss(n);
    printf("%.5f\n",ans);
    return 0;
}