6 TIN的构建
6.1 需求规格说明
【问题描述】
不规则三角网(TIN)是一系列互不交叉、互不重叠的连接在一起的三角形来表示物体/地形表面,能够很好的描述和维护空间关系,在地理信息系统里有广泛应用。
给定一系列点,利用这些点作为不规则三角网中三角形的节点,来构建TIN。TIN对三角形的几何形状有三个基本要求:
1、三角形的格网唯一;
2、每一个三角形都尽可能接近正三角形;
3、三角形边长之和最小,保证最近的点形成三角形。
在线参考资料:
不规则三角网的算法设计与实…_不规则三角网为什么进行lop优化-优快云博客
【测试数据】
TIN_data.csv为洪山区POI数据,包括POI的序号、经纬度等信息。
【基本要求】(60%)
(1)给出一系列点,通过这些点来构建TIN(数据量不小于500条)。
(2)输入点的序号和坐标,输出三角形边对应的点序号。
【提高要求】(40%)
(1)在数据量较大(数据量大于5000条)的时候,需要考虑更复杂的数据结构从而保证系统的稳定性和高效性。给出洪山区的所有POI,请利用这一相对较为庞大的数据,快速构建出不规则三角网,并进行可视化。
(2)可视化的过程可以采用第三方软件(比如QGIS)进行实现,可寻求学长和学姐的帮助。
6.2 总体分析与设计
(1)设计思想
①存储结构
②主要算法思想
Delaunay三角网是一种基于点集的三角网。Delaunay三角网为相互邻接且互不重叠的
三角形的集合,每一个三角形的外接圆内不包含其他点。Delaunay三角网由对应Voronoi多边形的点连接而成。Delaunay三角形有三个相邻点连接而成,这三个相邻顶点对应的Voronoi多边形有一个公共的顶点,此顶点是Delaunay三角形外接圆的圆心。
Delaunay三角网生成算法可以通过给定点集,构建满足Delaunay条件的三角网。基于凸包的Delaunay三角网生成算法,又称“增量法”,是一种广泛使用的Delaunay三角网生成算法。它的基本思路是:从已有的凸包开始,每次加入一个新的点,对其进行“翻转”操作,使其满足Delaunay条件,直到所有点都被加入到三角网中。三角网生成步骤如下:
1、构建点集的凸包,将凸包的所有边加入到三角网中,选取最短边作为基线。
2、遍历点集中的每个点。
3、对于加入的每个点,在初始基线右侧的离散点中查找距此基线距离最短的点,作为第三点,生成三角形,即找到所有与之相邻的三角形(每个三角形包含三个顶点),并计算它们的外接圆。
4、如果该点在某个三角形的外接圆内,则需要进行“翻转”操作,即删除该三角形,并将该点与与之相邻的三角形的未加入点构成的三角形加入到三角网中。
5、重复步骤3和4,直到所有点都被加入到三角网中。
简单来说,核心思想是构建凸包,通过三角剖分的方式将新的点加入到已有的三角形中,生成TIN(不规则三角网)。这个过程涉及到一系列的辅助函数,用于计算角度、距离等信息,以确定如何构建新的三角形。我们需要生成一定数量的随机点,这些点代表地图上的离散数据点。将这些点从经纬度坐标转换为屏幕坐标后,根据屏幕坐标将这些点分配到不同的区域中。完成分区匹配后,需要找到这些点的凸包。凸包是一种几何结构,它是由一组点的最外边界构成的。找到凸包后,我们将凸包的顶点按照极角从小到大进行排序。然后选择凸包上的两个点作为初始边,将凸包上的其他点逐一与初始边的中点连接,形成初始三角形。将这些三角形加入到所设置的三角形列表和边列表中。接下来进行三角剖分,直到所有点都被包含在三角形中。在每一步中,选择一个未使用的边作为“基线”,寻找下一个离中点最近的点,形成一个新的三角形。不断更新三角形列表和边列表,如果存在符合条件的点,则继续剖分;否则,结束算法。在剖分完成后,使用findEdgeMin函数找到构建的凸包中相邻两点距离最小的一对点,并将它们放在凸包的开始和结束位置,优化初始凸包,确保凸包的
边缘不过长。
(2)设计表示
TIN构建的UML类图如图6.2-1所示。
图6.2-1 程序UML类图
在TIN构建程序中,主要的类和函数包括:
BuildTIN 类:负责生成随机点、读取文件中的点、构建TIN以及管理点和三角形的数据。
CreateTin 方法:实现了TIN构建的核心算法,通过迭代的方式添加点到现有的三角网中。
MinDistence 方法:寻找点集中距离最小的两个点,作为构建TIN的起点。
FindPoint3 方法:寻找与已选两点构成三角形的第三个点。
IsRepeat 方法:检查一条边是否已经被两个三角形共用。
IsConCurrent 方法:确保新添加的点与现有三角形的两个顶点在同一侧,以维持TIN的结构。
maxAngle 方法:在候选点中找到与现有边构成最大角度的点,以优化TIN的形状。
GetAngle 方法:计算三个点构成的角度。
这些方法和类共同工作,实现了从点集到TIN的转换。
(3)详细设计表示
TIN构建的程序流程图如图6.2-2所示。
图6.2-2 TIN构建流程图
该程序的流程步骤如下所示:
1.初始化:程序开始时,通过BuildTIN类的构造函数初始化点集。如果提供了文件名,则从文件中读取点集;如果没有提供文件名,则生成指定数量的随机点。
2.寻找最小距离点:使用MinDistence方法找到点集中距离最小的两个点,这两个点将作为构建TIN的起点。
3.寻找第三个点:使用FindPoint3方法找到与前两个点构成三角形的第三个点。
4.构建TIN:通过CreateTin方法开始构建TIN。对于每个新添加的点,检查其与现有三角形的关系,并根据需要进行“翻转”操作,以确保新点满足Delaunay条件。
5.检查边的共用:使用IsRepeat方法检查新添加的边是否已经被两个三角形共用,如果是,则不添加该边。
6.检查点的同侧性:使用IsConCurrent方法确保新添加的点与现有三角形的两个顶点在同一侧,以维持TIN的结构。
7.寻找最大角度点:使用maxAngle方法在候选点中找到与现有边构成最大角度的点,以优化TIN的形状。
8.计算角度:使用GetAngle方法计算三个点构成的角度,用于确定点的添加顺序。
9.输出TIN:构建完成后,输出TIN的详细信息,包括三角形的顶点序号。
6.3 编码
【问题1】:最初在添加边时由于没有人为规定边的指向,同一个三角形一直再重复构建(即找右边的点时总能找到已经构建过的点)
【解决方法】:拓展边时,需要将三角形的另外两条边放到栈中,如拓展P0P1,则要按照向量的形式将P0P2和P2P1放到栈中(以及存到结构体中),一定要按照向量,不能用凸包构建中的余弦夹角。
【问题2】:在生成随机点后进行分区时,考虑到点是由qrand()函数随机生成,会出现坐标重复、坐标多次被分配的问题以及随机数生成范围的问题。
【解决方式】:为了解决坐标重复问题,通过引入三重嵌套循环来检查坐标是否重复。外层循环遍历生成的坐标点列表PointList。中层循环遍历分区的行,内层循环遍历分区的列。在进行坐标分配时,通过条件判断和标记机制来避免重复分配。具体做法是,在分配坐标点时检查其横纵坐标,将其分配到相应的区域中,并设置一个重复坐标的标记(PointList[i].pflag),以避免重复分配。为了确保坐标只分配到一个区域,使用continue语句来跳过已经被分配的坐标点。对于随机数生成范围的问题,我们可以通过(qrand()/double(RAND_MAX))来生成一个范围在[0,1)之间的随机数。完成以上处理,我们可以根据具体需求对随机数生成过程进行控制,避免超出预定的范围,从而有效解决坐标重复和多次分配的问题,并对随机数生成范围进行可控,分区过程将更加准确和可靠。
【问题3】:给出的POI点数据量达到了15万条,如何提高计算效率?
【解决方式】:在构建TIN时,首先使用Graham扫描算法找初始凸包,将问题简化为处理凸包上的点集的三角剖分,降低计算规模。算法采用顺序构建凸包,排序点集,按顺序遍历构建凸包,减少比较次数。通过逐步添加点,实时动态更新边和三角形的信息及标记边的状态和使用flag属性,提高算法效率。更新边和三角形信息时,判断基线,找到离中点最近的点,更新ToUsePoleList,构建新边并加入lineList,更新flag值,标记当前边为已使用。根据边的状态,更新triList。最后通过循环迭代,生成三角形,实时更新边和三角形的信息。从而本题能够快速加载5000条以下数据。经较长时间等待,可完成五万条数据输出,显示在屏幕上。也可以数据输出至txt文件,使用XY转线功能,能够完成所有点的TIN构建。
6.4 程序及算法分析
①使用说明
打开程序后,即可出现初始化程序样式,如图6.4-1所示。
图6.4-1 初始化程序
这里我使用GetPoint函数将读取的文件默认写死路径,方便进行测试,可以直接点击“生成点”即可随机加载点,这里主要是使用srandPoint函数进行伪随机数的生成,如图6.4-2所示。
图6.4-2 随机生成点
点击“构建TIN”即可以进行不规则三角网的构建,如图6.4-3所示。
图6.4-3 TIN的构建
在任务要求中还要求对构建的三角网的点进行编号从而支持查询的操作,这里笔者的软件也同样支持查询的操作,输入想要查询的点号,再点击“查询”即可得到相应的点位情况,如图6.4-4所示。
图6.4-4 查询TIN
也可以点击“清楚查询”后再对其它的点进行查询操作,如图6.4-5所示。
图6.4-6 其它点位查询
此外,程序还可以实现大规模数据的加载和TIN的生成,比如这里加载5550个点,如图6.4-7所示。
图6.4-7 随机点生成
不难发现,这里生成的点的密集程度远超此前加载550个数据点的情况,同样点击“构建TIN”后程序即会调用BuildTIN算法进行TIN的构建,不过由于数据点较多,因此这里的构建时间会比较久,在五千多个数据点的密度下,大概需要用时35秒钟左右,TIN如图6.4-8所示。
图6.4-8 TIN的构建
生成完成后的不规则三角网同样支持查询功能,查询后的结果如图6.4-9所示。
图6.4-9 大规模不规则三角网的查询
6.5 小结
TIN(Triangulated Irregular Network),即不规则三角网。这是一种用于表示三维地形表面的数据结构,常用于地理信息系统(GIS)中。TIN是由一系列相连的三角形组成的网络,每个三角形的顶点都与一个地理位置相关联,并且每个顶点都有一个高度值。这些三角形共同构成了一个表面,可以用来表示地形的起伏变化。这种数据结构的一个主要优点是它可以准确地表示复杂的地形,如峡谷、山脉和河流等。TIN模型相比于栅格模型(例如数字高程模型,DEM),在表现地形细节和复杂地形特征方面有着优势,可以更精确地表示地形的起伏变化,而且数据量相对较小,因此在GIS中有着广泛的应用,如水文分析、地形分析、视域分析等。
编写此题代码时,我一步步明白了生成TIN的具体流程,TIN的构建不是一个简单工程,需要一步步详细地处理,这包括从随机生成点开始,经过坐标转换、凸包生成,最终完成TIN的构建。在这个过程中,我对点的状态、边的状态进行了合理的标记和控制,以保证算法的正确性。这个过程让我认识到,在实际编程过程中,对数据状态的管理和控制是非常重要的。只有每个数据状态搞准了,算法才能跑对。构建TIN的过程中,由于已经有了凸包的代码编写基础,所以省略这部分Graham算法的编写,直接引用过来了凸包中的计算方法,为避免不必要的冲突,我将代码重新复制,代码有一定的冗余,此处在未来我会继续完成优化。进而我学习了如何利用基于最大内角的策略来逐步构建TIN。这个过程中,我深入理解了三角形的拓扑关系,并学会了如何通过不断寻找最大内角的点,逐步完成了TIN构建。
6.6 附录
//TIN算法
//构造函数生成指定数量的不重复的随机点
BuildTIN::BuildTIN(int num)
{
srand(time(NULL));
for (int i = 0; i < num; i++)
{
int X = rand() % 942 + 10;
int Y = rand() % 942 + 10;
QPointF point(X, Y);
bool flag = false;
for (int i = 0; i < PointList.size(); i++)
{
if (PointList[i].x() == X && PointList[i].y() == Y)
{
flag = true;
break;
}
}
if (flag == true)
{
continue;
}
PointList.push_back(point);
}
MinDistence();//寻找距离最小的两点
FindPoint3();//寻找第三点
}
//构造函数
BuildTIN::BuildTIN(QString fileName)
{
QFile file(fileName);
file.open(QIODevice::ReadOnly);
int i = 0;
while (!file.atEnd())
{
QString line = file.readLine();
if (i == 0)
{
i++;
continue;
}
const QStringList columnStrings = line.split(QLatin1Char(','));
QPointF point((columnStrings[2].toDouble() - 114) * 2000 - 300, (columnStrings[3].toDouble() - 30) * 2000 - 550);
PointList.push_back(point);
i++;
}
MinDistence();//寻找距离最小的两点
FindPoint3();//寻找第三点
}
//创建TIN网
void BuildTIN::CreateTin()
{
t1.clear();
t2.clear();
t3.clear();
//填入距离最短的两点
t1.push_back(point1);
t2.push_back(point2);
//寻找距离第一条边距离最短的第三个点
t3.push_back(point3);
int k = 0;//记录扩展的三角形数
int L = 1;//记录已有三角形数
int p1, p2, p3;//记录三角形的三个顶点
vector<int> reserve;//可能存在的扩展点
while (k != L)
{
p1 = t1[k];
p2 = t2[k];
p3 = t3[k];
k++;
//第一条扩展边(p1p2)没有重复使用
if (IsRepeat(p1, p2, L))
{
reserve.clear();
for (int i = 0; i < PointList.size(); i++)
{
if (i != p1 && i != p2 && i != p3 && IsConCurrent(p1, p2, p3, i))
{
reserve.push_back(i);
}
}
if (reserve.size() > 0)
{
int max = maxAngle(reserve, p1, p2);
t1.push_back(p1);
t2.push_back(p2);
t3.push_back(max);
L++;
}
}
//第二条扩展边p1p3没有重复利用
if (IsRepeat(p1, p3, L))
{
reserve.clear();
for (int i = 0; i < PointList.size(); i++)
{
if (i != p1 && i != p3 && i != p2 && IsConCurrent(p1, p3, p2, i))
{
reserve.push_back(i);
}
}
if (reserve.size() > 0)
{
int max = maxAngle(reserve, p1, p3);
t1.push_back(p1);
t2.push_back(p3);
t3.push_back(max);
L++;
}
}
//第三条扩展边p2p3没有重复利用
if (IsRepeat(p2, p3, L))
{
reserve.clear();
for (int i = 0; i < PointList.size(); i++)
{
if (i != p2 && i != p3 && i != p1 && IsConCurrent(p2, p3, p1, i))
{
reserve.push_back(i);
}
}
if (reserve.size() > 0)
{
int max = maxAngle(reserve, p2, p3);
t1.push_back(p2);
t2.push_back(p3);
t3.push_back(max);
L++;
}
}
}
}
//获取最小距离的两点
void BuildTIN::MinDistence()
{
int min1 = 0, min2 = 0;//最小距离的两点
int min = INF;//最小距离
for (int i = 0; i < PointList.size() - 1; i++)
{
for (int j = i + 1; j < PointList.size(); j++)
{
double length = pow((PointList[i].x() - PointList[j].x()) * (PointList[i].x() - PointList[j].x()) + (PointList[i].y() - PointList[j].y()) * (PointList[i].y() - PointList[j].y()), 0.5);
if (length < min)
{
min = length;
min1 = i;
min2 = j;
}
}
}
this->point1 = min1;
this->point2 = min2;
}
//寻找第三边
void BuildTIN::FindPoint3()
{
double meanX = double(PointList[point1].x() + PointList[point2].x()) / 2;
double meanY = double(PointList[point1].y() + PointList[point2].y()) / 2;
double min = INF;
int min3 = 0;
for (int i = 0; i < PointList.size(); i++)
{
if (i != point1 && i != point2)
{
double length = pow((PointList[i].x() - meanX) * (PointList[i].x() - meanX) + (PointList[i].y() - meanY) * (PointList[i].y() - meanY), 0.5);
if (length < min)
{
min = length;
min3 = i;
}
}
}
point3 = min3;
}
//判断p1p2是否被两个三角形共用
bool BuildTIN::IsRepeat(int p1, int p2, int length)
{
int sum = 0;
for (int i = 0; i < length; i++)
{
if (p1 == t1[i] && p2 == t2[i] || p1 == t2[i] && p2 == t1[i] ||
p1 == t1[i] && p2 == t3[i] || p1 == t3[i] && p2 == t1[i] ||
p1 == t2[i] && p2 == t3[i] || p1 == t3[i] && p2 == t2[i])
{
sum++;
}
if (sum >= 2)
{
return false;
}
}
return true;
}
//判断p1p2p3和p1p2p4是否在同一侧,确保三角形中没有其他点
bool BuildTIN::IsConCurrent(int point1, int point2, int point3, int p)
{
double a = double(PointList[point2].y() - PointList[point1].y()) / (PointList[point2].x() - PointList[point1].x());
double b = double(PointList[point1].x() * PointList[point2].y() - PointList[point2].x() * PointList[point1].y()) / (PointList[point2].x() - PointList[point1].x());
double fxy1 = PointList[point3].y() - (a * PointList[point3].x() - b);
double fxy2 = PointList[p].y() - (a * PointList[p].x() - b);
//位于对侧
if (fxy1 * fxy2 <= 0)
{
return true;
}
return false;
}
//获取最大张角位置的点
int BuildTIN::maxAngle(vector<int>& reserve, int point1, int point2)
{
int max = 0; double maxAngle = INF;
for (int i = 0; i < reserve.size(); i++)
{
double angle = GetAngle(point1, point2, reserve[i]);
if (angle < maxAngle)
{
maxAngle = angle;
max = reserve[i];
}
}
return max;
}
//获取p1p3p2的张角
double BuildTIN::GetAngle(int p1, int p2, int p3)
{
double d1 = (PointList[p2].x() - PointList[p3].x()) * (PointList[p2].x() - PointList[p3].x()) + (PointList[p2].y() - PointList[p3].y()) * (PointList[p2].y() - PointList[p3].y());
double d2 = (PointList[p1].x() - PointList[p3].x()) * (PointList[p1].x() - PointList[p3].x()) + (PointList[p1].y() - PointList[p3].y()) * (PointList[p1].y() - PointList[p3].y());
double d3 = (PointList[p1].x() - PointList[p2].x()) * (PointList[p1].x() - PointList[p2].x()) + (PointList[p1].y() - PointList[p2].y()) * (PointList[p1].y() - PointList[p2].y());
double angle = (d1 + d2 - d3) / (2 * pow(d1, 0.5) * pow(d2, 0.5));
return angle;
}
//获取点数
int BuildTIN::GetPointNum()
{
return PointList.size();
}
//获取三角形个数
int BuildTIN::GetTrangleNum()
{
return t1.size();
}
//获取第三个三角形的三个顶点
void BuildTIN::GetTranglePoint(int i, QPointF& p1, QPointF& p2, QPointF& p3)
{
p1 = PointList[t1[i]];
p2 = PointList[t2[i]];
p3 = PointList[t3[i]];
}
void BuildTIN::GetTranglePoint(int i, int& p1, int& p2, int& p3)
{
p1 = t1[i];
p2 = t2[i];
p3 = t3[i];
}
vector<int> BuildTIN::GetTrangleByPoint(int point)
{
vector<int> tempList;
for (int i = 0; i < t1.size(); i++)
{
if (t1[i] == point || t2[i] == point || t3[i] == point)
{
tempList.push_back(i);
}
}
return tempList;
}
//获取点
QPointF BuildTIN::GetPoint(int i)
{
if (i >= PointList.size())
{
return PointList[PointList.size() - 1];
}
//QPoint point((PointList[i].x() * 111 - 12654) * 6, (PointList[i].y() * 111 - 3300) * 5);
return PointList[i];
}