一种含边拓扑信息的不规则三角网弧扫式构建方案

1.一种用于地理信息系统的含边拓扑信息的不规则三角网弧扫式构建方法，其特征在于，包括以下步骤：
(1)使用含顶点、有向边、三角形的3种结构体类型定义，并创建这三种结构体类型的数组容器对象，用于存放这三种结构体类型的变量；
(2)确定一个参考中心位置，计算所有离散点相对参考中心的距离和方位角，将待构网的离散点集按照该距离从小到大升序排序；
(3)建立一个存放三角网外边界边序列的双向循环链表；并建立根据该双向循环链表中边的始点的方位角存放结点的方位角存储桶；方位角存储桶中只存放属于双向循环链表结点的记录，当双向循环链表中加入一条边时，由该边生成的结点对象的记录要按照该边起始方位角存入相应的方位角存储桶；当双向循环链表中要删去一条边时，相应的方位角存储桶中的记录也要删去；
(4)在排序过的点集中取最初三个点按逆时针顺序连成首三角形，并将三条边的记录以同样的逆时针顺序存入该双向循环链表中，形成初始三角网外边界；
(5)从点集中按序取下一个点，按照该点的所属方位角，从对应的方位角存储桶开始快速找出以右侧面向当前点的边，都作为与当前点连成新三角形的基边；查找基边时需要借助两种条件判断，条件1为待判断边是否以右侧面向当前点，条件2为当前点的方位角是否介于待判断边起始方位角和终止方位角之间；首先根据当前点的方位角在相应的方位角存储桶找到一个边结点，以确定待查找的多条边在双向循环链表中的大体位置；查找的方法是从当前方位角存储桶中找出使条件2为真的一个结点；若对应的方位角存储桶为空，则向方位角更大的后续方位角存储桶中查找，直至找到一个符合条件的结点时结束循环；对新找到的不为空的桶中查找，若找到一个能使条件2为真的结点则将其返回，否则直接返回当前桶中最后一个结点；接着从上面返回的结点开始，在双向循环链表中以循环的方式向当前结点的左右两侧查找能使条件1为真的其它边结点，直至遇上不符合条件1的任意边结点后结束循环；
(6)将步骤(5)中找出的每条基边与当前点构建成为三角形，将其注册到三角形数组中，同时生成和注册新三角形的三条邻边，即当前基边的反向边、当前基边始点与当前点构成的边即右侧边以及当前基边终点与当前点构成的边即左侧边；其中没有反向边的右侧边和左侧边需要加入该双向循环链表来实现链表更新，检查双向循环链表中是否有其反向边时，只需要直接将右侧边与当前基边的前驱比较，左侧边与当前基边的后继比较即可完成检查；当未找到右侧边或左侧边的反向边时，在双向循环链表中，将右侧边作为当前基边的前驱插入，将左侧边作为当前基边的后继插入，最后在双向循环链表中删除当前基边的记录；当新生成的左侧边或右侧边在双向循环链表中已有反向边时，在互设各自的反向边下标后，需要将原双向循环链表中反向边的记录删去；
(7)对每个新生成的三角形都要与其所有邻接三角形进行检验是否符合Delaunay三角网最优准则，否则交换三角形的对角边，并以递归的方式对交换后得到的三角形进行扩散式LOP优化；
(8)重复第(5)~(7)步，直至点集中所有的点都被处理过。
2.根据权利要求1所述的方法，其特征在于：顶点、有向边和三角形这三种类型对象之间的指针联系用它们在数组中的位置下标表示，其中有向边的定义设置有指向其左邻三角形和其反向边的数组下标索引，三角形的定义设置有指向其三条邻边的数组下标索引并使三边按照逆时针顺序排列。
3.根据权利要求1所述的方法，其特征在于：构网前需要计算所有离散点相对于所选参考中心的距离及方位角，方位角是参考中心-离散点构成的向量逆时针离开x轴正向的夹角。
4.根据权利要求1所述的方法，其特征在于：构网过程中要保持双向循环链表中的边按逆时针顺序相接，即一条边的始点与其前驱边的终点是同一顶点，该边的终点与其后继边的始点是同一顶点。
5.根据权利要求1所述的方法，其特征在于：方位角存储桶的建立是将方位角按等间隔划分，每个间隔建立一个对应的存放双向循环链表结点的方位角存储桶，每个方位角存储桶使用链表或动态数组类型，全部方位角存储桶用一个数组表示。
6.根据权利要求1所述的方法，其特征在于：在优化过程中交换两个三角形对角边时，修改两个三角形及两条反向边对象的属性，以及设置与其它4条邻边有关的拓扑属性。

