测绘程序设计大作业——TIN三角网生成+等高线生成

Posted 正义的伙伴啊

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了测绘程序设计大作业——TIN三角网生成+等高线生成相关的知识,希望对你有一定的参考价值。

文章目录

测绘程序设计大作业——TIN三角网生成+等高线生成

学校布置的测绘程序设计大作业,花了一个星期从算法学习到数据结构设计完成到最后出成果,虽然以后不搞测绘,但是还是记录一下在学习过程中的心得。

图形库的选取

测绘的最终目的在我看来是将数据可视化,所以选取一个更加直观的图形显示的方法是这个作业最核心的基础。本来我是不打算做这个作业的,因为C/C++在图像显示方面我不太了解,直到我知道了easyX图形库的存在,才让这次作业成为了可能。

easyX官网上拥有对该图形库里包含的函数接口的介绍,简单易上手是其最突出的特点,我的代码中就是使用这个库将数据可视化的

数据的读取

数据读取我使用的还是C语言的文件读取接口fopenfread将坐标数据从txt格式的文件中读取出来,这里将读取方式都封装getpoint()函数中,vector<point> v作为输出参数将数据读取出来到main函数

用到的函数:

  • strtok:读取的是字符串,所以一切都要用字符分割函数来进行操作,这里的分隔符有:','用来分割同一行不同数据数据、'\\n'用来分割行与行
  • atof:将字符串转换成浮点数类型(但是会有精度的损失)

注意:C语言读取文件中会把文件最后一行读取两次,所以读取完成后vector<point>最后两个元素是相同的,这时候要pop掉一个

代码如下:

void get_point(std::vector<Point>& v)

	char* buffer = new char[1024 * 8];
	FILE* fd = fopen("data.txt", "r");
	int sz = fread(buffer, 1, 8 * 1024, fd);
	int num;
	double x;
	double y;
	double z;
	char* name;

	char* p = strtok(buffer, ",\\n");
	sz -= strlen(p);
	num = atof(p);

	p = strtok(NULL, ",\\n");
	sz -= strlen(p);
	name = p;

	p = strtok(NULL, ",\\n");
	sz -= strlen(p);
	x = atof(p);

	p = strtok(NULL, ",\\n");
	sz -= strlen(p);
	y = atof(p);

	p = strtok(NULL, ",\\n");
	sz -= strlen(p);
	z = atof(p);
	v.push_back(Point(num, name, x , y, z));
	z -= 5;

	while (sz > 0)
	
		char* p = strtok(NULL, ",\\n");
		sz -= strlen(p);
		num = atof(p);

		p = strtok(NULL, ",\\n");
		sz -= strlen(p);
		name = p;

		p = strtok(NULL, ",\\n");
		sz -= strlen(p);
		x = atof(p);

		p = strtok(NULL, ",\\n");
		sz -= strlen(p);
		y = atof(p);

		p = strtok(NULL, ",\\n");
		sz -= strlen(p);
		z = atof(p);

		v.push_back(Point(num, name, x , y , z));
		z -= 5;
	
	v.pop_back();

坐标显示转换问题

从文件里面读取了坐标之后,如何将其转换到屏幕的坐标系上又是一个问题,首先我们要了解的是屏幕上的坐标系是是y轴朝下,x轴朝右,坐标原点在窗口的左上角。而且我们在一开始的时候需要初始化窗口的大小(也就是图幅的大小),所以变换后的测绘坐标一定要在这个范围内,且点与点之间的相对关系不能改变。

-首先我们先求出测绘坐标系中的坐标的范围,也就是确定一个相框将所有点框住,接下来只需要将这个相框与屏幕的相框重合即可
实际上我在处理的时候稍微的将这个相框扩大了一点,因为后来在屏幕上打的点实际上是一个有半径的圆,为了让圆显示完整(点显示完整)可以稍微将相框扩大

  • 接下来我们将这个相框经过平移和放缩之后就可以与屏幕相框重合

代码

