三维计算几何之三维凸包增量法-蒲公英云

三维凸包

三维凸包就是将凸包放在三维中求。在三维空间中有一堆点，现求一个多面体将所有点全部包住的最小凸多面体。这里可以类比一下二维的凸包。
这个凸多面体有一个性质：他是包住所有点的多面体中表面积最小的。

三维凸包的求法

三维凸包的求法可以类比二维凸包，但三维凸包的核心还是暴力，暴力的去求。这里主要用一个增量法的思想。
这里要清楚组成三维凸包的不是二维凸包中的一条条线段，在这是一个个平面，由这些平面组成的凸多面体就是三维凸包。
和二维一样，我们先确定一个平面，然后不断地往里加点，每次加一个点，更新一下凸包，直到所有点都加进来了，整个的凸包也就求完了。这就是增量法的思想，不断地一个一个加点。

然后问题就是如何更新凸包呢，使得加入一个点后，更新为一个新的多面体，是i+1个点的凸包。

如上图：每当新加一个点pr时，将这个点想成一个光源（太阳），向四周散发光线（如图二），照到原先的凸多面体CH上，所有能被照到的平面是要删去的面，所有照不到的面是要留下来的面。然后看所有边，每条边都在两个平面上（两平面的交线），当光线照到他所在平面时，代表该边被照到一次，然后我们找出所有被照到一次的边，说明他们是能照到和不能照到的分界线，所以我们新加一个该边和新加入的pr点所构成的一个平面（如图三），就是新凸多面体。

大体的步骤就是如上所说的，但还有很多细节需要注意：存每个平面时都只存三角形平面，即存三个顶点，但若存在n边形的平面怎么办，我们可以把他划分成n-2个三角形，然后存下来，因为任何的多边形都可以通过三角剖分得到多个三角形。
其次如何判断一个点能照到的平面呢？我们可以判断这个点在平面上方还是下方（通过和法向量的点积判断，具体在基础知识中），在上面说明能照到，在下面说明照不到，所以这里的平面是有方向的，因此我们存平面时一定要严格按照一个方向存三个点，这里统一规定逆时针存。
然后就是判断每条边被照到了几次，我们加个bool二维数组，g[a][b]为1的话说明a到b的边被照到一次，所以当g[a][b]为1，g[b][a]为0时说明ab这条边只被照到了一次。（因为都是逆时针存点，所以这个面是ab边时，邻面就是ba边）
还有就是如何只存三角形平面呢？我们给每个点都加一个微小扰动，就是让他的坐标发生一个很细微很细微的变化，使得不存在四点共面的情况，然后任意三个点就可以构成一个唯一平面了。（注意加了微小扰动之后就不能用a和b的绝对值差eps内就算相等了，必须严格大于等于小于判断）

下面具体可以看代码模板了解详细步骤

代码模板

#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
const int N = 210;
const double eps = 1e-12;    //精度要求高一点
int n,m;    //n是总点数，m是凸包中平面的数量
bool g[N][N];    //g用来判断一条边被照到几次
double rand_eps()    //用rand函数来生成一个非常小的随机数
{ 
    return ((double)rand() / RAND_MAX - 0.5) * eps;    //用rand生成一个-0.5到0.5之间的数，再乘eps，就得到了一个非常小的随机数
}
struct Point{     //定义点的结构体
    double x,y,z;    //xyz三个坐标
    void shake()    //微小扰动，给每个坐标都加一个极小的随机数
    { 
        x += rand_eps(), y += rand_eps(), z += rand_eps();
    }
    Point operator-(Point t)    //重载一下减号运算符
    { 
        return { x-t.x, y-t.y, z-t.z};
    }
    Point operator*(Point t)    //叉乘
    {     
        return { y*t.z - t.y*z, t.x*z - x*t.z, x*t.y - y*t.x};
    }
    double operator&(Point t)    //点积
    { 
        return x*t.x + y*t.y + z*t.z;
    }
    double len()        //求向量的模长，三个坐标也能存向量
    { 
        return sqrt(x*x + y*y + z*z);
    }
}p[N];    //p来存所有点
struct Plane{     //定义平面的结构体
    int v[3];    //三个顶点
    Point norm()    //求法向量
    { 
        return (p[v[1]] - p[v[0]]) * (p[v[2]] - p[v[0]]);    //平面中两向量的叉积
    }
    bool above(Point t )    //判断一个点是否在平面上方
    { 
        return ((t-p[v[0]]) & norm()) >= 0;    //用向量和法向量的点积判断
    }
    double area()        //求一平面的面积
    { 
        return norm().len() / 2;    //法向量的模长除以2
    }
}plane[N],tp[N];    //plane存凸包上的平面，tp用来更新凸包，每次凸包上要留的平面和新加的平面都存进tp，要删的平面不存，最后将tp在复制给plane，就实现了凸包的更新
void convex()
{ 
    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 fg = plane[j].above(p[i]);    //判断这个点是否在该平面上方
            if(!fg)        //如果是下方的话，说明照不到
                tp[cnt++] = plane[j];    //存进tp数组
            for(int k = 0;k < 3;k++)    //循环该平面的三条边
                g[plane[j].v[k]][plane[j].v[(k+1)%3]] = fg;    //ab边照不照得到情况赋值给g[a][b]
        }
        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])        //判断该边是否被照到了一次，即是否是交界线的边
                    tp[cnt++] = { a,b,i};    //若是，加新平面abi，ab一定是逆时针的，i在后面
            }
        }
        m = cnt;    //将tp再赋值给plane
        for(int j = 0;j < m;j++)
            plane[j] = tp[j];
    }
}
int main()
{ 
    cin >> n;
    for(int i = 0;i < n;i++)
    { 
        cin >> p[i].x >> p[i].y >> p[i].z;    //输入n个点
        p[i].shake();    //微小扰动
    }
    convex();    //求三维凸包
    double ans = 0;        //求面积和
    for(int i = 0;i < m;i++)    //循环最终凸包上的m个平面
        ans += plane[i].area();    //将平面的面积加和
    printf("%.6lf\n",ans);    //输出答案
    return 0;
}