CN105160706B - 一种单机多核环境下约束地形并行构建方法 - Google Patents
一种单机多核环境下约束地形并行构建方法 Download PDFInfo
- Publication number
- CN105160706B CN105160706B CN201510299433.6A CN201510299433A CN105160706B CN 105160706 B CN105160706 B CN 105160706B CN 201510299433 A CN201510299433 A CN 201510299433A CN 105160706 B CN105160706 B CN 105160706B
- Authority
- CN
- China
- Prior art keywords
- edge
- point
- parallel
- constraint
- algorithm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Image Generation (AREA)
Abstract
本发明公开了一种单机多核环境下约束地形并行构建方法，包括以下步骤：基于四方边缘结构的CD‑TIN数据结构及主要函数设计；基于负载平衡的数据划分策略；基于分治算法的子网并行构建；子网并行合并；格网索引构建；约束点地物并行插入；约束线地物并行插入；约束面地物并行插入。充分发挥并行计算技术为普通用户对地形高效构建的需求，给普通用户带来更强的体验。在现有技术基础上，设计一种合适的约束地形内存描述数据结构，并构建以之对应的数据库存储模型，进一步提出了基于负载平衡策略离散点数据划分策略；选择分治算法作为地形构建的基础算法，在空间索引技术支持下，制定了约束点地物、约束线地物和约束面地物的并行插入规则，并设计了并行插入算法。
Description
技术领域
本发明属于测绘、地理信息技术领域，涉及一种单机多核环境下约束地形并行构建方法。
背景技术
三角剖分可追溯到二十世纪三十年代，于1934年由俄国著名数学家Delaunay在解决数值分析问题提出的，然而当时的生成算法并不成熟，后经过学者们的努力探索研究，在七十年代后期得到了较大的发展和应用，如今已有多种Delaunay三角剖分算法被人们接受和采用。Tsai根据实现过程将剖分算法分成了三类，即分治算法、逐点插入法和三角网生长法。除此之外，此后的研究大多是基于这三类算法在两种方式上进行改进，一是对算法实现过程的某方面继续改进和完善，二是研究两种算法同时参与D-TIN构建的合成算法，当然还有学者提出了新的算法如：基于扫描线的剖分法、基于凸壳的剖分法、基于遗传算法的剖分法等。
Shamos和Hoey在1975首次提出了分治算法，算法设计的目的主要是为了生成Voronoi图，而Lewis和Robinson首次将分治思想应用于D-TIN构建中。Lawson在1977年首次提出了逐点插入法及其改进算法，Lee和schachter等人先后进行了改进和完善，这些改进的用于实现D-TIN构建的逐点插入法中，主要差别在于设置初始多边形和建立初始三角网的方法不同。Green和Sibson在1978年首次提出了三角网生长算法，并实现了一个生成Dirichlet多边形图的生长算法，McCullagh和Ross为了减少搜索第三点的时间，将点集进行了分块和排序从而达到缩短时间目的。Brassel,Reif,Maus等人进行了类似的改进，主要在改进搜寻“第三点”上入手。A.Mirante和N.Weingarten在1982年提出了辐射扫描法，然而这种方法产生的三角网并不是真正的D-TIN，主要因其不一定能满足外接圆性质。StevenFortune于1987年在Algorithmica上提出了V-图生成的平面扫描算法思想。J.R.Shewchuk等人应用扫描线算法实现了D-TIN构建。Marcelo Kallmann和Hanspeter Bieri等提出了在约束Delaunay三角网中进行约束点和线的插入和删除算法，并解决了自相交或重复的点自动检测，得到一个完全动态的约束Delaunay三角剖分，同时，在可视化、重构、地理信息系统等应用方面进行了讨论。
针对D-TIN串行算法在国内也有很多的研究成果，其中的研究涉及到无约束域和约束域的剖分算法改进、合并算法设计、在数据划分上改进、在数据结构上改进等。芮一康等提出了一种新的基于扫描线算法和分治算法的合成算法。算法兼顾了空间与时间性能，稳定性较高，运行效率和鲁棒性更优。刘云等在分别对比分治算法和逐点插入法后，针对两种算法存在的缺点提出了一些改进方案，主要是对逐点插入法中点的查询和三角形定位查询方面进行了改进，在点定位方面使用了网格存储空间数据，在三角形定位方面采用ClockWise检测的方向搜索技术。徐旭等在逐点插入算法的基础上，采用格网的方式同三角网结合建立“网格一点一三角形”之间的索引关系，从而大大提高了待插入点在三角形所在位置的定位效率。袁正午等针对定位待插点所在三角形效率不高的现状，提出基于对待插点集反复收集分配来来完成快速定位，并在数据结构和实现方式上进行了改进，使得总体时间复杂度为O(NlogN)。李立新等通过对已有算法的分析和借鉴，提出了两种多对角线交换算法并证明了循环算法不会存在死循环现象，并认为循环算法具有编程简单和运算速度快的特点。刘少华等提出了一种CD-TIN快速内插多边形算法，算法流程是先将多边形的边作为约束数据入网，然后将多边形内部多余边进行清空处理，且在影响区域及多边形内部三角形确定上给出了解决办法。贾晓林等针对约束线端点不在已有CD-TIN中、约束线通过三角形顶点、约束线影响域内有其他约束边三种特殊情况下的处理方法，并且提出了边的四相关局部优化算法。马洪滨等对于约束多边形的嵌入不在限定于规则化多边形，研究并给出了适用于任意多边形(含岛屿)的三角剖分算法。宋晓眉等深入探索了CD-TIN两步法，提出了采用递归割耳法解决影响区域包含悬挂点和影响区域为凹多边形。颜林等针对当前构建约束Delaunay三角网的算法在影响域为凹多边形时进行对角线交换可能失效的情况下，提出若为凹多边形则处理下一条相交边，处理所有相交边一遍后必会出现新的两个相邻三角形形成凸多边形，重新交换直到相交边处理完为止。张咏等继续拓展了当前约束多边形嵌入问题，研究和改进了线段相交判断和首三角形确定问题的相关算法，并给出了存在重复点情况下的解决方案。杨琦明等在D-TIN中嵌入约束边时，引入了Qi算法，减少了算法的运算次数，加快了嵌入的执行速度。任振娜在基于一次性构建CD-TIN算法上，从点数据方面入手，实现了快速、高效地插入点达到局部动态修改的目的。曾闽山等基于格网划分的自适应分割-合并算法，对已经通过块分割的格网数据重新排序再分割，然后以分割的逆序方式合并各子网。李小丽等利用线性四叉树方式将数据划分为包含点数相当的格网块，接着按照自下而上的方式将各格网块内子网进行无缝合并将，同时为了避免出现过小锐角的情况，加入了约束角来进行优化。刘永和等提出了一种构建三角形与边之间关系的数据结构，其中三角形包含了三边的对象指针，而对于边保存了左右邻接三角形的指针。李翔等提出了使用四方边缘结构来快速构建D-TIN，并将此数据结构应用于具体工程实践，证明了该结构的实用性和可靠性。谢增广提出了利用双链接边表(DCEL)结构结合分治算法实现了平面点集的三角剖分，并对特殊情形也做了处理。
D-TIN并行算法的研究开始于20世纪80年代末，基于分治算法的二维点集并行算法就是Davy首次提出的，D-TIN并行算法研究开始主要侧重并行索引机制的研究，然而较少涉及数据的并行划分，数据划分是并行算法中的一个不可缺少的步骤，每个核分配的数据若有较大的差异，负载将达不到均衡，从而严重影响算法性能。国内针对D-TIN并行算法的研究起步较晚，且主要以集群的分布式并行环境为主，对于单机多核方面的研究并不多。
在国外研究成果有：Kolingerová等介绍了一种基于Delaunay三角剖分的随机增量插入的并行算法，通过在多处理器工作站上测试证明算法可行、简单、易修改，并且对于约束三角剖分或者四面体都试用。Kohout等继续研究了随机增量插入法，由于它的简单性和稳定性被应用于E3展会，并基于这种方法研究了一种新的并行算法，在计算机体系结构的多核处理器和共享存储器环境下进行了测试。Kohout在共享内存系统的环境下，又利用逐点插入算法设计了基于外接圆准则的D-TIN并行算法。Foteinos等提出了三维D-TIN构建的并行算法，支持全动态并行插入和删除点，在共享内存下充分利用自定义内存管理和轻度锁机制减少了线程的通信和同步成本。
在国内研究成果有：张三元等人最初采用分治方法，即通过对散乱点集的分类，在不同处理器上分别求取分凸包，最后构造各分凸包的公共凸包。易法令等在分布式环境中提出了一种基于网格的D-TIN生成算法，算法较好的保证了负载均衡，并解决了四点共圆的不唯一性和边界处理的任意性问题。张三元提出的凸包并行化构造算法虽然使计算效率有所提高，但在删除内点上需要消耗大量时间，郝小柱认为通过采用并行的方法删除内点，对保留点排序形成简单有序多边形，最后对该多边形求取凸包，可以大大减小了时间复杂度。范刚龙等在分布式多处理机环境下提出对数据点按x坐标进行排序，并将排序后的数据按给定的阈值点数划分，通过分配给各个线程生成子三角网，相邻三角网两两归并；最后获得归并结果。申永源等在无约束域情况下提出了用串行算法的分治特征应用于并行生成泰森多边形，与串行算法相比提高了计算速度，并减少了执行时间。马劲松等针对线要素简化存在大量计算、难以做到实时的缺点，运用多核并行技术实现了Douglas-Peucker算法，并在多核处理器的计算机上进行实验，验证了并行算法的效率与实时性。张晓蒙等在共享内存模型下通过OpenMP采用任务并行的策略研究了并行Delaunay网格生成算法，实验结果表明与串行算法生成的网格在质量上区别较小，满足了大规模网格生成的需求。齐琳等针对传统的D-TIN并行算法尚未给出划分结果均衡、划分效率高效的理想解决方案情形下，在传统D-TIN并行算法规则条带划分方法的基础上，提出采用动态条带实现针对集聚分布点集数据的均衡、高效划分方法。
(1)与本发明相关的现有技术一为：
步骤1，将离散数据点在x、y投影平面上，沿先x方向排序，后y方向排序，将排序结果存入数组V[0…n-1]中；
步骤2，根据分布式环境中计算结点个数k及各计算结点的内存和计算能力的大小，将V[0…n-1]中的n个数据点分成m个对应长度的段V[S0，S1，…，Sm-1]；
步骤3，开辟一数组T[0…m-1]用于记录生成的初始子三角网，然后以每次k个段为单位，依次将Si，Si+1，…，Si+k-1分配给相应的结点，由该结点调用Delaunay三角网构建过程进行构网，并将构建的子三角网依次存入T[i…i+k-1]，换出到外存；循环上述过程，最终将在外存中形成m个初始归并段T[0…m-1]；
步骤4，用一链表T’记录生成的下一轮子三角网，根据内存的大小，依次从外存(T[0…m-1])给各结点调入相邻的若干个子三角网，并由该结点调用子网归并算法将这些子三角网归并为一个三角网，并按分配顺序依次将各结点合并所得的三角网插入链表尾部，换出到外存；然后再对后面的若干个相邻子三角网进行归并，重复上述过程直至m个三角网归并完毕为止；
步骤5，对链表T’重复步骤4中的过程，对相邻的子三角网进行下一轮归并，如此经过多轮的归并，最终将形成一个三角网。
(2)与本发明相关的现有技术二为：
步骤1，把离散点集V所在区域D划分成一定大小的正方形网格，点集V以网格为单元进行存储；
步骤2，主控处理机求点集V的凸包；
步骤3，按一定的间隔把凸包的边传给其它m-1个处理机，并相应建立边表和总边表，把点集V传送给其它m-1个处理机；
步骤4，建立总边表Gross_EdgeTable、新边表队列New_EdgeQueue及每个处理机边表Each_Table；
步骤5，每个处理机按规则从Each_Table中取出一条有向边，寻找其右边符合Delaunay三角新规则的点P，如果该点存在，删除Each_Table及Gross_EdgeTable中原来的边，连接点P构成新的边new_edge，并加到New_EdgeQueue；
步骤6，主控机将每个处理机生成的New_EdgeQueue进行合并，从而形成一个完整三角网。
现有技术存在以下不足：
1.现有技术并行构建建立在集群模式下，对于普通用户来说，代价太大，增加经济成本，个人难于承担，集群模式主要是为了满足那些有需求的部门和科研单位，而普通户面向的主要是目前流行的多核独立CPU的计算机。
2.现有的三角网并行构建技术主要面向的是无约束D-TIN的构建，而对于具有约束的D-TIN在并行环境下的构建需要考虑的问题更多，算法更复杂，因而较难实现。
3.数据划分过程中没有对每个处理机的计算能力，实现自适应负载划分。
4.在并行计算算法设计中没有考虑运用数据库实现数据的存储，数据都在内存中计算，内存消耗较大，而且难以支持海量数据的并行构建。
发明内容
本发明的目的是提供一种单机多核环境下约束地形并行构建方法，针对当前多核CPU计算机的计算性能越来越强大，以及内外存储容量越来越大情况，为了充分挖掘现有资源优势，提出了实现在该环境下的约束地形的并行构建算法，充分发挥并行计算技术为普通用户对地形高效构建的需求，给普通用户带来更强的体验。在现有技术基础上，设计一种合适的约束地形内存描述数据结构，并构建以之对应的数据库存储模型，进一步提出了基于负载平衡策略离散点数据划分策略；选择分治算法作为地形构建的基础算法，在空间索引技术支持下，制定了约束点地物、约束线地物和约束面地物的并行插入规则，并设计了并行插入算法。本发明构建的多核环境下约束地形并行构建算法，将促进空间矢量算子在单机环境下的并行化发展。
本发明所采用的技术方案是，一种单机多核环境下约束地形并行构建方法，包括以下步骤：
步骤1，基于四方边缘结构的CD-TIN数据结构及主要函数设计；
步骤2，基于负载平衡的数据划分策略；
步骤3，基于分治算法的子网并行构建；
步骤4，子网并行合并；
步骤5，格网索引构建；
步骤6，约束点地物并行插入；
步骤7，约束线地物并行插入；
步骤8，约束面地物并行插入。
步骤1的具体步骤为：
CCW代表逆时针方向，CW代表顺时针方向，该数据结构定义如下：
(1)假设初始边为e，则e有一个与其方向相反的伴生边e.twin；e.org为e的起始点坐标；e.dest为e的终点坐标；e.Rot和e.InvRot为e的对偶边；
(2)包含“next”的边是指以逆时针方向围绕邻面或顶点的下一条边；e.Lnext指逆时针围绕e边左面的下一条边，与e边同左面；e.Dnext指逆时针围绕e边终点的下一条边，与e边同终点；e.Onext指逆时针围绕e边起点的下一条边，与e边同起点；e.Rnext指逆时针围绕e边右面的下一条边，与e边同右面；
(3)包含“prev”的边是指以顺时针方向围绕邻面或顶点的下一条边；e.Dprev指顺时针围绕e边终点的下一条边，与e边同终点；e.Rprev指顺时针围绕e边右面的下一条边，与e边同右面；e.Lprev指顺时针围绕e边左面的下一条边，与e边同左面；e.Oprev指顺时针围绕e边起点的下一条边，与e边同起点；
为实现该数据结构，定义了三个数据类：节点类—VertexClass、边类—EdgeClass和三角形类—TriangleClass，以及一个操作类：UtilityClass；VertexClass记录节点三维坐标、关联边和布尔值的约束属性，EdgeClass记录四方边缘结构的各拓扑边、左右三角边形和布尔值的约束属性；UtilityClass中定义了CD-TIN的系列操作函数接口，
其中包含四个主要含数据接口：makeEdge、splice、connected和swap，函数功能描述如下：
(1)void splice(Edge*a,Edge*b)，该函数作用是对a、b边进行连接或者拆分；若将a、b边所在的两个具有共同顶点的独立子三角网合并，该函数将两个子三角网在共同顶点处关联在一起；同时，合并a和b左侧的面；相反，若要在共同顶点处沿a和b边拆分，该函数将在共同顶点处将三角网拆分成两个独立子三角网，同时，分裂a和b左侧的面；
(2)Edge*connected(Edge*a,Edge*b)，该函数作用是连接a的终点，b的起点生成新边c，在此过程中会建立a,b,c之间的四方边缘结构关系，next指边的Lnext，prev指边的Lprev；
(3)void swap(Edge*e)，该函数作用是用来交换对角线边，当用局部优化法LOP检测三角形不满足Delaunay法则时，则需要交换对角线边，将与e边相邻的两三角形组成的四边形另两顶点相连，形成新的边取代e，该边分别以e边起点和终点逆时针方向查找，第一点为新边的起点，第二点为终点。
步骤2的具体步骤为：
采用坐标排序法作为统一的数据划分策略，(1)用总点数N除于计算线程数K，假设每个线程平均分配点数为m；(2)将离散点集先X轴方向，后Y轴方向进行排序；(3)根据每个线程的点数m，将点号区间[i*m,(i+1)*m-1]，i为线程序号，i＝0,1,2,…,K，分配给每个线程，如果存在没有分配完的剩余点号，则全部给最后一个线程，所有线程共享同一份数据，每个线程计算时，根据点号区间从新排序点集中读取对应点。
步骤3的具体步骤为：以左右合并情况为例，具体实现过程如下：
(1)定义数据结构MaxEdge用于存储线程内部和线程间子网凸包的最左上/下和最右上/下边；
(2)定义左右凸包合并时的边指针参数：ldo、ldi、rdo和rdi，其中，ldo—指针为左凸包上以最左点为始点的逆时针方向的边；ldi—指针为左凸包上以最右点为始点的顺时针方向的边；rdo—指针为右凸包上以最右点为始点的逆时针方向的边；rdi—指针为右凸包上以最左点为始点的顺时针方向的边，该四个指针为全局变量，为了解决并行时数据共享冲突问题，在OpenMP并行模块头部设置指针变量为线程私有：#pragma ompthreadprivate(ldo,ldi,rdo,rdi)；
(3)分治算法程序模块，分治算法使用递归方式，将离散点集按子块数量在空间上分成左右均等的两个子集L、R，子块经过递归继续划分为更小子块直到划分数目size小于4，并且最小只会出现点集size为2和3的情况，此时，直接两两相连形成线段或三角形，递归下降过程结束之后，开始递归上升；在上升过程，size最小为4，左右子集L、R剖分形成左右凸包后，将凸包数据以MaxEdge结构方式保存在leftRet、rightRet，从左右凸包中可以获取指针ldi、ldo、rdi、rdo；
(4)多核并行构建程序模块。
步骤4的具体步骤为：
(1)获取左右凸包下公切线，
以左凸包最右边ldi和右凸包最左边rdi为起始边参数，通过判断各自起点org在对方边的哪侧，并结合数据结构快速找到下公切线；为了快速获取点在哪侧，使用向量叉乘k分量结果进行判断；假设有两个向量和这两个向量有共同点A，对这两个向量作叉乘运算：
假设A、B、C三点都为二维平面上的点则c1和c2都为零，那么有：
叉乘的k分量具有以下性质：
性质1：若(a1b2-b1a2)＞0,则在的顺时针方向，那么C点在的左边；
性质2：若(a1b2-b1a2)＜0,则在的逆时针方向，那么C点在的右边；
性质3：若(a1b2-b1a2)＝0,则与共线，但可能同向也可能反向，如果C点在A和B的最小外包MBR内，则同向且C点落在线段AB上；根据上述向量叉乘k分量性质，获取下公切线伪代码；
(2)获取最终候选点并连接合并，
假设下公切线记为baseLine，以其终点为起点、左候选邻居点为终点的边记为LEdge；以baseLine共起点、右候选邻居点为终点的边记为REdge，候选点判断过程如下：
①重新获取左右边：LEdge和REdge如果LEdge终点在baseLine右边，则判断baseLine终点、起点及LEdge终点三点构成的外接圆是否包含LEdge.Onext终点，若包含，则LEdge＝LEdge.Onext，重复①；如果REdge终点在baseLine右边，则判断baseLine终点、起点及REdge终点三点构成的外接圆是否包含REdge.Oprev终点，若包含，则REdge＝REdge.Oprev，重复①；
②获取最终候选点；如果LEdge终点不在baseLine右边，或者如果REdge终点在baseLine右边，且LEdge起点、终点及REdge起点三点构成的外接圆包含REdge终点，则由baseLine终点和REdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine；如果以上条件不满足，则从左凸包获取候选点，则由baseLine起点和LEdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine；最后，基于新下公切线，并回到①，重新获取LEdge、REdge；
③退出算法；如果LEdge和REdge的终点都不在baseLine右边，则此时baseLine为上公切线，算法结束；
(3)相邻子网两两递归并行合并，
根据并行构建子网的凸包找到相邻子网两两合并，并得到新的子网及其凸包；然后，通过新子网的凸包找到相邻子网进行合并，如此下去直到子网数为1，则合并结束。
步骤5的具体步骤为：
以D-TIN中所有三角形最小外包矩形MBR的长作为样本，并对该样本进行排序，得到一组从小到大的新样本，利用正态分布检验算法求出样本的均值μ和方差σ，接着以μ+σ为网格单元行高的初始值，每次增加0.1σ，增加m次，m由r决定，对于那些不服从正态分布要求的数据参照正态分布的“3σ准则”引进一个比例系数r，0<r≤1，r的默认值为[0.90，0.95]，最终以μ+σ+m*0.1σ作为网格单元的行高，从而得到网格划分的行数M，再以D-TIN中所有三角形最小外包矩形MBR的宽作为样本，确定网格划分的列数N。
步骤6的具体步骤为：
结合格网索引和四方边缘结构，加快点状约束地物影响区域的搜索速度，对正参与约束地物影响区域重构的三角形在内存中标记为：已占用，避免其它地物同时调用此三角形，从而避免算法混乱，计算完成后将参与计算的三角形从内外存删除，其它地物将从最新三角网搜索影响区域；
约束点并行插入时，对不同情形，采取不同并行策略：
(1)并行插入点的影响区域不重叠，则两线程同时执行对影响区域重新构网，完成后删除影响区域三角形，并更新三角网和格网索引；
(2)并行插入点的影响区域存在重叠，则不可以同时进行影响区域重新构网，需要待到约束点地物CP1完成计算后，约束点地物CP2再从新的三角网中，查找影响区域；对上述情形，CP2释放所占线程，进入下一约束点地物的插入，CP2则仍保留在未插入约束点地物集合中。
步骤7的具体步骤为：
(1)通过约束线地物外包矩形找到其所在范围的格网索引；
(2)根据格网索引初步找到约束线地物的影响区域三角形；
(3)从(2)的结果之中，逐一找到约束线地物顶点与相邻点间线段的影响区域三角形，并将找到的三角形标记为：已占用，此时存在以下两种情形：
(ⅰ)如果前一点与后一点的影响区域三角形相邻，则此两点的影响区域三角形也即为此两点间线段的影响区域；
(ⅱ)如果相邻两点影响区域三角形不相邻，则通过四方边缘数据的拓扑关系可以快速找到顶点V2影响区域三角形的相邻三角形C；然后，将线段V2V3与三角形C的其余两条边e1和e2作跨立实验，判断线段V2V3穿越哪条边；其次，根据穿越边找到下一影响区域三角形，如此下去，直到搜索到顶点V3所影响的三角形E为止。
约束线地物并行插入时可能存在的情形处理：
(1)约束区域不重叠时，约束线地物CPl1和CPl2各自所在影响区域则调用约束线地物插入算法并行执行局部重新构网；
(2)约束区域重叠时，约束线地物CPl1影响区域先行重新构网；则CPl2所占线程将其先跳过，同时，输入下一约束边线地物进行局部重新构网处理；
(3)最后，再从未处理的约束线地物堆栈对未插入的约束线地物进行并行插入，并对(1)、(2)的情形作与(1)、(2)的同样处理，如此往复循环，直到所有约束线地物插入为止。
步骤8的具体步骤为：
约束面地物的影响区域搜索则分解为两步：面地物边界影响区域搜索和面地物内部影响区域搜索，操作如下：
(1)面地物边界影响区域搜索调用步骤7的算法；
(2)面地物内部影响区域搜索；首先，得到面地物所在范围格网索引单元；然后，根据格网索引单元找到对应三角形，并排出边界影响区域三角形；最后，通过判断其余三角形任一顶点是否在面地物内，如果在，则整个三角形都包含在面地物内。
本发明的有益效果是并行算法确实可行，并在很大程度上提高了D-TIN构网性能。通过并行可以使整个效率明显提升，并且随着线程数的增加，执行时间逐渐缩短，整个过程主要是顶点嵌入和边嵌入最费时，而内部三角形清空占据总时间的极少部分，由40个多边形变为120个时所消耗的时间基本上不变，在此过程使用并行的加速比也是最小的，并行效率一般。
附图说明
图1是四方边缘结构图。
图2是共同顶点处的连接与拆分图。
图3是connected功能示意图。
图4是swap功能示意图。
图5是左右凸包指针。
图6是下公切线获取方法。
图7是左右凸包合并过程图。
图8是约束点地物插入情形图。
图9是约束线地物影响区域搜索图。
图10是约束线地物并行插入可能情形图。
图11是约束面地物影响区域搜索图。
图12是约束面地物并行插入可能情形图。
图13是约束面地物内部三角形清理图。
图14是并行与串行算法的耗时对比图。
图15是并行算法中各过程占用总耗时比例图。
图16是并行算法的加速比。
图17是并行算法的并行效率图。
图18是总体技术路线图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明涉及具有约束的Delaunay不规则三角网(Constrained TriangulatedIrregular Network，以下简称CD-TIN)数据结构设计、不带约束的不规则三角网(Triangulated Irregular Network，以下简称D-TIN)的单机多核并行构建算法等；在此基础上分别设计点状地物、线状地物、面状地物作为约束条件的单机多核并行嵌入D-TIN的算法。
单机多核环境下的约束地形并行构建关键在于约束D-TIN(CD-TIN)的并行构建，通常分两大步执行：首先，构建无约束D-TIN；然后，向其中插入约束点地物、约束线地物或约束面地物。目前D-TIN的构建方法主要有：分治法、逐点插入法和三角形生长法，属于算法密集型算子。根据CD-TIN的构建特点可以采用数据并行和任务并行的混合方式实施。CD-TIN可以分解成D-TIN构建和加入约束后的更新两步执行任务。D-TIN构建采用分治法，然后相邻子D-TIN逐一合并，最后逐一并行插入约束要素。
总体技术路线如图18所示。本发明所采取的方法在Visual Studio 2010平台下，利用OpenMP作为并行编程模型和Qt5Add-in插件设计并行程序界面，并采用C++编程语言实现了单机多核环境下的D-TIN并行构建算法。
本发明的单机多核环境下约束地形并行构建方法，分为非约束地形并行构建和约束地物并行插入两个过程。
步骤1，基于四方边缘结构(quad-edge structure)的CD-TIN数据结构及主要函数设计。
四方边缘结构描述是由Guibas和Stolfi在1985年提出的，它能完整的表达拓扑关系，所有边都可以通过任意边抵达，对于非邻边也能快速找到。该结构保存了边的对偶信息，使得Delaunay三角与Voronoi图之间的相互转换也非常容易。它的插入、删除操作也较方便实现，比其它结构更完善、灵活；在查询方面，也更方便、快速。开源的跨平台计算机视觉库OpenCV使用该结构来实现三角网的剖分。该数据结构如图1所示，CCW代表逆时针方向，CW代表顺时针方向。该数据结构定义如下：
(1)假设初始边为e，则e有一个与其方向相反的伴生边e.twin；e.org为e的起始点坐标；e.dest为e的终点坐标；e.Rot和e.InvRot为e的对偶边。
(2)包含“next”的边是指以逆时针方向围绕邻面或顶点的下一条边。e.Lnext指逆时针围绕e边左面的下一条边，与e边同左面。e.Dnext指逆时针围绕e边终点的下一条边，与e边同终点。e.Onext指逆时针围绕e边起点的下一条边，与e边同起点。e.Rnext指逆时针围绕e边右面的下一条边，与e边同右面。
(3)包含“prev”的边是指以顺时针方向围绕邻面或顶点的下一条边。e.Dprev指顺时针围绕e边终点的下一条边，与e边同终点。e.Rprev指顺时针围绕e边右面的下一条边，与e边同右面。e.Lprev指顺时针围绕e边左面的下一条边，与e边同左面。e.Oprev指顺时针围绕e边起点的下一条边，与e边同起点。
为实现该数据结构，定义了三个数据类：节点类—VertexClass、边类—EdgeClass和三角形类—TriangleClass，以及一个操作类：UtilityClass。VertexClass记录节点三维坐标、关联边和布尔值的约束属性，EdgeClass记录四方边缘结构的各拓扑边、左右三角边形和布尔值的约束属性等。UtilityClass中定义了CD-TIN的系列操作函数接口，其中包含四个主要含数据接口：makeEdge、splice、connected和swap，函数功能描述如下：
(1)void splice(Edge*a,Edge*b)该函数作用是对a、b边进行连接或者拆分。若将a、b边所在的两个具有共同顶点的独立子三角网合并，如图2(I)所示，该函数将两个子网在共同顶点处关联在一起；同时，合并a和b左侧的面，产生如图2(II)所示的结果。相反，若要在共同点处沿a和b边拆分，如图2(II)所示，该函数将在共同顶点处将三角网拆分成两个独立子网，同时，分裂a和b左侧的面，产生如图2(I)所示的结果。
(2)Edge*connected(Edge*a,Edge*b)该函数作用是连接a的终点，b的起点生成新边c，在此过程中会建立a，b，c之间的四方边缘结构关系，如图3所示，next指边的Lnext，prev指边的Lprev。
(3)void swap(Edge*e)该函数作用是用来交换对角线边，当用局部优化法LOP检测三角形不满足Delaunay法则时，如图4所示，则需要交换边线，将与e边相邻的两三角形组成的四边形(凸四边形)另两顶点相连，形成新的边取代e，该边分别以e边起点和终点逆时针方向查找，第一点为新边的起点，第二点为终点。
步骤2，基于负载平衡的数据划分策略。
数据划分结果的均衡性和划分算法的高效性是提高D-TIN并行算法加速比及并行算法效率达到最优的重要前提。目前针对D-TIN并行的数据划分方法主要有3种：①坐标排序的划分方法；②规则矩形划分方法；③基于投影分割的划分方法。从并行效率考虑，数据划分策略要考虑三个主要因素：首先，要考虑数据划分方法的适宜性；其次，数据划分算法必须简单高效，较少的占用总耗时；最后，划分的结果具有负载平衡最优，从而最大程度地减少计算线程之间的等待时间。本方法为了兼顾分治算法的效率，采用坐标排序法作为统一的数据划分策略，该方法简单高效，兼顾空间。具体划分思路：(1)用总点数N除于计算线程数K，假设每个线程平均分配点数为m；(2)将离散点集先X轴方向，后Y轴方向进行排序；(3)根据每个线程的点数m，将点号区间[i*m,(i+1)*m-1](i为线程序号，i＝0,1,2,…,K)分配给每个线程，如果存在没有分配完的剩余点号，则全部给最后一个线程。所有线程共享同一份数据，每个线程计算时，根据点号区间从新排序点集中读取对应点。
步骤3，基于分治算法的子网并行构建。
对逐点插入算法、生长算法和分治算法的时间复杂度进行了深入比较分析，从比较结果来看，分治算法的效率最高，最坏情况下的时间复杂度都能达到O(NlogN)，最好情况下其性能甚至能达到O(NloglogN)；而生长算法的效率最差，性能最高也只能达到O(N3/2)。分治算法不足在于由于递归的原因，剖分时消耗内存较大，但如今人们对时间的要求越来越高，并且随着计算机的发展，硬件将会更加强大，因而内存的影响将会得到缓解，也可以通过改进分治算法尽量使内存消耗最小化。分治算法最关键之处是要解决点集的划分，根据点集分割方法的不同，可以分为单向划分、双向划分、四叉树划分等。单向划分容易形成狭长三角形，增加交换边的次数；四叉树划分由于离散点的空间分布不均匀，可能出现叶子节点拥有1个或0个点的情况。本文在已经排序的点集基础上，将每个计算线程的数据实现基于自适应格网划分，根据分治算法要求，所有子集中的点数少于4个。分治算法关键步骤是子网凸包构建，凸包是线程内部和线程间子网合并重要依据，它随着子网的合并而不断更新。两子网的合并是从两个子网的合并起始线开始至终止线结束，伴随着三角形删除、生成的过程。相邻两子网合并由于位置的不同可分为左右和上下合并两种情况，原理相同，下文图7以左右合并情况为例。具体实现过程如下：
(1)定义数据结构MaxEdge用于存储线程内部和线程间子网凸包的最左上/下和最右上/下边。
(2)定义左右凸包合并时的边指针参数：ldo、ldi、rdo和rdi，其中，ldo—指针为左凸包上以最左点为始点的逆时针方向的边；ldi—指针为左凸包上以最右点为始点的顺时针方向的边；rdo—指针为右凸包上以最右点为始点的逆时针方向的边；rdi—指针为右凸包上以最左点为始点的顺时针方向的边，如图5所示。该四个指针为全局变量，为了解决并行时数据共享冲突问题，在OpenMP并行模块头部设置指针变量为线程私有：#pragma ompthreadprivate(ldo,ldi,rdo,rdi)。
(3)分治算法程序模块
分治算法使用递归方式，将离散点集按子块数量在空间上分成左右均等的两个子集L、R，子块经过递归继续划分为更小子块直到划分数目size小于4，并且最小只会出现点集size为2和3的情况，此时，直接两两相连形成线段或三角形，递归下降过程结束之后，开始递归上升。在上升过程，size最小为4，左右子集L、R剖分形成左右凸包后，将凸包数据以MaxEdge结构方式保存在leftRet、rightRet，从左右凸包中可以获取指针ldi、ldo、rdi、rdo，分治算法伪代码如下：
(4)多核并行构建程序模块
步骤4，子网并行合并。
子网合并是分治算法下三角网剖分的核心，也是并行计算结果完整性和准确性的保证。多核并行构建子网间的合并与串行算法内部子网合并原理具有一致性，都是根据左右凸包实现相邻子网的合并，关键之处在于提取左右凸包下公切线及与其构成Delaunay三角形的候选点。具体实现过程如下：
(1)获取左右凸包下公切线
以左凸包最右边ldi和右凸包最左边rdi为起始边参数，通过判断各自起点org在对方边的哪侧，如图6(a)所示，并结合数据结构快速找到下公切线，如图6(b)所示。为了快速获取点在哪侧，通常使用向量叉乘k分量结果进行判断。假设有两个向量和这两个向量有共同点A，对这两个向量作叉乘运算：
假设A、B、C三点都为二维平面上的点则c1和c2都为零，那么有：
叉乘的K分量在判断点与线段的位置关系具有十分重要的作用，它具有以下性质：
性质1：若(a1b2-b1a2)＞0,则在的顺时针方向，那么C点在的左边；
性质2：若(a1b2-b1a2)＜0,则在的逆时针方向，那么C点在的右边；
性质3：若(a1b2-b1a2)＝0,则与共线，但可能同向也可能反向，如果C点在A和B的最小外包MBR内，则同向且C点落在线段AB上。
根据上述向量叉乘k分量性质，获取下公切线伪代码如下：
(2)获取最终候选点并连接合并
如图7是左右凸包合并过程图，合并过程是通过下公切线往上获取构建Delaunay三角形候选点(也可以通过上公切线往下获取候选点)。当满足相邻三点构成的外接圆不包含其他候选点时，则该候选点为最终候选点，连接后作为新的下公切线，所有的候选点从下公切线右边的点集中选取；当左右候选点都在下公切线的左边，且当前下公切线为最顶边，即上公切线，则合并结束。假设下公切线记为baseLine，以其终点为起点、左候选邻居点为终点的边记为LEdge；以baseLine共起点、右候选邻居点为终点的边记为REdge，如图7(a)所示。候选点判断过程如下：
①重新获取左右边：LEdge和REdge如果LEdge终点在baseLine右边，则判断baseLine终点、起点及LEdge终点三点构成的外接圆是否包含LEdge.oNext终点，若包含，则LEdge＝LEdge.oNext，重复①；如果REdge终点在baseLine右边，则判断baseLine终点、起点及REdge终点三点构成的外接圆是否包含REdge.Oprev终点，若包含，则REdge＝REdge.Oprev，重复①。如图7(a)所示的情形，无需更新左右边。
②获取最终候选点：如果LEdge终点不在baseLine右边，或者如果REdge终点在baseLine右边，且LEdge起点、终点及REdge起点三点构成的外接圆包含REdge终点，如图7(a)所示，则由baseLine终点和REdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine；如果以上条件不满足，则从左凸包获取候选点，如图7(b)所示，则由baseLine起点和LEdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine。最后，基于新下公切线，并回到①，重新获取LEdge、REdge。
③退出算法：如果LEdge和REdge的终点都不在baseLine右边，则此时baseLine为上公切线，如图7(e)所示，算法结束。
(3)相邻子网两两递归并行合并
在分治算法中，子网合并在数据量比较大时，会严重影响算法的性能，应尽量避免数据分割的块数过多，通常根据计算线程数自适应划分，并行块数越多，最终比串行多出更多合并操作。子网并行合并过程如图8所示，根据并行构建子网的凸包找到相邻子网两两合并，并得到新的子网及其凸包；然后，通过新子网的凸包找到相邻子网进行合并，如此下去直到子网数为1，则合并结束。
步骤5，格网索引构建。
格网索引作用是为了当约束地物插入时快速找到其对D-TIN影响区域。格网单元大小的确定是构建格网索引的一个关键环节，太小会使得目标对象跨越格网就太多，检索的格网单元数就增多；而太大又会使得格网单元中落入对象数会增多，进行二次空间过滤的对象数就越多，两种情况都会影响空间检索效率。目前，确定格网单元大小比较常用方法是对目标对象MBR的长和宽做正态分布计算，确立划分格网最小单元的最佳行高和列宽。它的主要思路是：以D-TIN中所有三角形MBR的长作为样本，并对该样本进行排序，得到一组从小到大的新样本，利用正态分布检验算法求出样本的均值μ和方差σ，接着以μ+σ为网格单元行高的初始值，每次增加0.1σ，增加m次(m由r决定，对于那些不服从正态分布要求的数据参照正态分布的“3σ准则”引进一个比例系数r(0<r≤1)，r的默认值为[0.90，0.95])，最终以μ+σ+m*0.1σ作为网格单元的行高，从而得到网格划分的行数M。再以D-TIN中所有三角形MBR的宽作为样本，确定网格划分的列数N。
步骤6，约束点地物并行插入。
本方法结合格网索引和四方边缘结构，加快点状约束地物影响区域的搜索速度。本方法对已正参与约束地物影响区域重构的三角形在内存中标记为：已占用，避免其它地物同时调用此三角形，从而避免算法混乱，计算完成后将参与计算的三角形从内外存删除，其它地物将从最新三角网搜索影响区域。约束点并行插入时，可能存在如图8所示的情形，对不同情形，采取不同并行策略：
(1)对于图8(a)所示的情形，并行插入点的影响区域不重叠，则两线程同时执行对影响区域重新构网，完成后删除影响区域三角形，并更新三角网和格网索引；
(2)对于图8(b)、(c)所示的情形，并行插入点的影响区域存在重叠，则不可以同时进行影响区域重新构网，需要待到约束点地物CP1完成计算后，约束点地物CP2再从新的三角网中，查找影响区域。对上述情形，CP2释放所占线程，进入下一约束点地物的插入，CP2则仍保留在未插入约束点地物集合中。
步骤7，约束线地物并行插入。
约束线地物的影响区域搜索比约束点如地物复杂得多，如图9所示的约束线地物影响区域搜索过程：
(1)通过约束线地物外包矩形找到其所在范围的格网索引；
(2)根据格网索引初步找到约束线地物的影响区域三角形；
(3)从(2)的结果之中，逐一找到约束线地物顶点与相邻点间线段的影响区域三角形，并将找到的三角形标记为：已占用，此时存在以下两种情形：
(ⅰ)如果前一点与后一点的影响区域三角形相邻，则此两点的影响区域三角形也即为此两点间线段的影响区域，如图9中约束线CPl的顶点V1和V2；
(ⅱ)如果相邻两点影响区域三角形不相邻，如图9中约束线CPl的顶点V2和V3，则通过四方边缘数据的拓扑关系可以快速找到顶点V2影响区域三角形的相邻三角形C；然后，将线段V2V3与三角形C的其余两条边e1和e2作跨立实验，判断线段V2V3穿越哪条边；其次，根据穿越边找到下一影响区域三角形，如此下去，直到搜索到顶点V3所影响的三角形E为止。
约束线地物并行插入时可能存在的情形处理：
(1)约束区域不重叠时，如图10(a)所示，约束线段CPl1和CPl2各自所在影响区域可以调用约束线地物插入算法并行执行局部重新构网；
(2)约束区域重叠时，如图10(b)所示，此时，约束线地物CPl1影响区域先行重新构网；则CPl2所占线程将其先跳过，同时，输入下一约束边线地物进行局部重新构网处理；
(3)最后，再从未处理的约束线地物堆栈对未插入的约束线地物进行并行插入，并对(1)、(2)的情形作与(1)、(2)的同样处理，如此往复循环，直到所有约束线地物插入为止。
步骤8，约束面地物并行插入。
约束面地物的影响区域搜索可以分解为两步：面地物边界影响区域搜索和面地物内部影响区域搜索，如图11所示。操作如下：
(1)面地物边界影响区域搜索可以调用步骤7的算法；
(2)面地物内部影响区域搜索。首先，得到面地物所在范围格网索引单元；然后，根据格网索引单元找到对应三角形，并排除边界影响区域三角形；最后，通过判断其余三角形任一顶点是否在面地物内，如果在，则整个三角形都包含在面地物内。
约束面地物并行插入时可能存在的情形处理：
(1)约束区域不重叠时，如图12(a)所示，约束面polygon1和polygon2各自所在影响区域可以调用约束面地物插入算法并行执行局部重新构网；
(2)约束区域重叠时，如图12(b)所示，此时，约束面地物polygon1影响区域先行重新构网；则polygon2所占线程将其先跳过，同时，输入下一约束面地物进行局部重新构网处理；
(3)最后，再从未处理的约束面地物堆栈对未插入的约束面地物进行并行插入，并对(1)、(2)的情形作与(1)、(2)的同样处理，如此往复循环，直到所有约束面地物插入为止。
约束面地物形嵌入后，需要清空约束面内部三角形，即删除内部边，如图13所示。删除边后，四方边缘结构需要重建，整个清理过程可以并行完成。如何确定多边形包含的所有内部边是首要解决的问题，遍历所有三角形显然是不可取的，可以通过从嵌入边开始结合拓扑关系来判断，四方边缘结构保存了三角形之间的完整拓扑关系，正好可以满足要求，在上一步过程中已经将嵌入边加入到链表中，直接获取每个多边形地物的第一条嵌入边beginEdge即可，在边结构中加入一个是否为删除边的标志isMustDelete并假定一个存储多边形内部边的链表为InnerEdges。操作过程如下：
(1)先判断嵌入边的lnext、lprev是否为约束多边形边，如果不是则将lnext、lnext.twin、lprev、lprev.twin的isMustDelete设置为真，并在链表InnerEdges中加入lnext、lprev这两条边，然后循环获取链表数据，从中取出一条边firstEdge。
(2)如果firstEdge.twin.lnext不为多边形的边，则再判断其isMustDelete的真假，如果为假，将isMustDelete设置为真，并在链表InnerEdges中加入这条边，其反向边firstEdge.twin.lnext.twin的isMustDelete也设置为真，同理对于firstEdg.twin.lPrev也用样的方式判断，直到循环结束则链表InnerEdges中的所有边都是要删除的边。
删除的同时还需要重建结构关系，对于链表的删除来说，通过把删除边从链表中分离出来，改变链表中指针的指向，并用delete删除指针内存地址，同时设置为null来达到要求。
本发明专利所采取的方法在Visual Studio 2010平台下，利用OpenMP作为并行编程模型和Qt5Add-in插件设计并行程序界面，并采用C++编程语言实现了单机多核环境下的D-TIN并行构建算法。实验离散点个数(单位：万个)分别为：0.1、0.5、1、5、10、50和100，分别在两台操作系统均为Windows7的计算机上运行，硬件环境分别为：(1)联想启天M8250：CPU为Q9500 2.83GHz、内存4G、4核4线程；(2)兼容机：CPU为Xeon(至强)E3 3.30GHz、内存8G、4核8线程。测试指标主要包括执行时间、加速比和并行效率。
(1)无约束地形多核并行构建算法测试：
①时间性能
D-TIN串行和并行时间对比是基于所用剖分算法和数据结构等相同的情况下进行的，串行执行时间是指在单核单线下程序从开始执行到整个任务执行完所耗费的总时间；并行执行时间是从开始执行到最后一个线程执行完毕所耗费的总时间。
从图14可以看出当离散点数小于1万时，串行和并行构建D-TIN所花费的时间相当；当离散点数大于1万时，串行和并行耗时发生明显变化，串行所耗费的时间增长速度比并行快得多，同时，随着计算线程数的增加，并行所耗时间明显下降。从实验结果可知，串行耗时是四线程的两倍多，为八线程的四倍多。由此，可见发明设计的并行算法确实可行，并在很大程度上提高了D-TIN构网性能。
D-TIN并行构建的时间主要由数据划分排序时间、子网凸包构建时间、子网合并时间三部分组成，前者属于串行计算，后两者属于并行计算。
由图15可知，在并行算法中数据排序和子网合并这两个过程占用总时间非常少，主要时间集中在子网构建的时间开销上，无论四线程还是八线程，所占据的平均时间比例都超过90％，可见子网构建是整个并行算法中最费时的部分，需要对基础算法进一步优化来降低耗时。从图15中，还可以发现数据排序时间随着核数增多比例变大，原因是总耗时随线程数增多而减少，而数据排序为串行执行，耗时随线程数增多基本保持不变，因此，所占总耗时比增大；然而，子网合并总耗时比没有变大，反而减小，主要原因是子网合并也是并行的。
②加速比
加速比是指串行时间与并行时间的比值，它主要用来描述并行算法的加速性能。为了充分利用线程，本方法以计算机核数或一次最大可支持的逻辑线程数作为数据划分依据。由阿姆达尔定律可知，当串行比例越小，加速比越大，且最大不会超过核数，由图16可知四核四线程满足定律，而四核八线程并不满足该定律，因为四核八线程采用超线程技术，最大可支持八逻辑线程，相当于八核，但未必能达到真八核的性能效果。
③并行效率
并行效率是指加速比与核数即处理器数的比值，它反映了每个处理器被利用的程度。从图17可知处理器的平均利用率大致为50％左右，四核四线程的计算机处理器利用率略高些。
(2)约束地形多核并行插入算法测试
为了测试约束地物嵌入过程使用并行效果，选用了16978个离散点，先通过这些离散点生成无约束域的D-三角网，然后在其中插入多边形，分别选取多边形数为40、80和120个作为对比，其中测试的80个多边形是在不改变前40个多边形的位置和几何形状的基础上继续添加40个多边形得到的，而测试的120个多边形是不沿用原多边形重新获取的多边形，这些多边形包括非凸多边形，整个过程分为顶点并行嵌入、约束边并行嵌入和内部三角形并行清空。测试计算机环境同上，同样以执行时间、加速比和并行效率作为评价性能的指标具体结果如表1所示。
①约束点地物并行嵌入
表1不同约束点并行嵌入性能
②约束边并行嵌入
表2不同多边形数下边并行嵌入性能
③约束面地物内部三角形并行清空
表3不同多边形数下内部三角形并行清空性能
从以上表可以看出通过并行可以使整个效率明显提升，并且随着线程数的增加，执行时间逐渐缩短，整个过程主要是顶点嵌入和边嵌入最费时，而内部三角形清空占据总时间的极少部分，由40个多边形变为120个时所消耗的时间基本上不变，在此过程使用并行的加速比也是最小的，并行效率一般。
Claims (9)
1.一种单机多核环境下约束地形并行构建方法，其特征在于，包括以下步骤：
步骤1，基于四方边缘结构的CD-TIN数据结构及主要函数设计；
步骤2，基于负载平衡的数据划分策略；
步骤3，基于分治算法的子三角网并行构建；
步骤4，子三角网并行合并；
步骤5，格网索引构建；
步骤6，约束点地物并行插入；
步骤7，约束线地物并行插入；
步骤8，约束面地物并行插入；
所述步骤1的具体步骤为：
CCW代表逆时针方向，CW代表顺时针方向，该数据结构定义如下：
(1)假设初始边为e，则e有一个与其方向相反的伴生边e.twin；e.org为e的起始点坐标；e.dest为e的终点坐标；e.Rot和e.InvRot为e的对偶边；
(2)包含“next”的边是指以逆时针方向围绕邻面或顶点的下一条边；e.Lnext指逆时针围绕e边左面的下一条边，与e边同左面；e.Dnext指逆时针围绕e边终点的下一条边，与e边同终点；e.Onext指逆时针围绕e边起点的下一条边，与e边同起点；e.Rnext指逆时针围绕e边右面的下一条边，与e边同右面；
(3)包含“prev”的边是指以顺时针方向围绕邻面或顶点的下一条边；e.Dprev指顺时针围绕e边终点的下一条边，与e边同终点；e.Rprev指顺时针围绕e边右面的下一条边，与e边同右面；e.Lprev指顺时针围绕e边左面的下一条边，与e边同左面；e.Oprev指顺时针围绕e边起点的下一条边，与e边同起点；
为实现该数据结构，定义了三个数据类：节点类—VertexClass、边类—EdgeClass和三角形类—TriangleClass，以及一个操作类：UtilityClass；VertexClass记录节点三维坐标、关联边和布尔值的约束属性，EdgeClass记录四方边缘结构的各拓扑边、左右三角边形和布尔值的约束属性；UtilityClass中定义了CD-TIN的系列操作函数接口，
其中包含四个主要含数据接口：makeEdge、splice、connected和swap，函数功能描述如下：
(1)void splice(Edge*a,Edge*b)，该函数作用是对a、b边进行连接或者拆分；若将a、b边所在的两个具有共同顶点的独立子三角网合并，该函数将两个子三角网在共同顶点处关联在一起；同时，合并a和b左侧的面；相反，若要在共同顶点处沿a和b边拆分，该函数将在共同顶点处将三角网拆分成两个独立子三角网，同时，分裂a和b左侧的面；
(2)Edge*connected(Edge*a,Edge*b)，该函数作用是连接a的终点，b的起点生成新边c，在此过程中会建立a,b,c之间的四方边缘结构关系,next指边的Lnext，prev指边的Lprev；
(3)void swap(Edge*e)，该函数作用是用来交换对角线边，当用局部优化法LOP检测三角形不满足Delaunay法则时，则需要交换对角线边，将与e边相邻的两三角形组成的四边形另两顶点相连，形成新的边取代e，该边分别以e边起点和终点逆时针方向查找，第一点为新边的起点，第二点为终点。
2.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤2的具体步骤为：
采用坐标排序法作为统一的数据划分策略，(1)用总点数N除于计算线程数K，假设每个线程平均分配点数为m；(2)将离散点集先X轴方向，后Y轴方向进行排序；(3)根据每个线程的点数m，将点号区间[i*m,(i+1)*m-1]分配给每个线程，i为线程序号，i＝0,1,2,…,K，如果存在没有分配完的剩余点号，则全部给最后一个线程，所有线程共享同一份数据，每个线程计算时，根据点号区间从新排序点集中读取对应点。
3.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤3的具体步骤为：以左右合并情况为例，具体实现过程如下：
(1)定义数据结构MaxEdge用于存储线程内部和线程间子三角网凸包的最左上/下和最右上/下边；
(2)定义左右凸包合并时的边指针参数：ldo、ldi、rdo和rdi，其中，ldo—指针为左凸包上以最左点为始点的逆时针方向的边；ldi—指针为左凸包上以最右点为始点的顺时针方向的边；rdo—指针为右凸包上以最右点为始点的逆时针方向的边；rdi—指针为右凸包上以最左点为始点的顺时针方向的边，该四个指针为全局变量，为了解决并行时数据共享冲突问题，在OpenMP并行模块头部设置指针变量为线程私有：#pragma ompthreadprivate(ldo,ldi,rdo,rdi)；
(3)分治算法程序模块，分治算法使用递归方式，将离散点集按子块数量在空间上分成左右均等的两个子集L、R，子块经过递归继续划分为更小子块直到划分数目size小于4，并且最小只会出现点集size为2和3的情况，此时，直接两两相连形成线段或三角形，递归下降过程结束之后，开始递归上升；在上升过程，size最小为4，左右子集L、R剖分形成左右凸包后，将凸包数据以MaxEdge结构方式保存在leftRet、rightRet，从左右凸包中可以获取指针ldi、ldo、rdi、rdo；
(4)多核并行构建程序模块。
4.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤4的具体步骤为：
(1)获取左右凸包下公切线，
以左凸包最右边ldi和右凸包最左边rdi为起始边参数，通过判断各自起点org在对方边的哪侧，并结合数据结构快速找到下公切线；为了快速获取点在哪侧，使用向量叉乘k分量结果进行判断；假设有两个向量和这两个向量有共同点A，对这两个向量作叉乘运算：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mover>
<mrow>
<mi>A</mi>
<mi>B</mi>
</mrow>
<mo>&RightArrow;</mo>
</mover>
<mo>&times;</mo>
<mover>
<mrow>
<mi>A</mi>
<mi>C</mi>
</mrow>
<mo>&RightArrow;</mo>
</mover>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<mi>i</mi>
<mo>,</mo>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
<mi>j</mi>
<mo>,</mo>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mn>2</mn>
</msub>
<mi>i</mi>
<mo>,</mo>
<msub>
<mi>b</mi>
<mn>2</mn>
</msub>
<mi>j</mi>
<mo>,</mo>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfenced open = "|" close = "|">
<mtable>
<mtr>
<mtd>
<mi>i</mi>
</mtd>
<mtd>
<mi>j</mi>
</mtd>
<mtd>
<mi>k</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>a</mi>
<mn>2</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>b</mi>
<mn>2</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<msub>
<mi>b</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mi>i</mi>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<msub>
<mi>c</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<msub>
<mi>a</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mi>j</mi>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<msub>
<mi>b</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
<msub>
<mi>a</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mi>k</mi>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
假设A、B、C三点都为二维平面上的点则c1和c2都为零，那么有：
<mrow>
<mover>
<mrow>
<mi>A</mi>
<mi>B</mi>
</mrow>
<mo>&RightArrow;</mo>
</mover>
<mo>&times;</mo>
<mover>
<mrow>
<mi>A</mi>
<mi>C</mi>
</mrow>
<mo>&RightArrow;</mo>
</mover>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>a</mi>
<mn>1</mn>
</msub>
<msub>
<mi>b</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>b</mi>
<mn>1</mn>
</msub>
<msub>
<mi>a</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mi>k</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
叉乘的k分量具有以下性质：
性质1：若(a1b2-b1a2)＞0,则在的顺时针方向，那么C点在的左边；
性质2：若(a1b2-b1a2)＜0,则在的逆时针方向，那么C点在的右边；
性质3：若(a1b2-b1a2)＝0,则与共线，但可能同向也可能反向，如果C点在A和B的最小外包MBR内，则同向且C点落在线段AB上；根据上述向量叉乘k分量性质，获取下公切线伪代码；
(2)获取最终候选点并连接合并，
假设下公切线记为baseLine，以其终点为起点、左候选邻居点为终点的边记为LEdge；以baseLine共起点、右候选邻居点为终点的边记为REdge，候选点判断过程如下：
①重新获取左右边：LEdge和REdge如果LEdge终点在baseLine右边，则判断baseLine终点、起点及LEdge终点三点构成的外接圆是否包含LEdge.Onext终点，若包含，则LEdge＝LEdge.Onext，重复①；如果REdge终点在baseLine右边，则判断baseLine终点、起点及REdge终点三点构成的外接圆是否包含REdge.Oprev终点，若包含，则REdge＝REdge.Oprev，重复①；
②获取最终候选点：如果LEdge终点不在baseLine右边，或者如果REdge终点在baseLine右边，且LEdge起点、终点及REdge起点三点构成的外接圆包含REdge终点，则由baseLine终点和REdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine；如果以上条件不满足，则从左凸包获取候选点，则由baseLine起点和LEdge终点组成Delaunay三角形的一条边unionLine，然后更新下公切线，令baseLine＝unionLine；最后，基于新下公切线，并回到①，重新获取LEdge、REdge；
③退出算法：如果LEdge和REdge的终点都不在baseLine右边，则此时baseLine为上公切线，算法结束；
(3)相邻子三角网两两递归并行合并，
根据并行构建子三角网的凸包找到相邻子三角网两两合并，并得到新的子三角网及其凸包；然后，通过新子三角网的凸包找到相邻子三角网进行合并，如此下去直到子三角网数为1，则合并结束。
5.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤5的具体步骤为：
以D-TIN中所有三角形最小外包矩形MBR的长作为样本，并对该样本进行排序，得到一组从小到大的新样本，利用正态分布检验算法求出样本的均值μ和方差σ，接着以μ+σ为网格单元行高的初始值，每次增加0.1σ，增加m次，m由r决定，对于那些不服从正态分布要求的数据参照正态分布的“3σ准则”引进一个比例系数r，0<r≤1，r的默认值为[0.90，0.95]，最终以μ+σ+m*0.1σ作为网格单元的行高，从而得到网格划分的行数M，再以D-TIN中所有三角形最小外包矩形MBR的宽作为样本，确定网格划分的列数N。
6.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤6的具体步骤为：
结合格网索引和四方边缘结构，加快点状约束地物影响区域的搜索速度，对正参与约束地物影响区域重构的三角形在内存中标记为：已占用，避免其它地物同时调用此三角形，从而避免算法混乱，计算完成后将参与计算的三角形从内外存删除，其它地物将从最新三角网搜索影响区域；
约束点并行插入时，对不同情形，采取不同并行策略：
(1)并行插入点的影响区域不重叠，则两线程同时执行对影响区域重新构网，完成后删除影响区域三角形，并更新三角网和格网索引；
(2)并行插入点的影响区域存在重叠，则不可以同时进行影响区域重新构网，需要待到约束点地物CP1完成计算后，约束点地物CP2再从新的三角网中，查找影响区域；对上述情形，CP2释放所占线程，进入下一约束点地物的插入，CP2则仍保留在未插入约束点地物集合中。
7.根据权利要求1所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤7的具体步骤为：
(1)通过约束线地物外包矩形找到其所在范围的格网索引；
(2)根据格网索引初步找到约束线地物的影响区域三角形；
(3)从(2)的结果之中，逐一找到约束线地物顶点与相邻点间线段的影响区域三角形，并将找到的三角形标记为：已占用，此时存在以下两种情形：
(ⅰ)如果前一点与后一点的影响区域三角形相邻，则此两点的影响区域三角形也即为此两点间线段的影响区域；
(ⅱ)如果相邻两点影响区域三角形不相邻，则通过四方边缘数据的拓扑关系可以快速找到顶点V2影响区域三角形的相邻三角形C；然后，将线段V2V3与三角形C的其余两条边e1和e2作跨立实验，判断线段V2V3穿越哪条边；其次，根据穿越边找到下一影响区域三角形，如此下去，直到搜索到顶点V3所影响的三角形E为止。
8.根据权利要求7所述的一种单机多核环境下约束地形并行构建方法，其特征在于，约束线地物并行插入时可能存在的情形处理：
(1)约束区域不重叠时，约束线地物CPl1和CPl2各自所在影响区域则调用约束线地物插入算法并行执行局部重新构网；
(2)约束区域重叠时，约束线地物CPl1影响区域先行重新构网；则CPl2所占线程将其先跳过，同时，输入下一约束线地物进行局部重新构网处理；
(3)最后，再从未处理的约束线地物堆栈对未插入的约束线地物进行并行插入，并对(1)、(2)的情形作与(1)、(2)的同样处理，如此往复循环，直到所有约束线地物插入为止。
9.根据权利要求7所述的一种单机多核环境下约束地形并行构建方法，其特征在于，所述步骤8的具体步骤为：
约束面地物的影响区域搜索则分解为两步：面地物边界影响区域搜索和面地物内部影响区域搜索，操作如下：
(1)面地物边界影响区域搜索调用步骤7的算法；
(2)面地物内部影响区域搜索；首先，得到面地物所在范围格网索引单元；然后，根据格网索引单元找到对应三角形，并排除边界影响区域三角形；最后，通过判断其余三角形任一顶点是否在面地物内，如果在，则整个三角形都包含在面地物内。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510299433.6A CN105160706B (zh) | 2015-06-03 | 2015-06-03 | 一种单机多核环境下约束地形并行构建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510299433.6A CN105160706B (zh) | 2015-06-03 | 2015-06-03 | 一种单机多核环境下约束地形并行构建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105160706A CN105160706A (zh) | 2015-12-16 |
| CN105160706B true CN105160706B (zh) | 2017-12-19 |
Family
ID=54801548
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510299433.6A Expired - Fee Related CN105160706B (zh) | 2015-06-03 | 2015-06-03 | 一种单机多核环境下约束地形并行构建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105160706B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107045512B (zh) * | 2016-02-05 | 2020-11-24 | 北京京东尚科信息技术有限公司 | 一种数据交换方法及系统 |
| CN106294985B (zh) * | 2016-08-08 | 2019-06-11 | 福建农林大学 | 一种高效分布式与并行Delaunay三角形构建方法 |
| CN106485766A (zh) * | 2016-10-21 | 2017-03-08 | 西南大学 | 一种约束Delaunay三角网的并行构建方法 |
| CN110533739B (zh) * | 2019-08-15 | 2024-02-23 | 深圳供电局有限公司 | 一种地图网格的成图方法 |
| CN110633149B (zh) * | 2019-09-10 | 2021-06-04 | 中国人民解放军国防科技大学 | 均衡非结构网格单元计算量的并行负载均衡方法 |
| CN112344863B (zh) * | 2020-09-11 | 2022-08-09 | 湖北三江航天江北机械工程有限公司 | 一种自由曲面回转体工件壁厚检测方法 |
| CN112256816A (zh) * | 2020-11-03 | 2021-01-22 | 亿景智联(北京)科技有限公司 | 一种基于分治网格的空间大数据算法 |
| CN114396962A (zh) * | 2022-02-25 | 2022-04-26 | 北京世纪高通科技有限公司 | 可达域的生成方法、装置、设备及存储介质 |
| CN116385878B (zh) * | 2023-03-29 | 2023-10-20 | 阿里巴巴（中国）有限公司 | 道路中心线提取方法、分布式系统、服务器和存储介质 |
| CN116841742B (zh) * | 2023-07-03 | 2024-05-03 | 劳弗尔视觉科技有限公司 | 一种用于计算海量数据的流式处理方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101604018A (zh) * | 2009-07-24 | 2009-12-16 | 中国测绘科学研究院 | 高分辨率遥感影像数据处理方法及其系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5161936B2 (ja) * | 2010-08-11 | 2013-03-13 | 株式会社パスコ | データ解析装置、データ解析方法、及びプログラム |
-
2015
- 2015-06-03 CN CN201510299433.6A patent/CN105160706B/zh not_active Expired - Fee Related
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101604018A (zh) * | 2009-07-24 | 2009-12-16 | 中国测绘科学研究院 | 高分辨率遥感影像数据处理方法及其系统 |
Non-Patent Citations (1)
| Title |
|---|
| 单机多核环境下的TIN地形并行构建关键技术研究;熊证;《中国优秀硕士学位论文全文数据库》;20150131;摘要、正文第13-16,22-24,28,33-40,46-47,62页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105160706A (zh) | 2015-12-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105160706B (zh) | 一种单机多核环境下约束地形并行构建方法 | |
| CN102682103B (zh) | 一种面向海量激光雷达点云模型的三维空间索引方法 | |
| CN113628314B (zh) | 一种虚幻引擎中摄影测量模型的可视化方法、装置和设备 | |
| Puppo et al. | Parallel terrain triangulation | |
| CN112181991B (zh) | 基于快速构建kd树的地球模拟系统网格重映射方法 | |
| CN110147377A (zh) | 大规模空间数据环境下基于二级索引的通用查询算法 | |
| CN113268557B (zh) | 一种适应显示导向型可视化分析的快速的空间索引方法 | |
| CN110175175A (zh) | 一种基于spark的分布式空间二级索引与范围查询算法 | |
| CN111260784A (zh) | 一种城市三维空间网格压缩编码方法、装置及终端设备 | |
| CN106485766A (zh) | 一种约束Delaunay三角网的并行构建方法 | |
| Jing et al. | An improved distributed storage and query for remote sensing data | |
| Su et al. | An adaptive and rapid 3D Delaunay triangulation for randomly distributed point cloud data | |
| CN117036567A (zh) | 一种三维场景模型渲染方法 | |
| CN113419861A (zh) | 一种面向gpu卡群的图遍历混合负载均衡方法 | |
| CN113722415B (zh) | 点云数据的处理方法、装置、电子设备及存储介质 | |
| Dou et al. | A fine-granularity scheduling algorithm for parallel XDraw viewshed analysis | |
| CN113342999B (zh) | 一种基于多层跳序树结构的变分辨率点云简化方法 | |
| Nguyen et al. | B-EagleV: visualization of big point cloud datasets in civil engineering using a distributed computing solution | |
| CN117237503B (zh) | 一种地理要素数据加速渲染及装置 | |
| WO2021175110A1 (zh) | 区块链方法及系统、电子设备及计算机可读存储介质 | |
| CN110008597A (zh) | 基于并行计算框架的建筑信息模型三角剖分方法及装置 | |
| CN117608476A (zh) | 矢量数据分块存储方法、装置、电子设备及介质 | |
| Dou et al. | An equal‐area triangulated partition method for parallel Xdraw viewshed analysis | |
| Jones et al. | The implicit triangulated irregular network and multiscale spatial databases | |
| CN115409907A (zh) | 点云模型的压缩方法、装置及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20171219 |