Delaunay三角网之逐点插入法

秒速五厘米 2022-12-17 12:59 238阅读 0赞

一直想实现一下Delaunary三角网，趁着自己正好闲了下来就实现了一下。

### 文章目录 ###

*  一、简介
 *  二、逐点插入法
 *   *  2.1算法步骤
     *  2.2代码实现
 *  三、小结

# 一、简介 #

> Delaunary三角网作为一种主要的DTM表示方法，其具有非常广泛的用途。那么，一个点集的三角网可以有很多种形式，为什么Delaunary三角网会成为宠儿呢？这主要是由于其具有以下几个性质：  
> 1、任何点集的Delaunary三角网都是唯一的。  
> 2、Delaunary三角网力求最佳的三角形几何形状，每个三角形尽量接近等边三角形。  
> 3、保证最邻近的点构成三角形，即三角形的边长之和最小。  
> 也由于以上这几个性质，使得Delaunary三角网在所有可能的三角网中，对地形拟合方面表现的最为出色，因此其常常会被用到TIN(不规则三角网)的生成。除此之外，Delaunary三角网为相互邻接且互不重叠的三角形集合，每个三角形的外接圆不包含其他的点，这一基本法则也叫做空圆法则。

> 目前， Delaunary三角网主流的生成算法主要有三类：分治算法、逐点插入法以及生长算法，在此基础上衍生出了许多其他三角网生成算法，感兴趣的同学可以查找相关的文献进行学习。这几种主流的三角网生成算法都有其优劣势，如分治算法是一种典型的“以空间换时间”的算法，其时间效率较高，但是其占用的内存空间也会较多；而逐点插入法则正好相反，其时间消耗较差，但是其占用的内存空间会比较少，除此之外呢，逐点插入法的实现过程相对简单，易于被程序实现。

# 二、逐点插入法 #

## 2.1算法步骤 ##

> 1、构建平面点集 p p p的超三角形(也可以是其他超多边形)，并以此作为该点集的凸闭包。  
> 2、将点集 p p p中的一点插入三角网中，在插入的过程中要注意满足上述的Delaunary三角网的空圆法则，如果三角网中存在三角形 T i T\_i Ti的外接圆包含插入点即将其删除。  
> 3、以插入点和三角形 T i T\_i Ti的非公共边为基础构建新的三角网，以此循环迭代直至点集 p p p为空时则算法停止。

## 2.2代码实现 ##

代码可能写的有些粗糙，因为本来就是为了了解这个算法，所以如果有地方出现错误还请见谅\*~\*。  
Point2.hpp

#ifndef POINT2
    #define POINT2
    
    #include <math.h>
    
    namespace DT { 
    
    	//因为并没有用的指针动态分配内存，所以这里使用class与struct没有什么区别
    	template<typename T>
    	class Point2
    	{ 
    	public:
    		T x;
    		T y;
    		Point2() = default;		//使用默认的构造函数
    		Point2(const Point2<T> &point) = default;		//使用默认的拷贝构造函数
    		Point2(const T point_x, const T point_y) //一般构造函数
    			:x(point_x), y(point_y)
    		{ }
    
    		T Distance2(const Point2<T> &point) const { 		//两点之间距离的平方
    			const T dx = this->x - point.x;
    			const T dy = this->y - point.y;
    			return dx * dx + dy *dy;
    		}
    
    		T Distance(const Point2<T> &point)const { 
    			return sqrt(Distance2(point));
    		}
    
    		T Norm2() const { 		//点到原点的距离平方
    			return x*x + y*y;
    		}
    
    		Point2 &operator= (const Point2<T>&) = default;	//使用默认的拷贝构造函数
    		Point2 &operator= (Point2&&) = default;
    		bool operator == (const Point2<T> &point) const { 	//重载==符号，用于判断两个点是否相同
    			return (this->x == point.x) && (this->y == point.y);
    		}
    	};
    
    }
    
    #endif // !POINT2

Edge.hpp

