一种基于六边形网格的球面体可视化算法*

2021-01-08张继凯张然赵君顾兰君

内蒙古科技大学学报 2020年4期

张继凯，张然，赵君，顾兰君

(内蒙古科技大学信息工程学院，内蒙古包头 014010)

近年来，随着高分卫星的应用和遥感技术的飞速发展，地理信息相关的遥感数据量呈爆发式增长.因遥感数据的地理跨度较大，地表弧度带来的影响无法忽略不计，导致传统体可视化中的立方体法无法适用；另外，传统的全球范围数据可视化方式多采用经纬度网格进行数据存储，这就造成数据分布在全球各处密度不同，不仅会造成各纬度的可视化细节表现差异，更是会带来严重的极点问题和数据存储冗余现象.国内外研究人员对此提出过多种不同的球面网格剖分方式，其中六边形离散网格因其良好的数据分布特性，成为地理信息领域和数据可视化领域近年来的一个研究热点.

本文在对球面进行体可视化的操作过程中将六边形离散网格引入其中，使得球面各处得数据密度可以保持一致.同时，为了简化光线投射算法的插值过程，设计一种基于三角形单元的编码方式，以纹理作为其存储载体，使光线投射过程的数据插值点更容易查找，更适合硬件加速实现.

1 相关工作

随着硬件计算能力地提升和图形显示手段地发展，对海量数据进行可视化成为数据分析和知识发现的有效手段之一，而其中三维遥感数据的体可视化是信息含量最大难度最高的可视化形式.

直接体可视化技术经过多年的发展，已经产生诸如抛雪球算法[1，2]、频域算法[3]、错切-变形算法[4]以及光线投射算法[5]等多类算法.其中，因算法原理相对简单，可视化效果相对更符合实际数据特征等优势，光线投射算法成为当前体可视化中应用最为广泛的一类算法.

在球面离散网格的研究上，赵学胜等[6，7]于2009年提出一种基于菱形块四叉树节点的邻近搜索算法，使用三角形单元作为顶层节点，而其余各层使用四边形单元且层次越低越接近于矩形，但因各纬度带的剖分差异性造成复杂的拓扑结构.针对六边形离散网格，贲进、童晓冲和张永生等人进行了一系列探索[8-11]，他们通过分析各种六边形剖分方法，实现一种新的全球六边形离散网格剖分及编码算法，并证明了其在计算几何领域有较好地应用潜力[12,13].但因该方法的编码规则比较复杂且涉及大量的归属性判断，使得很难将其应用于数据可视化尤其是三维数据场可视化中.

XIE J等[15]于2013年提出了基于六边形网格的光线投射算法，才使得球面离散网络应用于体可视化领域的探索有所突破.该算法通过采用六边形网格实现了直接体绘制，但是由于需提前按照六边形网格相互之间的对应关系将其顶点、边以及网格单元全都连接起来，然后实时的依据该相互之间的对应关系检索相邻的网格单元，从而实现光线投射，但该算法复杂程度高.为了提高渲染速度，他们之后还研究使用多GPU并行处理的方式对全球体数据进行体可视化[16].

2 球面六边形网格数据预处理

2.1 数据组织方式

通过合理的方式构造和储存六边形网格信息，从而应用于光线投射算法.图1所示为六边形网格的构造组成，六边形网格单元的每个顶点使用灰色数字进行编码，六边形的中心则用黑色圆点表示，并用黑色数字对其中心进行编码，由图可知每个六边形网格的中心都与其6个顶点相对应，而每个六边形的顶点又与其周围相邻的3个六边形中心相对应，网格数据根据该结构进行构造.

为了满足不同分辨率数据下的可视化需求，依据文献[15]提出的球面离散网格生成算法对不同分层的六边形网格数据进行构造生成.对于每个六边形网格单元，使用XML的文件形式构造保存网格数据，主要包含单元网格编码、六边形网格各个顶点及其中心点的经度坐标数据与纬度坐标数据.表1展示了不同分层的网格构造过程，其中在当前分辨率下，每个六边形网格单元的面积用表示，地球的表面积使用表示，网格分层层数用N表示.

表1 网格剖分说明表

2.2 数据项处理

在对体数据进行可视化时，若通过具有相同分辨率的第一层网格中心为原点按不同步长向径向进行缩放，用同样的层次数且分辨率也一致的六边形网格作为数据载体，以充分体现体数据的多层次特性，即可得到具有不同高度但分辨率相同的多层次网格.