//先找到x,y,z方向上的最大值和最小值
	double max_x = mypoint[0].x, max_y = mypoint[0].y, min_x = mypoint[0].x, min_y = mypoint[0].y, min_z = mypoint[0].z, max_z = mypoint[0].z;
	for (auto& e : mypoint)
	
		if (e.x < min_x)
			min_x = e.x;
		else if (e.y < min_y)
			min_y = e.y;
		else if (e.x > max_x)
			max_x = e.x;
		else if (e.y > max_y)
			max_y = e.y;
		else if (e.z < min_z)
			min_z = e.z;
		else if (e.z > max_z)
			max_z = e.z;
	
	//让整个范围与坐标系空出 整个图幅的1/50 这样边界点能完整显示
	min_x -= (max_x - min_x) / 70;
	min_y -= (max_y - min_y) / 70;
	max_x += (max_x - min_x) / 70;
	max_y += (max_y - min_y) / 70;

	double scaling_x = (max_x - min_x) / WIDTH;   //这里求出x方向上的放缩因子
	double scaling_y = (max_y - min_y) / HEIGHT;  //求出y方向上的放缩因子

Delaunay三角网的递归生成算法

什么是Delaunay三角网

空泛的定义:

  1. 三角网的格网唯一;
  2. 最佳三角形形状,尽量接近正三角形
  3. 三角形边长之和最小,保证最近的点形成三角形

剖分准则:(划分三角形的核心依据)

  • 空接外接圆准则

    在三角网中过每个三角形的外界圆均不包含点集中的其余任意一点

  • 张角最大准则

    一点到基边的张角最大

如何生成

数据结构的定义

一个TIN三角网的基本组成单位是: 点,线,面(三角形),所以必须定义这三个数据结构用作为描述一个三角网的基础,其中一些数据结构中包含的参数会在算法中解释
点的定义

class Point

public:
	Point(int _num = 0, char* _p = nullptr, double _x = 0, double _y = 0, double _z = 0)
		:num(_num)
		, x(_x)
		, y(_y)
		, z(_z)
	
		if(_p)
		strcpy(p, _p);
	

	void operator=(const Point& m)
	
		num = m.num;
		if(m.p)
			strcpy(p, m.p);
		x = m.x;
		y = m.y;
		z = m.z;
	


	bool operator == (const Point& m)
	
		return (x == m.x) && (y == m.y);
	


	bool operator !=(const Point& m)
	
		return !((x == m.x) && (y == m.y));
	


	int num;        //序号
	char p[20];     //数据点名称
	double x;       //x坐标
	double y;       //y坐标
	double z;       //z坐标
;

线的定义

class Line

public:
	Line(Point m = Point(), Point n = Point(), Point z = Point(),int flag=0)   //所有线的生成都用这个定义
		:left_point(m)
		,right_point(n)
		,third(z)
		,flag(0)
	

	Line(const Line& l)
		:left_point(l.left_point)
		, right_point(l.right_point)
		, third(l.third)
		, flag(l.flag)
		



	void make_line(Point m = Point(), Point n = Point())   //所有线的生成都用这个定义
	
		left_point = m;
		right_point = n;
	


	bool operator == (const Line & l2)
	
		if ((this->left_point == l2.left_point && this->right_point == l2.right_point) || (this->left_point == l2.right_point && this->right_point == l2.left_point))
			return true;
		else
			return false;
	
	bool operator != (const Line& l2)
	
		return !(*this == l2);
	


	Point left_point;      
	Point right_point;        //一条直线的 两个顶点
	Point third;             // 如果有一条直线,那么一定存在一个三角形包含该直线(第一个三角形除外),这个点表示该三角形中与直线相对的点
	int flag;                //标识一条直线被多少个三角形共有
;

三角形定义

class Triangle

