首页 > ACM题库 > HDU-杭电 > HDU 4449-Building Design-计算几何-[解题报告]HOJ
2015
07-16

HDU 4449-Building Design-计算几何-[解题报告]HOJ

Building Design

问题描述 :

You are a notable architect.
Recently, a company invites you to design their new building for them. It is not an easy task, because they make some strange claims to the new building:
1.Its shape should be the same as a convex polyhedron which they have given to you, though you can rotate it arbitrarily.
2.One face of the building should cling to the ground. (Of course! Have you ever seen a building only touch the ground by an edge, or a point, or even suspend in the air?)
3.They want the building to be as striking as possible. We call the highest point of the building as the “peak”. The higher the peak is the more striking the building will be.
4.If there are many designs meet all claims above, the occupied land of the building should be less. In another word, the building’s vertical projection area on the ground should be as small as possible.
Now, give you the relative positions of vertices of the building, please design the building and tell them H �C the height of the peak, as well as S – the projection area of the building in your best design.

输入:

There are several test cases in the input.
Each test case begins with one integer n (1 ≤ n ≤ 50), indicating the number of vertices of the building.
Then n lines follow. Each line contains three integers x, y, z (-10000 ≤ x, y, z ≤ 10000), separated by spaces, indicating the relative positions of one vertex of the building.
The input ends with n = 0.

输出:

There are several test cases in the input.
Each test case begins with one integer n (1 ≤ n ≤ 50), indicating the number of vertices of the building.
Then n lines follow. Each line contains three integers x, y, z (-10000 ≤ x, y, z ≤ 10000), separated by spaces, indicating the relative positions of one vertex of the building.
The input ends with n = 0.

样例输入:

6
1 0 0
-1 0 0
0 1 0
0 -1 0
0 0 1
0 0 -1
0

样例输出:

1.155 1.732

模板题~我表示直接复制粘贴了三维凸包的模板和平面旋转的模板~模板是别人的自己没改过~所以代码看起来比较别扭~

平面旋转的部分我是从http://hi.baidu.com/gdtangwu/item/ff5e56147a614b4de75e0696这里学来的~不过这个旋转有一个比较纠结的地方就是两个平面的法向量夹角是PI的时候会出现除0错误~所以旋转的时候要 if(sign(ang) != 0 && sign(ang-PI)!=0) 这样判断而不能只是sign(ang)
!= 0…………我因为这个点WA了5次吧大约~自己出了一组数据才发现问题的~

下面是我直接贴模板的代码~乱死了没救了~

