三角剖分算法

目录

一、引言

1、背景

 2、研究意义

二、散乱点集剖分算法

1、Lawson算法

2、Bowyer算法

3、算法的优化 

4、程序实现 

三、三角形剖分算法 

1、实现思想

2、代码实现

---

 

一、引言

1、背景

        本算法来源于美术馆问题，美术馆问题如下：假如有一个凹四边形的美术馆，我们需要在美术馆当中安排尽可能少的警卫，使得他们能观察到美术馆的所有角落，对于一个形状复杂的美术馆（如下图示），至少需要安排多少个警卫呢？这就是著名的美术馆问题（Art Gallery Problem）：给定一个多边形，确定最少需要在多边形内放置多少名警卫，才能让他们的视野覆盖整个多边形。

 

 

        多边形中的视野问题不好处理，但三角形中的视野问题是很简单的：站在其中的任意一个位置都能观察到这个三角形中的每一个点（究其原因，是因为三角形永远是一个凸图形）。我们可以作出多边形的若干条对角线，把整个多边形划分成一个个的三角形。在每一个三角形里面放一个警卫，就能看守整个多边形了。不过，把警卫放在每个三角形的正中间也实在是有些亏。为什么不把警卫放在多边形的对角线上呢？多边形的每一条对角线都是两个三角形的公共边，站在它上面的每一个警卫都能同时看守两个三角形。这样一来，我们便能大大减少所需警卫的数目。其实，把警卫放在边上还是很亏。为什么不把警卫放在多边形的顶点上呢？这样的话，每个警卫都能看遍顶点交汇于此的所有三角形了。

        综上所述，只需要将复杂的美术馆划分为一个个的三角形就可以了，但是问题又来了，如何将多边形划分为三角形？答案是三角剖分算法。

        另外三角剖分也有优劣之分，一般来说，剖分出来的三角形越匀称越好，匀称的三角形在图形图像处理方面效果更好，如下图所示，都是对同一个六边形进行三角剖分，但是右边的效果更好。

        因此，为了得到更好的三角剖分，提出了Delaunay三角网，Delaunay三角网是一系列相连的但不重叠的三角形的集合, 而且这些三角形的外接圆不包含这个面域的其他任何点，其具有下面两个性质：

1. 空接圆性质：每个德洛内(Delaunay) 三角形的外接圆不包含面内的其他任何点
2. 每两个相邻的三角形构成的凸四边形的对角线，在相互交换后，六个内角的最小角不再增大。

Delaunay三角剖分可以保证离散数据点集的三角剖分结果的唯一性，它是离散点集三角剖分中一种优越的三角剖分方法，由于Delaunay三角网是很强健的一种三角网格，所以Delaunay三角剖分也是实际运用中得到广泛应用的一种三角剖分方法。

Delaunay三角的性质保证了Delaunay三角网是一种比较优的三角剖分，但是如何实现剖分出来的三角形是Delaunay三角网呢？本文便介绍实现Delaunay三角剖分的算法。

 2、研究意义

        数据可视化技术伴随着计算机技术的发展日益走向成熟。现代的数据可视化技术指的是运用计算机图形学和图像处理等技术，将数据换为图形或图像在屏幕上显示出来，并进行交互式处理的相关理论、方法和技术。它涉及到的领域广泛，如计算机辅助设计、计算机视觉、人机交互技术、计算机图形学、图像处理等多个领域。基于三维空间的数据可视化方法在计算机图形技术发展的影响下也越来越普及，现有的许多二维可视化技术也开始逐步推向三维可视化。

        Delaunay 三角剖分是实现数据可视化的一种行之有效的方法和工具，三维空间Delaunay三角剖分理论在三角剖分技术由二维逐渐向三维甚至更高维推进的影响下，也得到了进一步的丰富和完善，并在医学可视化，GIS，计算机图形图像处理、三维有限元方法的预处理等多个领域有着重要的意义和广泛的应用。

