目录

一、概述

1.1 三角网的介绍

1.2 Delaunay三角形

二、三角网生长算法

2.1 建立第一个三角形

2.2 扩展TIN

三、各部分代码实现

3.1 数据结构

3.2 算法过程

3.3 代码不足

四、完整算法代码

五、相关参考


一、概述

1.1 三角网的介绍

        三角网是由一系列连续三角形构成的网状的平面控制图形,是三角测量中布设连续三角形的两种主要扩展形式,同时向各方向扩展而构成网状.适用于地势起伏大,通视条件比较好的场地。

        三角网是实现地形三维可视化数字地面模型(Digital Terrain Model,简称DTM)是一种很有效的途径。DTM主要是由栅格和不规则三角网(Triangulated Irregular Network,简称TIN)两种数据格式来表示,相比于栅格TIN具有许多优点,几乎能适用于任何复杂的地形,所以TIN是DTM常采用的一种格式。

1.2 Delaunay三角形

        Delaunay三角网的定义: 它是一系列相连的但不重叠的三角形的集合, 而且这些三角形的外接圆不包含这个面域的其他任何点。它具有两个特有的性质:

(1) 每个Delaunay三角形的外接圆不包含面内的其他任何点, 称之为Delaunay三角网的空外接圆性质, 这个特征已经作为创建Delaunay三角网的一项判别标准;

(2) 它的另一个性质最大最小角性质: 每两个相邻的三角形构成的凸四边形的对角线,在相互交换后,六个内角的最小角不再增大。

Delaunay三角网
Delaunay三角网

【举例】

        生成下图1、2、3、4、5、6、7、8点的泰森多边形

 

  1. 把点两两连起来,两两连接有很多种选择方式!于是使用了Delaunay的构造方法对其进行约束,得到了一个唯一的三角网(白色的部分,即为Delaunay三角网)
  2. 对三角形的每条边做中垂线,相交的点成了泰森多边形的顶点,构成了泰森多边形(紫色的多边形)

二、三角网生长算法

        Delaunay三角网有一个特性,每个三角网形成的外接圆都不包含其他参考点。利用这一个性质,我们可以直接构成Delaunay三角网:

2.1 建立第一个三角形

1、判断用来建立TIN的总脚点数,小于3则报错退出。
2、第一点的选择:
链表的第一节点,命名为Pt1;
3、第二点的选择:
A.非Pt1点; B.距Pt1最近命名为Pt2
4、第三点的选择
A.非Pt1,Pt2点;
B.与Pt1,Pt2点组成的三角形的外接圆内无其他节点;
C.与Pt1,Pt2组成的三角形中的角Pt1Pt3Pt2最大。命名为Pt3
5、生成三边,加入边表
6、生成第一个三角形,组建三角形表

2.2 扩展TIN

1、从边表头取一边,要求:该边flag标志为假(只在一个三角形中)
2、从点链表中搜索一点,要求:
A、与边中的Pixel3在边的异侧;
B、该点与边组成的三角形的外接圆内无其他点
C、满足上面两条件的点中角Pt1Pt3Pt2最大的点为Pt3。
3、判断新生成的边,如果边表中没有,则加入边表尾,设定标志flag为假,如果有,则设定该边flag为真。
4、将生成的三角形假如三角形表
5、设定选中的边的标志flag为真
6、转至1,直至边表边的标志flag全部为真。

注:如果需要进行限制三角剖分,则可利用重心法取出不要的三角形,必要时,可对边界进行限制,不让生成的边与边界相交

三、各部分代码实现

3.1 数据结构

struct Pixel    //点数据
{
	double x, y, z;
	int Num;
	bool flag;
};
struct List //数据链表
{
	Pixel *pixel;
	List *next;
};
struct Line //三角形边
{
	Pixel *pixel1;  //三角形边一端点
	Pixel *pixel2;  //三角形边另一端点
	Pixel *pixel3;  //三角形边所对顶点
	bool flag;
};
struct Linelist //三角形边表
{
	Line *line;
	Linelist *next;
};
struct Triangle //三角形表
{
	Line *line1;
	Line *line2;
	Line *line3;
	Triangle *next;
};

3.2 算法过程

此部分对原文代码进行注释,对每一部分进行剖析,方便进行理解。

