文章目录
测绘程序设计大作业——TIN三角网生成+等高线生成
学校布置的测绘程序设计大作业，花了一个星期从算法学习到数据结构设计完成到最后出成果，虽然以后不搞测绘，但是还是记录一下在学习过程中的心得。
图形库的选取
测绘的最终目的在我看来是将数据可视化，所以选取一个更加直观的图形显示的方法是这个作业最核心的基础。本来我是不打算做这个作业的，因为C/C++在图像显示方面我不太了解，直到我知道了easyX图形库的存在，才让这次作业成为了可能。
easyX官网上拥有对该图形库里包含的函数接口的介绍，简单易上手是其最突出的特点，我的代码中就是使用这个库将数据可视化的
数据的读取
数据读取我使用的还是C语言的文件读取接口fopen和fread将坐标数据从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三角网
空泛的定义：
- 三角网的格网唯一；
- 最佳三角形形状，尽量接近正三角形
- 三角形边长之和最小，保证最近的点形成三角形
剖分准则：（划分三角形的核心依据）
如何生成
数据结构的定义
一个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] == root.third)
continue;
}
if (p3 == Point(-1, NULL, min_x, min_y, 0))
{
p3 = mypoint[j];
}
else
{
if (JudgePoint(root.left_point, root.right_point, p3, mypoint[j]))
{
p3 = mypoint[j];
}
}
}
if (p3 != Point(-1, NULL, min_x, min_y, 0))
{
newpoint = p3;
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);
}
}
对递归生成的反思
虽然我们按照递归的思路生成了三角网，但是众所周知递归是一种很挫的方法，特别是在处理一些计算量比较大的运算的时候，由于每递归一次就要函数压栈一次，可能会导致：时间复杂度过大而导致无法运算出结果、甚至栈溢出。所以在此基础上对递归这种算法进行改进，上面递归的方法实际上是树的深度优先遍历，我们可以使用树的层序遍历（广度优先遍历）用空间换时间，来达到数据量大的时候的计算。
优化代码
这里就不赘述层序遍历如何实现，直接放代码了：
//层序遍历生成三角网
//生成第一个三角形
double r, a, b;
Point p3 = mypoint[0];
for (int i = 0; i < mypoint.size(); i++)
{
if (i != pt1 && i != pt2)
{
if (JudgePoint(MIN_line.left_point, MIN_line.right_point, p3, mypoint[i]))
{
p3 = mypoint[i];
}
}
}
myline.push_back(make_line(MIN_line.left_point, MIN_line.right_point, p3));
//line_to_gragh(MIN_line, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,RED);
myline.push_back(make_line(p3,mypoint[pt1],mypoint[pt2]));
//line_to_gragh(make_line(p3, mypoint[pt1], mypoint[pt2]), Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,RED);
myline.push_back(make_line(p3, mypoint[pt2], mypoint[pt1]));
//line_to_gragh(make_line(p3, mypoint[pt2], mypoint[pt1]), Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,RED);
mytriangle.push_back(Triangle(&myline[0], &myline[1], &myline[2]));
//开始循环遍历
Point p1;
Point p2; //用来记录 遍历时线两个端点的临时变量
for (int i = 0; i < myline.size(); i++)
{
/*setlinecolor(GREEN);
line_to_gragh(myline[i], Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,GREEN);
Sleep(TIME);*/
int ret = 0;
if (myline[i].flag > 0)
continue;
p1 = myline[i].left_point;
p2 = myline[i].right_point;
p3 = Point(-1, NULL, min_x, min_y, 0);
for (int j = 0; j < mypoint.size(); j++)
{
if (p1 == mypoint[j] || p2 == mypoint[j])
{
continue;
}
if (is_same_side(myline[i], mypoint[j],min_x,min_y))
{
continue;
}
if (p3 == Point(-1, NULL, 0, 0, 0))
p3 = mypoint[j];
else
{
if (JudgePoint(p1, p2, p3, mypoint[j]))
p3 = mypoint[j];
}
}
if (p3 != Point(0, NULL, min_x, min_y, 0))
{
Line ln13 = make_line(p3, p1, p2);
Line ln23 = make_line(p3, p2, p1);
Line* tri_l1 = nullptr; //要插入三角形的边的指针
Line* tri_l2 = nullptr;
int ret1 = 0;
if ((ret1=is_in_myline(ln13, myline))<0 )
{
myline.push_back(ln13);
/* line_to_gragh(ln13, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,RED);
Sleep(TIME);*/
tri_l1 = &myline[myline.size() - 1]; //说明ln13不在myline中属于新生成的边，直接指向myline数组新插入的最后一个元素
}
else
{
tri_l1 = &myline[ret1]; //ln13在myline中已经存在了，指针指向ret1下标处的myline数组的值
myline[ret1].flag=2;
}
int ret2 = 0;
if ((ret2 = is_in_myline(ln23, myline))<0)
{
myline.push_back(ln23);
tri_l2 = &myline[myline.size() - 1];
/*
line_to_gragh(ln23, Point(0, NULL, min_x, min_y, 0), scaling_x, scaling_y,RED);
Sleep(TIME);*/
}
else
{
tri_l2 = &myline[ret2];
myline[ret2].flag=2;
}
mytriangle.push_back(Triangle(&myline[i], tri_l1, tri_l2));
}
if (p3 != Point(0, NULL, min_x, min_y, 0))
myline[i].flag = 2;
else
myline[i].flag = 1;
}
这里还是要解释一下，关于生成过程中一些的细节：
- 类
class triangle中存的是类Line的地址，也就是vector<Line> myline里面元素的地址，由于vector容器的增容机制，我们要提前将vector的空间开好（reserve接口），不让会出现野指针的情况，至于为什么要用指针，是因为前面的Line的flag会被后面生成的过程中修改，这个修改需要找到对应的Triangle进行更新，效率上就会变低，而且存储指针能大大的减少空间 - 在
class Line中我们注意到我们有一个int flag成员变量，用来存储该直线被多少个三角形共用，这个虽然在三角网生成中好像无关紧要但是在后面的等高线生成中起着至关重要的作用，我们可以将所有flag为1的直线生成为绿色，flag为2的直线生成为红色：
优化后的生成图👇
等高线生成
等高线我的做法是：先找到z坐标的最大值和坐标的最小值，然后按照一定比率（代码中的宏定义：RATE），将高差分成若干段，每一段取一次等高线
关于等高线生成，我是这样理解的：两点生成一条直线，三点生成一个平面，三个数据结构：Point、Line、Triangle分别对应的几何意义是：点、线、面。而等高线生成可以看成的高线所在的登高面与Triangle所形成的的面立体几何关系：只可能是 相离、相交，而交线就是等高线。
除此之外交线还有如下性质：
- 性质一：一个三角形最多 要么有两条直线与等高面相交 要么没有直线与等高面相交
- 性质二：如果一个三角形有一条边与等高面相交，那么他的另外两条边一定有一条与等高面相交
等高线数据结构建立
等高线也需要数据结构的描述，所以也需要建立相应的数据结构，由于和和之前的一样这里就不赘述了代码如下
代码
等高线按照形状分成：
开曲线
等高线理论上在图幅无穷大、点无穷多的时候是闭合的，但是图幅显然不是无穷大的，所以就会出现不闭合曲线。开曲线的开始于三角网中的边界边，结束语三角网的的边界边。这里的边界边也就是上面的只被一个三角形占用的边，也就是flag==1的边，这里知道为啥在生成三角网的时候就标定flag，不然你这地方会很麻烦。
思路：
- 遍历所有直线，看是否存在
flag==1且与登高面相交的直线，如果不存在那等高线类型就是闭曲线，如果找到了就命名为first，并且找到包含直线first的三角形t，执行下面循环 - 根据性质二找到三角形t中另外一条直线next，并将first和next与登高面相交的点记录下来
- 判断next的flag值 ——如果为1 ，就跳出循环——如果为2，继续下面步骤
- 找到除了三角形t之外包含直线next的三角形n（下一个三角形）
- 将next的值给first，再把n的值给t
- 如果first的值为2，继续执行第二步，如果不是跳出循环
闭曲线
思路：
- 从开曲线判断的过来，找到任意一条直线图等高面相交的直线begin，和他所在的三角形t，用一个临时变量first来记录begin的直线
- -根据性质二找到三角形t中另外一条直线next，并将first和next与登高面相交的点记录下来
- 判断next的是否与begin的值相同 ——如果相同 ，就跳出循环——如果不相同，继续下面步骤
- 找到除了三角形t之外包含直线next的三角形n（下一个三角形）
- 将next的值给first，再把n的值给t
- 如果first的值与begin不相同，继续执行第二步，如果相同跳出循环
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