六边形单元都是由6个顶点和1个中心点组成，将体数据与六边形网格互相关联时，可以在顶点和中心点间进行选择.若选择前者，则图2所示的每个黑色网格代表一个“数据单元”，任一层六边形单元在位置上都与其他层单元在径向上呈对应关系，且相邻2层的“数据单元”可以形成1个六棱柱结构，如图3(a)所示；若选择中心点绑定，每3个相邻六边形的中心点相连得到1个三角形“数据单元”，图2的红色4虚线所示，相邻两层的三角形“数据单元”之间可构成三棱柱，如图3(b)示.

基于六边形网格的光线投射插值算法会涉及2步耗时量的工作：一是确定采样点所在的棱柱，二是计算采样点的数据值.前者需要获取采样点的经纬度以及高度信息，后者需要计算采样点在棱柱内部所处的位置，并根据棱柱顶点数据进行插值计算.鉴于中心绑定的计算量比顶点绑定要少，因此选择六边形中心绑定数据.在减少计算时间的同时，为进一步减少存储空间，在纵向上利用棱柱体形成的各层数据的一一映射关系，在横向上利用每个六边形中心与一系列处于同一经纬度的数据形成的关联关系，可得到其他层的数据，所以只需存储最底层网格数据.

关于空间某点所属三棱柱的确定，首先根据三棱柱的投影三角形与坐标点的经度、维度进行关联，得到该点底部三角形的位置，而后通过点的高度信息向上正方向映射从而确定该点所在网格单元的层级，得到该点所属三棱柱的信息.假设在各种条件都满足的条件下，三角形的网格单元与经度信息和纬度信息理应具有线性转换关系，也就是说即使经度值与纬度值具有非常高的精度，也都可以通过相关的转换或检索计算，最终获得相应的三角形网格单元.利用体数据分布的离散性，构造一张精度极高的“检索表”，该表可以对经纬度到三角形单元的映射进行详细记录，进而有效地减少巨大的计算量，如对空间中某点所属的三角形单元进行检索时，通过“检索表”快速缩小查找范围，再经少量计算即可完成转换.该“检索表”的原理如图4(a)所示，其具有个大小均匀的单元，为最大经度，最小经度，为最大纬度，为最小纬度，每个单元都一个有独一无二的经纬度范围，并且固定经纬度范围的单元中存放了所属的三角形单元编号.在查找某一经纬度所属三角形单元时，只需查找其所在单元格即可.

如图4(a)所示，图中所属三角形的边缘处用蓝色网格单元表示，当划分该单元格所属的三角形时，由于其或许横跨多个三角形，所以存储于该蓝色网格单元中的编码信息也许存在稍许误差.对此，如图4(b)所示，将单元格的密度增大，减小误差，从而也使得处于三角形边缘的单元格区域大大减小.然而，如果单纯的增大划分单元格的密度，则会占用“检索表”大量的存储空间，由此需要找到一个合适划分的密度，既不会让误差增大，也无需占用大量存储空间.通过多次对比实验表明，当将“检索表”的单元格精度设置为(单位：度)时，最终可视化呈现的效果最佳，的计算公式如下：

(1)

式中：dc为每2个相邻的六边形中心的距离；β为纬度值.式(1)显示纬度与划分的网格精度成反比.

“检索表”初始化之前必须对其精度进行确认，而每张“检索表”的单元格的标准均设置为其中心点的经纬度坐标位置，因为已知每个三角形单元的所有顶点的坐标，由此可运用重心法计算得到各个三角形和每个“检索表”单元的从属关系，每个单元存储属于其三角形的编码.为了在可视化进行中可以快速读取，通过二进制或文本的结构形式将该表按照“行优先”的存储方式置于外存中，该表的生成均在预处理阶段完成.

重心法的基本原理是:已知三角形的所有顶点均处于同一个平面，若将该三角形的任一顶点设为原点坐标，那么由该点指向其他2点的方向即为2个坐标轴的方向.如图5所示，将点A定为原点，则AB方向为坐标轴V的正坐标方向A方向为坐标轴U的正坐标方向.

因此，对于同一平面内的任意一点，均可通过式(2)表示得出：

p=A+(C-A)·x+(B-A)·y.

(2)

可见，如果x或y为负数，则点是向U或V轴的负方向移动，就处于三角形外部，因此为了保证点处于三角形内部，需要同时满足如下条件：x≥0，y≥0且x+y≤1.其中，当x=0且y=0时，p点与A点重合；同理，当x=0且y=1时，p点B与点重合；当x=1且y=0时，p点与C点重合.

