HDOJ4449解题报告

题目地址：

　　http://acm.hdu.edu.cn/showproblem.php?pid=4449

题目概述：

　　在空间中给出一个凸多面体，现在让你旋转这个多面体使得它有一面与平面z=0重合，此时多面体有一个在平面z=0上的投影以及多面体上的点到平面z=0的最大距离，求旋转过程中所有这些距离中最大的值，距离相同的投影面积应该尽可能的小。

大致思路：

　　计算几何显然啦，首先来一个三维凸包找到这个凸多面体，然后就要开始旋转了。

　　小知识：把平面P：ax+by+cz=d 旋转、平移到z=0，方法是设一个向量V为向量（a,b,c）与（0,0,1）的叉积，则P绕着V做逆时针旋转，角度为前面两个向量的夹角，然后对这个凸多面体的所有点做同样的操作即可。

　　这样问题就只剩下求投影面积了，既然你会三维凸包，那底面显然直接上二维凸包求一个面积即可。

复杂度分析：

　　计算几何大部分题时间都不是大问题，代码复杂度才是最大的！！！！

代码：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <ctime>
#include <map>
#include <assert.h>
#include <stack>
#include <set>
#include <bitset>
#include <queue>
#include <iomanip>
#include <cstring>
#include <algorithm>
using namespace std;

#define sacnf scanf
#define scnaf scanf
#define scnafi scanfi
#define pb push_back
#define Len size()
#define gchar getchar()
#define FOR(i,j,k) for(int (i)=(j);(i)<=(k);++(i))
#define mes(a,b) memset((a),(b),sizeof(a))
#define scanfi(a) scanf("%d",&(a))
#define scanfti(a,b) scanf("%d%d",&(a),&(b))
#define println(a) printf("%d
",(a))
#define printIs(a) printf((a)?"Yes
":"No
")
#define print_b printf("
")
#define IsNum(x) '0'<=(x)&&(x)<='9'
#define lt dir*2
#define rt dir*2+1

#define maxn 60
#define maxm 26
#define inf 1061109567
//const long long inf=1e15;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define E 2.718281828459
const double PI=acos(-1.0);
#define mod 1000
//const long long mod=1000000007;
#define MAXNUM 1000000000

typedef long long ll;
typedef unsigned long long ulld;
ll Abs(ll x) {return (x<0)?-x:x;}
int Read() {char ch;int res=0;while(ch=getchar(),!(IsNum(ch)));while(IsNum(ch)) res=res*10+ch-'0',ch=getchar();return res;}

int sgn(double x)
{
    if(x<-eps) return -1;
    if(x>eps) return 1;
    return 0;
}

struct vec
{
    double x,y;
    vec() {x=y=0;}
    vec(double _x,double _y) {x=_x;y=_y;}

    vec operator - (vec v) {return vec(x-v.x,y-v.y);}
};

double cross(vec a,vec b) {return a.x*b.y-a.y*b.x;}

bool cmpXY(vec a,vec b)
{
    if(sgn(a.x-b.x)) return a.x<b.x;
    return a.y<b.y;
}

int convex_hull(vec* v,int n,vec *z)
{
    sort(v,v+n,cmpXY);
    int m=0;
    FOR(i,0,n-1)
    {
        while(m>1&&cross(z[m-1]-z[m-2],v[i]-z[m-2])<=0) --m;
        z[m++]=v[i];
    }
    int k=m;
    for(int i=n-2;i>=0;--i)
    {
        while(m>k&&cross(z[m-1]-z[m-2],v[i]-z[m-2])<=0) --m;
        z[m++]=v[i];
    }
    if(n>1) --m;
    return m;
}

#define PR 1e-8

struct TPoint
{
    double x,y,z;
    TPoint() {}
    TPoint(double _x,double _y,double _z):x(_x),y(_y),z(_z) {}
    TPoint operator + (const TPoint p) {return TPoint(x+p.x,y+p.y,z+p.z);}
    TPoint operator - (const TPoint p) {return TPoint(x-p.x,y-p.y,z-p.z);}
    TPoint operator * (const TPoint p) {return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);}  //叉积
    TPoint operator * (double p) {return TPoint(x*p,y*p,z*p);}
    TPoint operator / (double p) {return TPoint(x/p,y/p,z/p);}
    double operator ^ (const TPoint p) {return x*p.x+y*p.y+z*p.z;}
} center;