首先是构建第一个三角形

void BuildFirstTri(const List* list,Triangle* &tglnode,Linelist* &linehead, Linelist* &linelast,Line* &ln1, Line* &ln2, Line* &ln3,Pixel* &pt1, Pixel* &pt2, Pixel* &pt3)
{
	pt1 = list->pixel;//数据链中的首点
	//找离pt1最近点pt2
	pt2 = list->next->pixel;//数据链表中的第二点

	List* node;
	node = list->next->next;
	while (node != NULL)
	{
		if (JudgeDis(pt1, pt2, node->pixel))
			pt2 = node->pixel;
		node = node->next;
	}

	//找第三个点
	//重新初始化node为链表第一点的next
	node = list->next;
	getPoint3(node, pt3, pt1, pt2);

	//根据三点构建三角形
	Linelist* linenode;
	//构建边ln1
	linenode = new Linelist();
	BuildLine(pt1, pt2, pt3, linenode);
	ln1 = linenode->line;

	//设置边链表的头尾指针
	linehead = linelast = linenode;

	//构建边ln2
	linenode = new Linelist();
	BuildLine(pt2, pt3, pt1, linenode);
	ln2 = linenode->line;
	//设置边链表的尾指针
	AddLine(linelast, linenode);

	//构建边ln3
	linenode = new Linelist();
	BuildLine(pt3, pt1, pt2, linenode);
	ln3 = linenode->line;
	//设置边链表的尾指针
	AddLine(linelast, linenode);

	//first Triangle

	//构建三角形
	tglnode = new Triangle;
	tglnode->line1 = ln1;
	tglnode->line2 = ln2;
	tglnode->line3 = ln3;
	tglnode->next = NULL;
}

完成第一个三角形的构建之后,将三角形链表进行更新,把第一个三角形加入其中。

//第一个三角形构建完毕!更新三角形链表
	tglhead = tgllast = tglnode;

再将点集、三角网指针、边的集合、点的集合传入TIN的构建函数之后,以第一个三角形构建的方式构建其他三角形,知道所有点均被遍历到结束,完成三角网的构建。

void BuildTIN(List* &list,Linelist* &linehead, Linelist* &linelast,Line* &ln1, Line* &ln2, Line* &ln3,Pixel* &pt1, Pixel* &pt2, Pixel* &pt3,bool &flag, Triangle* &tglnode, Triangle* &tgllast)
{
	Linelist *linetmp, *linetemp;
	List *pixeltmp;
	double x1, y1, x2, y2, x3, y3;
	linetmp = linehead;//三角网边链表的头指针
	while (linetmp != NULL)
	{
		if (linetmp->line->flag == true)//该边已经在三角网中
		{
			linetmp = linetmp->next;
			continue;
		}
		//获取当前边、当前边的结点pt1 pt2 及该边所属三角形中的第三个点(x3,y3)
		ln1 = linetmp->line;//初始化三角形边pt1 pt2 --ln1

		pt1 = linetmp->line->pixel1;
		pt2 = linetmp->line->pixel2;
		//当前三角形的三个结点坐标
		x1 = linetmp->line->pixel1->x;
		y1 = linetmp->line->pixel1->y;
		x2 = linetmp->line->pixel2->x;
		y2 = linetmp->line->pixel2->y;
		x3 = linetmp->line->pixel3->x;
		y3 = linetmp->line->pixel3->y;

		pixeltmp = list;//获取结点链表头指针
		//寻找下一结点pt3
		//getPoint3(pixeltmp, pt3, pt1, pt2);//有新增条件,不可复用
		pt3 = NULL;

		while (pixeltmp != NULL)
		{
			if (pixeltmp->pixel == pt1 || pixeltmp->pixel == pt2)
			{
				pixeltmp = pixeltmp->next;
				continue;
			}
			//判断该结点是否与该边所属三角形另一结点同侧--新增条件
			if (JudgeSameSide(x1, x2, y1, y2, x3, y3, pixeltmp->pixel))
			{
				pixeltmp = pixeltmp->next;
				continue;
			}

			//为pt3赋值
			if (pt3 == NULL)pt3 = pixeltmp->pixel;
			else//遍历所有结点,寻找最佳pt3
			{
				if (JudgePoint(pt1, pt2, pt3, pixeltmp->pixel))
					pt3 = pixeltmp->pixel;
			}
			pixeltmp = pixeltmp->next;
		}

		//初始化三角形边ln2
		if (pt3 != NULL)
		{
			linetemp = linehead;
			flag = false;
			//遍历所有边,将pt1 pt3所在边ln2 flag置为true
			setLineFlag(linetemp, pt1, pt2, pt3, ln2, linelast, flag);


			//初始化三角形边ln3
			linetemp = linehead;
			flag = false;
			//遍历所有边,将pt2 pt3所在边ln3 flag置为true
			setLineFlag(linetemp, pt2, pt1, pt3, ln3,linelast, flag);

			//构建三角形
			tglnode = new Triangle;
			tglnode->line1 = ln1;
			tglnode->line2 = ln2;
			tglnode->line3 = ln3;
			tglnode->next = NULL;
			//更新三角形链表
			tgllast->next = tglnode;
			tgllast = tglnode;
		}
		//当前边的flag置为true(遍历构建三角网完成)
		linetmp->line->flag = true;
		linetmp = linetmp->next;
	}
}

