基于Lawson 插入法的Delaunay三角网格构造算法C++实现

sankingvenice

已于 2023-04-12 20:04:26 修改

阅读量707

点赞数 2

分类专栏：点云文章标签：算法 c++

于 2023-04-04 22:16:35 首次发布

本文链接：https://blog.csdn.net/baidu_34931359/article/details/129962825

版权

点云专栏收录该内容

7 篇文章 4 订阅

订阅专栏

文章目录

一、原理与伪代码
二、C++函数实现
实现效果
参考资料：

一、原理与伪代码

本文的实现参考计算几何经典教材《Computational Geometry-Third Edition》邓俊辉译。
算法伪代码如下图所示：
在这里插入图片描述

首先构造一个足够大的三角形p-1p-2p-3，将整个点集P包围起来；

随后按随机的次序逐一插入点集P各点Pr，在插入的同时更新当前点集的三角剖分。考虑引入点Pr时的情况。首先，我们要在当前的三角剖分中，确定pr落在哪个三角形之内（具体方法，我们也将在后面介绍）。然后，我们将pr与该三角形的三个顶点分别联接起来，生成三条边。倘若Pr碰巧落在三角剖分的某条边PiPj上，我们就需要找到与PiPj关联的那两个三角形，然后将pr与对顶的那两个顶点Pk分别联接起来，生成新的边。
在引入点Pr 之后，我们需要针对每一条可能的非法边，调用一次子函数LEGALIZEEDGE进行Delaunay空圆检测。若该边不合法，则通过Flip函数实现边翻转操作，同时对翻转后新生成的边也调用子函数LEGALIZEEDGE进行空圆检测，对于所有可能的新非法边，LEGALIZEEDGE都要递归地调用自己，逐一进行核查，如此持续直到将所有的非法边转换为合法边。LEGALIZEEDGE的伪代码如下：
在这里插入图片描述
一个例子如下：

二、C++函数实现

1.构建超级三角形

//构建超级三角形
Triangle  Delaunay2D::CreateSuperTriangle(std::vector<Point>& points)
{
	Triangle triangle;
	if (points.size() < 3)
		return triangle;

	Point v0, v1, v2;
	CalcSuperTriangle(points, v0, v1, v2);

	triangle = Triangle (v0, v1, v2);
	return triangle;
}

// 计算距离最远的三个点，左下角、右下角和上方的中心位置，作为超级三角形的顶点
void Delaunay2D::CalcSuperTriangle(const std::vector<Point>& points, Point& v0, Point& v1, Point& v2) 
{
	double xmin, ymin, xmax, ymax;
	GetBoundingBox(points, xmin, ymin, xmax, ymax);

	double dx = xmax - xmin;
	double dy = ymax - ymin;
	double dmax = std::max(dx, dy);
	double xmid = xmin + dx * 0.5;
	double ymid = ymin + dy * 0.5;

	v0 = Point(xmid - 2 * dmax, ymid - dmax);
	v1 = Point(xmid, ymid + 2 * dmax);
	v2 = Point(xmid + 2 * dmax, ymid - dmax);

}

// 计算点集中x、y坐标的最小值和最大值
void Delaunay2D::GetBoundingBox(const std::vector<Point>& points, double& xmin, double& ymin, double& xmax, double& ymax)
{
	if (0==points.size())
	{
		return;
	}
	xmin = xmax = points[0].x;
	ymin = ymax = points[0].y;
	for (int i = 1; i < points.size(); i++) 
	{
		if (points[i].x < xmin) xmin = points[i].x;
		if (points[i].y < ymin) ymin = points[i].y;
		if (points[i].x > xmax) xmax = points[i].x;
		if (points[i].y > ymax) ymax = points[i].y;
	}
}

2.逐点插入

