hdu4449Building Design（三维凸包+平面旋转）

链接

看了几小时也没看懂代码表示的何意。。无奈下来问问考研舍友。

还是考研舍友比较靠谱，分分钟解决了我的疑问。

可能三维的东西在纸面上真的不好表示，网上没有形象的题解，只有简单"明了"的讲解。

这题说起来很简单，求下三维凸包，枚举每一个面，进行坐标旋转，使得当前面作为xoy面时的其他坐标，然后求下投影面的凸包的面积。

为什么要旋转面而不直接算点到面的距离，是因为投影的面积没有办法算。

面旋转时是怎么旋转的，首先求得当前面的法向量p1，再求得它与向量e(0,0,1)的法向量pp,所有的点都是绕pp这个向量旋转的，并且旋转的角度是p1与e的夹角。

如果还有跟我一样不懂空间叉积代表的是什么的同学向后看---》法向量 = 平面上两条不平行的向量的叉积。

这样求出每个点的都是相对于当前平面的点，只取点的x,y坐标即每个点的投影点-->求二维凸包-->求面积。

#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>>2;
//三维点
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=0,double y=0):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*8]; //凸包表面的三角形

    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 = 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)
    {
        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 = 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((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((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 (point3of(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 && point3of(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 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(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;
    }

} hull;
point3 get_point(point3 st,point3 ed,point3 tp)//tp在直线st-en上的垂足
{
    double t1=(tp-st)^(ed-st);
    double t2=(ed-st)^(ed-st);
    double t=t1/t2;
    point3 tt = (ed-st)*t;
    point3 ans=st + tt;
    return ans;
}
point3 rotate(point3 st,point3 ed,point3 tp,double A)//将点tp绕st-ed逆时针旋转A度  从ed往st看为逆时针
{
    point3 root=get_point(st,ed,tp);
    point3 e=(ed-st)/hull.vlen(ed-st);
    point3 r=tp-root;
    point3 vec=e*r;
    point3 ans=r*cos(A)+vec*sin(A)+root;
    return ans;
}
int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    return x<0?-1:1;
}
double cross(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}
double mul(point p0,point p1,point p2)
{
    return cross(p1-p0,p2-p0);
}
double dis(point a)
{
    return sqrt(a.x*a.x+a.y*a.y);
}
bool cmp(point a,point b)
{
    if(dcmp(mul(p[0],a,b))==0)
        return dis(a-p[0])<dis(b-p[0]);
    else
        return dcmp(mul(p[0],a,b))>0;
}
double Polyarea(int n,point p[])
{
    double area = 0;
    for(int i = 1 ; i < n-1 ; i++)
        area+=cross(p[i]-p[0],p[i+1]-p[0]);
    return fabs(area)/2;
}
int Graham(int n)
{
    int i,k = 0,top;
    point tmp;
    for(i = 0 ; i < n; i++)
    {
        if(p[i].y<p[k].y||(p[i].y==p[k].y&&p[i].x<p[k].x))
            k = i;
    }
    if(k!=0)
    {
        tmp = p[0];
        p[0] = p[k];
        p[k] = tmp;
    }
    sort(p+1,p+n,cmp);
    ch[0] = p[0];
    ch[1] = p[1];
    top = 1;
    for(i = 2; i < n ; i++)
    {
        while(top>0&&dcmp(mul(ch[top-1],ch[top],p[i]))<0)
            top--;
        top++;
        ch[top] = p[i];
    }
    return top;
}
void solve()
{
    int i,j;
    double h = 0,sa = INF;
    int cnt = hull.cnt,n = hull.n;
    for(i = 0; i < cnt ; i++)
    {
        for(j = 0 ; j < n; j++)
            rp[j] = hull.P[j];//rp数组为待旋转点

        point3 p1 = (rp[hull.F[i].b]-rp[hull.F[i].a])*(rp[hull.F[i].c]-rp[hull.F[i].a]);//平面的法向量，注意法向量的方向。
        point3 e = point3(0,0,1);//z轴上的向量
        point3 vec = p1*e;//垂直于p1和e的向量，即所有点将要绕其旋转的向量

        double A = p1^e/hull.vlen(p1);
        A = acos(A); //p1与e的夹角

        if(dcmp(A)!=0&&dcmp(A-pi)!=0)
        {
            point3 p0 = point3(0,0,0);
            for(j = 0 ; j < n; j++)
                rp[j] = rotate(p0,vec,rp[j],A);//绕直线p0-vec旋转
        }
        double tt = rp[hull.F[i].a].z;
        for(j = 0 ; j < n; j++)
            rp[j].z-=tt;
        double th = 0,ts;
        for(j = 0 ; j < n; j++)
        {
            th = max(th,hull.ptoface(hull.P[j],i));
        }
        for(j = 0 ; j< n; j++)
        {
            p[j].x = rp[j].x;
            p[j].y = rp[j].y;
        }
        int m = Graham(n);
        ch[++m] = ch[0];

        ts = Polyarea(m,ch);//cout<<ts<<endl;
        if(dcmp(th-h)>0||(dcmp(th-h)==0&&dcmp(ts-sa)<0))
        {
            h = th;
            sa = ts;
        }
    }
    printf("%.3f %.3f
",h,sa);
}
int main()
{
    int n,i;
    while(scanf("%d",&n)&&n)
    {
        hull.n = n;
        for(i = 0 ; i < n; i++)
        {
            scanf("%lf%lf%lf",&pp[i].x,&pp[i].y,&pp[i].z);
            hull.P[i] = pp[i];
        }
        hull.construct();
        solve();
    }
}