宗 真,袁林旺,2,罗 文,俞肇元,2,胡 勇,3
1.南京师范大学 虚拟地理环境教育部重点实验室,江苏 南京210023;2.南京师范大学 江苏省大规模复杂系统数值模拟重点实验室,江苏南京210023;3.南京师范大学计算机科学与技术学院,江苏南京210023
三角网表达及求交是三维表面模型建模的基础,也是地质工程、科学可视化以及模型分析的关键算法之一[1,2]。GIS研究对象向三维以及时空维的扩展,不仅导致了数据量的激增[3],也导致了对象表达结构的复杂性。GIS对象表达结构的复杂性,不仅使得在空间关系计算过程中对象类型及求交条件判断模式的复杂,也增加了对诸如空间、边界等复杂的空间约束[4]判断的难度,因而极大地影响了现有矢量GIS的效率。伴随着GIS研究对象向三维以及时空维的扩展,三角网求交的复杂度也随之增加。在三维动态场景中,三角网中对象的维度、形态和结构均存在显著的变化[5]。基于欧氏几何或计算几何发展而来的三角网求交算法需要对其相对位置关系进行预判断,且在处理不同维度三角网时,需要对线、面进行分别处理,导致了其对动态变化的三角网求交支撑不足[6-7]。以文献[8]算法为代表的标量法在平面距离计算及交线位置判断上的计算相对复杂,在大规模环境中时间耗费较多。文献[9]的矢量法利用行列式符号的几何意义来判断三角形的关系而后求交,计算过程中必须根据相关情况进行复杂的判断,导致程序结构复杂,适应性不强。
几何代数是以代数语言进行几何对象表达、构造与运算的数学系统,被认为是有效连接几何与代数,实现几何计算的数学语言[10]。共形几何代数(conformal geometric algebra,CGA)基于协变视角引入了基准原点和无穷远点,使得分别基于外积和内积的形体表达和几何度量都具有明确的几何意义[11]。通过定义几何代数运算空间,可实现高维空间中不依赖于坐标的几何计算,前期研究表明几何代数可用于多维空间对象的统一表达、建模与复杂对象间的空间关系分析[12-14],且基于其设计的算法具有较强的高维拓展性[15]。本文针对三角网求交中三角形空间关系判定及点、线、面等不同几何对象统一求交问题,从代数表达与关系运算相结合的角度,基于几何代数中交并关系运算的meet算子构建空间关系判断算子,为三角网表达结构的动态构建及各维度组分的统一求解提供基础,进而实现多维统一三角网求交的几何代数算法。
基于内积和外积进行代数空间构造及几何对象表达,基于几何代数多重向量(multivector)可实现多维几何对象的整合表达。在几何代数中,基于外积表达的几何对象维度与Grassmann维度具有一致性,可以实现不同维度对象的构造与相互转化[10]。多重向量则为不同类型、不同维度几何对象的描述和存储提供了高维可拓展性。几何代数自身的坐标无关性与维度无关性使其所表达的几何特征是几何对象自身内蕴的几何特征,计算获得的几何关系也是各几何对象间不依赖坐标的相对关系[11,16]。这为几何对象的自适应表达、关系判断与关系计算提供了有利前提。CGA表达中包含有如方向、位置、值大小等几何语义信息,可通过直观的算子运算直接判断对象间相离、相交、包含等拓扑关系。表1给出了CGA中常用blade的分类及其所包含的几何特征及其计算方法[15]。
表1 blade分类及其几何特征Tab.1 Blades and their geometric characteristics
根据CGA中blade的表达,并结合前期研究中所构造的基于CGA的三维空间数据模型[12],可以构造三角网的表达模式如图1所示。首先构建多重向量集合MVi。在MVi集合中任意一个对象均为构成该三角网的特定三角形Ti,采用外积构造点、线、面对象,并利用其边界范围将其约束为点、线段与三角形。利用CGA中几何对象维度与Grassmann维度的一致性,实现点、线段与三角形的层次组织。在三角网的表达中,点、线段与三角形均用blade表达,且依据其层次关系组合形成最终的三角网多重向量表达。不仅在维度结构上具有简明性,也为三角网的统一表达与空间关系计算提供了独立、完备的数据组织与数据结构基础。
几何代数算子的多维统一性使其可以直接用于处理多维对象,而无需根据对象维度进行分别处理[17],从而可构建简明、清晰、多维统一的三角网快速求交算法。基于CGA几何对象表达的多维统一性和meet算子构建三角网自适应求交算法(图1),通过对表征两三角网的多重向量直接应用meet算子实现对多重向量中所包含的各几何对象的分别求交。对于基本几何对象间交并关系的求解,如点、线、平面、圆环、球面等一般blade,可直接基于meet算子进行运算;对于求解含固定边界几何对象间交并关系,如多边形、多面体对象,除了先基于meet算子进行求交判断外,还需借助对象表达中所蕴含的丰富的几何属性特征,并定义一系列具有统一结构的关系判断算子,去除无效结果。最后将meet求解结果重新投影回欧氏空间,可设计相应的语法解析器对运算结果加以解析,并转换成对应几何对象,形成三角网求交的最终算法流程。
图1 基于几何代数三角网求交流程Fig.1 Calculation flow of triangulation intersection based on geometric algebra
不同几何对象的相交关系会随着对象类型和位置关系的变化而变化。传统欧氏几何中需要对不同维度对象各种相交情况分别进行处理,导致时变对象的求交计算较为复杂,对动态空间关系计算的支撑不足。在几何代数中,不同几何对象间所共有的最大子空间可以通过meet算子获得。由于meet运算的多维统一性,使其可以同时处理点、线、面以及其他各类blade表达的几何对象间的相交运算[17]。meet运算的结果为多重向量,且其运算结果的几何对象类型仅与参与meet运算的几何对象相关,而与维度和坐标系统无关。meet算子计算结果所具有的自适应性,有助于实现利用单一算子实现不同几何对象相交关系的统一求解。
给定任意两bladeA和B,两者间的meet运算M=A∩B满足下式
引入可容纳A、B的最小子空间伪标量IAB,则meet算子可利用A、B间的对偶运算求解
其一般形式为
式中,β为标量。对于blade对象间的求交,可通过meet算子M的符号和它的平方M2的符号来作初步的相交关系判断。Eduardo等基于meet算子详细推导了线-线,线-面,面-面之间相交关系的计算方法与判断规则[18]:当M2>0时,A与B至少相交于两点;当M2=0时,A与B至少相切于一点;当M2<0时,A与B不交。各类型几何对象基于meet算子的相交关系的判断如表2。meet算子表达形式上的统一性及其计算结果的自适应性使其可适用于动态场景中变化对象的自适应求交,并能保证动态对象条件下,不同维度对象间空间关系计算的统一性与一致性。
表2 基于meet算子的基本几何对象相交关系求解Tab.2 The solutions of intersection relations among basic geometric objects based on meet operators
meet算子对于求解A=a1∧a2∧…∧ak形式的blade之间的交集简单有效,但对于由多个blade的加和构成的k-blade间的meet运算,需要基于join(并)和meet的关系对其进行拆分和变换[10,19]。鉴于 meet运算的可拆分性,可以通过将blade对象分解成线性无关的几个支撑向量以简化运算。对于k-bladeB,其因子分解就是要找到k个线性不相关的支撑向量fi,使得B=βf1∧f2∧…∧fk。使用投影算子qi=(pi」B-1)」B搜寻B中k个相互独立的分解因子。设计meet算子快速求解算法,其基本流程为:
(1)先将“测试向量”pi投影到二维平面B上;
(2)用pi」B-1计算出pi在B中投影的正交补余;
(3)据公式qi=(pi」B-1)」B求解pi」B-1顺时针旋转90度的正交投影,qi即为所求支撑向量;
(4)循环直至找到k个支撑向量。
由于对象分解过程具有可加和性和可逆性,且分解得到的qi线性不相关,可用于对象间meet的快速运算。为快速求得meet算子结果的维度,可引入delta积运算,其几何意义为两对象间的并集减交集[19-20],存在
给定初始交积M=1、并积J=In,对于kbladeA与B,计算其delta积在空间IAB中的补集S=(AΔB)*,并将S分解为线性无关的k个部分Si;对于每个Si计算其投影pi,如果pi非零,则通过M=M∧pi扩充meet,直到M的维度与所求得meet的维度相等。
在复杂三角网对象求交中,对象间空间关系判断是进行求交对象筛选和处理的基础。在三角形求交过程中,主要需考虑不同线段、线段与三角形以及点与三角形之间的关系判断。为此,基于三角形及其相关几何对象的CGA表达,通过构造相应的判断算子进行上述对象的空间关系判断。其计算过程如图2所示。
图2 线段、三角形、点的位置关系判断Fig.2 Location relations judgment between segments,triangles and points
计算线段R1及其端点p1与线段R2两端点所构成向量间的外积2-bladeB1=v3∧v1、B2=v3∧v2,同理可求得由线段R2及R1端点所构成向量间外积B3、B4。构建平面上线段关系判断算子:Oseg-seg=(dir(B1)= =dir(B2))‖(dir(B3)==dir(B4)),当Oseg-seg为true时,线段不交,否则线段相交。
构建由线段R的顶点q1到三角形的三个顶点所构成向量的外积3-bladeTr1=v1∧v2∧v3,同理可得到由顶点q2与三角形顶点所构成向量的外积Tr2。构建线段与三角形的关系判断算子:Oseg-tri=(dir(Tr1)==dir(Tr2)),当Oseg-tri为true时,二者不交,否则线段与三角形相交。
构建由三角形各顶点指向待求点的向量,按一定的方向计算三角形边同各向量的外积:B1=v12∧v10、B2=v23∧v20、B3=v31∧v30。构建点与三角形的关系判断算子:Opt-tri=(dir(B1)≥0&&dir(B2)≥0&&dir(B3)≥0),当Opt-tri为true时,点在三角形内部或边界,否则点在三角形外部。
不失一般性,据文献[17]构建的基于几何代数的对象表达模型,给定两三角形T1、T2,其共形几何代数表达为
上述三角形的多重向量表达中,已同时包含有各三角形中线段和面对象的几何代数表达,其meet的求解过程为
式中,S1∩l21+S1∩l22+S1∩l23由于与S2∩l11+S2∩l12+S2∩l13具有等同性,二者取一即可,并根据三角形与直线的构建关系可将上式改写成如下形式
式中,“*”仅为一种连接关系,不参与运算,表示当“*”左侧满足一定要求时,才进行其右侧运算,即在求解过程中,可利用面—面、线—面、线—线对象meet算子结果的符号对其相交关系进行初步判定,从而降低算法的计算复杂度。
基于上述形式化表达,利用三角形求交的空间关系判定方法,构造相应的判定算子对其空间关系进行判定,进而获得三角形在共面、异面等情况下,其中各要素的相交情况,然后将求交结果逆向回溯,由交点构建线段,即为所求交线结果。
其算法流程如图3所示。首先利用meet算子进行三角形所在平面间相交关系的预判断,即M=(T1∩T2)。由于当M2<0时,T1与T2不交,因此,可筛选出所有M2<0(即不满足相交条件)的三角形,对所有可能满足相交条件的三角形进行空间关系的判断。构造由三角形T1任一顶点A1、B1、C1到三角形T2顶点A2、B2、C2所构成向量的外积3-bladeVol=A2∧B2∧C2值是否为零,筛选出共面相交的情况。对于异面相交的三角形,计算三角形各边与另一三角形平面的meet。共面相交的三角形,则需要遍历计算两三角形边与边的meet,并通过交点、交线与三角形空间位置关系判断算子筛选出落在三角形内部或边界上的对象。
图3 三角形求交流程示意Fig.3 The triangle intersection process
对上述算法流程的求解步骤为:
输入:待求交的两三角网T1、T2。
输出:两三角网交线、交点等几何对象Vout。
(1)对T1、T2进行几何代数编码,形成多重向量表达MVT1、MVT2;
(2)顺序提取MVT1、MVT2中三角形Π1、Π2,对其所在平面求meet,并进行三角形相交与共面的判断;
(3)对Π1、Π2中各边求交,判断其相交情况;
(4)对所求交点、交线与三角形关系进行判断,当交点、交线在三角形内时,保留结果;
(5)依次遍历三角网中各三角形,直至所有三角形求交完成;
(6)解析计算结果,转化为对应欧氏几何对象输出或可视化。
以南极冰盖时空演化[21]的三维表面模型求交为例进行算法验证,选用33 800Ka BP、33 550Ka BP、33 400Ka BP、33 300Ka BP、33 200Ka BP 5个时间节点的南极冰盖三角网曲面模拟数据作为案例数据,数据结点及三角网个数信息见表3。利用几何代数对三角形的统一表达与运算,结合相关的运算与判断算子,设计三角网求交的数据存储结构和算法。首先对三角网格曲面数据进行共形空间转换与维度映射,实现三角格网的几何代数表达与存储;然后利用空间三角形的几何代数求交方法,基于meet算子和关系判断算子Oseg-seg、Oseg-tri、Opt-tri进行相交结果的筛选与位置关系判断,并获取相应的交点与交线段,最后求得两个时期冰盖表面的曲面交线。结果显示本文算法提取出的交线能够清晰反映某一时刻冰盖曲面相对原始冰盖变化的空间结构与位置,可为地质地形建模与分析提供可行的思路与方法(图4)。
表3 三角网数据的基本信息Tab.3 General information of the triangulation data
图4 本文算法求得的南极冰盖表面三角网交线结果Fig.4 The triangulation intersection of Antarctic ice sheet based on CGA algorithm
对上述三角网交线结果进行分析易知其中含有大量的重复交线,这些冗余交线会影响后期交线拓扑的构造效率。冗余交线主要产生于如下3种情况:①退化交线(交点),例如当一三角形顶点与另一三角形相切(接)时,交线段退化为交点,此交点必为其邻接三角形交线的顶点;②切交线,例如当三角形一边切另一三角形面时,该三角形的邻接三角形必求解出相同的切交线;③接交线,例如当两个三角形的边部分或整体重合时,该求解结果将会在其邻接三角形中重复出现。情况①中的点由于不参与拓扑构建可直接删除,情况②与情况③中需要保留一条交线用于构造交线拓扑。基于上述考虑,可分别为线线关系判断算子Oseg-seg和线面关系判断算子Oseg-tri添加判断条件Bi==0,i=1,…,4和Oseg-segi==true,i=1,2,3,标记出线线相接和线面相接的情况,对退化交线加以筛除,式中Bi为线段端点构成的2-blade,segi指代三角形tri的第i条边。
分别运用Möller标量法、Guigue&Devillers矢量法和本文算法(CGA算法)按时间序对上述三角网数据两两求交,统计出交线数目及各算法运行时间,对本文算法的准确性和效率加以验证。最终交线的统计结果如表4所示。由于基于CGA的三角形求交算法去除了一定量的冗余交线,其求得的交线总数均小于另两算法,减去冗余交线后,三算法求得的有效交线数目一致,表明了本算法求解结果的准确性与有效性。从计算效率上看(图5),本方法较 Möller算法略快,但慢于Guigue&Devillers算法。其主要原因在于,本文实现的CGA算法是将几何对象由三维欧氏空间投影到五维空间后,使用数值运算而非逻辑位运算进行几何代数计算的求解,导致了内积、外积以及meet算子等基本几何代数运算仍然依赖于加、减、乘、除等数值运算。通过基于逻辑位运算优化的几何代数引擎的利用[22],可进一步大幅提升运算效率。而几何代数对象表达与运算的统一性,及其算法的结构简明性与自适应性,也为基于几何代数的并行计算方法提供了便利[23]。
表4 三角网相交结果统计Tab.4 The statistical results of the triangulation intersection
图5 三角形求交算法效率对比Fig.5 Efficiency comparisons of the three triangulation intersection algorithms
三角形是三角网格曲面、四面体格网、单纯形剖分等三维矢量模型构建的基本单元,三角形求交运算是上述模型空间拓扑关系运算的基础。空间对象的相交、距离关系的判断是空间拓扑分析的基础,几何代数通过定义相应的运算空间实现多维对象的统一表达与计算。本文利用几何代数多维统一表达与运算的优势,设计了三角形等基本对象的统一表达;基于meet算子维度统一与对象无关特性,构建统一的三角形求交运算流程;并设计相应的关系判断算子实现计算过程中对象的筛选与判断。
在三角形的求交过程中存在三角形退化为线段、三角形相接、三角形相切等多种临界与奇异情况,在所选的两类传统对比算法中,均存在人为设定的相交准则,使得计算流程存在一定的定制化,而对不同维度、不同对象缺乏自适应性。本文算法基于统一的多重向量表达及灵活的关系判断算子可根据应用需要合理配制相交准则,对动态场景下的时变对象间相交关系的计算具有更强的自适应性。基于南极冰盖曲面网格数据的案例表明,本文提出的基于几何代数的三角形快速求交方法,可有效支撑三角网的求交运算,且逻辑结构清晰,运算简洁高效,可为其他复杂几何对象间建模表达与空间分析统一求解提供借鉴。
支持复杂对象与现象的动态化表达、建模与模拟是GIS发展的必然趋势[24]。建立几何—拓扑—语义一致的且可有效支撑地学分析的地理格网与空间分析模型是实现上述目标的重要途径[25]。基于几何代数的三角网参数化表达在融入对象的几何、拓扑特征的同时可直接支撑几何运算,并具有对象和维度的自适应性。除meet算子外,大部分的几何代数算子均具有类似的多维统一和自适应特性,因此,基于几何代数算子设计对象和维度无关的空间关系计算算法,将可为多维统一分析与计算提供了新的理论与方法参考。未来主要的研究拓展包括:①研究可支撑复杂对象的几何代数表达方法,实现地理场景中不规则边界形体数据的语义结构与拓扑特性分析;②基于几何代数算法流程的一致性与逻辑结构的简洁性,构建并行化分析算法以支撑大规模复杂场景及海量数据分析;③将算法向复杂、高维对象扩展,设计几何与语义统一的空间关系计算算法,建立具有以维度融合为特征的GIS空间分析体系。
[1] AI Tinghua,LIU Yaolin,HUANG Yafeng.The Hierarchical Watershed Partitioning and Generalization of River Network[J].Acta Geodaetica et Cartographica Sinica,2007,36(2):231-236,243.(艾廷华,刘耀林,黄亚锋.河网汇水区域的层次化剖分与地图综合[J].测绘学报,2007,36(2):231-236,243.)
[2] ZHANG Yao,FAN Hong,HUANG Wang.The Method of Generating Contour Tree Based on Contour Delaunay Triangulation[J].Acta Geodaetica et Cartographica Sinica,2012,41(3):461-467.(张尧,樊红,黄旺.基于Delaunay三角网的等高线树生成方法[J].测绘学报,2012,41(3):461-467.)
[3] WANG Yuan,LIU Jianyong,JIANG Nan,et al.A Dynamic Triangulation Algorithm for the View-dependent and Realtime LoD Model of Terrain[J].Acta Geodaetica et Cartographica Sinica,2003,32(1):47-52.(王源,刘建永,江南,等.视点相关实时LoD地形模型动态构网算法[J].测绘学报,2003,32(1):47-52.)
[4] LIU Qiliang,DENG Min,SHI Yan,et al.A Novel Spatial Clustering Method Based on Multi-constraints[J].Acta Geodaetica et Cartographica Sinica,2011,40(4):509-516.(刘启亮,邓敏,石岩,等.一种基于多约束的空间聚类方法[J].测绘学报,2011,40(4):509-516.)
[5] SUN Dianzhu,LI Xincheng,TIAN Zhongchao,et al.Accelerated Boolean Operations on Triangular Mesh Models Based on Dynamic Spatial Indexing[J].Journal of Computer-aided Design & Computer Graphics,2009,21(9):1232-1237.(孙殿柱,李心成,田中朝,等.基于动态空间索引结构的三角网格模型布尔运算[J].计算机辅助设计与图形学学报,2009,21(9):1232-1237.)
[6] JIANG Qianping,TANG Jie,YUAN Chunfeng.Fast Triangle Mesh Surface Intersection Algorithm Based on Uniform Grid[J].Computer Engineering,2008,34(21):172-174.(蒋钱平,唐杰,袁春风,等.基于平均单元格的三角网曲面快速求交算法[J].计算机工,2008,34(21):172-174.)
[7] YIN Changlin,YU Dingquan.Rapid Topological Searchingbased Intersection Algorithm of Triangulated Networks[J].Computer Engineering and Applications,2006,42(36):209-211.(尹长林,喻定权.一种基于拓扑搜索的三角网求交算法[J].计算机工程与应用,2006,42(36):209-211.)
[8] MÖLLER T.A Fast Triangle-triangle Intersection Test[J].Journal of Graphics Tools,1997,2(2):25-30.
[9] GUIGUE P,DEVILLERS O.Fast and Robust Triangle:Triangle Overlap Test Using Orientation Predicates[J].Journal of Graphics Tools,2003,8(1):25-42.
[10] DORST L,FONTIJNE D,MANN S.Geometric Algebra for Computer Science[C]∥ An Object-oriented Approach to Geometry.San Fransisco:Morgan Kaufmann Publishers,2007.
[11] LI Hongbo.Conformal Geometric Algebra:A New Framework for Computational Geometry[J].Journal of Computer Aided Design & Computer Graphics,2005,17(11):2383-2393.(李洪波.共形几何代数:几何代数的新理论和计算框架[J].计算机辅助设计与图形学学报,2005,17(11):2383-2393.)
[12] YUAN L W,YU Z Y,CHEN S F,et al.CAUSTA:Clifford Algebra-based Unified Spatio-temporal Analysis [J].Transactions in GIS,2010,14(s1):59-83.
[13] LASENBY J.New Geometric Methods for Computer Vision:An Application to Structure and Motion Estimation[J].Inter-national Journal of Computer Vision,1998,26(3):191-213.
[14] YUAN L W,YU Z Y,LUO W,et al.A 3DGIS Spatial Data Model Based on Conformal Geometric Algebra[J].Science China:Earth Sciences,2011,54(1):101-112.
[15] YI Lin,YUAN Linwang,Yu Zhaoyuan,et al.Cliffrod Algebra-based Voronoi Algorithm[J].Geography and Geo-Information Science,2011,27(5):37-41.(易琳,袁林旺,俞肇元,等.Voronoi生成的Clifford代数实现方法[J].地理与地理信息科学,2011,27(5):37-41.)
[16] PERWASS C.Geometric Algebra with Applications in Engineering[M].Heidelberg:Springer-Verlag,2009.
[17] YUAN L W,LüG N,LUO W,et al.Geometric Algebra Method for Multidimensionally:Unified GIS Computation[J].China Science Bulletin,2012,57(7):802-811.
[18] EDUARDO R.Operaciones de Cómputo Gráfico en el Espacio Geométrico Conforme 5DUsando GPU[D].Caracas:Universidad Simón Bolívar,2011.
[19] FONTIJNE D,DORST L.Efficient Algorithms for Factorization and Join of Blades[C]∥Proceedings of the 3rd International Conference on Applications of Geometric Algebras in Computer Science and Engineering(AGACSE 2008).London:Grimma,2010.
[20] BOUMA T A,DORST L,PIJLS H G J.Geometric Algebra for Subspace Operations[J].Acta Applicandae Mathematicae,2002,73:285-300.
[21] DECONTO R.POLLARD D.Rapid Cenozoic Glaciation of Antarctica Induced by Declining Atmospheric CO2[J].Nature,2003,421:245-249.
[22] FRANCHINI S,GENTILE A,SORBELLO F,et al.Fixedsize Quadruples for a New,Hardware-oriented Representation of the 4DClifford Algebra[J].Advances in Applied Clifford Algebras,2011,21(2):315-340.
[23] WAREHAM R J.Computer Graphics Using Conformal Geometric Algebra[D].Cambridge:University of Cambridge,2006.
[24] GOODCHILD M F,GUO H,ANNONI A,et al.Nextgeneration Digital Earth[J].Proceedings of the National Academy of Sciences,2012,109(28):11088-11094.
[25] CAMELLI F,LIEN J M,SHEN D,et al.Generating Seamless Surfaces for Transport and Dispersion Modeling in GIS[J].Geoinformatica,2012,16(2):307-327.