一种面向不等尺寸分块海量数据集的并行体绘制算法

朱小烨 薛 健 吕 科

(中国科学院大学工程科学学院 北京 100049)

0 引 言

在科学计算可视化领域，体绘制技术是使用计算机图形和图像处理策略将获取的三维体数据转换为在屏幕上显示和交互处理的图形和图像的一种可视化技术。体绘制被广泛用于医学和工业成像领域，因为它提供了高质量的渲染图像和方便的并行处理模式。研究人员将光线提前终止、空体素跳过等方法应用于基于GPU(图形处理器)的体绘制中，进一步提高了绘制速度。

迅速发展的成像技术和应用带来了爆炸性的数据量增长，然而传统的方法受到内存和外存访问速度和数据存取的限制，导致渲染速度相当慢。海量数据是指数据量大到无法完全加载到计算机内存中使用的数据，这种应对海量数据的研究方向也反过来促进了三维体数据可视化技术的新发展。在过去的二十年中，不断有人提出算法来克服这个问题。基于分块的体绘制方法是将原始体数据划分成便于装入内存的子数据块并依次绘制，最后将分块绘制结果合成最终结果。该方法也可以与硬件加速技术一起使用，以实现更快的渲染速度。

本文提出并实现了一种基于非均匀尺寸分块策略的体绘制方法，通过分块间的遮挡关系构建全局遮挡关系图，对分块进行分组并行渲染。与现有方法的对比实验表明，本文提出的针对Out-of-core数据的体绘制方法在海量数据集上达到了比现有方法更快的渲染速度。

1 相关工作

三维可视化的研究始于20世纪70年代中期，并且随着计算机断层扫描(CT)技术的诞生逐渐发展起来。20世纪80年代，各种影像技术不断出现，如磁共振成像、超声等，极大地促进了三维体可视化技术的发展。Levoy等[1]首次提出了光线投射算法，不久之后又相继出现了许多经典的体绘制算法，如Splatting、ShearWarp、CellProjection以及基于3D纹理的方法[2-5]。

然而，随着近几年图像采集硬件的不断更新，采集图像的分辨率也越来越高，体数据的规模也变得越来越大，促使体绘制技术的方向朝着Out-of-core技术方向发展。针对海量数据的体绘制技术应用十分广泛，其中一个重要的应用领域是工业超声无损检测领域。超声无损检测的原理是通过超声探头发射超声波针对样品进行扫描，并利用超声波与样品的回波进行成像，其优势是能够穿透检测物体内部，并且可以检测内部缺陷、观察内部结构等。由高频超声探头扫描得到的超声图像具有高分辨率的特点(一百兆的超声探头的扫描精度最高可达到一百微米)，当扫描的目标器件尺寸达到一定规模时，就会形成海量的超声体数据，利用传统的体绘制方法很难在有限内存的机器上得到渲染结果。

目前针对海量数据体绘制的研究存在两个主要的研究方向。一方面，通过经典GPGPU语言，基于CPU-GPU组合的异构并行计算平台在PC系统中变得可行，文献[6-7]描述了使用GPU加速渲染的经典方法；另一方面，传统的体绘制方法对海量数据的效果不佳，是因为它们的设计初衷并没有考虑到数据量超过内部存储器大小的情景。基于此，文献[8]提出了一种基于GPU的用于并行架构的体绘制算法，可以直接渲染不规则的海量体数据。文献[9]提出了一种从三角网格构造稀疏体素八叉树的算法，只使用核心算法所需的一小部分内存，但运行时间与上述针对海量数据的体绘制方法大致相同。文献[10]采用了新颖的SparseLeap方法使得光线投射阶段跳过了大量的空体素区域，可用于海量数据中含有大量背景区域的体绘制。文献[12]的方法将整个切片序列分为一个个由相邻两张切片组成的Slab，按光线投射方向依次导入所有Slab，在当前Slab中推进在该Slab范围内的光线，最终得到光线投射的绘制结果。基于分块策略的体绘制方法也有了长足的进步，例如：文献[13]提出的自适应分块方法使得每个分块数据能够单独装入纹理内存，并按照从后往前的顺序绘制，这种方法能够最大限度地将背景与非背景部分分离。文献[14]使用等尺寸分块策略在单机平台上实现了体绘制的硬件加速算法。文献[15]使用统一的海量数据可视化框架来处理Out-of-core数据集，并提出了一种基于半自适应分块策略的快速体绘制方法。