public:
	Triangle(Line* x = nullptr, Line* y = nullptr, Line* z = nullptr)
		:root(x)
		, left(y)
		,right(z)
	

	bool operator ==(const Triangle& x)
	
		if (root == x.root)
		
			if (left == x.left)
			
				if (right == x.right)
					return true;
			
			else if (left == x.right)
			
				if (right == x.left)
					return true;
			
		
		else if (root == x.right)
		
			if (left == x.left)
			
				if (right == x.root)
					return true;
			
			else if (left == x.root)
			
				if (right == x.left)
					return true;
			

		
		else if (root == x.left)
		
			if (left == x.root)
			
				if (right == x.right)
					return true;
			
			else if (left == x.right)
			
				if (right == x.root)
					return true;
			
		
		return false;
	


	bool operator !=(const Triangle& x)
	
		bool ret = (*this == x);
		return !ret;
	


	bool is_have_line(Line& l)
	
		return (*root == l) || (*left == l) || (*right == l);
	





	Line* root;
	Line* left;
	Line* right;
;

递归生成算法

递归生成的算法核心的思想是:
将三角形看成由一条边和一个点组成,初始条件下给定一条基线边(一般是离得最近的两个点连成的线),然后去找符合三角网定义的点,找到的点与 基线边的两点形成新的直线,在以新形成的边为基线向下生成。其实就是一颗二叉树,如果按照深度优先遍历:必须是前序遍历,按照根 左 右的顺序向下遍历。一直遍历到 以该节点为基线边时,再也找不到合适的点来形成三角形为止(这也是叶子节点的条件),所以这棵树在宏观上看是一边遍历一边向下生成。

宏观上了解之后,但是还是有一些非常的重要的小细节(当时就是这些小细节没想清楚卡了我很长时间):
给定一条基线边寻找的待定点形成三角形,待定点要满足那些条件?

  • 首先不能和基线边所在的三角形中与基线边相对点在基线边的同一侧(第一条基线边除外,因为第一条基线边没有相对点),因为除第一条基线边时任意给定的之外,其他所有的基线边都是由三角形生成的,所以一定含有相对的点,如果要找的相对点与目标点同侧,那么三角网就会相交。这也是为什么要在Line类中加入成员变量Point third来记录相对点!
  • 其次待定点不能是基线边端点 和 相对点这三个点,除此以外点集里面的所有点都有可能是待定点
  • 如果找到了符合条件的待定点判断是否满足 空接外接圆准则张角最大准则,满足的话就继续向下递归,不满足这条边就是叶子节点退出该层递归回到上一层

空接外接圆准则

生成过程视频👇
生成视频😋

//以最小外接圆内没有其他点为条件生成三角网
bool form_trangle_by_circle(Line& root, Point& newpoint ,std::vector<Point>& mypoint, std::vector<Line>& myline, std::vector<Triangle>& mytriangle, double& min_x, double& min_y)

	Point circle;
	double r;
	for (int i = 0; i < mypoint.size(); i++)
	
		if(mypoint[i]!=root.left_point && mypoint[i]!=root.right_point && (myline.size()==1 || (is_same_side(root, mypoint[i], min_x, min_y) == false && mypoint[i]!=root.third) ))  //与直线建立三角形的点必须是没有用过的
		
			//这里计算外接圆圆心的时候需要进行一个坐标平移,以(min_x,min_y)为原点的坐标系建立,否则计算外接圆的时候数据可能会溢出
			CircleCenter(root.left_point.x - min_x, root.left_point.y - min_y, root.right_point.x - min_x, root.right_point.y - min_y, mypoint[i].x - min_x, mypoint[i].y - min_y, &(circle.x), &(circle.y), &r);
			int flag = 1;
			Point p3 = Point(-1, NULL, min_x, min_y, 0);
			for (int j = 0; j < mypoint.size(); j++) //这里不用做检查 因为即使这个点是三角形的三点,会在圆的边界上
			
				if (is_in_circle(mypoint[j], circle.x, circle.y, r, min_x, min_y))   //所有点都必须不在圆内
				
					flag = 0;
					break;
				
			
			if (flag)
			
				newpoint = mypoint[i];
				return true;
			
		
	

	return false;