主体代码

//以下是程序中关于建网的部分:
//功能:   用给定的数据链表数据,组建Delaunay不规则三角网
//输入参数:数据链表list;
//输出参数:不规则三角网首三角形地址

Triangle* CreateDelaunayTIN(List *list)
{
	//组建第一个三角形
	
	//List *node;
	Pixel *pt1, *pt2, *pt3;//三角形三个顶点
	bool flag;
	Line *ln1, *ln2, *ln3;//三角形的三条边
	Linelist *linehead, *linelast;// *linenode ;//三角网的边链表
	Triangle *tglhead, *tglnode, *tgllast;//三角网的三角形链表

	BuildFirstTri(list,tglnode,linehead,linelast,ln1,ln2,ln3,pt1,pt2,pt3);

	//第一个三角形构建完毕!更新三角形链表
	tglhead = tgllast = tglnode;

	BuildTIN(list,linehead,linelast,ln1,ln2,ln3,pt1,pt2,pt3,flag,tglnode,tgllast);
	
	Triangle *TIN;
	//创建三角网
	TIN = tglhead;

	return TIN;
}

3.3 代码不足

1.通过对结点链表的一次次遍历,效率较低,还需优化。此外,算法大多数代码重复,需要对函数进行封装,提高代码复用效率。

2.对于前两个结点的确定,寻找与当前结点最近结点算法需优化;对于已确定三角形两顶点,寻找三角形第三个顶点方法需进行优化(减少遍历次数)

3.对于第一个点的确定,最后找位于点集中心的点而非随机取点

四、完整算法代码

此处代码为原始博客内容Delaunay三角剖分—三角网生长算法_sinat_27434117的博客-CSDN博客_三角网生长算法

此处为对该代码进行封装后,并完成测试的完整代码:

#include <iostream>
using namespace std;
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <list>

struct Pixel    //脚点数据
{
	double x, y, z;
	int Num;
	bool flag;
};
struct List //数据链表
{
	Pixel *pixel;
	List *next;
};
struct Line //三角形边
{
	Pixel *pixel1;  //三角形边一端点
	Pixel *pixel2;  //三角形边另一端点
	Pixel *pixel3;  //三角形边所对顶点
	bool flag;
};
struct Linelist //三角形边表
{
	Line *line;
	Linelist *next;
	Linelist() { line = new Line; }
};
struct Triangle //三角形表
{
	Line *line1;
	Line *line2;
	Line *line3;
	Triangle *next;
};

bool JudgeDis(Pixel* pt1, Pixel* pt2, Pixel *node) {//比较链表中node与pt2到pt1的距离
	return (pt1->x - node->x)*(pt1->x - node->x) + (pt1->y - node->y)*(pt1->y - node->y) < (pt1->x - pt2->x)*(pt1->x - pt2->x) + (pt1->y - pt2->y)*(pt1->y - pt2->y);
}

