Luogu 4724 三维凸包

Posted 2021-02-09 jklover

Luogu 4724 三维凸包

• 增量法,维护当前凸包,每次加入一个点 (P) ,视其为点光源,将可见面删去,新增由"晨昏线"(分割棱)与 (P) 构成的平面.
• 注意每个平面表面积为其三个端点算出的法向量模长一半.
• 参考讲解.

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pii pair<int,int>
inline int read()
{
    int x=0;
    bool pos=1;
    char ch=getchar();
    for(;!isdigit(ch);ch=getchar())
        if(ch=='-')
            pos=0;
    for(;isdigit(ch);ch=getchar())
        x=x*10+ch-'0';
    return pos?x:-x;
}
const int MAXN=2e3+10;
const double eps=1e-9;
double Rand()
{
    return rand()/(double)RAND_MAX;
}
int dcmp(double x)
{
    return fabs(x)<eps?0:(x>0?1:-1);
}
struct point{
    double x,y,z;
    point(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
    void randtrans()
        {
            x+=(Rand()-0.5)*eps;
            y+=(Rand()-0.5)*eps;
            z+=(Rand()-0.5)*eps;
        }
};
point A[MAXN];
struct Vector{
    double x,y,z;
    Vector(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
    double modulus()
        {
            return sqrt(x*x+y*y+z*z);
        } 
};
Vector vec(point a,point b)
{
    return Vector(b.x-a.x,b.y-a.y,b.z-a.z);
}
Vector cross(Vector a,Vector b)
{
    return Vector(a.y*b.z-b.y*a.z,-(a.x*b.z-b.x*a.z),a.x*b.y-b.x*a.y);
}
double cdot(Vector a,Vector b)
{
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
struct Plane{
    int v[3];
    Vector Normal()
        {
            return cross(vec(A[v[0]],A[v[1]]),vec(A[v[0]],A[v[2]]));
        }
    double area()
        {
            return Normal().modulus()/2.0;
        }
};
Plane getplane(int a,int b,int c)
{
    Plane res;
    res.v[0]=a;
    res.v[1]=b;
    res.v[2]=c;
    return res;
}
bool insight(Plane a,point p)
{
    Vector w=a.Normal();
    return dcmp(cdot(w,vec(A[a.v[0]],p)))>0;
}
int n;
Plane f[MAXN],C[MAXN];
int vis[MAXN][MAXN];
double ConvexHullSurfaceArea()
{
    int cnt=0;
    f[++cnt]=(Plane){1,2,3};
    f[++cnt]=(Plane){3,2,1};
    for(int i=4,cc=0;i<=n;i++)
    {
        for(int j=1,v;j<=cnt;j++)
        {
            if(!(v=insight(f[j],A[i]))) 
                C[++cc]=f[j];
            for(int k=0;k<3;k++) 
                vis[f[j].v[k]][f[j].v[(k+1)%3]]=v;
        }
        for(int j=1;j<=cnt;j++)
            for(int k=0;k<3;k++)
            {
                int x=f[j].v[k],y=f[j].v[(k+1)%3];
                if(vis[x][y]&&!vis[y][x]) 
                    C[++cc]=(Plane){x,y,i};
            }
        for(int j=1;j<=cc;j++) 
            f[j]=C[j];
        cnt=cc;cc=0;
    }
    double ans=0;
    for(int i=1;i<=cnt;++i)
        ans+=f[i].area();
    return ans;
}
int main()
{
    n=read();
    for(int i=1;i<=n;++i)
        {
            scanf("%lf%lf%lf",&A[i].x,&A[i].y,&A[i].z);
            A[i].randtrans();
        }
    double ans=ConvexHullSurfaceArea();
    printf("%.3lf
",ans);
    return 0;
}