一、研究背景

在GIS场景中，有很多面要素类型的数据，如地层面、水位面等，这些面的数据格式基本都是三角网，而基于这些数据生成等高线是用户常用的需求，网上搜索了一番也没有找到可直接用于生成的库，为此自研了基于python语言的三角网生成等高线的算法。

二、等高线的两种情况

使用过等高线的用户应该知道，在一个有限范围的平面图上，等高线是一系列高程相同的点的连线，由于范围有限，这些连线会出现闭合和不闭合两种情况。

三、算法思考

3.1 从几何的角度看等高线

1. 从宏观上看：等高线可以看成是具有相同高程间隔的一系列水平面与目标面（如地层面、水位面等）求交得到的交线集合。且不闭合的线的端点在面的边缘。
2. 从细节上看：等高线是具有相同间隔的z平面集与三角网的三角形的三条边的交点的连线的集合。且不闭合的线的端点在三角形的非共有的边（该边不同时属于两个三角形）上。

3.2 基本准则

根据上述的思考，建立如下算法准则

1. 如果三角形的某条边与某z平面相交，则有且仅有另一条边与该z平面相交
2. 如果等高线的起点为三角形的非共有的边，则终点必然为另一个三角形的非共有的边
3. 如果一个点为闭合等高线上的点，则不论沿哪个方向搜索，都能遍历完所有的点并回到自身原点
4. 非闭合等高线的点从端点开始可遍历完线上所有的点，且不会回到自身原点

3.3 实现算法

根据上述准则，建立如下某z平面与三角网求交线的算法

1. 先从三角网数据中解析出所有不重复的线段，并记录每条线段所属的三角形索引号
2. 计算每条线段与z平面的交点（插值计算），如果有则保存点坐标并记录该点所属的线段索引号
3. 为每个点设置是否被使用的标记
4. 先寻找属于非闭合等高线的点：循环遍历交点集，判断点是否被使用，已使用则跳过，未使用则根据点所属的线索引号找到线所属的三角形列表，判断列表大小，如果为2，则说明该线段为共有边，此时跳过该点；如果所属三角形列表大小为1，则说明为非共有边，可作为非闭合等高线的起点开始搜索。根据准则1，可以寻找当前点所属的线所属的三角形（因为是非共有边，可以确定该三角形）的另外一个交点，如此往复，可以得到一条非闭合等高线线的点集，这些点集都标记为已使用。一条线遍历完后，重复上述流程，直至找不到非共有边上的点。这样我们就得到了所有非闭合的等高线。
5. 然后寻找闭合等高线的点：循环遍历交点集，判断点是否被使用，已使用则跳过，未使用则根据准则1，可以寻找当前点所属的线所属的三角形（因为是共有边，有2个三角形，选任意一个即可）的另外一个交点，如此往复，可以得到一条闭合等高线线的点集，这些点集都标记为已使用。一条线遍历完后，重复上述流程，直至找不到未标记的点。这样我们就得到了所有闭合的等高线。

四 算法编写

1. 构建Tin类来承载三角网的数据结构以及相关计算方法，三角网的数据源选用市面上常用的stl文件（测试文件百度网盘下载链接: https://pan.baidu.com/s/1T1DRpI3nExAcryXt_izB9g 提取码: sxvg）
2. 由于用到了stl文件，需安装stl解析库 pip install numpy-stl

# pip install numpy-stl
# file: tin.py
from stl import mesh