【2013-01-28更新代码】

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include <algorithm>
using namespace std;
#define INF (1<<30)
#define PI acos(-1)
#define SET(a,b) memset(a,b,sizeof(a))
#define M 10010
//---------------
#define N 505
#define EPS 1e-8
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;
    }
};
struct _3DCH
{
    struct fac
    {
        int a,b,c;
        bool ok;
    };
    int n;
    pt P[N];
    int cnt;
    fac  F[N*8];
    int to[N][N];
    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));
    }
    double volume(pt a,pt b,pt c,pt d)
    {
        return (b-a)*(c-a)^(d-a);
    }
    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;
    }
    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((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(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 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])));
    }
    void get_panel(double &a,double &b,double &c,double &d,pt p1,pt p2,pt p3)
    {
        a=(p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y);
        b=(p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z);
        c=(p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x);
        d=-(a*p1.x+b*p1.y+c*p1.z);
    }
    inline double dist(pt p1, pt p2)
    {
        return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));
    }
}hull;
//----------------------------------------------
inline int sign(double d)
{
    if(d>EPS) return 1;
    if(d<-EPS) return -1;
    return 0;
}
struct point
{
    double x,y,z;
    point(double x=0,double y=0,double z=0): x(x),y(y),z(z) {}
    point operator - (point p)
    {
        return point(x-p.x,y-p.y,z-p.z);
    }
    point operator + (point p)
    {
        return point(x+p.x,y+p.y,z+p.z);
    }
    point operator / (double len)
    {
        return point(x/len,y/len,z/len);
    }
    point operator * (double len)
    {
        return point(x*len,y*len,z*len);
    }
    double operator ^ (point p)//点积
    {
        return x*p.x+y*p.y+z*p.z;
    }
    point operator * (point p)//叉积
    {
        return point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
    }
    double getlen()//长度
    {
        return sqrt(x*x+y*y+z*z);
    }
}con[N],ps[N],org;
double a,b,c,d;
int n;
inline point get_point(point u,point v,point p)//p在uv上的垂足
{
    double a,b,t;
    a=(v.x-u.x)*(v.x-u.x)+(v.y-u.y)*(v.y-u.y)+(v.z-u.z)*(v.z-u.z);
    b=(p.x-u.x)*(v.x-u.x)+(p.y-u.y)*(v.y-u.y)+(p.z-u.z)*(v.z-u.z);
    t=b/a;
    point ans;
    ans.x=v.x*t+(1-t)*u.x;
    ans.y=v.y*t+(1-t)*u.y;
    ans.z=v.z*t+(1-t)*u.z;
    return ans;
}
inline double dist(point a,point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
point rotate(point u,point v,point p,double ang)//点p沿uv右旋ang角度
{
    point root=get_point(u,v,p),e,r;
    point ans;
    e=(v-u)/dist(u,v);
    r=p-root;
    e=e*r;
    ans=r*cos(ang)+e*sin(ang)+root;
    return ans;
}
double inter_pro(point u1,point v1,point u2,point v2)
{
    return (v1.x-u1.x)*(v2.y-u2.y)-(v1.y-u1.y)*(v2.x-u2.x);
}
bool cmp(point a,point b)
{
    return a.y<b.y || a.y==b.y && a.x<b.x;
}
void Graham(point* pol, int n, point* con, int& len){  //st为起点的下标
 sort(pol, pol+ n, cmp);
 con[0] = pol[0];
 con[1] = pol[1];
 int top = 1;
 for(int i = 2; i < n; i++){
  while(top > 0 && inter_pro(con[top - 1], con[top], con[top - 1], pol[i]) <= 0) top--;
  con[++top] = pol[i];
 }
 int tmp = top;  //注意要赋值给tmp!
 for(int i = n - 2; i >= 0; i--){
  while(top > tmp && inter_pro(con[top - 1], con[top], con[top - 1], pol[i]) <= 0) top--;
  con[++top] = pol[i]; 
 }
 len = top;
}
double polyArea(point* ps, int n){
    //printf("n = %d\n",n);
 ps[n] = ps[0];
 int i;
 double ans=0;
 for(i = 0; i < n; i++){
  ans += (ps[i].x*ps[i+1].y-ps[i].y*ps[i+1].x);
 }
 return fabs(ans/2.0);
}
double solve()
{
    point tp(0,0,0),end(0,0,1),vec;
    double ang;
    int i,cn;
    if(sign(a)) tp.x=d/a;
    else if(sign(b)) tp.y=d/b;
    else if(sign(c)) tp.z=d/c;
    ps[n+1]=tp;
    vec=(point(a,b,c))*(end);
    if(sign(vec.x)==0) vec.x=0;
    if(sign(vec.y)==0) vec.y=0;
    if(sign(vec.z)==0) vec.z=0;
    ang=(a*end.x+b*end.y+c*end.z)/(point(a,b,c).getlen());
    ang=acos(ang);
    if(sign(ang)!=-0 && sign(ang-PI)!=0)
        for(i=0;i<n;i++)
            ps[i]=rotate(org,vec,ps[i],ang);
    for(int i=0;i<n;i++) ps[i].z=0;
    Graham(ps,n,con,cn);
    double ans=fabs(polyArea(con,cn));
    return ans;
}
int main()
{
    while (scanf("%d",&hull.n) && hull.n)
    {
        for(int i=0;i<hull.n;i++)
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        if(n==1)
        {
            printf("0.000 0.000\n");
            continue;
        }
        if(n==2)
        {
            printf("%.3lf 0.000\n",hull.dist(hull.P[0],hull.P[1]));
            continue;
        }
        hull.construct();
        n=hull.n;
        double ans1=0,ans2=1e100;
        if(hull.cnt==0)
        {
            double max0=0;
            for(int i=0;i<hull.n;i++)
            for(int j=0;j<hull.n;j++)
                max0=max(max0,hull.dist(hull.P[i],hull.P[j]));
            printf("%.3lf 0.000\n",max0);
            continue;
        }
        n=hull.n;
        for(int i=0;i<hull.cnt;i++)
        {
            for(int j=0;j<hull.n;j++)
                ps[j].x=hull.P[j].x,
                ps[j].y=hull.P[j].y,
                ps[j].z=hull.P[j].z;
            hull.get_panel(a,b,c,d,hull.P[hull.F[i].a],hull.P[hull.F[i].b],hull.P[hull.F[i].c]);
            d=-d;
            double dist=0;
            for(int j=0;j<hull.n;j++)
                dist=max(dist,hull.ptoface(hull.P[j],i));
            double ans=solve();
            if(dist>ans1) ans1=dist,ans2=ans;
            else if(fabs(dist-ans1)<EPS) if(ans2>ans) ans1=dist,ans2=ans;
        }
        printf("%.3f %.3f\n",ans1,ans2);
    }
    return 0;
}

版权声明:本文为博主原创文章,未经博主允许不得转载。

参考:http://blog.csdn.net/d891320478/article/details/8127374