二、散乱点集剖分算法

        实现对散乱点集的三角剖分主要是通过逐点插入法来实现的，其主要算法有Lawson算法和Bowyer算法，接下来主要讲解Bowyer算法。

1、Lawson算法

        最早的逐点插入算法是Lawson提出来的，它的时间复杂度为N4/3，该算法的大致思想是:首先定义一个包含所有数据点集的初始多边形，然后将数据点逐个插入到多边形中，最后所有的数据点都被三角形覆盖形成一个三角网络。逐点插入算法的关键在于进行每一次插入点的时候，对新形成的三角形要进行空圆检测，同时用LOP进行局部优化处理。

        该算法基于散点的构网算法理论严密、唯一性好，网格满足空圆特性，是一种比较理想的算法。其逐点插入的构网过程中如果遇到非Delaunay边时，可以通过删除调整，进而构造形成新的Delaunay边。在完成构网后，增加新点时，无需对所有的点进行重新构网，只需对新点的影响三角形范围进行局部构网。但是在实际应用当中，点集数量较大，分布复杂时构网速度较慢，而且对于区域为非凸，存在内环时，则会产生非法三角形。

2、Bowyer算法

        Bowyer逐点插入法是一个很经典的实现方法，网上的资料大多数都是采用的这种算法。

        其基本思想是：先由给定的点集生成一初始网格，再根据 Delaunay 剖分原理，逐次向网格内加点 ,并重新连接生成新网格，直到所有点添加完毕。首先是假定已经生产了连接若干个顶点的Delaunay三角网格，如果需要加入一个新的节点，找出所有外接圆包含新加入节点的三角形，并将这些三角形删除，形成一个空腔加入一个新的节点，找出所有外接圆包含新加入节点的三角形，并将这些三角形删除，形成一个空腔，空腔的节点与新加入的节点连接，形成新的 Delaunay 三角形网格，不断循环直到遍历完所有的点。

        在其论文中使用伪代码展示了其算法的实现过程，如下：

subroutine triangulate
input : vertex list 　　输入：顶点链表
output : triangle list　　输出：三角形链表
initialize the triangle list　　初始化三角形链表
determine the supertriangle　　确定超级三角形
add supertriangle vertices to the end of the vertex list　　 将超级三角形的顶点加入到顶点链表中
add the supertriangle to the triangle list　　将超级三角形加入到三角形链表中
for each sample point in the vertex list　　对顶点链表中的每一个点
initialize the edge buffer　　初始化边数组
for each triangle currently in the triangle list　　对于三角形链表中的每一个三角形
calculate the triangle circumcircle center and radius　　计算出外接圆圆心和半径
if the point lies in the triangle circumcircle then　　如果这个点在三角形的外接圆内
add the three triangle edges to the edge buffer　　把这个三角形的三条边加入到边数组中
remove the triangle from the triangle list　　从三角形链表中将这个三角形删除
endif
endfor
delete all doubly specified edges from the edge buffer　　把边数组中所有重复的边都删除，注意这里是把所有的重复边都删除，比如有边e1,e2,e2,e3，删除后应该只剩e1和e3
this leaves the edges of the enclosing polygon only　　这步会得到一个闭合的多边形
add to the triangle list all triangles formed between the point 用这个点和边数组中的每一条边都组合成一个三角形，加入到三角形链表中
and the edges of the enclosing polygon
endfor
remove any triangles from the triangle list that use the supertriangle vertices从三角形链表中删除使用超级三角形顶点的三角形
remove the supertriangle vertices from the vertex list将超级三角形的顶点从顶点链表中删除
end

        接下来使用一个实例解释此算法。

以四个点A、B、C、D举例，首先建立一个超级三角形PQR，这个三角形要把所有点都包含进去，如下图所示：

接着分析A点，因为A点在三角形PQR的外接圆内部，所以利用A点将三角形PQR分拆成三个子三角形。

