1. 概述

1.1. 多边形分类

需要首先明确的是多边形的分类，第一种是最简单的凸多边形：

凸多边形的每个内角都是锐角或钝角，这种多边形最普通也最常见。如果至少存在一个角是优角(大于180度小于360度)，那么就是凹多边形了：

以上多边形有一个共同特征就是由单个环线的边界组成。如果存在一个外环和多个内环组成多边形，那么就是带洞多变形了：

如上图所示的多边形是由一个外环和两个内环组成的，两个内环造成了外环多边形的孔洞，也就是带洞多边形了。

1.2. 三角剖分

三角剖分也叫做三角化，或者分格化（tessellation/triangulation），将复杂的多边形剖分成多个三角形。这在图形学上有非常多的好处，便于绘制和计算。这类算法往往与Delaunay三角网算法相关，多边形的边界作为Delaunay三角网的边界约束，从而得到比较好的三角网。

2. 详论

我曾经在《通过CGAL将一个多边形剖分成Delaunay三角网》这篇文章中，通过CGAL实现了一个多边形的剖分，不过这个文章介绍的算法内容不支持凹多边形和带洞多边形（这也是很多其他算法实现的问题）。所以我继续翻了CGAL的官方文档，在《2D Triangulation》这一章中确实介绍了带洞多边形的三角剖分的案例。由于带洞多边形最复杂，那么我通过这个案例，来实现一下带洞多边形的三角剖分。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Polygon_2.h>
#include <iostream>#include <gdal_priv.h>
#include <ogrsf_frmts.h>struct FaceInfo2
{FaceInfo2() {}int nesting_level;bool in_domain() {return nesting_level % 2 == 1;}
};typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>    Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb>        Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>               TDS;
typedef CGAL::Exact_predicates_tag                                Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag>  CDT;
typedef CDT::Point                                                Point;
typedef CGAL::Polygon_2<K>                                        Polygon_2;
typedef CDT::Face_handle                                          Face_handle;void
mark_domains(CDT& ct,Face_handle start,int index,std::list<CDT::Edge>& border)
{if (start->info().nesting_level != -1) {return;}std::list<Face_handle> queue;queue.push_back(start);while (!queue.empty()) {Face_handle fh = queue.front();queue.pop_front();if (fh->info().nesting_level == -1) {fh->info().nesting_level = index;for (int i = 0; i < 3; i++) {CDT::Edge e(fh, i);Face_handle n = fh->neighbor(i);if (n->info().nesting_level == -1) {if (ct.is_constrained(e)) border.push_back(e);else queue.push_back(n);}}}}
}
//explore set of facets connected with non constrained edges,
//and attribute to each such set a nesting level.
//We start from facets incident to the infinite vertex, with a nesting
//level of 0. Then we recursively consider the non-explored facets incident
//to constrained edges bounding the former set and increase the nesting level by 1.
//Facets in the domain are those with an odd nesting level.
void
mark_domains(CDT& cdt)
{for (CDT::Face_handle f : cdt.all_face_handles()) {f->info().nesting_level = -1;}std::list<CDT::Edge> border;mark_domains(cdt, cdt.infinite_face(), 0, border);while (!border.empty()) {CDT::Edge e = border.front();border.pop_front();Face_handle n = e.first->neighbor(e.second);if (n->info().nesting_level == -1) {mark_domains(cdt, n, e.first->info().nesting_level + 1, border);}}
}int main()
{//创建三个不相交的嵌套多边形Polygon_2 polygon1;polygon1.push_back(Point(-0.558868038740926, -0.38960351089588));polygon1.push_back(Point(2.77833686440678, 5.37465950363197));polygon1.push_back(Point(6.97052814769976, 8.07751967312349));polygon1.push_back(Point(13.9207400121065, 5.65046156174335));polygon1.push_back(Point(15.5755523607748,-1.98925544794189));polygon1.push_back(Point(6.36376361985472, -6.18144673123487));  Polygon_2 polygon2;polygon2.push_back(Point(2.17935556413387, 1.4555590039808));polygon2.push_back(Point(3.75630057749723, 4.02942327866582));polygon2.push_back(Point(5.58700685737883, 4.71820385921534));polygon2.push_back(Point(6.54767450919789, 1.76369768475295));polygon2.push_back(Point(5.71388749063795, -0.900795613688593));polygon2.push_back(Point(3.21252643495814, -0.320769861646896));Polygon_2 polygon3;polygon3.push_back(Point(7.74397762278389, 0.821155837685192));polygon3.push_back(Point(9.13966458863422, 4.24693293568146));polygon3.push_back(Point(10.1909612642098, 1.83620090375816));polygon3.push_back(Point(12.1485481773505, 4.84508449247446));polygon3.push_back(Point(11.4416417920497, -2.29648257953892));polygon3.push_back(Point(10.1547096547072, 0.712401009177374));//将多边形插入受约束的三角剖分CDT cdt;cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);cdt.insert_constraint(polygon2.vertices_begin(), polygon2.vertices_end(), true);cdt.insert_constraint(polygon3.vertices_begin(), polygon3.vertices_end(), true);//标记由多边形界定的域内的面mark_domains(cdt);//遍历所有的面int count = 0;for (Face_handle f : cdt.finite_face_handles()){if (f->info().in_domain()) ++count;}std::cout << "There are " << count << " facets in the domain." << std::endl;//将结果输出成shp文件，方便查看{GDALAllRegister();GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");if (!driver){printf("Get Driver ESRI Shapefile Error！\n");return false;}const char *filePath = "D:/test.shp";GDALDataset* dataset = driver->Create(filePath, 0, 0, 0, GDT_Unknown, NULL);OGRLayer* poLayer = dataset->CreateLayer("test", NULL, wkbPolygon, NULL);//创建面要素for (Face_handle f : cdt.finite_face_handles()){if (f->info().in_domain()){OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());OGRLinearRing ogrring;for (int i = 0; i < 3; i++){            ogrring.setPoint(i, f->vertex(i)->point().x(), f->vertex(i)->point().y());}ogrring.closeRings();OGRPolygon polygon;polygon.addRing(&ogrring);poFeature->SetGeometry(&polygon);if (poLayer->CreateFeature(poFeature) != OGRERR_NONE){printf("Failed to create feature in shapefile.\n");return false;} }}//释放GDALClose(dataset);dataset = nullptr;}return 0;
}

