【模板】【计几】三维凸包

两道模板题

求三维凸包表面积：

洛谷p4724

求三维凸包表面多少个多边形：

hdu3662

 

#include<bits/stdc++.h>
using namespace std;
#define eps 1e-8
const int N = 2e3+9;
struct Point3{
    double x,y,z;
    Point3(){}
    Point3(double _x,double _y,double _z) : x(_x),y(_y),z(_z){}
    Point3 operator + (Point3 b){ return Point3(x+b.x,y+b.y,z+b.z); }
    Point3 operator - (Point3 b){ return Point3(x-b.x,y-b.y,z-b.z); }
    Point3 operator * (double k){ return Point3(x*k,y*k,z*k); }
};
#define Vector3 Point3
double dot(Vector3 a,Vector3 b){return a.x*b.x + a.y*b.y + a.z*b.z;}
Point3 cross(Vector3 a,Vector3 b){
    return Point3(a.y*b.z-a.z*b.y , a.z*b.x-a.x*b.z , a.x*b.y-a.y*b.x);
}
double len(Vector3 a){return sqrt(dot(a,a)); }
//面积的2倍
double area2(Point3 a,Point3 b,Point3 c){ return len(cross(b-a,c-a)); }
//四面体体积*6
double volume4(Point3 a,Point3 b,Point3 c,Point3 d){
    return dot( cross(b-a,c-a) , d-a );
}
struct CH3D{
    struct face{
        //三个点编号，逆时针
        int a,b,c;
        //是否在最终凸包上
        bool ok;
    };
    //n点数(0开始),cnt面数(1开始),face是凸包表面，都是三角形,g[x][y]记录x到y这个向量右边的面编号
    int n,cnt; Point3 P[N]; face F[2*N]; int g[N][N];
    //判断点和面的法向量是否同向
    double dblcmp(Point3 p,face f){
        Point3 m = P[f.b] - P[f.a];
        Point3 n = P[f.c] - P[f.a];
        Point3 t = p - P[f.a];
        return dot(cross(m,n),t);
    }
    //一直往下删除
    void dfs(int p,int a,int b){
        int f = g[a][b];
        if(F[f].ok == 0 ) return;
        if(dblcmp(P[p],F[f]) > eps){
            //如果p能看到面f，那么把f删除
            F[f].ok = 0;
            dfs(p,F[f].b,F[f].a);
            dfs(p,F[f].c,F[f].b);
            dfs(p,F[f].a,F[f].c);
        }
        else{
            //看不到则add一个面
            face add;
            add.a = b; add.b = a;
            add.c = p; add.ok = 1;
            ++cnt;
            g[p][b] = g[a][p] = g[b][a] = cnt;
            F[cnt] = add;
        }

    }
    void del(int p,int now){
        //删除凸包面
        F[now].ok = 0;
        dfs(p,F[now].b,F[now].a);
        dfs(p,F[now].c,F[now].b);
        dfs(p,F[now].a,F[now].c);
    }
    void create(){
        cnt = 0;
        memset(g,0,sizeof(g));
        if(n<4) return;
        bool find = 0;
        //前两个点不共点
        for(int i = 1;i<n;++i){
            if(len(P[0] - P[i]) > eps){
                swap(P[1],P[i]);
                find = 1;
                break;
            }
        }
        if(!find) return;
        find = 0;
        //前三个点不共线
        for(int i = 2;i<n;++i){
            if(len( cross(P[0] - P[1] , P[1]- P[i]) ) > eps){
                swap(P[2],P[i]);
                find = 1;
            }
        }
        if(!find ) return;
        find = 0;
        //前四个点不共面
        for(int i = 3;i<n;++i){
            if( fabs( dot( cross(P[0]-P[1],P[1]-P[2]),P[0]-P[i] )) > eps){
                swap(P[3],P[i]);
                find = 1;
                break;
            }
        }
        if(!find) return;
        for(int i = 0;i<4;++i){
            face add;
            add.a = (i+1)%4;
            add.b = (i+2)%4;
            add.c = (i+3)%4;
            add.ok = 1;
            //保证逆时针，就是法向量朝外，这样新点才可看到
            if(dblcmp(P[i],add) > 0 ) swap(add.b,add.c);
            ++cnt;
            g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = cnt;
            F[cnt] = add;
        }
        for(int i = 4;i<n;++i){
            for(int j = 1;j<=cnt;++j){
                //i点，j面
                if(F[j].ok && dblcmp(P[i],F[j]) > eps){
                    del(i,j);
                    break;
                }
            }
        }
        int tem = cnt;
        cnt = 0;
        for(int i = 1;i<=tem;++i){
            if(F[i].ok) F[++cnt] = F[i];
        }
    }
    //凸包表面积
    double area(){
        double res = 0;
        for(int i = 1;i<=cnt;++i){
            res += area2(P[F[i].a],P[F[i].b],P[F[i].c]);
        }
        return res/2.0;
    }
    //凸包体积
    double volume(){
        double res = 0;
        Point3 tem(0,0,0);
        for(int i = 1;i<=cnt;++i) res += volume4(tem,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(res/6.0);
    }
    //判断两面是否同一个面
    bool same(int s,int t){
        Point3 a = P[F[s].a];
        Point3 b = P[F[s].b];
        Point3 c = P[F[s].c];
        return fabs( volume4(a,b,c,P[F[t].a])) < eps &&
               fabs( volume4(a,b,c,P[F[t].b])) < eps &&
               fabs( volume4(a,b,c,P[F[t].c])) < eps;
    }
    //表面多边形个数
    int polygon(){
        int res = 0;
        for(int i = 1;i<=cnt;++i){
            bool ok = 1;
            for(int j = 1;j<i;++j){
                if(same(i,j)){
                    ok = 0;
                    break;
                }
            }
            res += ok;
        }
        return res;
    }
};
CH3D hull;
int main(){
    while(~scanf("%d",&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);
        }
        hull.create();
        printf("%.3f\n",hull.area());
        // printf("%d\n",hull.polygon());
    }
    return 0;
}