再考虑Ｂ点，它只在三角形AQR的内部，再将三角形AQR分拆成三个子三角形。

接着是C点，这时我们有5个三角形，对这5个三角形每一个检查C点在不在它的外接圆内。经过检测，发现它在三角形APR和三角形ABR的外接圆内。

删除三角形APR和三角形ABR的公共边AR，得到四边形ABPR，并将C点与每一个边组成一个三角形。

最后对D点进行分析，它在三角形ABC和三角形BCR的外接圆内，所以应删除公共边BC，再用D点与这两个三角形的其他边形成子三角形。

最后把含有超级三角形的顶点的三角形全部删除，就得到这四个点的三角剖分

3、算法的优化 

        通过上面的讲解我们发现Bowyer算法每增加一个点，已经产生的左右三角形均需要计算三角形的圆形和半径，随着点集中点数的增加，三角形的增加是非常快的，这样会造成时间复杂度的极度增大，所以我们通过查阅资料对算法进行了优化。

        优化的思想是对所有点的坐标进行排序，这样的话判断三角形和插入点的位置时有三种情况，一种是该点在三角形外接圆的外部，则可以判断该三角形为Delaunay三角形，之后加入新的点不需要再次判断，后面介绍为什么不需要判断；第二种情况是该点在外接圆的外部但不是外接圆的右侧，这时候跳过，该三角形不确定是否为Delaunay三角形；第三种情况是该点在三角形外接圆内部，则需要删除此三角形，并生成新的三个三角形。

        接下来解释一下为什么在右侧之后不需要进行判断，由于已经对点的x坐标进行了排序，所以处理的时候是从左到右进行处理的，遍历下一个点时，该点在外接圆的右侧，则表示以后所有的点都在该外接圆的右侧，则保证了Delaunay三角形的空圆特性，如下图所示：

 

 

如果点1在外侧但不是右侧，后面的点不确定是否在外接圆内，点2可能出现下图所示的特殊情况：

 

 

算法的伪代码如下所示：

input: 顶点列表(vertices)　　　　　　　　　　　　　　　　　　　　  　//vertices为外部生成的随机或乱序顶点列表

output:已确定的三角形列表(triangles)

　　　　初始化顶点列表

　　　　创建索引列表(indices = new Array(vertices.length))　　　　//indices数组中的值为0,1,2,3,......,vertices.length-1

　　　　基于vertices中的顶点x坐标对indices进行sort　　  　　　　　  //sort后的indices值顺序为顶点坐标x从小到大排序（也可对y坐标，本例中针对x坐标）

　　　　确定超级三角形

　　　　将超级三角形保存至未确定三角形列表（temp triangles）

　　　　将超级三角形push到triangles列表

　　　　遍历基于indices顺序的vertices中每一个点　　　　　　　　　  　//基于indices后，则顶点则是由x从小到大出现

　　　　　　初始化边缓存数组（edge buffer）

　　　　　　遍历temp triangles中的每一个三角形

　　　　　　　　计算该三角形的圆心和半径

　　　　　　　　如果该点在外接圆的右侧

　　　　　　　　　　则该三角形为Delaunay三角形，保存到triangles

　　　　　　　　　　并在temp里去除掉

　　　　　　　　　　跳过

　　　　　　　　如果该点在外接圆外（即也不是外接圆右侧）

　　　　　　　　　　则该三角形为不确定        　　　　　　　　　     //后面会在问题中讨论

　　　　　　　　　　跳过

　　　　　　　　如果该点在外接圆内

　　　　　　　　　　则该三角形不为Delaunay三角形

　　　　　　　　　　将三边保存至edge buffer

　　　　　　　　　　在temp中去除掉该三角形

　　　　　　对edge buffer进行去重

　　　　　　将edge buffer中的边与当前的点进行组合成若干三角形并保存至temp triangles中

　　　　将triangles与temp triangles进行合并

　　　　除去与超级三角形有关的三角形