struct fac
{
    int a,b,c;
    bool ok;
};

struct T3dhull
{
    int n;//初始点的个数
    TPoint ply[maxn];//初始点
    int trianglecnt;//凸包上三角形数
    fac tri[maxn*8];//凸包三角形
    int vis[maxn][maxn];//点i到点j是属于哪个面
    double dist(TPoint a) {return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}
    double area(TPoint a,TPoint b,TPoint c) {return dist((b-a)*(c-a));}//三角形面积*2
    double volume(TPoint a ,TPoint b,TPoint c,TPoint d) {return (b-a)*(c-a)^(d-a);}//四面体有向体积*6
    double ptoplane(TPoint &p,fac &f) //正：点在面同向
    {
        TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a];
        return (m*n)^t;
    }
    void deal(int p,int a,int b)
    {
        int f=vis[a][b];
        fac add;
        if(tri[f].ok)
        {
            if(ptoplane(ply[p],tri[f])>PR) dfs(p,f);
            else
            {
                add.a=b,add.b=a,add.c=p,add.ok=1;
                vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt;
                tri[trianglecnt++]=add;
            }
        }
    }
    void dfs(int p,int cnt) //维护凸包，如果点p在凸包外更新凸包
    {
        tri[cnt].ok=0;
        deal(p,tri[cnt].b,tri[cnt].a);
        deal(p,tri[cnt].c,tri[cnt].b);
        deal(p,tri[cnt].a,tri[cnt].c);
    }
    bool same(int s,int e) //判断两个面是否为同一面
    {
        TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c];
        return fabs(volume(a,b,c,ply[tri[e].a]))<PR
             &&fabs(volume(a,b,c,ply[tri[e].b]))<PR
             &&fabs(volume(a,b,c,ply[tri[e].c]))<PR;
    }
    void construct() //构建凸包
    {
        trianglecnt=0;
        if(n<4) return;
        bool tmp=1;
        FOR(i,1,n-1) //前两点不共点
        {
            if(dist(ply[0]-ply[i])>PR)
            {
                swap(ply[1],ply[i]);
                tmp=0;break;
            }
        }
        if(tmp) return;
        tmp=1;
        FOR(i,2,n-1) //前三点不共线
        {
            if(dist((ply[0]-ply[1])*(ply[1]-ply[i]))>PR)
            {
                swap(ply[2],ply[i]);
                tmp=0;break;
            }
        }
        if(tmp) return;
        tmp=1;
        FOR(i,3,n-1) //前四点不共面
        {
            if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)
            {
                swap(ply[3],ply[i]);
                tmp=0;break;
            }
        }
        if(tmp) return;
        fac add;
        FOR(i,0,3) //构建初始四面体
        {
            add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1;
            if(ptoplane(ply[i],add)>0) swap(add.b,add.c);
            vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;
            tri[trianglecnt++]=add;
        }
        FOR(i,4,n-1) //构建更新凸包
        {
            FOR(j,0,trianglecnt-1)
            {
                if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>PR)
                {
                    dfs(i,j);break;
                }
            }
        }
        int cnt=trianglecnt;
        trianglecnt=0;
        FOR(i,0,cnt-1)
            if(tri[i].ok) tri[trianglecnt++]=tri[i];
    }
    double area() //表面积
    {
        double ret=0;
        FOR(i,0,trianglecnt-1) ret+=area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
        return ret/2.0;
    }
    double volume() //体积
    {
        TPoint p(0,0,0);
        double ret=0;
        FOR(i,0,trianglecnt-1) ret+=volume(p,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
        return fabs(ret/6);
    }
    int facetri() {return trianglecnt;} //表面三角形数
    int facepolygon() //表面多边形数
    {
        int ans=0,j,k;
        FOR(i,0,trianglecnt-1)
        {
            for(j=0,k=1;j<i;j++)
                if(same(i,j)) {k=0;break;}
            ans+=k;
        }
        return ans;
    }
    TPoint barycenter() //重心
    {
        TPoint ans(0,0,0),o(0,0,0);
        double all=0;
        FOR(i,0,trianglecnt-1)
        {
            double vol=volume(o,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
            ans=ans+(o+ply[tri[i].a]+ply[tri[i].b]+ply[tri[i].c])/4.0*vol;
            all+=vol;
        }
        return ans/all;
    }
    double ptoface(TPoint p,int i) //点到面的距离
    {
        return fabs(volume(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c],p)
                    /area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]));
    }
} a;