基于以上研究，在等尺寸分块策略下，文献[16]根据分块之间的并行关系对分块进行并行分组，并通过并行绘制分组实现了良好的体绘制加速效果。将文献[16]的算法与上述算法进行比较，可以得出这样的结论：以分块方式处理海量数据仍有改进的空间。例如，文献[16]使用的等尺寸分块策略无法有效地区分背景区域。自然地，可以考虑对不等尺寸分块的大规模数据集使用Out-of-core技术实现快速的体绘制。我们期望这种方法可用于处理任意尺寸分块的数据集。分块良好的数据集能够有效地区分背景区域且不会产生过多的由于尺寸太小而影响内存和外部数据的交换效率的“碎片”分块，例如经过文献[15]的半自适应分块方法得到的体数据。本文工作就是针对这种非等尺寸大小的分块数据集，实现类似于Yao的并行加速体绘制方法，并将其应用于超声无损检测设备扫描得到的海量超声数据的体绘制当中。

2 基于不等尺寸分块的并行体绘制算法

首先，原始的海量体数据在外存中被任意划分成不等大小的数据作为算法的输入，我们可以将它们想象成不等尺寸的体数据分块，它们在所占的存储空间中各自具有相对的空间位置坐标。

算法的整个流程如图1所示，一次完整的渲染流程(初次渲染)大致可分为以下几步：1) 设计数据结构来存储每一个子块信息，根据它们之间的空间坐标信息建立不等尺寸分块之间的相邻关系，并将相邻关系存储在分块数据结构中；2) 根据三维标量场数据体绘制的传递函数以及当前视点位置和不等尺寸分块之间的相邻关系，遍历不等尺寸分块中的非空分块，建立三维标量场数据中非空分块的全局遮挡关系图；3) 对绘制的全局遮挡关系图进行拓扑排序，得到三维标量场数据中非空分块的并行分组队列；4) 依次对并行分组队列中每一并行分组的非空分块进行并行体绘制渲染，得到三维标量场数据的并行体绘制渲染结果。

图1 算法总体流程图

2.1 建立分块数据结构

为了建立分块间的遮挡关系从而得到分组并行的绘制体数据，首先需要得到分块间的相邻关系并将其存储于分块数据结构中。分块的数据结构设计如下：

typedef struct _cuboid_block_node{

unsigned long fidx;

//分块编号

int size[3];

//分块尺寸

float bbox[6];

//分块包围盒的边界

vector adjBbox[6];

//分块包围盒各面上的相邻分块

} CuboidBlockNode;

其中，已知信息为分块编号fidx，分块尺寸size[3]，分块包围盒的边界bbox[6]。根据已知的分块编号及其包围盒的边界判断两两分块是否相邻的准则如下：

假设当前不等尺寸分块为B，其6个边界值分别为B.bbox[i]，其中，i表示边界编号，编号顺序为左右下上前后。如果有不等尺寸分块B′，其边界值满足：

B.bbox[2]B′.bbox[2]

B.bbox[4]B′.bbox[4]

B.bbox[1]-B′.bbox[0]<1

则不等尺寸分块B′是不等尺寸分块B的相邻不等尺寸分块，其接触面为不等尺寸分块B的左边界。将不等尺寸分块B′的编号放入不等尺寸分块B左边界面的相邻块节点编号数组B.adjBbox[0]中。判断分块B的包围盒上其余5个面的准则以此类推。

根据上述准则，建立所有两两分块间的相邻关系，步骤如下：1) 对于每一不等尺寸分块，根据建立的数据结构，顺序遍历所有其他不等尺寸分块，建立不等尺寸分块之间的相邻关系；2) 根据不等尺寸分块之间的相邻关系，将所有不等尺寸分块的相邻块节点编号数组存入三维标量场数据中不等尺寸分块的数据结构。完成上述步骤后，不等尺寸分块的数据结构应包括：分块编号(外存中的数据块序号)、分块尺寸、分块的包围盒边界和分块包围盒6个面上的相邻分块编号。