vector<Triangle> Delaunay2D::GetDelaunay(vector<Point> &points)
{
Tins.push_back(superTriangle);
vector<Triangle> delTriangle; //要删除的三角形
vector<Triangle> newTriangle;//新生成的三角形，之后用于检测空圆
int numberRelatedTriangle = 0; //与该点相关的三角形数量
//离散点内插过程
for (int i = 0; i < points.size(); i++)
{
	newTriangle.clear();
	delTriangle.clear();
	if (i>2)
	{
		int a = 1;
	}
	for (int j = 0; j < Tins.size(); j++)
	{
		if (Tins[j].IsPointOnTriangle(points[i]) == true)
		{
			delTriangle.push_back(Tins[j]);
			numberRelatedTriangle =2;
			continue;
		}
		if (Tins[j].IsPointInTriangle(points[i]) == true)
		{
			delTriangle.push_back(Tins[j]);
			numberRelatedTriangle = 1;
		}
	}

	vector<Edge> borderLines;//离散点的边界边
	if (0 == numberRelatedTriangle) //点不在所有的三角形中  构建的三角形存在异常
	{
		break;
	}
	else if (1== numberRelatedTriangle) //在三角形内部，删除所在三角形。连接检索点P与所在三角形的三个边，生成三个新三角形。
	{
		Tins.erase(std::remove(Tins.begin(), Tins.end(), delTriangle[0]), Tins.end());
		borderLines.push_back(delTriangle[0].edges[0]);
		borderLines.push_back(delTriangle[0].edges[1]);
		borderLines.push_back(delTriangle[0].edges[2]);
		for (int j = 0; j < borderLines.size(); j++)
		{
			Tins.push_back(Triangle(borderLines[j].a, borderLines[j].b, points[i]));
			//Tins.back().outVerticePoints(); //debug
		}
		//检查新生成的三角形（的各边）是否符合Delaunay准则
		for (int j = 0; j < borderLines.size(); j++)
		{
			LegalizeEdge(points[i], borderLines[j], Tins);
		}
	}
	else //numberPointOnTriangle =2 在三角形边e上，去掉重复的边e，确定与e关联的两个三角形，并删掉。连接检索点P与四条边界边，生成四个新三角形
	{
		Tins.erase(std::remove(Tins.begin(), Tins.end(), delTriangle[0]), Tins.end());
		Tins.erase(std::remove(Tins.begin(), Tins.end(), delTriangle[1]), Tins.end());

		int index = -1;
		for (int m = 0; m < 3; m++)
		{
			if (delTriangle[1].IsEdgeOnTriangle(delTriangle[0].edges[m]) == 0) //不重合
			{
				borderLines.push_back(delTriangle[0].edges[m]);
			}
			else//重合，返回边的index
			{
				index = delTriangle[1].IsEdgeOnTriangle(delTriangle[0].edges[m]) - 1;
			}
		}

		for (int m = 0; m < 3; m++)
		{
			if (m != index)
				borderLines.push_back(delTriangle[1].edges[m]);
		}

		for (int j = 0; j < borderLines.size(); j++)
		{
			newTriangle.push_back(Triangle(borderLines[j].a, borderLines[j].b, points[i]));
			Tins.push_back(Triangle(borderLines[j].a, borderLines[j].b, points[i]));
			//newTriangle.back().outVerticePoints();
		}

		//检查新生成的三角形（的各边）是否符合Delaunay准则
		for (int j = 0; j < borderLines.size(); j++)
		{
			LegalizeEdge(points[i], borderLines[j], Tins);
		}
	}
}
return Tins;
}

3.合法性检测

bool Delaunay2D::IsIllegal(Point p1, Point p2, Point p3, Point pk) 
{
	Point center = CircumCenter(p1, p2, p3);
	long double radius = sqrt(pow(p1.x - center.x, 2) + pow(p1.y - center.y, 2));
	long double dist = sqrt(pow(pk.x - center.x, 2) + pow(pk.y - center.y, 2));

	if (dist > radius) //合法的三角形
	{
		return false;
	}	
	else
	{
		return true;//不合法的三角形
	}	
}