#ifndef EDGE_H
    #define EDGE_H
    
    #include"Point2.hpp"
    
    namespace DT { 
    
    	template<typename T>
    	struct Edge
    	{ 
    		size_t p1;		//边的端点的索引
    		size_t p2;
    		bool isBad = false;
    
    		Edge(const size_t p1, const size_t p2)
    			:p1(p1), p2(p2)
    		{ }
    	};
    }
    
    #endif // !EDGE_H

Triangle.hpp

#ifndef TRIANGLE_H
    #define TRIANGLE_H
    
    #include"Edge.hpp"
    
    namespace DT { 
    
    	template<typename T>
    	struct Triangle { 
    		size_t a;		//三角形端点的索引
    		size_t b;
    		size_t c;
    		bool isBad = false;
    
    		Triangle() = default;
    		Triangle(const Triangle&) = default;
    		Triangle(Triangle&&) = default;
    		Triangle(const size_t a, const size_t b, const size_t c)
    			:a(a), b(b), c(c)
    		{ }
    
    		Triangle &operator=(const Triangle&) = default;
    		Triangle &operator=(Triangle&&) = default;
    	
    	};
    
    }
    
    #endif // !TRIANGLE_H

Delaunay.hpp

#ifndef DELAUNARY_H
    #define DELAUNARY_H
    
    #include"Point2.hpp"
    #include"Edge.hpp"
    #include"Triangle.hpp"
    
    #include<vector>
    #include<algorithm>
    
    namespace DT { 
    
    	template<typename T>
    	class Delaunay
    	{ 
    		using Type = T;
    		using Point2Type = Point2<Type>;
    		using EdgeType = Edge<Type>;
    		using TriangleType = Triangle<Type>;
    
    		std::vector<TriangleType> _triangles;		//三角形面片集合
    		std::vector<EdgeType> _edges;		//边集合
    		std::vector<Point2Type> _points;	//点集合
    
    		//判断点、线、面相同的重载函数
    		bool almost_equal(const Point2<T> &p1, const Point2<T> &p2) const		//点
    		{ 
    			return std::abs(p1.x - p2.x) < 0.0001f&&
    				std::abs(p1.y - p2.y) < 0.00001f;
    		}
    
    		bool almost_equal(const Edge<T> &e1, const Edge<T> &e2)const		//线
    		{ 
    			return	(almost_equal(_points[e1.p1], _points[e2.p1]) && almost_equal(_points[e1.p2], _points[e2.p2])) ||
    				(almost_equal(_points[e1.p1], _points[e2.p2]) && almost_equal(_points[e1.p2], _points[e2.p1]));
    		}
    
    		bool almost_equal(const Triangle<T> &t1, const Triangle<T> &t2)const		//面
    		{ 
    			return	(almost_equal(_points[t1.a], _points[t2.a]) || almost_equal(_points[t1.a], _points[t2.b]) || almost_equal(_points[t1.a], _points[t2.c])) &&
    				(almost_equal(_points[t1.b], _points[t2.a]) || almost_equal(_points[t1.b], _points[t2.b]) || almost_equal(_points[t1.b], _points[t2.c])) &&
    				(almost_equal(_points[t1.c], _points[t2.a]) || almost_equal(_points[t1.c], _points[t2.b]) || almost_equal(_points[t1.c], _points[t2.c]));
    		}
    
    		bool ContainVertex(const Triangle<T>& t, const Point2<T> &p) { 		//判断三角形中是否包含这个点
    			return almost_equal(_points[t.a], p) || almost_equal(_points[t.b], p) || almost_equal(_points[t.c], p);
    		}
    
    		bool CircumCircleContains(const Triangle<T>& t,const Point2<T> &p) { 		//判断一个点是否在三角形内
    			const T aa = _points[t.a].Norm2();		//点p1的平方和
    			const T bb = _points[t.b].Norm2();
    			const T cc = _points[t.c].Norm2();
    
    			const T ax = _points[t.a].x;		//a点的x坐标
    			const T ay = _points[t.a].y;
    			const T bx = _points[t.b].x;
    			const T by = _points[t.b].y;
    			const T cx = _points[t.c].x;
    			const T cy = _points[t.c].y;
    
    			//求取圆心坐标
    			const T circum_x = (aa * (cy - by) + bb * (ay - cy) + cc * (by - ay)) / (ax * (cy - by) + bx * (ay - cy) + cx * (by - ay));
    			const T circum_y = (aa * (cx - bx) + bb * (ax - cx) + cc * (bx - ax)) / (ay * (cx - bx) + by * (ax - cx) + cy * (bx - ax));
    
    			const Point2<T> circumCenter(circum_x / 2, circum_y / 2);		//圆心
    			const T circumRadius = _points[t.a].Distance2(circumCenter);		//a点到圆心的距离
    			const T dist = p.Distance2(circumCenter);	//p点到圆心的距离
    
    			return dist <= circumRadius;
    		}
    
    		
    
    	public:
    		Delaunay() = default;
    		Delaunay(const Delaunay &) = delete;		//不使用默认的拷贝构造函数
    		Delaunay(Delaunay&&) = delete;
    
    		const std::vector<TriangleType>& triangulate(const std::vector<Point2Type>& points) { 		//对给出的点集进行三角化
    			//将传进来的点集保存至_points中
    			_points = points;
    
    			//构建一个超三角
    			T minX = _points[0].x;
    			T minY = _points[0].y;
    			T maxX = minX;
    			T maxY = minY;
    
    			for (std::size_t i = 0; i < _points.size(); ++i)		//挑选出x,y中的最大最小值
    			{ 
    				if (_points[i].x < minX) minX = _points[i].x;
    				if (_points[i].y < minY) minY = _points[i].y;
    				if (_points[i].x > maxX) maxX = _points[i].x;
    				if (_points[i].y > maxY) maxY = _points[i].y;
    			}
    
    			const T dx = maxX - minX;
    			const T dy = maxY - minY;
    			const T deltaMax = std::max(dx, dy);
    			const T midx = (minX + maxX) / 2;
    			const T midy = (minY + maxY) / 2;
    
    			const Point2Type p1(midx - 10 * deltaMax, midy - deltaMax);		//超三角的三个顶点
    			const Point2Type p2(midx, midy + 10 * deltaMax);
    			const Point2Type p3(midx + 10 * deltaMax, midy - deltaMax);
    
    			_points.push_back(p1);		//将超三角形的顶点添加到_points中，但是有感觉这种方法笨笨的，不过一时之间又想不到其他的方法
    			_points.push_back(p2);
    			_points.push_back(p3);
    			_triangles.push_back(TriangleType(_points.size()-3, _points.size() - 2, _points.size() - 1));		//将超三角添加至面集合中
    
    																//遍历点集中所有的点
    			for (size_t i=0;i<_points.size()-3;i++)
    			{ 
    				std::vector<EdgeType> polygon;		//存储包含插入点的三角形的边
    
    				//遍历三角形列表，判断三角形是否包含插入点
    				for (auto & t : _triangles)
    				{ 
    					if (CircumCircleContains(t,_points[i]))
    					{ 
    						t.isBad = true;
    
    						polygon.push_back(Edge<T>{ t.a, t.b});
    						polygon.push_back(Edge<T>{ t.b, t.c});
    						polygon.push_back(Edge<T>{ t.c, t.a});
    					}
    				}
    
    				//如果三角形中含有插入点，则将其从三角形列表中删除
    				_triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [](TriangleType &t) { 
    					return t.isBad;
    				}), end(_triangles));
    
    				//遍历边表，检查是否有相同的边
    				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.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, i));
    			}
    
    			//最后删除所有包含超三角三个顶点的三角形,类中的每一个函数都有一个this指针
    			_triangles.erase(std::remove_if(begin(_triangles), end(_triangles), [p1, p2, p3, this](TriangleType &t) { 
    				return ContainVertex(t,p1) || ContainVertex(t,p2) || ContainVertex(t,p3);
    			}), end(_triangles));
    
    
    			//添加三角形的边至_edges边表中
    			for (const auto t : _triangles)
    			{ 
    				_edges.push_back(Edge<T>{ t.a, t.b});		//使用初始化列表来实现
    				_edges.push_back(Edge<T>{ t.b, t.c});
    				_edges.push_back(Edge<T>{ t.c, t.a});
    			}
    
    			//将超三角形中的点弹出
    			_points.pop_back();
    			_points.pop_back();
    			_points.pop_back();
    
    			return _triangles;
    		}
    
    		const std::vector<TriangleType>& getTriangles() const { 
    			return _triangles;
    		}
    		const std::vector<EdgeType>& getEdges() const { 
    			return _edges;
    		}
    		const std::vector<Point2Type>& getVertices() const { 
    			return _points;
    		}
    
    		Delaunay& operator=(const Delaunay&) = delete;
    		Delaunay& operator=(Delaunay&&) = delete;
    	};
    
    }
    
    #endif // !DELAUNARY_H