一种含边拓扑信息的不规则三角网弧扫式构建方案\n技术领域\n[0001] 本发明涉及地理信息系统(GIS)以及激光雷达(LiDAR)数据处理所需的不规则三角网(TIN)的构建技术，尤其是含有边信息的不规则三角网的快速构建技术。\n背景技术\n[0002] 不规则三角网是GIS领域中一种重要的数字高程模型，由于它不受分辨率的限制，数据冗余度低，便于三维图形的绘制，在地形可视化、水文模拟模型领域具有广泛的应用前景。目前不规则三角网采用的数据结构多为基于三角形记录的结构，而不直接包含有关边的拓扑信息。这造成有关边的拓扑信息需要通过一定的搜索和计算来得到，影响拓扑信息获取的效率。而且，TIN的构建是以Delaunay三角网的原理在采样得到的平面上随机分布的散点基础上实现的。有关平面Delaunay三角网生成的方法较为丰富，但大多数主要为逐点插入法(增点发)、三角网扩张法、分治法(或分块法)以及扫描线法。现代测量手段不断发展，TIN要面对海量数据的应用，如由LiDAR获取的点云文件通常包含数百万个采样点。因此，Delaunay三角网生成方法的执行效率是十分重要的。对数十万至数百万个点的应用，传统的增点法、扩张法由于构网效率较低而几乎无法使用。目前，Fortune的扫描线法【 Fortune, S., 1987. A sweepline algorithm for Voronoi diagrams. Algorithmica, \n2(1): 153-174.】、分治法【Dwyer, R., 1987. A faster divide-and-conquer algorithm for constructing delaunay triangulations. Algorithmica, 2(1): 137-151】和Zalik的扫描线法【 Zalik, B., 2005. An efficient sweep-line Delaunay triangulation algorithm. Computer-Aided Design, 37(10): 1027-1038.】的执行效率最高，但这三种方法都各有不足之处，如Fortune使用抛物线边界来隔离已扫过的点，搜索较为复杂且计算量较大，Zalik法试图避免三角形生成过程中狭长三角形的生成以减少优化过程中的计算量，但考虑的问题分支变得复杂，编程较难实现；分治法的缺限是分块生成的各三角网的合并过程较为复杂，一般难以实现。刘永和等【刘永和, 王燕平, 齐永安, 2008. 一种快速生成平面Delaunay三角网的横向扩张法. 地球信息科学, 10(1): 20-25.】提出的横向扩张法也属于扫描线方法，构网效率较高，但该法有几个缺点：一是每次生成一条边后就得为查找其反向边而遍历整个三角网外侧边界上所有边；二是三角形优化采用的是在联网结束后的多次批量优化方式，而没有每构建一个三角形就进行递归式优化的方式，影响构网效率；三是每次向网络中增加一个点时需要从整个凸包边界中搜索边，仍有必要加速该搜索过程。此外，国内外都缺少有关使用扫描线方法生成三角网时高效维护边拓扑信息的方法。\n发明内容\n[0003] 本发明的目的是为解决现有Delaunay三角网数据结构不含有边信息以及现有扫描线方法存在的不足，而提出了一种针对含有边拓扑信息的Delaunay三角网的扫描线式构建方法。该方法不仅能够直接维护三角形与边之间的拓扑关系，还能够充分利用边信息来提高构建效率，使其成为最快的TIN构网方法之一。\n[0004] 本发明的技术方案包括以下步骤：\n[0005] (1)定义顶点、有向边、三角形3种数据类型，分别用三个数组存放；\n[0006] (2)确定一个参考中心位置，计算所有离散点相对参考中心的距离和方位角，将待构网的离散点集按照该距离从小到大升序排序；\n[0007] (3)建立一个存放三角网外边界边序列的双向循环链表，并建立一个根据链表中始点的方位角存放结点的方位角存储桶；\n[0008] (4)在排序过的点集中取最初三个点按逆时针顺序连成首三角形，并将三条边的记录以同样的逆时针顺序存入一个双向循环链表中，形成初始三角网外边界；\n[0009] (5)从点集中按序取下一个点，按照该点的所属方位角，从对应的方位存储桶开始快速找出以右侧面向当前点的边，都作为与当前点连成新三角形的基边；\n[0010] (6)将第(5)步中找出的所有基边与当前扫描过的点构建成为三角形，将其加入到三角形数组中，同时基边的反向边以及另外两条边(分别称为左侧边和右侧边)，更新外边界链表；\n[0011] (7)对每个新生成的三角形都要与其所有邻接三角形进行检验是否符合Delaunay三角网最优准则，否则交换三角形的对角边，并以递归的方式对交换后得到的三角形进行扩散式LOP优化；\n[0012] (8)重复进行第(5)~(7)步，直至点集中所有的点都被处理过。\n[0013] 本发明的有益效果是：\n[0014] (1)使用一种有向边数据结构，通过为它设置指向其反向边及其左侧三角形(即该边所属的三角形)的索引，增强了三角网的拓扑性能；并通过与三角形结构结合，能够隐式获得三角形之间的邻接关系。对三角形仅设置指向其三边的索引，省去了标识其顶点的信息，节省了空间占用和顶点信息维护时的执行时间。三角形顶点的信息可由其三条边的始点直接获得，这对三角网绘制等需要顶点信息的计算效率的影响很小。\n[0015] (2)顶点、有向边及三角形之间的相互索引通过组下标方式实现，其执行效率不低于常规内存指针方式，更重要的是它便于在磁盘上直接存储。\n[0016] (3)三角形三边(隐含了其三顶点)的存放以逆时针方式排列，能够保证三角形都位于属于它的每条边的左侧，而边的右侧为它所属三角形的外侧。因而三角网外围边界右侧的点必在已有三角网的外部，这便于在构网时只从已有边的右侧寻找待联入的点。\n[0017] (4)通过对点集按照与参考中心的距离进行排序，有助于以圆弧扫描方式从里到外逐个将点联入三角网。而且每次以扫描线方式找到下一个点时，就将所有符合联网条件的已有边界边与之建立三角形，从而使加入每个点后的三角网外边界为凸包，且凸包上各边以逆时针方向排列，这使凸包上任意边与其右侧的点连成的三角形不会与凸包边界存在交叉，因此避免了是否相交的判断过程，使流程得以简化。\n[0018] (5)建立相对参考中心空间方位的边界边结点存储桶，使查找边界边不需要遍历整个边界链表，加快了该查找过程。\n[0019] (6)边界边集合用链表存储，便于插入和删除。对每个新生成的边搜索链表中的反向边时，不需要遍历链表，而是通过直接检查当前边的前驱和后继即可，减少了运算量。\n[0020] (7)本发明中三角形结构没有邻接三角形的指针或索引，在LOP优化交换对角线时不需要管理涉及的两个三角形与周围4个邻接三角形的邻接关系，因此简化了LOP优化的步骤。\n[0021] (8)本发明方法是目前带边信息的Delaunay三角网的构网方法中效率最高的，且还具有简单、清晰、稳定的优势。生成的三角网带有完整的拓扑关系，可实现灵活的基于拓扑信息的搜索。\n附图说明\n[0022] 图1是圆弧扫描式不规则三角网构建及方位角的确定图示。\n[0023] 图2是构建首三角形时，前三个点排列不符合逆时针顺序时变换基边方向的示例。\n[0024] 图3是存放边界边的双向链表与方位角存储桶的数据结构与关系。\n[0025] 图4是向右扩展一个点时的示例。\n[0026] 图5是局部优化时交换三角形对角边的示例。\n具体实施方式\n[0027] (1)顶点、有向边及三角形数据结构定义\n[0028] 顶点的结构定义仅有横纵坐标值，无任何指针。有向边的结构定义则包含始点和终点、左侧邻接三角形，反向边等信息，这些信息全部采用在数组中的下标表示。为避免不必要的搜索，边结构还要包含本边在边数组中的下标。三角形的结构定义仅包含一个存放其三条邻边下标的数组和本三角形的下标，其中三条邻边的以逆时针顺序存放。两个相邻三角形之间的边为两条方向相反的有向边。\n[0029] struct DPoint: //顶点的定义\n[0030] double x;\n[0031] double y;\n[0032] 存放点的数组容器对象： vector Points; \n[0033] struct DEdge: //有向边的定义\n[0034] int beginPoint; //边的始点，用点在Points数组中的索引下标来表示[0035] int endpoint; //边的终点在Points中的下标\n[0036] int ltri; //本边的左邻三角形下标\n[0037] int ContraryEdge; //本边的反向边下标\n[0038] int Index; //本边的存放位置下标\n[0039] vector edgeList; //存放边的数组容器\n[0040] struct DTriangle: //三角形的定义\n[0041] int Edges[3]; //三角边的位置下标\n[0042] int index; //本三角形的存放位置下标\n[0043] vector triangleList; //存放三角形的数组容器[0044] list hull; //存放边界上边的链表 \n[0045] (2)参考距离的计算和排序\n[0046] 为点集选择一参考中心位置，可以是点集的位置重心，也可以是主观指定的大体上位于点集中心处的某一位置。先计算所有离散点与参考中心的距离，然后根据距离对点集按升序排序，以便使用圆弧扫描的方式将所有点按序加入三角网中。\n[0047] 如图1，圆内包含的点全部被三角化，随着圆弧的逐渐增大时每次扫描到一个点时，该点就要被联入三角网，这个扫描过程实际上是由对排序后的点集进行按序逐个处理来实现的。\n[0048] (3)方位角的确定和方位存储桶的建立\n[0049] 各点的方位角为按照逆时针方向离开x轴正方向的角度(0-360°) ，如图1中方位角α为向量OP 按逆时针方向离开x轴正向的夹角。方位角存储桶可按照5°或10°的间隔等分，如按10°间隔划分则有36个存储桶，用一个只有36个元素的数组存放，每个元素即存储桶可以是一个链表或动态数组。\n[0050] 方位存储桶中只存放属于边界边链表(hull)结点的记录，当边界边链表中加入一条边时，由该边生成的结点对象的记录要按照该边起始方位角存入相应的存储桶。当边界边链表中要删去一条边时，相应的存储桶中的记录也要删去。\n[0051] (4)边和三角形对象的创建与注册\n[0052] 每次通过指定始点索引和终点索引生成一条边后，要把该边在edgeList中注册，并把在数组中存放的下标赋值给边的index。边的ltri和ContraryEdge默认初始化为-1(表示空)。对三角形的处理与对边的处理类似。\n[0053] (5)首三角形的确定\n[0054] 本发明方法步骤的第四步是用最左侧的三个点p0,p1,p2构建首三角形，首先判断p2是否位于p0p1构成有向边的左侧，若是(符合逆时针排列规则)则直接生成首三角形Δp0p1p2及有向边p0p1、p1p2、p2p0，否则逆转前两个点的有向边为p1p0，生成首三角形Δp1p0p2及有向边p1p0、p0p2、p2p1，如附图2。\n[0055] 将首三角形在triangleList中注册后，新生成的三边的ltri设置为首三角形的索引，ContraryEdge保留为-1(空)。并将三边以逆时针顺序加入外边界链表hull中。\n[0056] (4)外边界边的快速查找\n[0057] 本发明的第五步中，每次扩张一个点p时，需要找出能与点p连成三角形的若干条边界边。查找边界边时需要借助两种条件判断：条件1为待判断边是否以右侧面向p点，条件2为点p的方位角是否介于待判断边起始方位角和终止方位角之间。\n[0058] 首先根据p 点的方位角在相应的存储桶中找到一个边结点，以确定待查找的多条边在链表中的大体位置。查找到的方法是从当前存储桶中找出使条件2为真的一个边结点。若对应的存储桶为空桶，则向方位角更大的后续存储桶中查找，直至找到一个符合条件的结点时结束循环。参照图3，假设点p的方位角为41.235°，将41.235取整再整除10，得到对应存储桶的位置是4，若第4个存储桶为空则接着检查第5个、第6个存储桶，以此类推，直至遇到不为空的桶结束。对新找到的不为空的桶中查找，若找到一个能使条件2为真的结点则将其返回，否则直接返回当前桶中最后一个结点。\n[0059] 接着从上面返回的链表边结点开始，在链表hull中以循环的方式向当前结点的左右两侧(即链表hull中的前驱和后继)查找能使条件1为真的其它边结点，直至遇上不符合条件1的任意边后结束循环。至此，查找过程全部结束。在查找过程中要用动态数组收集找到的所有符合条件1的边，在后续的步骤中将这些边与点p构建新三角形。\n[0060] 如图4，查找过程中，在存储桶中可能找到ab,bc,cd 三条边的结点的任意一个，比如首先从存储桶中找到了bc，接着使用两个循环分别借助结点的前驱和后继找到ab和cd，直至遇到ab 的前驱和cd 的后继，因它们不满足条件2，则分别结束两个循环。\n[0061] (4)向右扩张点\n[0062] 本发明方法步骤的第五步是从顶点数组中的下标为3的点开始，每次重复向外扩张一个点，直至所有点全部处理过。\n[0063] 每扩张一个点p时，先要判断hull链表中p位于哪些边的右侧，凡是符合该条件的边都要作为基边与点p构成一个新三角形。如图4(a)，e 为扫描到的新点，它位于ab,bc,cd 三条边界边的右侧，因此这三条边都要和e构建三角形。\n[0064] 本发明方法步骤的第六步中，每次生成一个新三角形时，需要在triangleList数组中注册该三角形。并生成和注册新三角形的三条邻边，同时要更新边界链表。这三条邻边分别是基边的反向边、新三角形的右侧边和左侧边。由于基边的始点总位于终点坐标之下，则始点与p构成的边即右侧边与基边的前驱相接，而其终点与p构成的边即左侧边则与基边的后继相接。\n[0065] 只有那些没有反向边的新边才能加入hull作为边界边，因此还需要检查hull中是否有其反向边。而右侧边只可能与原基边的前驱构成反向边，左侧边只可能与原基边的后继构成反向边，因此不需要遍历hull，只需要直接将右侧边与基边的前驱比较，左侧边与基边的后继比较即可完成检查。\n[0066] 当未找到右侧边或左侧边的反向边时，在边界链表hull中，将右侧边作为基边的前驱插入，将左侧边作为基边的后继插入，最后在hull中删除基边的记录。如图4(a)中，当ab 为基边时，则需要生成其反向边ba，以及右侧边ae、左侧边eb ；在hull中将ae 和eb分别作为基边ab 的前驱和后继插入，并删除基边ab，得到图4(b)的效果。\n[0067] 当新生成的左侧边(或右侧边)在hull中已有反向边时，则它与反向边均不再是外部边界，因此在互设各自的ContraryEdge后，需要将原hull中反向边的记录删去。如图\n4(b)中，当bc作为基边时，右侧边be与hull中原有的边eb 互为反向边，从hull中删除eb 的记录，而左侧边ec由于没有反向边，因此它作为bc 的后继被加入hull，最后bc的记录被从hull中删除，得到了图4(c)的效果。\n[0068] 在更新hull中结点的同时要维护结点在方位桶中的记录。\n[0069] (5)三角网的优化\n[0070] 本发明方法步骤的第七步中，当生成一个新三角形后，需要与其各个邻接三角形是否符合Delaunay三角网优化准则进行判断，若不符合就需要与其邻接三角形交换对角线。采用递归执行方式，将其发生交换后的三角形再执行同样的判断和优化。需要两个主要的子例程来执行：优化子例程和对角边交换子例程。\n[0071] 优化子例程中，给定参数是待优化的三角形索引，包括如下步骤：\n[0072] 获取当前三角形的备份对象，赋值给triangle；\n[0073] 令循环控制变量i=0；\n[0074] 若i<3，继续执行第 步，否则退出本例程；\n[0075] 获取triangle的第i个邻边和邻接三角形；\n[0076] 若当前邻接三角形存在，则执行第 步，否则执行i=i+1后执行第 步；\n[0077] 获得triangle中与当前邻接三角形不相邻的点，赋给ap；\n[0078] 若ap在当前邻接三角形的最小外接圆中，则调用对角边交换子例程完成对当前邻边两侧三角形对角边的交换；否则执行i=i+1后执行第 步；\n[0079] 以交换对角边后获得的两个新三角形的索引为参数，递归调用本子例程。\n[0080] 交换两个三角形对角边的示例如图5，给定参数为与待交换的两个三角形之间的一条邻边edge，具体步骤如下：\n[0081] 通过edge.ltri获取左侧三角形对象，赋给LTri；\n[0082] 若edge的反向边不为空，则通过edge反向边的ltri属性获得右侧三角形对象，赋值给RTri，否则退出本例程；\n[0083] 获取不在edge上的另外两个对角的顶点：令a=LTri中不在edge上的顶点，令b 为RTri中不在edge上的顶点；\n[0084] 创建新的三角形对象tri1和tri2，分别在数组中替换edge的左侧和右侧三角形对象；注意这时被替换掉的两个对象的信息仍可通过RTri和LTri获得。\n[0085] 创建b点为始点和a点为终点的新边ba，用它来替换edgeList中edge的记录，并创建ba 的反向边ab，用它来替换edgeList中edge的反向边的记录；\n[0086] ba的ltri设为tri2的记录，ab的ltri设为tri1的记录，并互设ContraryEdge属性；\n[0087] 为tri1设置三边：从原备份三角形RTri中找到以b点为始点的边作为第一边，从原备份三角形LTri中找到以a 点为终点的边作为第二边，以新边ab 作为第三边；\n[0088] 为tri2设置三边：以新边ba作为第一边,从原备份三角形LTri中找到的以a为始点的边为第二边,从原备份三角形RTri中找到的以b 为终点的边为第三边；\n[0089] 将tri1三边对象的ltri全部更新为tri1的记录，将tri2三边对象的ltri全部更新为tri2的记录；\n[0090] 退出本例程。\n[0091] (6)方法验证\n[0092] 本方法通过Visual C#进行了程序编写和测试验证。在内存为2G、CPU型号为Intel i5 M460及WindowsXP操作系统的笔记本电脑上，对10万个点构网只需1.1秒，100万个点构网只需14秒。优化后生成的三角网不发生任何构网错误和遗漏问题，且具有完善的拓扑关系。