class Tin(object):
    """三角网类：用于封装三角网的数据结构，以及外部三角网数据的读取、三角网生成等高线等功能
    """
    def __init__(self) -> None:
        self.__tin = None # 三角网格格式数据：由三个点坐标构成的三角形集合，如[[[0, 0, 0], [1, 0, 0], [1, 1, 0]], [[1, 0, 0], [2, 1, 0], [1, 1, 0]]]
        self.__box = None # 三角网的包围盒：由左下角点和右上角点组成，如[[0, 0, 0], [1, 1, 1]]

    def init_from_data(self, tin:list):
        """直接用三角网数据初始化
        :param list tin: 三角网格格式数据
        """
        self.__tin = tin
        
    def init_from_stl(self, path:str):
        """从stl文件读取三角网数据
        :param str path: stl文件路径
        :param bool return: True-读取成功/False-读取失败
        """
        try:
            data = mesh.Mesh.from_file(path)
            self.__tin = data.vectors.tolist()
            self.__box = [[data.x.min(),data.y.min(),data.z.min(),], 
                [data.x.max(),data.y.max(),data.z.max(),]]
               
            return True
        except Exception as e:
            print(repr(e))
            return False

    def get_tin(self):
        """获取三角网的数据
        :return list self.__tin: 三角网数据
        """
        return self.__tin

    def get_box(self):
        """获取三角网的包围盒
        :return list self.__box: 包围盒
        """
        return self.__box

    def generate_contour(self, z_interval:float):
        """生成等高线
        :param float z_interval: 等高线z间隔
        :return list contour: 等高线数据
        """

        # step1: 判断三角网是否有存在
        if self.__tin == None or len(self.__tin) == 0:
            print("Error: 三角网格不存在")
            return None

        # step2: 从tin中获取线的信息以及坐标的最值
        z_min, z_max, lines, lines_belongto_tri = self.__get_lines_from_tin()

        # step3: 获取等高线的z值列表
        z_set = self.__get_z_value_set_by_interval(z_min, z_max, z_interval)

        # step4: 使用每个z平面与三角网格面求交线
        contour = []
        for z in z_set:
    
            pts, pts_belongto_line = self.__get_inter_pts_between_lines_and_z_plane(lines, z)
            uses_pt = [False]*len(pts) # 点是否被使用,默认False
        
            polylines = []
        
            # 先获取非闭合的
            while True:
                start = -1
                for idx, use in enumerate(uses_pt):
                    if not use and len(lines_belongto_tri[pts_belongto_line[idx]]) == 1:
                        start = idx
                        break
                
                if start == -1:
                    break
        
                polyline = []
        
                self.__get_next_pt(start, uses_pt, pts_belongto_line, lines_belongto_tri, polyline)
        
                polylines.append(polyline)
        
            # 再获取闭合的
            while True:
                start = -1
                for idx, use in enumerate(uses_pt):
                    if not use:
                        start = idx
                        break
                
                if start == -1:
                    break
        
                polyline = []
        
                self.__get_next_pt(start, uses_pt, pts_belongto_line, lines_belongto_tri, polyline)
                polyline.append(polyline[0])
        
                polylines.append(polyline)
        
            # 将索引转成坐标
            for idx_line in polylines:
                contour_pts = []
                for idx_pt in idx_line:
                    contour_pts.append(pts[idx_pt])
            
                contour.append(contour_pts)
    
        return contour

    def __get_lines_from_tin(self):
        """从self.__tin中解析出所有不重复的线段以及点z坐标的最值
        :return float z_min: 三角网的最小z值
        :return float z_max: 三角网的最大z值
        :return list lines: 三角网中的线段(两个点)，每一项都是两个点的集合，如[[[],[]],[[],[]],[[],[]]...]
        :return list lines_belongto_tri: 每条线所属的三角形索引号，每一项都是一个列表，记录所属三角形索引号的集合
        """
        lines = []
        lines_belongto_tri = []
    
        z_min = 1e10
        z_max = -1e10

        # 遍历三角网
        for idxTri, tri in enumerate(self.__tin):
            # 一个三角形有三条线，默认设为未记录
            line1_exist = False
            line2_exist = False
            line3_exist = False

            # 判断线是否已被记录
            for idxLine, line in enumerate(lines):
                if (not line1_exist and tri[0] == line[0] and tri[1] == line[1]) or (tri[0] == line[1] and tri[1] == line[0]):
                    lines_belongto_tri[idxLine].append(idxTri)
                    line1_exist = True
    
                if (not line2_exist and tri[0] == line[0] and tri[2] == line[1]) or (tri[0] == line[1] and tri[2] == line[0]):
                    lines_belongto_tri[idxLine].append(idxTri)
                    line2_exist = True
    
                if (not line3_exist and tri[1] == line[0] and tri[2] == line[1]) or (tri[1] == line[1] and tri[2] == line[0]):
                    lines_belongto_tri[idxLine].append(idxTri)
                    line3_exist = True

            # 如果没有记录，则将线记录到lines，同时记录该线所属三角形的索引    
            if not line1_exist:       
                lines.append([tri[0], tri[1]])
                lines_belongto_tri.append([idxTri])
    
            if not line2_exist:  
                lines.append([tri[0], tri[2]])
                lines_belongto_tri.append([idxTri])
    
            if not line3_exist:       
                lines.append([tri[1], tri[2]])
                lines_belongto_tri.append([idxTri])
    
            # 获取Z的最值
            for pt in tri:
                if pt[2] < z_min:
                    z_min = pt[2]
    
                if pt[2] > z_max:
                    z_max = pt[2]
    
        return z_min, z_max, lines, lines_belongto_tri
    
    def __get_z_value_set_by_interval(self, z_min:float, z_max:float, z_interval:float):
        """根据阈值和步长得到z数组
        param: float z_min: z的最小值
        param: float z_max: z的最大值
        param: float z_interval: z的步长
        return: list z_set: z数值列表
        """
        z_set = []
        z_input = int(z_max)
        while True:
            if z_input < z_min:
                break
    
            z_set.append(z_input)
    
            z_input -= int(z_interval)
    
        return z_set

    def __get_inter_pts_between_lines_and_z_plane(self, lines:list, z:float):
        """获取线集和z平面的交点
        param: list lines: 线段集合
        param: float z: 高程值
        return: list pts: 交点集合
        return: list pts_belongto_line: 每个点所属的线段索引号
        """
        pts = []
        pts_belongto_line = []
        for idx, line in enumerate(lines):
            # 判断z是否在线段之内
            if line[0][2] <= line[1][2]:
                pt = self.__inter_by_z(line[0], line[1], z)
            else:
                pt = self.__inter_by_z(line[1], line[0], z)
    
            if pt != None:
                pts.append(pt)
                pts_belongto_line.append(idx)
    
        return pts, pts_belongto_line

    def __inter_by_z(self, pt_min:float, pt_max:float, z:float):
        """计算线段与z平面的交点
        param: list pt_min: 线段中z较小的点
        param: list pt_max: 线段中z较大的点
        param: float z: z值
        return: list pt: 交点
        """
        pt = None
        if pt_min[2] != pt_max[2] and z >= pt_min[2] and z <= pt_max[2]:
            rate = (z - pt_min[2])/(pt_max[2] - pt_min[2])
            x = pt_min[0] + rate*(pt_max[0] - pt_min[0])
            y = pt_min[1] + rate*(pt_max[1] - pt_min[1])
            pt = [x, y, z]
    
        return pt

    def __get_next_pt(self, start:int, uses_pt:list, pts_belongto_line:list, lines_belongto_tri:list, contour:list):
        """获取同一个三角形的其他边的交点
        param: int start: 交点的索引号
        param: list uses_pt: 交点是否被使用的标记
        param: list pts_belongto_line: 交点所属的线段索引号
        param: list lines_belongto_tri: 线段所属的三角形索引号
        param: list contour: 由点索引号组成的等高线数据
        """
        contour.append(start)
        uses_pt[start] = True
    
        for cur_tri_idx in lines_belongto_tri[pts_belongto_line[start]]:
            for idx_pt, use_pt in enumerate(uses_pt):
                if use_pt:
                    continue
        
                for idx_tri in lines_belongto_tri[pts_belongto_line[idx_pt]]:
                    if cur_tri_idx == idx_tri:
                        self.__get_next_pt(idx_pt, uses_pt, pts_belongto_line, lines_belongto_tri, contour)
                        return
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252