bool JudgePoint(Pixel* pt1, Pixel* pt2, Pixel* pt3, Pixel* node) {//比较链表中node pt3与pt1 pt2所成角的大小
	double dist11 = sqrt((pt1->x - node->x)*(pt1->x - node->x) + (pt1->y - node->y)*(pt1->y - node->y));
	double dist12 = sqrt((pt2->x - node->x)*(pt2->x - node->x) + (pt2->y - node->y)*(pt2->y - node->y));
	double dist12_3 = sqrt((pt1->x - pt2->x)*(pt1->x - pt2->x) + (pt1->y - pt2->y)*(pt1->y - pt2->y));
	double dist21 = sqrt((pt1->x - pt3->x)*(pt1->x - pt3->x) + (pt1->y - pt3->y)*(pt1->y - pt3->y));
	double dist22 = sqrt((pt3->x - pt2->x)*(pt3->x - pt2->x) + (pt3->y - pt2->y)*(pt3->y - pt2->y));
	return ((pow(dist11, 2) + pow(dist12, 2) - pow(dist12_3, 2)) / (2 * dist11*dist12) < (pow(dist21, 2) + pow(dist22, 2) - pow(dist12_3, 2)) / (2 * dist21*dist22)); //余弦判断角度
}

bool JudgeSameSide(double x1,double x2,double y1,double y2,double x3,double y3,Pixel* pixeltmp) {//判断该结点是否与该边所属三角形另一结点同侧
	//return ((pt2->y - pt1->y)*pixeltmp->x + (pt1->x - pt2->x)*pixeltmp->y + (pt2->x*pt1->y - pt1->x * pt2->y))*((pt2->y - pt1->y)*pt3->x + (pt1->x - pt2->x)*pt3->y + (pt2->x*pt1->y - pt1->x * pt2->y)) >= 0;
	return 	((y2 - y1)*pixeltmp->x + (x1 - x2)*pixeltmp -> y + (x2*y1 - x1 * y2))	*((y2 - y1)*x3 + (x1 - x2)*y3 + (x2*y1 - x1 * y2)) >= 0;
}

bool isInLineList(Pixel* pt1, Pixel* pt2,Line* linetemp) {//判断边pt1 pt2是否与linetemp相同
	return (pt1 == linetemp->pixel1 && pt2 == linetemp->pixel2) || (pt2 == linetemp->pixel1 && pt1 == linetemp->pixel2);
}

void BuildLine(Pixel* pt1,Pixel* pt2,Pixel* pt3, Linelist* &linenode) {//将当前边加入linenode边链表中
	linenode->line->pixel1 = pt1;
	linenode->line->pixel2 = pt2;
	linenode->line->pixel3 = pt3;
	linenode->line->flag = false;
	linenode->next = NULL;
}

void AddLine(Linelist* &linelast, Linelist* linenode) {//边链表中添加边
	linelast->next = linenode;
	linelast = linenode;
}

void getPoint3(List* &node, Pixel* &pt3, Pixel* pt1, Pixel* pt2) {//根据已知pt1 pt2 从node链表中寻找最佳pt3
	pt3 = NULL;
	while (node != NULL)
	{
		if (node->pixel == pt1 || node->pixel == pt2)
		{
			node = node->next;
			continue;
		}

		if (pt3 == NULL)
		{
			pt3 = node->pixel;
		}
		else
		{
			//比较当前结点pixel与pt3谁更合适--选取夹角最大的点
			if (JudgePoint(pt1, pt2, pt3, node->pixel))
				pt3 = node->pixel;
		}
		node = node->next;
	}
}

void setLineFlag(Linelist* &linetemp,Pixel* &pt1,Pixel* &pt2,Pixel* &pt3,Line* &ln,Linelist* &linelast,bool& flag)
{//遍历linetemp所有边,将pt1 pt3所在边ln的flag置为true,并更新三角网中的边链表linelast
	while (linetemp != NULL)
	{
		if (isInLineList(pt1, pt3, linetemp->line)) {
			linetemp->line->flag = true;
			flag = true;
			ln = linetemp->line;//初始化三角形边ln2
			break;
		}
		linetemp = linetemp->next;
	}
	if (!flag)//该边不在三角网的边链表中
	{

		Linelist* linenode = new Linelist();

		linenode->line->pixel1 = pt3;
		linenode->line->pixel2 = pt1;
		linenode->line->pixel3 = pt2;
		linenode->line->flag = false;
		linenode->next = NULL;

		//更新三角网的边链表
		linelast->next = linenode;
		linelast = linenode;

		ln = linenode->line;//初始化三角形边ln2

	}
}