end

接下来通过一个实例解释优化后的算法步骤：

如下图示，随机的三个点，仍然是使用超级三角形将其包括在内：

 

 

接下来就像是伪代码中描述的那样，对temp triangle中的的三角形遍历画外接圆，这时先对左边的第一个点进行判断，其在圆内，所以该三角形不为Delaunay三角形，将其三边保存至edge buffer中，temptriangle中删除该三角形。

 

 

将该点与edge buffer中的每一个边相连，组成三个三角形，加入到temp triangles中, 再将重复对temp triangles的遍历并画外接圆，这时使用的是第二个点来进行判断，有如下三种情况：

1. 该点在三角形1外接圆右侧，则表示左侧三角形为Delaunay三角形，将该三角形保存至triangles中
2. 该点在三角形2外接圆外侧，为不确定三角形，所以跳过（后面会讲到为什么要跳过该三角形），但并不在temp triangles中删除
3. 该点在三角形3外接圆内侧，则这时向清空后的edge buffer加入该三角形的三条边，并用该点写edge buffer中的三角边进行组合，组合成了三个三角形并加入到temp triangles中

一直持续上述的过程，知道所有点均完成了判断，不再一一详解，如下图所示：

4、程序实现 

 

        使用MFC完成了Bowyer算法，并将其以可视化的形式展示了出来，如下图所示：

 

 

        关键代码如下所示：

template <class T>
class Delaunay
{
public:
    using TriangleType = Triangle<T>;
    using EdgeType = Edge<T>;
    using VertexType = Vector2<T>;

    //Deluanay 三角剖分核心算法  ---  逐点插入法
    const std::vector<TriangleType>& triangulate(std::vector<VertexType>& vertices)
    {
        // 将点云拷贝一份
        _vertices = vertices;

        // 计算超级三角形的一些参数
        T minX = vertices[0].x;
        T minY = vertices[0].y;
        T maxX = minX;
        T maxY = minY;

        //计算超级三角形的上下左右边界
        for (std::size_t i = 0; i < vertices.size(); ++i)
        {
            if (vertices[i].x < minX) minX = vertices[i].x;
            if (vertices[i].y < minY) minY = vertices[i].y;
            if (vertices[i].x > maxX) maxX = vertices[i].x;
            if (vertices[i].y > maxY) maxY = vertices[i].y;
        }

        const T dx = maxX - minX;
        const T dy = maxY - minY;
        const T deltaMax = dx > dy ? dx : dy;
            //std::max(dx, dy);
        const T midx = half(minX + maxX);
        const T midy = half(minY + maxY);

        //尽可能扩大超级三角形范围，可以取比20更大的数字
        const VertexType p1(midx - 20 * deltaMax, midy - deltaMax);
        const VertexType p2(midx, midy + 20 * deltaMax);
        const VertexType p3(midx + 20 * deltaMax, midy - deltaMax);

        //std::cout << "Super triangle " << std::endl << Triangle(p1, p2, p3) << std::endl;

        // Create a list of triangles, and add the supertriangle in it
        //将超级三角形添加至 _triangles
        _triangles.push_back(TriangleType(p1, p2, p3));

        //开始依次遍历每个点
        for (auto p = begin(vertices); p != end(vertices); p++)
        {
            //std::cout << "Traitement du point " << *p << std::endl;
            //std::cout << "_triangles contains " << _triangles.size() << " elements" << std::endl;

            //构造变量用来存储临时新产生的边
            std::vector<EdgeType> polygon;

            //依次遍历 _triangles  中的每个三角形
            for (auto& t : _triangles)
            {
                //std::cout << "Processing " << std::endl << *t << std::endl;

                if (t.circumCircleContains(*p))  //如果包含点 p，那么就要产生新的3条边
                {
                    //std::cout << "Pushing bad triangle " << *t << std::endl;
                    t.isBad = true;  //flag 发生改变，准备接下来在 _triangles  中将其剔除
                    polygon.push_back(t.e1);
                    polygon.push_back(t.e2);
                    polygon.push_back(t.e3);
                }
                else
                {
                    //std::cout << " does not contains " << *p << " in his circum center" << std::endl;
                }
            }
            //更新 _triangles
            _triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [](TriangleType& t) {
                return t.isBad;
                }), end(_triangles));   //std::remove_if只移动不删除，erase将其删除，这里还用到了C++11的 lambda 表达式