为了求取已知直角坐标的点在三角形中的UV坐标，对式(2)做如下变换：

p-A=(C-A)·x+(B-A)·y.

(3)

令向量v0=C-A，v1=B-A，v2=p-A，则v2=v0·x+v1·y，有

v2·v0=(v0·x+v1·y)·v0

v2·v1=(v0·x+v1·y)·v1，

(4)

对其求解可得：

(5)

所以，将三角形网格的单元中心点与其顶点的位置坐标代入式(5)，若和值是满足条件，那么可以得出三角形与单元格之间满足相对应的从属关系.

3 基于六边形网格的体绘制

本文使用光线投射算法对体数据进行可视化的操作过程中，各光线的积分、插值等运算没有任何融合等现象存在.由于GPU有着强大的并行处理能力，利用该方法进行计算，并且每条光线独立积分运算.算法的大致流程图如图6所示.

球面光线投射算法的具体步骤如下：

3.1 球面光线插值

传统的光线投射算法将数据分别存放于六面体的6个面中，且所有数据均通过长、宽、高3个维度相互对应.当有屏幕像素随着视线穿过六面体时，六面体的表面就留下了2个交点，分别为输入点与输出点，如图7(a)所示，射线与六面体分别相交于和.光线投射算法的基本原理为，均匀的在和之间取若干个数据，通过颜色叠加和变化等运算获得屏幕上光线的颜色值.

在全球数据体可视化中，体数据范围一般分布在一定高度的球壳中，坐标系也从直角坐标系转换成球面坐标系，所以数据维度也就不存在上述所说的一一对应的关系，数据被压缩到地球表面的空心球壳中，光线在该球壳体中的投射过程如图7(b)所示.本文采用六边形网格作为数据载体，由1.2节获取光线中的采样点数据可知，求得各采样点的经纬高坐标才是关键，由于光线只在直角坐标系椽笔，所以均匀采样也要在该坐标系进行，因此球面光线投射中采样点经纬高坐标的计算步骤如下：

Step1 利用几何着色器和渲染到纹理技术将代理几何体前后面的坐标信息保存成纹理形式，三维位置信息中的XYZ分别保存至纹理的RGB通道中.

Step2 在Step1中的2个纹理中，按对应坐标查找光线的2个交点，将其坐标分别设为和.

Step3 按照设定步长从前到后选择采样点，假设采样步长为，对于第个采样点，其空间坐标计算式(6)所示，其中length()为求取向量长度的函数.

(6)

Step4 利用式(7)将采样点坐标从直角坐标(X，Y，Z)转化为地球球面经纬高坐标(lon，lat，alt).

lon=arctan(Y/X)

lat=arctan(U/V)

alt=p/cos(lat)-N.

(7)

式中：P，N，U，V的计算公式如下：

(8)

式中：re，rp分别为地球的赤道半径和2极半径，ds，theta的计算公式如下：

ds=(re2-rp2)/rp2

theta=arctan(Z·re/(p·rp)) .

(9)

2.2 采样点数据计算

通过在数据预处理阶段生成的三角网格编码的“检索表”，可以从采样点的经纬度直接得到相应的三角网格编码TCn，然后，根据1.1中的XML文件，查询到相应的3个六边形编码HCn1，HCn2和HCn3，并找到它们对应的数据链表，为了确定对应的三棱柱6个顶点的相关数据,应根据式(10)计算采样点的网格层次，即采样点相对于地面的高度h，体数据的空间高度hv，体数据的层数l，以及体数据所在棱镜的上表面索引index.

(10)

因此，可以从3个数据链接表中获得index和index-1层的6个三角形顶点a，b，c，a′，b′和c′，及其对应的6个数据值，并且它们与采样点的关系如图8所示.其中，O点是地球的中心点，p为采样点，O和P的连线在上下表面的交点分别是p′和p″.首先，采用重心插值法计算P′和P″的数据值，然后采用线性插值则可获得P的数据值.

数据值和颜色值之间的转换需要使用传递函数.在文中，线性映射函数被用作体积可视化的传递函数，即图9所示的颜色映射表的蓝色末端.对应于数据的最小值，红色端对应于数据的最大值,并根据数据大小通过线性插值得到其他数据的相应位置.相应的屏幕像素的最终颜色值可以通过累积光线上的采样点的颜色来获得.

2.3 绘制模式