void BuildFirstTri(const List* list,Triangle* &tglnode,Linelist* &linehead, Linelist* &linelast,Line* &ln1, Line* &ln2, Line* &ln3,Pixel* &pt1, Pixel* &pt2, Pixel* &pt3)
{
	pt1 = list->pixel;//数据链中的首点
	//找离pt1最近点pt2
	pt2 = list->next->pixel;//数据链表中的第二点

	List* node;
	node = list->next->next;
	while (node != NULL)
	{
		if (JudgeDis(pt1, pt2, node->pixel))
			pt2 = node->pixel;
		node = node->next;
	}

	//找第三个点
	//重新初始化node为链表第一点的next
	node = list->next;
	getPoint3(node, pt3, pt1, pt2);

	//根据三点构建三角形
	Linelist* linenode;
	//构建边ln1
	linenode = new Linelist();
	BuildLine(pt1, pt2, pt3, linenode);
	ln1 = linenode->line;

	//设置边链表的头尾指针
	linehead = linelast = linenode;

	//构建边ln2
	linenode = new Linelist();
	BuildLine(pt2, pt3, pt1, linenode);
	ln2 = linenode->line;
	//设置边链表的尾指针
	AddLine(linelast, linenode);

	//构建边ln3
	linenode = new Linelist();
	BuildLine(pt3, pt1, pt2, linenode);
	ln3 = linenode->line;
	//设置边链表的尾指针
	AddLine(linelast, linenode);

	//first Triangle

	//构建三角形
	tglnode = new Triangle;
	tglnode->line1 = ln1;
	tglnode->line2 = ln2;
	tglnode->line3 = ln3;
	tglnode->next = NULL;
}

void BuildTIN(List* &list,Linelist* &linehead, Linelist* &linelast,Line* &ln1, Line* &ln2, Line* &ln3,Pixel* &pt1, Pixel* &pt2, Pixel* &pt3,bool &flag, Triangle* &tglnode, Triangle* &tgllast)
{
	Linelist *linetmp, *linetemp;
	List *pixeltmp;
	double x1, y1, x2, y2, x3, y3;
	linetmp = linehead;//三角网边链表的头指针
	while (linetmp != NULL)
	{
		if (linetmp->line->flag == true)//该边已经在三角网中
		{
			linetmp = linetmp->next;
			continue;
		}
		//获取当前边、当前边的结点pt1 pt2 及该边所属三角形中的第三个点(x3,y3)
		ln1 = linetmp->line;//初始化三角形边pt1 pt2 --ln1

		pt1 = linetmp->line->pixel1;
		pt2 = linetmp->line->pixel2;
		//当前三角形的三个结点坐标
		x1 = linetmp->line->pixel1->x;
		y1 = linetmp->line->pixel1->y;
		x2 = linetmp->line->pixel2->x;
		y2 = linetmp->line->pixel2->y;
		x3 = linetmp->line->pixel3->x;
		y3 = linetmp->line->pixel3->y;

		pixeltmp = list;//获取结点链表头指针
		//寻找下一结点pt3
		//getPoint3(pixeltmp, pt3, pt1, pt2);//有新增条件,不可复用
		pt3 = NULL;

		while (pixeltmp != NULL)
		{
			if (pixeltmp->pixel == pt1 || pixeltmp->pixel == pt2)
			{
				pixeltmp = pixeltmp->next;
				continue;
			}
			//判断该结点是否与该边所属三角形另一结点同侧--新增条件
			if (JudgeSameSide(x1, x2, y1, y2, x3, y3, pixeltmp->pixel))
			{
				pixeltmp = pixeltmp->next;
				continue;
			}

			//为pt3赋值
			if (pt3 == NULL)pt3 = pixeltmp->pixel;
			else//遍历所有结点,寻找最佳pt3
			{
				if (JudgePoint(pt1, pt2, pt3, pixeltmp->pixel))
					pt3 = pixeltmp->pixel;
			}
			pixeltmp = pixeltmp->next;
		}

		//初始化三角形边ln2
		if (pt3 != NULL)
		{
			linetemp = linehead;
			flag = false;
			//遍历所有边,将pt1 pt3所在边ln2 flag置为true
			setLineFlag(linetemp, pt1, pt2, pt3, ln2, linelast, flag);


			//初始化三角形边ln3
			linetemp = linehead;
			flag = false;
			//遍历所有边,将pt2 pt3所在边ln3 flag置为true
			setLineFlag(linetemp, pt2, pt1, pt3, ln3,linelast, flag);

			//构建三角形
			tglnode = new Triangle;
			tglnode->line1 = ln1;
			tglnode->line2 = ln2;
			tglnode->line3 = ln3;
			tglnode->next = NULL;
			//更新三角形链表
			tgllast->next = tglnode;
			tgllast = tglnode;
		}
		//当前边的flag置为true(遍历构建三角网完成)
		linetmp->line->flag = true;
		linetmp = linetmp->next;
	}
}

