HDU 4449 Building Design 三维凸包+空间坐标转换＋二维凸包

来源：互联网发布：青岛seo顾问公司编辑：程序博客网时间：2024/06/05 17:44
题意：求点到面的最短距离的同时求出投影面积的最小套一下模版就好了
代码：
#include <iostream>#include <cstdio>#include <cmath>#include <algorithm>#define eps 1e-7using namespace std;const double inf=0x3f3f3f3f;const int MAXV=80;const double EPS = 1e-9;const double pi=acos(-1.0);//三维点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 + (const pt p1) {        return pt(x + p1.x, y + p1.y, z + p1.z);    }    pt operator *(const pt p){        return pt(y*p.z-z*p.y,z*p.x-x*p.z, x*p.y-y*p.x);    }    pt operator *(double d)    {        return pt(x*d,y*d,z*d);    }        pt operator / (double d)    {        return pt(x/d,y/d,z/d);    }    double operator ^ (pt p) {        return x*p.x+y*p.y+z*p.z;    //点乘    }    double len(){        return sqrt(x*x+y*y+z*z);    }};struct pp{    double x,y;    pp(){}    pp(double x,double y):x(x),y(y){}};struct _3DCH {    struct fac {        int a, b, c;    //表示凸包一个面上三个点的编号        bool ok;        //表示该面是否属于最终凸包中的面    };        int n;    //初始点数    pt P[MAXV];    //初始点        int cnt;    //凸包表面的三角形数    fac F[MAXV*8]; //凸包表面的三角形        int to[MAXV][MAXV];        pt Cross3(pt a,pt p){        return pt(a.y*p.z-a.z*p.y, a.z*p.x-a.x*p.z, a.x*p.y-a.y*p.x);    //叉乘    }    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(Cross3((b-a),(c-a)));    //三角形面积*2    }    double volume(pt a, pt b, pt c, pt d) {        return Cross3((b-a),(c-a))^(d-a);    //四面体有向体积*6    }    //三维点积    double Dot3( pt u, pt v )    {        return u.x * v.x + u.y * v.y + u.z * v.z;    }        //平面的法向量    pt pvec(pt a,pt b,pt c)    {        return (Cross3((a-b),(b-c)));    }    //点到面的距离    double Dis(pt a,pt b,pt c,pt d)    {        return fabs(pvec(a,b,c)^(d-a))/vlen(pvec(a,b,c));    }    //正：点在面同向    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 Cross3(m , n) ^ t;    }        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 = 1;                to[p][b] = to[a][p] = to[b][a] = cnt;                F[cnt++] = add;            }        }    }        void dfs(int p, int cur) {        F[cur].ok = 0;        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 = 0;        if (n < 4)            return;                /*********此段是为了保证前四个点不公面，若已保证，可去掉********/        bool sb = 1;        //使前两点不公点        for (int i = 1; i < n; i++) {            if (vlen(P[0] - P[i]) > eps) {                swap(P[1], P[i]);                sb = 0;                break;            }        }        if (sb)return;                sb = 1;        //使前三点不公线        for (int i = 2; i < n; i++) {            if (vlen(Cross3((P[0] - P[1]) , (P[1] - P[i]))) > eps) {                swap(P[2], P[i]);                sb = 0;                break;            }        }        if (sb)return;                sb = 1;        //使前四点不共面        for (int i = 3; i < n; i++) {            if (fabs(Cross3((P[0] - P[1]) , (P[1] - P[2])) ^ (P[0] - P[i])) > eps) {                swap(P[3], P[i]);                sb = 0;                break;            }        }        if (sb)return;        /*********此段是为了保证前四个点不公面********/                        fac add;        for (int i = 0; i < 4; i++) {            add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, add.ok = 1;            if (ptof(P[i], add) > 0)                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 = 4; i < n; i++) {            for (int j = 0; j < cnt; j++) {                if (F[j].ok && ptof(P[i], F[j]) > eps) {                    dfs(i, j);                    break;                }            }        }        int tmp = cnt;        cnt = 0;        for (int i = 0; i < tmp; i++) {            if (F[i].ok) {                F[cnt++] = F[i];            }        }    }        //表面积    double area() {        double ret = 0.0;        for (int i = 0; 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(0, 0, 0);        double ret = 0.0;        for (int i = 0; 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 = 0;        for (int i = 0; i < cnt; i++) {            bool nb = 1;            for (int j = 0; j < i; j++) {                if (same(i, j)) {                    nb = 0;                    break;                }            }            ans += nb;        }        return ans;    }    pt centroid(){        pt ans(0,0,0),o(0,0,0);        double all=0;        for(int i=0;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;    }    double res(){        pt a=centroid();        double _min=1e10;        for(int i=0;i<cnt;++i){            double now=Dis(P[F[i].a],P[F[i].b],P[F[i].c],a);            _min=min(_min,now);        }        return _min;    }    double ptoface(pt 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])));    }}lou;bool mult(pp sp,pp ep,pp op){    return (sp.x-op.x)*(ep.y-op.y)>=(ep.x-op.x)*(sp.y-op.y);}double Cross(pp a,pp b,pp c){    return (c.x-a.x)*(b.y-a.y) - (b.x-a.x)*(c.y-a.y);}bool cmp(pp a,pp b){    if(a.y==b.y)return a.x<b.x;    return a.y<b.y;}int n,res[60],top;pp ps[60];void Graham(){    int len;    n=lou.n;    top=1;    sort(ps,ps+n,cmp);    if(n==0)return;res[0]=0;    if(n==1)return;res[1]=1;    if(n==2)return;res[2]=2;    for(int i=2;i<n;i++){        while(top&&mult(ps[i],ps[res[top]],ps[res[top-1]]))            top--;        res[++top]=i;    }    len=top;    res[++top]=n-2;    for(int i=n-3;i>=0;i--){        while(top!=len&&mult(ps[i],ps[res[top]],ps[res[top-1]]))top--;        res[++top]=i;    }}inline pt get_point(pt st,pt ed,pt tp){    double t1=(tp-st)^(ed-st);    double t2=(ed-st)^(ed-st);    double t=t1/t2;    pt ans=st + ((ed-st)*t);    return ans;}inline double dist(pt st,pt ed){    return sqrt((ed-st)^(ed-st));}pp rotate(pt st,pt ed,pt tp,double A){    pt root=get_point(st,ed,tp);    pt e=(ed-st)/dist(ed,st);    pt r=tp-root;    pt vec=e*r;    pt ans=r*cos(A)+vec*sin(A)+root;    return pp(ans.x,ans.y);}int main(){    while (scanf("%d",&lou.n)&&lou.n) {        for (int i=0; i<lou.n; ++i)            scanf("%lf%lf%lf",&lou.P[i].x,&lou.P[i].y,&lou.P[i].z);        lou.construct();        double ansh=0,ansa=inf;        if(lou.n<=2)        {            printf("0.000 0.000\n");        }        else if(lou.n==3)        {            ansh=0;            ansa=(lou.P[1]-lou.P[0])^(lou.P[2]-lou.P[0]);            ansa/=2.0;            printf("%.3lf %.3lf\n",ansh,ansa);                    }        else {            for (int i=0; i<lou.cnt; ++i) {                pt p1=(lou.P[lou.F[i].b]-lou.P[lou.F[i].a])*(lou.P[lou.F[i].c]-lou.P[lou.F[i].a]);                pt e=pt(0,0,1);                pt vec=p1*e;                double A=p1^e/p1.len();                A=acos(A);                if(fabs(A-pi)>EPS&&fabs(A)>EPS){                    pt s=pt(0,0,0);                    for (int k=0; k<lou.n; ++k) ps[k]=rotate(s,vec,lou.P[k],A);                }                else                {                    for(int k=0; k<lou.n; ++k) ps[k].x=lou.P[k].x,ps[k].y=lou.P[k].y;                }                double h=0;                for (int j=0; j<lou.n; ++j)                    h=max(lou.ptoface(lou.P[j],i),h);                if (h<ansh)                    continue;                Graham();                double a=0;                for (int k=1; k<top-1; ++k){                    a+=fabs(Cross(ps[res[0]],ps[res[k]],ps[res[k+1]]));                }                a/=2.0;                if (fabs(h-ansh)<EPS) {                    ansa=min(ansa,a);                }                else                    ansa=a;                ansh=h;            }            printf("%.3lf %.3lf\n",ansh,ansa);        }    }}
0 0