由于体绘制过程中需要进行大量复杂的计算，如果对图像的每一帧都进行计算，将占用很大的计算资源.根据光线投射原理，当视觉参数不变，即视点、场景、窗口等参数一致时，光线投射产生完全相同的结果，即无需重复计算；除此之外，视点的视觉辨别力在移动的过程中会下降，从而对场景渲染精度的要求也会随之减小，所以，在视点移动的过程中，可以通过适当的降低光线的计算精度，从而使绘制速度有效提升.

为了有效利用绘制资源，避免额外的开销，提升预测的效率，设置了3种绘制方式：方式Ⅰ，零绘制，即将上一帧的绘制结果直接用于下一帧的绘制；方式Ⅱ，粗绘制，即光线投射计算时设置较大的步长进行采样；方式Ⅲ，精细绘制，即实施完整的光线投射计算，同时在计算过程中按规定步长采样.这3种方式的选择和调度过程如图10所示，其中阈值f是指当视觉参数保持f帧静止时设置精细计算,设置阈值的目的是为了预防交互过程中出现假静止，避免不必要的绘制.

4 实验结果分析

为了证明该算法的有效性，在windows7操作系统下，分别使用OSG三维图形引擎、GPU编程语言CUDA以及GLSL着色语言对算法进行了实验；实验所使用的硬件环境是Intel i7-4770k 3.5 GHz CPU，NVIDIA GTX780TI GPU，以及16G RAM；测试数据分为单层数据和多层体数据，其中前者为全球海洋某一时刻的温度场数据，后者数据是来自太平洋西部海域某时刻的海洋盐度数据，经度在(99.0°E-150.0°E)和(0.5°N-52°N)范围之内，分辨率为 409×497×51.接下来将从以下方面对实验的结果进行论述.

4.1 六边形网格的全球分布及其可视化

图11为单层六边形网格在全球的分布效果示意图，其中图11(a)，(b)，(c)分别为地球中部低纬度地区和两极地区的可视化结果，可见在地球球面的任意位置，六边形网格的分布都呈现比较均匀的特性.将单层六边形网格用于全球海洋温度场可视化的效果如图12所示.

4.2 渲染效果对比

图13是基于经纬度网格的体可视化效果，图14是基于六边形网格的体可视化效果，使用了不同的颜色代表的了该位置此刻的盐度强度，越接近蓝色说明盐度越低，相反越接近红色说明盐度越高.由图13(b)和图14(b)进行对比可以看出，基于经纬度网格的体可视化过程中会出现颗粒感，特别是在一些起伏较大的区域，在使用文中的算法之可以使可视化后的效果更加平滑，整体效果表现良好.

4.3 渲染帧率对比

文中在场景渲染效率方面基于六边形网格的体可视化方法与传统方法进行了比较.如图15所示，步长越短，采样的点就越多，积分的速度就越慢，导致光线投射的时间延长，但是最终绘制的可视化效果也会因此变佳，需要在最终呈现的可视化效果和性能权衡中取1个平衡值，将其用于实际应用中.

从图15中的4条曲线可以看出，基于经纬度的方法比基于六边形网格渲染帧率高，前者在渲染帧率上增加了19.7%，造成该现象的原因首先是传统经纬网格只需1次三线性插值即可获取数据值，而基于六边形的方法在平头截锥体中进行多次插值才能获得；其次另一个原因是六边形网格在光线投射过程中要进行几次的空间坐标与经纬度坐标的转换，这同时就会增加运算的时间.但在实际的应用中两种方法在绘制帧率上的差别几乎难以发觉.

4.4 数据量对比

六边形网格方法不仅能够有效提升绘制的效果，而且还能够减少体数据量，以本文为例，经纬度网格的分辨率是409×497×51，总共有10 366 923个标量值.如果按照六边形网格进行存储，仅需9 246 934个标量值就可以达到低纬数据密度，数据量减少了10.8%，如在全球范围内使用六边形网格存储将比经纬度存储减少33.7%的数据存储量.由于硬件的发展不均衡，目前显存价格是远远要比同等容量的内存要昂贵，所以缩小体数据的容量以减少显存占用明显具有重要的意义.

4.5 算法复杂度

文献[14]中提出的基于六边形网格的光线投射算法，使用了4张二维表将各个六边形单元的中心点、顶点还有边的对应关系进行了存储，用于紧邻搜索，因为要遍历网格单元，所以算法的复杂度是O(n).而本文算法实在数据预处理阶段就生成了一张经纬度-三角形网格的检索表，所以将算法的复杂度降为了O(1)，有效的减少了计算量.