五 算法测试

使用python绘制三角网以及等高线如下图所示，结果正常。

# pip install matplotlib
# file: test.py
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from tin import Tin


# 参数准备：stl文件路径以及等高线的z间隔
path = "test.stl"
z_interval = 5

# 构建Tin类
obj = Tin()

# 读取stl文件中的三角网数据
obj.init_from_stl(path)

# 根据z间隔计算等高线
contour = obj.generate_contour(z_interval)

# 获取包围盒（绘制的时候用于控制显示范围）
box = obj.get_box()

# 获取用于绘制三角网数据
triangles = obj.get_tin()

# 初始化三维绘制视图
ax = plt.gca(projection="3d")

# 绘制三角网
ax.add_collection(Poly3DCollection(triangles))

# 绘制等高线
for polyline in contour:
    polyline_x = []
    polyline_y = []
    polyline_z = []

    for point in polyline:
        polyline_x.append(point[0])
        polyline_y.append(point[1])
        polyline_z.append(point[2])
    
    ax.plot(polyline_x, polyline_y, polyline_z, color="red")

# 设置视图的视角范围
ax.set_xlim([box[0][0], box[1][0]])
ax.set_ylim([box[0][1], box[1][1]])
ax.set_zlim([box[0][2], box[1][2]])

# 显示视图
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53