            for (auto e1 = begin(polygon); e1 != end(polygon); ++e1)  //这个用来删除重复的边
            {
                for (auto e2 = e1 + 1; e2 != end(polygon); ++e2)
                {
                    if (almost_equal(*e1, *e2))
                    {
                        e1->isBad = true;
                        e2->isBad = true;
                    }
                }
            }

            //更新 polygon
            polygon.erase(std::remove_if(begin(polygon), end(polygon), [](EdgeType& e) {
                return e.isBad;
                }), end(polygon));

            for (const auto e : polygon)
                _triangles.push_back(TriangleType(e.p1, e.p2, *p));

        }

        //删除超级三角形
        _triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [p1, p2, p3](TriangleType& t) {
            return t.containsVertex(p1) || t.containsVertex(p2) || t.containsVertex(p3);
            }), end(_triangles));

        for (const auto t : _triangles)
        {
            _edges.push_back(t.e1);
            _edges.push_back(t.e2);
            _edges.push_back(t.e3);
        }

        return _triangles;
    }

    const std::vector<TriangleType>& getTriangles() const { return _triangles; }
    const std::vector<EdgeType>& getEdges() const { return _edges; }
    const std::vector<VertexType>& getVertices() const { return _vertices; }

private:
    std::vector<TriangleType> _triangles;
    std::vector<EdgeType> _edges;
    std::vector<VertexType> _vertices;
};

三、三角形剖分算法 

1、实现思想

上面的散乱点集三角剖分产生的Dalaunay三角网是一种比较优的三角形，接下来我们又实现了多边形的三角剖分，但是产生的不是比较优的三角剖分，和之前产生较优的Dalaunay三角网进行对比，接下来介绍我们的实现思路。

首先先考虑对于任意的多边形，是凸多边形容易处理还是凹多边形容易处理？显然，凸多边形更容易处理。对于凸多边形直接从一点出发即可全部实现三角剖分，对于凹多边形容易出现线相交的情况，如下图所示：

 

    （凸多边形三角剖分）               （凹多边形三角剖分）

    对于多边形三角剖分一个简单的实现思路是顺序选取输入的三个顶点划分为三角形，但是这种思路对于凸多边形正确，但是对于凹多边形会出现下面所示的情况：

 

 

    所以问题就转换成了如何识别上面的情况。我们的解决方案是规定的输入顺序和方向，首先我们了解两件事，一个是{i,j,k}能不能构成三角形与整个图形除ijk三点其他点的“大体”位置有关；另外一个是{i,j,k} 能否构成三角形隐含i先于j，j先于k。

    规定了输入的方向和顺序之后，情况只有两种：可以构成三角形和不可以构成三角形。接下来只需要识别下面这两种情况即可。

   为了识别上面的两种情况，这里需要一些数学知识，在数学中， 若{i,j,k}是右手系，那么向量(i,j)和向量(i,k)的向量积为正；若{i,j,k}是左手系，那么向量(i,j)和向量(i,k)的向量积为负。所以上述两种情况的判断就可以通过向量的计算来实现。

    但是如果仅仅是实现上述的条件是不够的，还有可能出现下面的这种特殊情况：

 

 

    对于这种特殊情况的解决方案只需要判断三角形ijk中没有其他的点即可。这时就可以将所有情况考虑到进行三角剖分了。