TPoint GetPoint(TPoint st,TPoint ed,TPoint tp)
{
    double t1=(tp-st)^(ed-st);
    double t2=(ed-st)^(ed-st);
    double t=t1/t2;
    return st+((ed-st)*t);
}

TPoint Rotate(TPoint st,TPoint ed,TPoint tp,double A)    //tp以ed-st为轴旋转A度
{
    TPoint root=GetPoint(st,ed,tp);
    TPoint e=(ed-st)/a.dist(ed-st);
    TPoint r=tp-root;TPoint vec=e*r;
    /*TPoint ans=r*cos(A)+vec*sin(A)+root;
    printf("%.f %.f %.f
",ans.x,ans.y,ans.z);*/
    return r*cos(A)+vec*sin(A)+root;
}

TPoint tmp[maxn];
vec tmp1[maxn],tmp2[maxn];

int main()
{
    //freopen("data.in","r",stdin);
    //freopen("std.out","w",stdout);
    //clock_t st=clock();
    int n;
    while(~scanfi(n))
    {
        if(n==0) break;
        int x,y,z;double ansH=0,ansS=0;
        FOR(i,0,n-1)
        {
            scnaf("%d%d%d",&x,&y,&z);
            a.ply[i]=(TPoint){x*1.0,y*1.0,z*1.0};
        }
        a.n=n;a.construct();
        FOR(i,0,a.trianglecnt-1)
        {
            FOR(j,0,n-1) tmp[j]=a.ply[j];
            TPoint p1=(tmp[a.tri[i].b]-tmp[a.tri[i].a])*(tmp[a.tri[i].c]-tmp[a.tri[i].a]);
            TPoint e(0,0,1);
            TPoint vec=p1*e;
            double A=p1^e/a.dist(p1);A=acos(A);

            if(sgn(A)!=0&&sgn(A-PI)!=0)
            {
                TPoint p0(0,0,0);
                FOR(j,0,n-1)
                    tmp[j]=Rotate(p0,vec,tmp[j],A);
            }

            double H=0;FOR(j,0,n-1) H=max(H,a.ptoface(a.ply[j],i)),tmp1[j].x=tmp[j].x,tmp1[j].y=tmp[j].y;

            sort(tmp1,tmp1+n,cmpXY);int tn=0;
            FOR(j,0,n-1)
            {
                int t=j+1;
                while(t<n&&sgn(tmp1[t].x-tmp1[j].x)==0&&sgn(tmp1[t].y-tmp1[j].y)==0) t++;
                tmp1[tn++]=tmp1[j];j=t-1;
            }
            int m=convex_hull(tmp1,tn,tmp2);

            double S=0;
            FOR(j,1,m-2) S+=cross(tmp2[j]-tmp2[0],tmp2[j+1]-tmp2[0]);
            S=fabs(S);

            if(sgn(ansH-H)<0)
            {
                ansH=H;ansS=S;
            }
            else if(sgn(ansH-H)==0)
            {
                if(sgn(ansS-S)>0) ansS=S;
            }
        }
        printf("%.3f %.3f
",ansH,ansS/2);
    }
    //clock_t ed=clock();
    //printf("

Time Used : %.5lf Ms.
",(double)(ed-st)/CLOCKS_PER_SEC);
    return 0;
}