void Delaunay2D::LegalizeEdge(Point &pr, Edge &edge, vector<Triangle> &triangle)
{
	vector<unsigned int> arr;
	Triangle borderTriangle;
	Point pk(0, 0);
	int triangleIndex = -1;
	for (unsigned int j = 0; j < triangle.size(); j++) 
	{
		int edgeIndex = triangle[j].IsEdgeOnTriangle(edge);
		if (triangle[j].IsPointAtTriangle(pr) && triangle[j].IsEdgeOnTriangle(edge))
		{
			continue;
		}
		if (0 == edgeIndex)
		{
			continue;
		}

		 if (1 == edgeIndex) //ab
		{
			pk = triangle[j].c;
		}
		else if (2 == edgeIndex) //bc
		{
			pk = triangle[j].a;
		}
		else // (3 == edgeIndex) //ac
		{ 
			pk = triangle[j].b;
		}
		triangleIndex = j;
		borderTriangle = triangle[j];
		break;
	}

	if (triangleIndex>=0) //存在相邻三角形
	{
		if (IsIllegal(pr, edge.a, edge.b, pk)) //不合法
		{
			vector<Triangle> triangleNew;
		    for (int i =0;i< triangle.size(); i++)
			{       
				int accTriangleRelated =0;
				if (triangle[i].IsEdgeOnTriangle(edge) !=0)
				{
					accTriangleRelated++;
				}
				else
				{
					triangleNew.push_back(triangle[i]);
				}
			}
			triangle.clear();
			triangle = triangleNew;

			triangle.push_back(Triangle(pr, edge.a, pk));
			triangle.push_back(Triangle(pr, edge.b, pk));

            Edge newEdge1(edge.a, pk);
			Edge newEdge2(pk, edge.b);

			LegalizeEdge(pr, newEdge1, triangle);
			LegalizeEdge(pr, newEdge2, triangle);
		}
	}
	
}

4.删除与超级三角形顶点相关三角形

vector<Triangle> Delaunay2D::RemoveSuperTriangles(vector<Triangle> &triangles)
{
	vector<Triangle> resultTriangles;
	for (vector<Triangle>::iterator i = triangles.begin(); i != triangles.end(); i++)
	{
		Triangle triangleTemp = *i;
		bool a = !triangleTemp.IsPointAtTriangle(superTriangle.a);
		bool b = !triangleTemp.IsPointAtTriangle(superTriangle.b);
		bool c = !triangleTemp.IsPointAtTriangle(superTriangle.c);
		if (a&&b&&c)
		{
			resultTriangles.push_back(*i);
		}
	}
	return resultTriangles;
}

实现效果

最后通过一组随机平面点集验证效果。
在这里插入图片描述
完整版代码已上传：
基于Lawson 插入法的Delaunay三角网格构造算法C++实现

转载请注明出处：https://blog.csdn.net/baidu_34931359/article/details/129962825?spm=1001.2014.3001.5501

参考资料：

Computational Geometry:Algorithms and Applications - Third Edition
M.de Berg等著；邓俊辉译计算几何：算法与应用

sankingvenice

关注

2
点赞
踩
9

收藏

觉得还不错? 一键收藏
打赏
0
评论
基于Lawson 插入法的Delaunay三角网格构造算法C++实现

若该边不合法，则通过Flip函数实现边翻转操作，同时对翻转后新生成的边也调用子函数LEGALIZEEDGE进行空圆检测，对于所有可能的新非法边，LEGALIZEEDGE都要递归地调用自己，逐一进行核查，如此持续直到将所有的非法边转换为合法边。考虑引入点Pr时的情况。首先，我们要在当前的三角剖分中，确定pr落在哪个三角形之内（具体方法，我们也将在后面介绍）。倘若Pr碰巧落在三角剖分的某条边PiPj上，我们就需要找到与PiPj关联的那两个三角形，然后将pr与对顶的那两个顶点Pk分别联接起来，生成新的边。
复制链接

扫一扫