//以下是程序中关于建网的部分:
//功能:   用给定的数据链表数据,组建Delaunay不规则三角网
//输入参数:数据链表list;
//输出参数:不规则三角网首三角形地址

Triangle* CreateDelaunayTIN(List *list)
{
	//组建第一个三角形
	
	//List *node;
	Pixel *pt1, *pt2, *pt3;//三角形三个顶点
	bool flag;
	Line *ln1, *ln2, *ln3;//三角形的三条边
	Linelist *linehead, *linelast;// *linenode ;//三角网的边链表
	Triangle *tglhead, *tglnode, *tgllast;//三角网的三角形链表

	BuildFirstTri(list,tglnode,linehead,linelast,ln1,ln2,ln3,pt1,pt2,pt3);

	//第一个三角形构建完毕!更新三角形链表
	tglhead = tgllast = tglnode;

	BuildTIN(list,linehead,linelast,ln1,ln2,ln3,pt1,pt2,pt3,flag,tglnode,tgllast);
	
	Triangle *TIN;
	//创建三角网
	TIN = tglhead;

	return TIN;
}

void printList(List* list) {
	cout << "List:\n";
	while (list != NULL) {
		cout << list->pixel->Num << ": " << list->pixel->x << "," << list->pixel->y << endl;
		list = list->next;
	}
}

int main()
{
    //此处已删除读入点集数据文件
    //pts信息:vector<MyPoint>
    //MyPoint数据结构如下
    //ptNum:点的序号 x:点的投影坐标X y:点的投影坐标Y z:点的高程Z
	List *ListHead, *ListNode, *ListLast;
	ListHead = ListNode =ListLast=NULL;
	for (int i = 0; i < pts.size(); i++) {
		Pixel* p = new Pixel();
		p->x = pts.at(i).x;
		p->y = pts.at(i).y;
		p->z = pts.at(i).z;
		p->Num = pts.at(i).ptNum;
		List* Node = new List();
		Node->pixel = p;
		Node->next = NULL;
		if (i == 0) { 
			ListHead = Node;
			ListNode = ListHead;
		}
		else {
			ListNode->next = Node;
			ListNode = ListNode->next;
			ListLast = ListNode;
		}
		
	}

	printList(ListHead);
	Triangle* tri = new Triangle();

	tri=CreateDelaunayTIN(ListHead);
	int i = 0;
	while (tri != NULL) {
		cout << i << ": " << " e1:  p1: " << tri->line1->pixel1->Num << " p2: " << tri->line1->pixel2->Num << " e2:  p1: " << tri->line2->pixel1->Num << " p2: " << tri->line2->pixel2->Num
			<< " e3:  p1: " << tri->line3->pixel1->Num << " p2: " << tri->line3->pixel2->Num <<"\t"<< endl;
		tri = tri->next;
		i++;
	}

	return 0;
}

输出结果

  

五、相关参考

参考内容:

三角网_百度百科

[GIS原理] 8.2 空间分析 - Delaunay三角网/Thiessen多边形(Voronoi图) - GIS开发者

Delaunay三角剖分—三角网生长算法_sinat_27434117的博客-CSDN博客_三角网生长算法

有问题联系作者。如有侵权,联系删除 

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