总体实现的步骤如下：当ijk能形成第三边在多边形内的三角形时，分割该三角形，且去除j点，进行下一个三角形的剖分，若ijk形成的三角形第三边在多边形外，则不能进行该三角剖分，此时ijk均向后移动一位，进行下一次剖分。若可以进行三角形的剖分，i回到起始点。

2、代码实现

对于多边形三角剖分算法，也使用MFC完成了程序实现，实现了可视化展示，完成如下：

关键代码如下：

bool CcgSXQFillPolyView::PT(int pNumbers, CPoint MPOlygon[], CDC* pDC)
{
	if (pNumbers <= 4)/// 如果当前多边形的边数小于3那么判定该多边形非法
		return FALSE;

	DWORD dwStart = GetTickCount();/// 获得开始运行时间
	CPen pNewPen;///定义画笔颜色为红色
	pNewPen.CreatePen(PS_SOLID, 1, RGB(255, 0, 0));
	CPen* poldPen = pDC->SelectObject(&pNewPen);
	/// 定义循环变量i,j,k
	int i = 0;
	int j = 1;
	int k = 2;
	while (pNumbers > 4) {/// 当多边形的边数大于3时进行
		bool cond1;/// 条件一，向量ij的和向量ik的夹角必须小于π
		double a1 = MPOlygon[j].x - MPOlygon[i].x;
		double a2 = MPOlygon[j].y - MPOlygon[i].y;
		double b1 = MPOlygon[k].x - MPOlygon[i].x;
		double b2 = MPOlygon[k].y - MPOlygon[i].y;
		cond1 = (a1*b2 <= a2*b1);
		bool cond2 = TRUE;/// 条件二，当前多边形中不能有其他点
		double tempx[4];/// 当前多边形上移半格
		double tempy[4];
		tempx[0] = MPOlygon[i].x;
		tempy[0] = MPOlygon[i].y+0.5;
		tempx[1] = MPOlygon[j].x;
		tempy[1] = MPOlygon[j].y+0.5;
		tempx[2] = MPOlygon[k].x;
		tempy[2] = MPOlygon[k].y+0.5;
		tempx[3] = tempx[0];
		tempy[3] = tempy[0];
		for (int p = 0; p < pNumbers; p++) {
			if (p != i && p != j && p != k) {/// 判断其他点是否在当前三角形内
				if (JudgeIn(MPOlygon[p].x, MPOlygon[p].y, tempx, tempy, 3 + 1)) {
					cond2 = FALSE;
					break;
				}	
			}
		}
		///cond1 = true;
		if (cond1 && cond2) {/// 如果同时满足条件一和条件二
			pDC->MoveTo(MPOlygon[i]);/// 从i点到k点画线
			pDC->LineTo(MPOlygon[k]);
			for (int t = j; t < pNumbers - 1; t++) {/// 更新当前多边形
				MPOlygon[t] = MPOlygon[t + 1];
			}
			pNumbers--;/// 更新多边形的边数
			i = 0;/// 重新赋值i，j，k
			j = 1;
			k = 2;
		}
		else {/// 否则
			i ++;/// i下移一点
			j = i + 1;
			k = j + 1;
		}
	}

	pDC->SelectObject(poldPen);/// 跟新pDC

	DWORD dwStop = GetTickCount();/// 获得结束时间
	DWORD dwInterval = dwStop - dwStart;/// 获得运行时间
	CcgSXQFillPolyDoc* pDoc = (CcgSXQFillPolyDoc*)GetDocument();/// 获得Doc指针
	pDoc->polygonalTriangualtionRuntime = (double)dwInterval;/// 将数据写入Doc类中
	CString str;/// 将数据写入polygonalTriangulationInformation，并将其格式化
	str.Format(_T("%.2lf"), pDoc->polygonalTriangualtionRuntime);
	pDoc->polygonalTriangulationInformation += _T("Runtime Consuming : ") + str + _T("ms\r\n");

	pDoc->UpdateAllViews(this);/// 更新视图
	return TRUE;
}