Skip to content

三维凸包 ​

cpp
namespace Geometry_3D{
    const double pi=acos(-1);
    const double eps=1e-12;
    double rand_eps(){
        return ((double)rand()/RAND_MAX-0.5)*eps;
    }
    struct Point_3D{
        double x,y,z;
        Point_3D(double x=0,double y=0,double z=0):x(x),y(y),z(z){}
        void shake(){
            x+=rand_eps(),y+=rand_eps(),z+=rand_eps();
        }
        double len(){
            return sqrt(x*x+y*y+z*z);
        }
    };
    typedef Point_3D Vector_3D;

    Vector_3D operator+(const Vector_3D &A, const Vector_3D &B) { 
        return Vector_3D(A.x + B.x, A.y + B.y, A.z + B.z); 
    }
    Vector_3D operator-(const Point_3D &A, const Point_3D &B) { 
        return Vector_3D(A.x - B.x, A.y - B.y, A.z - B.z); 
    }
    Vector_3D operator*(const Vector_3D &A, double p) { 
        return Vector_3D(A.x * p, A.y * p, A.z * p); 
    }
    Vector_3D operator/(const Vector_3D &A, double p) { 
        return Vector_3D(A.x / p, A.y / p, A.z / p); 
    }

    double Dot_3D(const Vector_3D &A, const Vector_3D &B) { 
        return A.x * B.x + A.y * B.y + A.z * B.z; 
    }
    Vector_3D Cross_3D(const Vector_3D &A, const Vector_3D &B) { 
        return Vector_3D(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x); 
    }


    int cmp(double x,double y){//比较函数:精度为eps
        if(fabs(x-y)<eps) return 0;//x与y相等
        if(x<y) return -1;//x小于y
        return 1;//x大于y
    }

    Point_3D q[N];

    struct Plane{
        int v[3];
        Vector_3D Normal(){
            return Cross_3D(q[v[1]] - q[v[0]], q[v[2]] - q[v[0]]);
        }
        double area(){
            return Normal().len()/2;
        }
        bool above(int i){
            return Dot_3D(q[i]-q[v[0]],Normal())>=0;
        }
    };

    Plane plane[N],np[N];
    int m=0;
    bool g[N][N];
    void get_convex_3d(int n){
        plane[m++]={0,1,2};
        plane[m++]={2,1,0};
        for(int i=3;i<n;i++){
            int cnt=0;
            for(int j=0;j<m;j++){
                bool t=plane[j].above(i);
                if(!t)np[cnt++]=plane[j];
                for(int k=0;k<3;k++){
                    g[plane[j].v[k]][plane[j].v[(k+1)%3]]=t;
                }
            }
            for(int j=0;j<m;j++){
                for(int k=0;k<3;k++){
                    int a=plane[j].v[k],b=plane[j].v[(k+1)%3];
                    if(g[a][b]&&!g[b][a]){
                        np[cnt++]={a,b,i};
                    }
                }
            }
            m=cnt;
            for(int j=0;j<m;j++)plane[j]=np[j];
        }
    }
}
using namespace Geometry_3D;
void solve(){
    int n;
    cin>>n;
    for(int i=0;i<n;i++){
        cin>>q[i].x>>q[i].y>>q[i].z;
        q[i].shake();
    }
    get_convex_3d(n);
    double res=0;
    for(int i=0;i<m;i++){
        res+=plane[i].area();
    }
    cout<<setprecision(3)<<fixed<<res<<"\n";
}