複雜多邊形的三角剖分

charlee44發表於2021-06-19

1. 概述

1.1. 多邊形分類

需要首先明確的是多邊形的分類,第一種是最簡單的凸多邊形:

imglink1

凸多邊形的每個內角都是銳角或鈍角,這種多邊形最普通也最常見。如果至少存在一個角是優角(大於180度小於360度),那麼就是凹多邊形了:

imglink2

以上多邊形有一個共同特徵就是由單個環線的邊界組成。如果存在一個外環和多個內環組成多邊形,那麼就是帶洞多變形了:

imglink3

如上圖所示的多邊形是由一個外環和兩個內環組成的,兩個內環造成了外環多邊形的孔洞,也就是帶洞多邊形了。

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檔案,疊加到原來的多邊形中:
imglink4
效果似乎不是很明顯,那麼我將原來的兩個內環範圍塗黑:
imglink5

說明這個演算法可以適配於凹多邊形以及帶洞多邊形的三角網剖分,不得不說CGAL這個庫真的非常強大。可惜就是這個庫太難以使用了,需要一定的計算幾何知識和Cpp高階特性的知識,才能運用自如,值得深入學習。

3. 參考

CGAL 5.2.2 - 2D Triangulation

相關文章