//
void TIN(Line* root,std::vector<Point>& mypoint,std::vector<Line>& myline,std::vector<Triangle>& mytriangle,double &min_x,double &min_y ,double scaling_x, double scaling_y)

	Point newpoint;
	if (form_trangle_by_angle(*root,newpoint,mypoint,myline,mytriangle, min_x, min_y))    //看是否还能形成三角形  能->返回真,并且将第三点的下标push到数组used_point中		          
	                                                                                                  //不能->返回假
		//根据line类型的参数寻找第三个点
		Line newleftline;
		Line newrightline;
		Line* l = nullptr, * r = nullptr;

		//line_to_gragh(*root, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y, GREEN);//解除注释查看生成过程
		//Sleep(TIME);


		newleftline=make_line((*root).left_point, newpoint,(*root).right_point);    //新形成的线的 左右点 一定要确定清楚 
		//line_to_gragh(newleftline, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y, RED);//解除注释查看生成过程
		//Sleep(TIME);
		
		newrightline=make_line((*root).right_point, newpoint, (*root).left_point);    //生成三角形的右边
		//line_to_gragh(newrightline, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y, RED);//解除注释查看生成过程
		//Sleep(TIME);
		
		//一线一点连成三角形,分别将线和三角形录入相应的数组

		if (std::find(myline.begin(), myline.end(), newleftline) == myline.end())
		
			myline.push_back(newleftline);
			l = &myline[myline.size() - 1];
			TIN(l, mypoint, myline, mytriangle, min_x, min_y, scaling_x, scaling_y);

		
		if (std::find(myline.begin(), myline.end(), newrightline) == myline.end())
		
			myline.push_back(newrightline);
			r = &myline[myline.size() - 1];
			TIN(r, mypoint, myline, mytriangle, min_x, min_y, scaling_x, scaling_y);
		
		Triangle newT(root,r,l);
		mytriangle.push_back(newT);

		//递归  根左右 前序遍历
		//TIN(left...............);
		//TIN(right..............);
	
	else
	
		//line_to_gragh(*root, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y, GREEN);//解除注释查看生成过程
		//Sleep(TIME);
	

张角最大准则

生成过程视频👇
生成视频

//以最大角为标准生成三角形
bool form_trangle_by_angle(Line& root, Point& newpoint, std::vector<Point>& mypoint, std::vector<Line>& myline, std::vector<Triangle>& mytriangle, double& min_x, double& min_y)

	double r;
	//for (int i = 0; i < mypoint.size(); i++)
	

		//这里计算外接圆圆心的时候需要进行一个坐标平移,以(min_x,min_y)为原点的坐标系建立,否则计算外接圆的时候数据可能会溢出
		int flag = 1;
		Point p3 = Point(-1, NULL, min_x, min_y, 0);
		for (int j = 0; j < mypoint.size(); j++) //这里不用做检查 因为即使这个点是三角形的三点,会在圆的边界上
		
			if (mypoint[j] == root.left_point || mypoint[j] == root.right_point)  //与直线建立三角形的点必须是没有用过的
			
				continue;
			

			if (myline.size() != 1)//如果不是第一条基线边,就要满足待定点和相对点不同侧
			
				//if((is_same_side(root, mypoint[i], min_x, min_y)) || (mypoint[i] == root.third))
				if (is_same_side(root, mypoint[j], min_x, min_y))
					continue;
				if (mypoint[j] 

以上是关于测绘程序设计大作业——TIN三角网生成+等高线生成的主要内容,如果未能解决你的问题,请参考以下文章

测绘程序设计大作业——TIN三角网生成+等高线生成

如何从ArcGIS中导出已经生成好的tin图层为tin文件?

请问Arcgis中由高程点生成TIN数据时,为啥失败?

gis10.0 生成tin,怎样剔除异常高程值?

arcgis怎么做点要素的核密度分析

arcgis9.2 中的3D analyst ——tin 模块不可用 怎么办?