在完成建立分块数据结构之后，为了加速体绘制，还需要剔除不需绘制的背景区域，即经传递函数映射为完全透明的数据块。这些块称为空块，不需要进入后续的绘制流程，反之则为非空块，需要进行后续处理和绘制。若在渲染过程中由于交互操作，传递函数发生改变，则需要重新对所有数据块进行空块和非空块的分类。

2.2 建立全局遮挡关系图

得到分块间的相邻关系之后就可以根据当前视点位置建立非空分块的全局遮挡关系图，整个过程的流程图如图2所示。其大致步骤如下：首先任取一个未访问的非空分块B0，在全局遮挡关系图中建立一个对应结点，并将其加入队列。然后从队列头取出非空分块B，依次检查非空分块B的6个边界面，从未访问的边界面开始依次检查与非空分块B的未访问边界面邻接的所有非空分块。若还有未访问的非空分块，则取出与非空分块B的未访问边界面邻接的未访问非空分块Bi，并将其加入队列，并标记非空分块Bi为“已访问”。在全局遮挡关系图中创建一个对应非空分块Bi的结点，并根据当前视点位置，判断非空分块Bi与非空分块B之间的遮挡关系。若非空分块B遮挡非空分块Bi，则在全局遮挡关系图中创建一条非空分块B到非空分块Bi的边；若非空分块B没有遮挡非空分块Bi，则在全局遮挡关系图中创建一条非空分块Bi到非空分块B的边。若没有未访问的非空分块，则返回继续检查下一边界面。最后判断是否还有未访问的非空分块，若还有未访问的非空分块，则将其加入队列重复以上步骤；若没有未访问的非空分块，则得到最终的全局遮挡关系图。

图2 建立全局遮挡关系流程图

相邻两个分块的遮挡关系判断准则如下：如图3所示，从当前视点到当前非空分块1的中心构成向量I1，垂直于非空分块1与非空分块2的接触面且指向非空分块2的向量为I2。当I1·I2>0时，非空分块1遮挡非空分块2；否则，非空分块1不遮挡非空分块2。

图3 当前视点位置下的相邻分块遮挡关系判断示意图

2.3 获取分块并行分组并绘制

对上述过程得到的分块的全局遮挡关系图进行拓扑排序所采用的方法与文献[16]的方法类似，用以得到分块的并行分组及其绘制顺序。

得到并行分组后就可以按绘制顺序依次对每个并行分组中的分块数据进行并行体绘制渲染。假设当前要绘制的分组中有n个非空分块，其具体过程为：1) 绘制当前并行分组中n个非空分块的后表面，得到每一非空分块的光线出点，将其作为三通道颜色值记录在一二维纹理中，同时，将每一非空分块的分块编号作为颜色值写入另一个二维纹理中；2) 绘制当前并行分组中n个非空分块的前表面，在片段着色器中并行处理每一条光线，根据当前片段着色器处理的光线入点在最终投影图像中的位置，访问两个二维纹理，得到当前片段着色器处理的光线出点和非空分块的分块编号；3) 根据非空分块的分块编号，在包含当前并行分组中n个非空分块的子三维纹理中找到当前光线所在非空分块的位置，得到当前光线对应的分块数据场；4) 根据当前光线对应的分块数据场，从当前非空分块的光线入点和光线出点之间按照预设的步长进行采样，并获取每一采样点的颜色值与不透明度值；5) 根据每一采样点的颜色值与不透明度值，依次对并行分组队列中每一并行分组的非空分块进行并行体绘制渲染，得到三维标量场数据的并行体绘制渲染结果。其中，由于需要在不同批次的分块绘制间保留每条投射光线的累积不透明度，需采用如下采样合成公式：

Cn+1=Cn+rn×α(sn)×C(sn)

(1)

式中：αn+1和Cn+1分别表示累加合成的不透明度和颜色值；α(sn)和C(sn)分别表示光线sn处采样所得三维纹理的不透明度和颜色值；rn=1-α(sn)表示剩余不透明度，将被记录在渲染结果纹理中，传递到下一批分块绘制流程。

