基于C++的最小凸多面体生成和可视化

0 前言

最小凸多面体的生产，主要是在给定一堆三维点的情况下，自动生成能包围这些三维点的最小凸多面体，也就是三维凸包的建立。同时，对于随机给出的一些空间点，判断是否在上述凸多面体内。
本文在对前人的代码进行说明的基础上，增加了交互和可视化的功能。（可在浏览完整代码的基础上再看详细说明）

参考程序：https://www.cnblogs.com/mypsq/p/4348140.html

1 整体思路

这里采用增量法进行求解，主要思路为：首先选择四个不共面的点形成一个四面体，然后每次增加一个新的点，分成两种情况：当新点在当前凸包内部时，忽略该点；当新点在当前凸包外部时，找到从该点可以“看见”的面，将这些面删除，然后对于这些面的每条边进行判断，看该点能否看见这些边的另一侧的面，若不行，则将该边与这个点相连构成新的面。
判断一个点P在凸包内还是凸包外，这里利用有向体积的方法。在存储面时，保证面的法线方向朝向凸包外部，若存在某一平面和点P所组成的四面体的有向体积为正，则P点在凸包外部，并且此面可被P点“看见”。

增量法详见：https://www.cnblogs.com/-sunshine/archive/2012/08/25/2656794.html

2 具体方案说明

注：该部分仅为程序说明，完整程序见后

2.1 点及其向量运算方法的定义

为方便表示，对最小凸多面体生成过程中的点定义为node类。由于三维点的坐标是向量的形式，向量之间的运算与实数有所不同，为了方便之后点与点之间的运算，对所用到的向量运算方法进行对应运算符的重载，包括向量的相加、相减、叉乘、点乘、乘以和除以实数等运算。

struct node
{
   
    double x,y,z,dis;
    node(){
   }
    node(double xx,double yy,double zz):x(xx),y(yy),z(zz){
   }
	
	//运算符重载
    node operator +(const node p)	//向量间相加操作
    {
   
        return node(x+p.x,y+p.y,z+p.z);
    }
    node operator -(const node p)	//向量间相减操作
    {
   
        return node(x-p.x,y-p.y,z-p.z);
    }
    node operator *(const node p)	//向量间叉乘操作
    {
   
        return node(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
    }
    node operator *(const double p)	//向量乘以一个数
    {
   
        return node(x*p,y*p,z*p);
    }
    node operator /(const double p)	//向量除以一个数
    {
   
        return node(x/p,y/p,z/p);
    }
    double operator ^(const node p)	//向量间点乘操作
    {
   
        return x*p.x+y*p.y+z*p.z;
    }
};

2.2 面和其它变量的定义

将最小凸多面体的面定义为face类，其中a，b，c表示一个面上三个点的编号，ok用来判断该面是否为凸多面体上的点。

除此之外，对所需要用到的一些重要变量进行全局定义，方便各种函数的调用，其中p[M]为给定的三维点，pm为所需判断的点，f[M*8]为凸包的三角形面。

struct face
    {
   
        int a,b,c;  //凸包一个面上的三个点的编号
        int ok;     //该面是否是最终凸包中的面
    };
node p[M];//初始点
node pj[M];
face f[M*8];//凸包三角形
vector<int> v;  

2.3 三维凸包的构建函数

定义一个创建三维凸包的类，在类中定义所需用到的一些变量和函数，并在最后定义变量“hull”。

struct Convex_hull//三维凸包
{
   
	int n;		//初始点数
	int m;		//判断点数
	int cnt;		//凸包三角形数      
	int to[M][M];	//点i到j是属于哪个面
	double len(node p);		//向量的长度、模	
	double area(node a, node b, node c);	//三个点的面积*2	
	double volume(node a, node b, node c, node d);	//四面体体积*6，即混合积	
	double ptof(node q, face f);		//若为正：点与面同向	
	void deal(int q, int a, int b);
	void dfs(int q, int cur);	//维护凸包，若点q在凸包外则更新凸包	
	int same(int s, int t);		//判断两个三角形是否共面
	void make();		//构建3D凸包
	
}hull;

2.3.1 功能函数的定义

首先对凸包类中常用的计算方法进行函数的编写，主要为向量的长度计算len函数、三个点的面积计算aera函数、四面体的体积计算volume函数。

	//向量的长度、模
    double len(node p)		
    {
   
        return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
    }
	//三个点的面积*2
    double area(node a,node b,node c)	
    {
   
        return len((b-a)*(c-a));
    }
	//四面体体积*6，即混合积
    double volume(node a,node b,node c,node d)
    {
   
        return (b-a)*(c-a)^(d-a);   
    }

定义ptof函数，通过有限体积法用来判断点与面的方向。该函数的输入为所需判断的三维点和面，先得到面上的两个向量m、n，连接三维点和面上任意一点为向量t，通过m与n的叉乘得到面的法向量，再将其与向量t点乘判断两者的方向关系。如果函数返回的值为正，则点与面的法向量同向，由于在存储面时，保证面的法线方向朝向凸包外部，因此可以知道该点在凸包外部，并且此面可被P点看见。

