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

本文介绍了一个计算点到面最短距离及投影面积的算法,并详细展示了三维凸包构建的过程,包括如何计算体积、表面积等关键指标。

题意:求点到面的最短距离的同时求出投影面积的最小   套一下模版 就好了

代码:

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#define eps 1e-7
using 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);
        }
    }
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值