在代码的最后，我将生成的三角网输出成shp文件，叠加到原来的多边形中：

效果似乎不是很明显，那么我将原来的两个内环范围涂黑：

说明这个算法可以适配于凹多边形以及带洞多边形的三角网剖分，不得不说CGAL这个库真的非常强大。可惜就是这个库太难以使用了，需要一定的计算几何知识和Cpp高级特性的知识，才能运用自如，值得深入学习。

3. 参考

CGAL 5.2.2 - 2D Triangulation

复杂多边形的三角剖分相关推荐

leetcode1039. 多边形三角剖分的最低得分（动态规划）
给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], -, A[N-1]. 假设您将多边形剖分为 N-2 个三角形.对于每个三角形,该三角形的值是顶点标记的乘积,三 ...
LeetCode 1039. 多边形三角剖分的最低得分（区间DP）
文章目录 1. 题目 2. 解题 1. 题目给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], ..., A[N-1]. 假设您将多边形剖分为 N-2 个三角形 ...
leetcode - 1039. 多边形三角剖分的最低得分
给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], -, A[N-1]. 假设您将多边形剖分为 N-2 个三角形.对于每个三角形,该三角形的值是顶点标记的乘积,三 ...
5047. 多边形三角剖分的最低得分
5047. 多边形三角剖分的最低得分转载请注明出处附上我的博客链接四元君题目难度 Medium 题目描述给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i] ...
【1039】多边形三角剖分的最低得分
给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], -, A[N-1].假设您将多边形剖分为N-2 个三角形.对于每个三角形,该三角形的值是顶点标记的乘积,三角剖 ...
【Leetcode 1039】多边形三角剖分的最低得分
问题描述给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], ..., A[N-1].假设您将多边形剖分为 N-2 个三角形.对于每个三角形,该三角形的值是顶点标 ...
动态规划--凸多边形最优三角剖分
算法设计与实现王晓东题目描述: 用多边形顶点的逆时针序列表示凸多边形,即P={v0,v1,-,vn-1}表示具有n条边的凸多边形. 给定凸多边形P,以及定义在由多边形的边和弦组成的三角形上的权函 ...
leetcode-多边形三角剖分的最低得分
给定 N,想象一个凸 N 边多边形,其顶点按顺时针顺序依次标记为 A[0], A[i], -, A[N-1]. 假设您将多边形剖分为 N-2 个三角形.对于每个三角形,该三角形的值是顶点标记的乘积,三 ...
动态规划---例题5.凸多边形最优三角剖分问题
一.题目描述通常,用多边形顶点的序列来表示一个凸多边形,即P=<v0 ,v1 ,- ,vn-1>表示具有n条边v0v1,v1v2,- ,vn-1vn的一个凸多边形,其中,约定v0 = v ...

复杂多边形的三角剖分

文章目录