	double ptof(node q,face f)     //若为正：点与面同向
    {
      
        node m=p[f.b]-p[f.a];
        node n=p[f.c]-p[f.a];
        node t=q-p[f.a];
        return m*n^t;
    }  

定义deal和dfs函数，用来更新凸包。当判断得到点q在凸包外并能“看见”面cur时，调用dfs函数，先删除该面，然后调用deal函数，在deal函数中，输入为三维点q和面上的其中两点，因此需调用三次。首先通过面上的两点a、b组成的边找到与当前面共边的另一个面f [fa]，如果这个面为凸包表面，判断点q能否“看到”该面，若可以，则再次调用dfs()函数将该面删除，若不可，将该点q与输入的另外两点a、b组成一个新的三角形作为面。因为每个三角形的的三边是按照逆时针记录的，所以把边反过来后对应的就是与ab边共线的另一个面，因此to[a][b]上的即为与当前面cur共边的另一个面

 	void deal(int q,int a,int b)
    {
   
        int fa=to[a][b];      //与当前面cnt共边的另一个面
        face add;  
        if(f[fa].ok)         //若fa面目前是凸包的表面则继续
        {
   
            if(ptof(p[q],f[fa])>eps)   //若点q能看到fa面继续深搜fa的三条边，更新新的凸包面
                dfs(q,fa);
            else            //当q点可以看到cnt面的同时看不到a，b共边的fa面，则p和a，b点组成一个新的表面三角形
            {
   
                add.a=b;
                add.b=a;
                add.c=q;
                add.ok=1;
                to[b][a]=to[a][q]=to[q][b]=cnt;
                f[cnt++]=add;
            }
        }
    }
	//因为每个三角形的的三边是按照逆时针记录的，所以把边反过来后对应的就是与ab边共线的另一个面
	void dfs(int q,int cur)//维护凸包，若点q在凸包外则更新凸包
    {
   
        f[cur].ok=0;//删除当前面，因为此时它在更大的凸包内部

		//下面把边反过来(先b,后a)，以便在deal()中判断与当前面(cnt)共边(ab)的那个面。
		/*即判断与当头面(cnt)相邻的3个面(它们与当前面的共边是反向的，
			如下图中(1)的法线朝外(即逆时针)的面130和312,它们共边13，但一个方向是13,另一个方向是31)*/
		deal(q,f[cur].b,f[cur].a);
        deal(q,f[cur].c,f[cur].b);
        deal(q,f[cur].a,f[cur].c);
    }

定义same函数，来判断两个三角形是否共面

	//判断两个三角形是否共面
	int same(int s,int t)
    {
   
        node a=p[f[s].a];
        node b=p[f[s].b];
        node c=p[f[s].c];
        if(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)
            return 1;
        return 0;
    }

2.3.2 三维凸包构建函数

1. Part1：三维凸包的构建函数为make函数，按上述思路，首先选择四个不共面的点形成四面体，然后对每个新点进行判断，删除新点可以“看见”的面，进行深度搜索判断。
   初始四面体的构建首先需要找到四个不共面的点，先找两个不同点p[0]、p[1]，然后寻找与它们不共线的第三个点p[2]，再找与前三个点不共面的第四个点p[3]，将这四个点作为初始四面体的四个点。
2. Part2：在找到四个初始点之后，创建初始四面体，在四个点中选取三个点，共创建四个面。如选择p[1]、p[2]、p[3]三个点，则将这三个点的编号1、2、3放入面f[0]中。为方便判断点在凸包内部还是外部，在存储面时，需要保证面的法线方向朝向凸包外部，因此通过ptof函数判断另一个点p[0]与该面的朝向关系，若法线朝内，则需要调换面f[0]中其中两个编号的顺序。同样的方法创建f[0]、f[1]、f[2]、f[3]四个面。
3. Part3：初始四面体创建完成后，对每个点判断是在当前的三维凸包内部还是外部，使用倍增法更新凸包。其中i表示当前点，j表示当前面，通过ptof函数判断点能否“看到”这个面，如果可以，调用dfs函数更新凸包的面，然后对下一个面进行判断；如果该点对凸包上所有面都不可见，则说明该点在凸包内部，开始进行下一个点的判断。当所有点都判断完成后，即可得到这些三维点构成的最小凸多面体。
4. Part4：删除不是凸包上的面，只保存最外层的凸包，将构成凸包的面储存至f中，将面上的三维点编号储存至v中，将v中的点进行排序、去重，最后输出构成凸包的三维点坐标。

    void make()//构建3D凸包
    {
   
		/*Part1*/
        cnt=0;
        if(n<4)
            return;
        int sb=1;              //立flag
        for(int i=1;i<n