3 实验结果与应用

为了验证所提出的算法，设计实验分别对不同数据量级的三个数据集进行渲染。三个数据集均为超声无损检测设备采集获得的超声数据，第1个数据集是一枚普通芯片的扫描数据(记为D1)，大小为113.7MB。第2和第3两个数据集(分别记为D2、D3)是扫描用于探伤测试的芯片的结果，大小分别为9.55GB和25.48GB。所有测试均在配备Intel Core i5-65003.20GHz CPU，8GB物理内存和GeForce GTX1060显卡的普通PC机上运行。D2和D3数据集大小已超出实验机器的内存容量限制，在实验环境下可认为是Out-of-core数据。图4为本文方法分别在三个数据集上关闭光照与开启光照的绘制结果。由图4可知，本文提出的方法在渲染普通数据集与Out-of-core数据集时均表现良好。

(a) D1数据集(左：关闭光照，右：开启光照)

(b) D2数据集(左：关闭光照，右：开启光照)

(c) D3数据集(左：关闭光照，右：开启光照)图4 本文算法在D1、D2、D3数据集上的部分绘制结果

为了测试新算法的效率，分别将它与当前的几种分块体绘制算法作对比，包括：A1，一种类似于Novins等[12]提出的，但通过OpenMP进行并行加速的算法；A2，一种基于半自适应分块策略且串行渲染每个分块的体绘制算法；A3，一种基于等尺寸分块策略且串行渲染每个分块的算法；A4，文献[16]提出的基于等尺寸分块策略且分组并行渲染分块的算法。A5为本文所提出的方法。

表1为上述5种算法在不同数据集(D1-D3)上的渲染速度表现。渲染速度是指对应算法在各自三种不同的数据集上渲染一帧所需的时间，每组数据都经过10组以上的实验结果取平均值得到。

表1 渲染速度的对比

续表1

表2为上述5种算法在不同数据集(D1-D3)上的内存使用情况。缓冲数据占用内存是指对应算法存放临时分块数据所需的内存缓冲区大小；辅助数据占用内存是指算法执行时用到的辅助数据结构占用的内存，如A5算法所用到的分块数据结构与遮挡关系图等；C是指常量，即该项数据与数据集规模大小无关。

表2 内存消耗的对比

由表1可以看出，A5算法在渲染速度方面有良好的性能表现。当数据量足够大时，与采取分块的串行渲染方法A2和A3相比，无论划分策略如何，A5算法的渲染速度都快50%以上。另外，A5算法也总是比采取等尺寸分块策略的A4快，且随着数据量的增加，A5算法的速度优势变得更加明显。总之，无论是对于可以一次性装入内存的普通数据集还是Out-of-core数据集，A5算法均有最快的渲染速度。

就内存消耗方面而言，由表2可以看出，A5算法无论是缓冲数据还是辅助数据占用内存都略大于A1-A3算法，但是与原始数据集的大小相比两种数据占用的空间都是微不足道的。因此，就运行时的内存消耗而言，A5算法也是可以接受的。

本文提出的算法已经成功应用于海量超声数据的处理和可视化系统中。图5为用实验室研制的超声无损检测设备扫描一个约5 cm×5 cm×8 mm大小的芯片得到的一组渲染图片，渲染芯片所用的超声数据是由超声信号频率为50 MHz的超声探头采集得到。

图5 本文算法在超声数据集上的部分绘制结果(右下为超声数据采集设备)

4 结 语

本文提出了一种针对Out-of-core数据的并行体绘制算法。该算法可用于处理采用非均匀大小子块划分的Out-of-core数据集，在最大化背景区域的同时控制分块的数量并减少内部和外部数据交换开销；在丢弃不需绘制的空子数据块之后，根据当前视点下分块间的遮挡关系构造全局遮挡关系图，进一步对有并行关系的分块进行分组并行绘制，在GPU上并行渲染非空子块，以此提升渲染速度。实验表明，本文算法可以在获得高分辨率的渲染结果的同时取得比现有方法更快的渲染速度。

本文算法目前在普通单机环境下运行良好，未来可以扩展到分布式环境下，在众多计算节点上同时绘制分组并行的非空分块，进一步提高渲染效率。