複雜多邊形的三角剖分
阿新 • • 發佈：2021-06-19
目錄
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高階特性的知識，才能運用自如，值得深入學習。