main.cpp

//基础头文件
    #include <iostream>
    #include <vector>
    #include <iterator>
    #include <algorithm>
    #include <array>
    #include <random>
    #include <chrono>
    
    //Delaunay相关头文件
    #include"Delaunay.hpp"
    
    using namespace DT;
    
    int main() { 
    	int numberPoints = 3;
    
    	//随机产生点,按照自己熟悉的随机函数产生点
    	/*std::default_random_engine eng(std::random_device{}()); std::uniform_real_distribution<double> dist_w(0, 800); std::uniform_real_distribution<double> dist_h(0, 600); std::cout << "Generating " << numberPoints << " random points" << std::endl;*/
    
    	std::vector<Point2<double>> points;
    	//for (int i = 0; i < numberPoints; ++i) { 
    	// points.push_back(Point2<double>{dist_w(eng), dist_h(eng)}); //随机产生点，并添加至points中
    	//}
    
    	points.push_back(Point2<double>{ 0, 0});		//点号0,1,2,3
    	points.push_back(Point2<double>{ 0, 1});
    	points.push_back(Point2<double>{ 1, 1});
    	points.push_back(Point2<double>{ 1, 0});
    
    	//构建Delaunay三角网
    	Delaunay<double> triangulation;
    	const auto start = std::chrono::high_resolution_clock::now();		//计时器
    	const std::vector<Triangle<double>> triangles = triangulation.triangulate(points);
    	const auto end = std::chrono::high_resolution_clock::now();
    	const std::chrono::duration<double> diff = end - start;
    	std::cout << triangles.size() << " triangles generated in " << diff.count()
    		<< "s\n";
    	const std::vector<Edge<double>> edges = triangulation.getEdges();
    
    	int num = 0;
    	for (size_t i = 0; i <edges.size(); i++)
    	{ 
    		if (i%3==0)
    		{ 
    			std::cout << "\n"
    				<<"第"<< ++num<<"个三角形：\n";
    		}
    		std::cout << edges[i].p1<< "\t" << edges[i].p2<< std::endl;
    	}
    
    	std::cin.get();
    	return 0;
    }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70_pic_center]

# 三、小结 #

> 逐点插入法作为一种经典的凸闭包收缩算法，其思想是：首先找到包含数据区域的最小凸包边形，并从该多边形开始从外到内形成Delaunary三角网。因此其每次插入一个新的点就会删除相应的三角形来构建性的三角网，这个过程中常常伴随着大量的查询计算过程，这也导致了其在面对大量数据的情况下会无能为力。但是其思想非常的简单巧妙，是比较适合通过程序进行实现的，不过由于自己对C++可视化工具学习的还没有到家，所以就没有可视化三角网。等有机会学的差不多了在补上吧，如果上面的代码出现错误，还望大家能及时指正\*~\*。

> 参考资料：《地理信息系统原理与方法》

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70_pic_center]: /images/20221122/dd56b31ee4da4b5b9b272e5c1861a2ed.png