一种近似无损的三维体数据压缩域体绘制方法

摘要

本发明公开了一种近似无损的三维体数据压缩域体绘制方法：将体数据分割成大小为4*4*4的块后逐块拉普拉斯分解为三个层次，在对第一、第二层次的分解数据分别构建矢量量化编码表后进行矢量量化，获取第一、第二级索引体数据，再以两级索引体数据和两级矢量量化编码表恢复体数据，将原始体数据与恢复体数据作差得到稀疏的残差体数据；使用完美空间哈希压缩该稀疏残差体数据；以两级索引体数据、两级矢量量化编码表和完美空间哈希的结果作为体绘制阶段的输入，在GPU里进行实时解码和绘制。本发明使解码数据更接近原始数据；解决了有限的存储资源、有限的CPU及GPU传输带宽与不断增加的体数据尺度及领域专家对高压缩质量的实时压缩域体绘制的需求之间的矛盾。

说明书

技术领域

本发明涉及一种三维体数据压缩域体绘制方法，具体涉及一种基于层次矢 量量化和完美空间哈希技术的近似无损的压缩域体绘制方法。

背景技术

体绘制技术是可视化研究中的关键技术之一，在计算流体力学、地球空间、 医疗科学等各个领域中得到了广泛应用，其研究对象是三维体数据(以下简称 体数据)。体数据产生于科学计算的各个领域，在现今的领域需求以及数据获取 途径下，体数据的尺度越来越大、分辨率越来越高，并且经常是海量的、时变 的、多属性的，太字节(TB)级甚至拍字节(PB)级的体数据已屡见不鲜。利 用可视化技术分析如此大尺度的数据充满挑战性。近年来，商业图形硬件发展 迅速，利用GPU加速体数据的直接体可视化已经成为分析体数据的通用方法。

然而，一般计算机的GPU显存无法一次性装载大尺度的体数据；即使可以 装载，GPU和CPU之间的数据传输速度过慢，制约了硬件加速直接体绘制算法 的应用，同时也阻碍用户通过可视化探求复杂大尺度数据中蕴藏的未知知识和 现象。例如，大小为1024³的单变量单帧浮点型体数据，其大小为 32bit×1024×1024×1024＝4GB，又如帧数为100帧，每帧的分辨率为512³的浮点型 单变量时变体数据，其数据量大小为32bit×512×512×512×100＝50GB，这些体数据 的数据量远超普通计算机的GPU显存，甚至内存。压缩域体绘制(Compression  domain volume rendering,CDVR)融合了体数据压缩和绘制，能有效节省存储和 带宽，是解决这一矛盾的重要而有效的方法。

目前已有的CDVR方法或者数据压缩方法中，压缩比高的通常数据恢复质 量不好，依据恢复数据所得到的体绘制结果图像的信噪比(Signal to Noise Ratio， SNR)或者峰值信噪比(Peak Signal to Noise Ratio，PSNR)较低,而均方差(Mean  Squared Error，MSE)则较高。数据恢复质量好的(如一些无损压缩算法)压缩 比太低，或者实现复杂、解码速度慢，不适合在GPU中实现，无法达到实时 CDVR。矢量量化(Vector Quantization,VQ)作为一种有损的数据压缩策略，常 用于图像和语音等多媒体信号的处理中。作为一种非对称的压缩方法，其压缩 相对复杂而解压缩异常简单的特性使得其成为CDVR的一类重要方法。层次矢 量量化(Hierarchical Vector Quantization,HVQ)是一种新的层次化的矢量量化 方法，较传统的矢量量化压缩方法，能极大地改善压缩效率和绘制质量。完美 空间哈希(Perfect Spatial Hashing,PSH)技术使用无冲突的静态哈希方法表达稀 疏数据，其实质在于充分利用数据的空间连续性，将一个稀疏的数据转换为紧 凑的哈希表进行表达，因而可以作为一种高效的稀疏数据压缩方法。

发明内容

本发明的目的在于提供一种近似无损的三维体数据压缩域体绘制方法，以 克服现有技术的全部或部分缺陷。

为实现上述目的，本发明所采取的技术方案是：本发明近似无损的三维体 数据压缩域体绘制方法包括如下步骤：

(1)将体数据值类型为浮点型、无符号字符型或者短整型的三维原始体数 据分割成大小为4*4*4的数据块，再对分割完的数据块逐个进行拉普拉斯分解， 由此将每个数据块分解为三个层次的数据：其中，每一个数据块的第一层次分 解数据是一个大小为4*4*4的64维矢量，所有数据块的第一层次分解数据构成 第一层次分解数据矢量集合；每一个数据块的第二层次分解数据是一个大小为 2*2*2的8维矢量，所有数据块的第二层次分解数据构成第二层次分解数据矢量 集合；每一个数据块的第三层次分解数据是该数据块的体数据值的算术平均值， 所述第三层次分解数据的大小为1*1*1，所有数据块的第三层次分解数据构成第 三层次分解数据集合；

(2)对所述第一层次分解数据矢量集合进行矢量量化，得到第一级矢量量 化编码表和第一级索引体数据，所述第一级索引体数据中记录了每一个4*4*4 数据块的第一层次分解数据与第一级矢量量化编码表中码元矢量的对应关系； 构造第一级偏移地址表，由此确定第一级矢量量化编码表中的码元矢量的元素 与每个数据块的第一层次分解数据中的体素的对应关系：

对所述第二层次分解数据矢量集合进行矢量量化，得到第二级矢量量化编 码表和第二级索引体数据，所述第二级索引体数据中记录了每一个4*4*4数据 块的第二层次分解数据与第二级矢量量化编码表中码元矢量的对应关系；构造 第二级偏移地址表，由此确定第二级矢量量化编码表中的码元矢量的元素与每 个数据块的第二层次分解数据的体素的对应关系：

(3)利用第一级矢量量化编码表、第一级索引体数据、第一级偏移地址表 以及第二级矢量量化编码表、第二级索引体数据、第二级偏移地址表分别进行 矢量量化解码而相应得到所有数据块的第一层次分解数据和第二层次分解数据 的解码恢复值，再逐体素地对所有数据块的每个体素所对应的第一层次分解数 据和第二层次分解数据的解码恢复值以及第三层次分解数据相加得到解码体数 据，将原始体数据与解码体数据逐体素作差而得到稀疏的残差体数据；

(4)使用完美空间哈希压缩所述残差体数据而得到偏移表数据、位置表数 据和哈希表数据；

(5)将第一级索引体数据、第二级索引体数据、第一级矢量量化编码表、 第二级矢量量化编码表、一级偏移地址表、二级偏移地址表、偏移表数据、位 置表数据和哈希表数据传入到GPU中并在GPU中进行实时解码和绘制。

进一步地，在所述步骤(2)中，按照公式(1)所示的方法构造第一级偏 移地址表：

value1[i,j,k]＝((i/2*4+j/2*2+k/2)*8+i％2*4+j％2*2+k％2)/64.0   (1)

公式(1)中，i,j,k分别表示待解码体素在其所属的4*4*4数据块内的三个 方向上的位置，i,j,k的取值范围皆为[0,3]；value1[i,j,k]表示位置为i,j,k的待 解码体素在第一级矢量量化编码表中它所对应的码元矢量内部的偏移地址。

进一步地，在所述步骤(2)中，按照公式(2)所示的方法构造第二级偏 移地址表：

value2[i,j,k]＝(k/2+j/2*2+i/2*4)/8.0     (2)

公式(2)中，i,j,k分别表示待解码体素在其所属的4*4*4数据块内的三个 方向上的位置，i,j,k的取值范围皆为[0,3]；value2[i,j,k]表示位置为i,j,k的待 解码体素在第二级矢量量化编码表中它所对应的码元矢量内部的偏移地址。

由上可见，本发明近似无损的三维体数据压缩域体绘制方法的主要过程为： 首先使用层次矢量量化技术压缩原始体数据，再利用完美空间哈希压缩层次矢 量量化后的残差数据，解码时把HVQ解码数据值加回PSH的解码数据值以达到 近似无损的效果。

与传统的压缩域体绘制技术相比，本发明的有益效果在于：提出一套基于 HVQ和PSH的体数据压缩域体绘制方法，对体数据进行高压缩比的近似无损压 缩，解决了有限的存储资源、有限的CPU与GPU传输带宽和不断增加的体数 据尺度以及领域专家对于高压缩质量的实时压缩域体绘制的需求之间的矛盾。 结合HVQ以及PSH技术的特点，在绘制时利用GPU进行实时的解压缩，直接 在压缩域进行体绘制，充分地利用了GPU的优势，解码端实现简洁快速，可满 足实时渲染的性能要求；同时针对领域专家、用户需要更精确的解码数据和可 视化效果这一实际问题，提供了近似无损的压缩效果，可明显提高解码数据的 质量，使得解码数据更接近原始数据，从而提高用户对于三维数据场体绘制效 果的感知，为用户准确地分析、观察体数据、发现数据特征、探索未知现象营 造更好的客观环境。

附图说明

图1是本发明方法的流水线的示意图【包含本发明的步骤(1)-(5)】。

图2是本发明的步骤(1)和(2)对原始体数据进行HVQ压缩的示意图。

图3是本发明的步骤(4)对稀疏残差体数据进行PSH压缩的示意图。

图4是本发明的步骤(5)对HVQ压缩数据进行HVQ解码的示意图。

图5是对Hurricane Isabel Cloud数据使用DVR的绘制结果。

图6是对Hurricane Isabel Cloud数据使用基于HVQ的CDVR的绘制结果。

图7是对Hurricane Isabel Cloud数据使用本发明的方法的绘制结果。

具体实施方式

下面结合附图对本发明的近似无损三维体数据压缩域体绘制方法作进一步 说明，具体步骤如下：

步骤(1)：首先，将体数据值类型为浮点型、无符号字符型或者短整型的 三维原始体数据分割成大小为4*4*4的数据块(以下简称一级数据块)，然后对 每一个一级数据块进行拉普拉斯分解，将块内数据分解为三个频率波段，进行 分层表示，三个层次分解数据的大小分别为4*4*4、2*2*2和1*1*1，如图2所 示。

以下以一个具体实例来详细说明对一级数据块进行拉普拉斯分解的过程。 首先，计算一级数据块的体数据值的算术平均值M₀，并将所述一级数据块分为 8个大小为2*2*2的二级数据块(以下简称二级数据块)，再分别计算各个二级 数据块的体数据值的算术平均值M₁～M₈。一级数据块的64个体素的体数据值与 体素各自所在的二级数据块的体数据值的算术平均值M_i(i∈[1,8])之差构成了64 维的第一层次分解数据，M₁～M₈与M₀之差构成了8维的第二层次分解数据，M₀构成第三层次分解数据。最终，对原始数据逐块进行拉普拉斯分解以后，整个 体数据转变为3个层次的表达，使得每个一级数据块可以用分解完成的三层数 据表示。

所有数据块的第一层次分解数据构成第一层次分解数据矢量集合，所有数 据块的第二层次分解数据构成第二层次分解数据矢量集合，所有数据块的第三 层次分解数据构成第三层次分解数据集合。

步骤(2)：使用矢量量化技术对第一层次分解数据矢量集合和第二层次分 解数据矢量集合进行压缩：即分别构建第一级矢量量化编码表和第二级矢量量 化编码表，以矢量维度64和8分别对第一层次分解数据矢量集合和第二层次分 解数据矢量集合进行矢量量化，获取对应的第一级索引体数据和第二级索引体 数据(如图2所示)。具体过程如下：

矢量量化是一种简单的、有损的数据压缩方法，其主要过程包含三个阶段： 码表设计、编码和解码，其中码表设计最为重要。矢量量化的本质在于通过构 造一个包含W个B维矢量的码表，将B维的输入矢量X用码表中的某一个码元 矢量进行近似表示，并用所述码元矢量在码表中的索引值代替X进行存储，当 W较小时，每一个索引只需要很少的比特位数，从而达到压缩的效果。解码则 是利用所述索引值在码表中查询获取相应的码元矢量以重构输入矢量X。对于 第一和第二层次的分解数据矢量集合，用矢量量化进行压缩，设定码表的大小 为256，索引数据类型为无符号字符型(unsigned char)，每个索引大小为8比特， 矢量的数据维度分别为64和8。第三层的数据是每个一级数据块的体数据值的 算术平均值，直接作为对应层次的索引数据。如图2中部所示，一级数据块最 终可以用两个索引值和一个算术平均值替代，两个索引值所索引的码元矢量维 度分别是64和8。在具体的实现过程中，采用基于主成分分析(Principal  Component Analysis,PCA)的分裂策略训练初始码表，然后进行矢量量化。

主成分分析是多元统计分析中用来分析数据的一种方法，借助变量的协方 差矩阵对信息进行处理，提取较少数量的特征对样本进行描述以达到降低特征 空间维数的目的。通过基于主成分分析的分裂策略获得初始码表的过程如下：

1)将当前所有矢量归到一个类D中，计算类D的中心矢量Y(即矢量均值)， 将类D中所有矢量与Y的欧氏距离的总和记为dis。构建一个双向链表，链表中 每个节点代表一个类D，节点按dis值大小的降序排列。

2)每次从链表头中取一个节点，得到一个类D₀，计算其协方差矩阵M。

3)计算协方差矩阵M的特征矢量，取其最大特征矢量e_max。

4)对D₀中的每个矢量X，计算<X,e_max>的值，若小于0，则将X归入左子 节点类D_l中，否则X归入右子节点类D_r中。

5)计算类D_l、D_r的中心矢量Y_l、Y_r及各自类内的所有矢量距离其类中心 矢量的欧式距离之和dis_l、dis_r，删除链表头节点，并将D_l、D_r插入到链表中。

6)重复2)—6)的操作直到链表大小达到码表大小。取链表中每个类的中 心矢量Y作为初始码表的码元矢量，中心矢量的集合便构成了初始码表。

得到初始码表之后，对码表进行精炼，使码表更加精确。码表精炼的核心 为K均值聚类(K-means)算法，其过程描述如下：

1.使用PCA生成的初始码表的码元矢量作为K-means每个类的初始中心， 码表大小即为K的取值，亦即K-means分类的个数。

2.对每一个类，计算类内的每一个矢量到其中心矢量的距离之和dist_old；对 所有的矢量，计算其和每一个类的中心矢量的欧式距离，将其归类到距离最近 的中心矢量所在的类。

3.对当前所有的聚类，重新计算聚类的中心矢量。

4.计算所有矢量到其所在类的中心矢量的距离之和dist_new，与上一次距离和 dist_old做差，得到dist_change，当dist_change/dist_old绝对值小于一定的阈值，则退出得 到最终码表，否则回到步骤2。

至此，得到第一级矢量量化编码表和第一级索引体数据以及第二级矢量量 化编码表和第二级索引体数据，所述第一级索引体数据中记录了每一个一级数 据块的第一层次分解数据与第一级矢量量化编码表中码元矢量的对应关系，所 述第二级索引体数据中记录了每一个一级数据块的第二层次分解数据与第二级 矢量量化编码表中码元矢量的对应关系。构造第一级偏移地址表，由此确定第 一级矢量量化编码表中的码元矢量的元素与每个一级数据块的第一层次分解数 据中的体素的对应关系。构造第二级偏移地址表，由此确定第二级矢量量化编 码表中的码元矢量的元素与每个一级数据块的第二层次分解数据的体素的对应 关系。一级偏移地址和二级偏移地址的计算分别如公式(1)和(2)：

value1[i,j,k]＝((i/2*4+j/2*2+k/2)*8+i％2*4+j％2*2+k％2)/64.0   (1)

value2[i,j,k]＝(k/2+j/2*2+i/2*4)/8.0     (2)

其中，i,j,k分别表示待解码体素在其所属的一级数据块内的三个方向上的位 置，取值范围皆为[0,3]；value1[i,j,k]表示位置为i,j,k的待解码体素在第一级矢 量量化编码表中它所对应的码元矢量内部的偏移地址；value2[i,j,k]表示位置为 i,j,k的待解码体素在第二级矢量量化编码表中它所对应的码元矢量内部的偏移 地址。

步骤(3)：构建完两级矢量量化编码表和两级索引体数据以后，进行一次 HVQ解码以获取恢复体数据，以便将原始体数据与恢复体数据作差得到稀疏的 残差体数据，进而将该残差体数据作为PSH的输入数据进行压缩。HVQ解码过 程如下：

对原始体数据的每一个体素，明确其在原始体数据中所属的一级数据块的 编号以及在该一级数据块中所属的二级数据块的编号，从而即可在第一级索引 体数据和第二级索引体数据中找到该体素在第一级矢量量化编码表和第二级矢 量量化编码表中的索引值以及该体素所在的一级数据块的体数据值的算术平均 值，再利用所述两个索引值在对应的第一级矢量量化编码表和第二级矢量量化 编码表中查询到所对应的两个码元矢量，将所述两个码元矢量中与待恢复体素 所对应的元素取出相加，再将相加的结果与该体素所在的一级数据块的体数据 值的算术平均值相加，得到恢复体素值，从而完成HVQ解码。还原得到的恢复 体素值与原始体数据对应体素值之差作为该体素的残差值。将整个残差体数据 保存下来以作为完美空间哈希方法的输入数据。

步骤(4)：使用完美空间哈希压缩所述残差体数据而得到偏移表数据、位 置表数据和哈希表数据。具体如下：

PSH利用数据的空间连续性，将稀疏的数据转换为紧凑的哈希表，因此它 是一种压缩率接近稀疏数据集密度的压缩方法。将步骤(3)中得到的稀疏残差 数据用PSH进行压缩，构建哈希表的具体过程如图3所示。

图3中，p表示有效数据的体素(即有残差值的体素)在残差体数据中的空 间位置，D(p)表示p上的残差数据值，S表示所有有效体素空间位置的集合，U 表示整个三维残差体数据空间域，s表示哈希表中的位置，h₀和h₁是两个非完美 哈希函数，实际的哈希操作为模寻址，模操作的对象分别为哈希表的大小ots以 及偏移表的大小pts，即h₀(p)＝p mod ots,q＝h₁(p)＝p mod pts，mod是模操作。 最终的完美哈希函数可表示为

h(p)＝h₀(p)+φ(q)＝h₀(p)+φ(h₁(p))          (3)

公式(3)中，φ(p)表示位置p处的哈希值，h(p)表示位置p处的最终的完 美哈希值。在具体实现中，用三个三维纹理来实现偏移表数据、位置表数据和 哈希表数据。第一个为三通道的三维偏移表纹理，每个元素分别表示三个方向 的偏移值，它将非完美哈希函数h₀转换成完美哈希函数；第二个为三通道的三 维位置标记表纹理，表中元素存放哈希表中每个元素在原始数据中的位置；第 三个是单通道的三维哈希表纹理，表中元素存放最终所要查找的数据值。

步骤(5)：在本步骤中，将第一级索引体数据、第二级索引体数据、第一 级矢量量化编码表、第二级矢量量化编码表、一级偏移地址表、二级偏移地址 表、偏移表数据、位置表数据和哈希表数据传入到GPU中并在GPU中进行实时 解码和绘制。具体如下：

通过步骤(1)-(4)已经得到一系列的压缩数据，其中，将第一级索引体 数据、第二级索引体数据和算术平均值数据构造成一个三通道三维索引纹理， 将第一级矢量量化编码表、第二级矢量量化编码表分别构造成一个二维单通道 码表纹理，将偏移表数据构造成一个三通道三维偏移表纹理，将位置表数据构 造成一个三通道三维位置标记表纹理，将哈希表数据构造成一个单通道三维哈 希表纹理，将一级偏移地址表和二级偏移地址表构造成一个大小为4*4*4(即和 一级数据块的大小相同)的双通道三维寻址纹理提供实际待解码体素与所索引 矢量具体元素的对应关系，上述纹理数据作为体绘制阶段的输入数据，在GPU 里进行实时解码和标准光线投射体绘制(以下简称体绘制)，解码和体绘制是紧 耦合的关系。

此外，步骤(5)中所提及的所有纹理的滤波器模式设定为GL_NEAREST， 以避免插值变化。双通道三维寻址纹理的地址取样属性GL_WRAP取为 GL_REPEAT，将该纹理的数据块重复平铺于原始体数据空间。

解码分成两个步骤。第一步是HVQ的GPU解码，如图4所示，在体绘制 过程中，通过当前体素的空间位置确认其在三通道三维索引纹理里的纹理坐标， 然后获得当前体素所属的一级数据块的体数据值的算术平均值以及第一级码表 索引值和第二级码表索引值，再从双通道三维寻址纹理中获取第一级偏移地址 和第二级偏移地址。使用第一级码表索引值和第一级偏移地址构造第一个纹理 坐标，使用第二级码表索引值和第二级偏移地址构造第二个纹理坐标，通过上 述两个纹理坐标可分别确定当前体素在两个矢量量化编码表中的码元矢量里的 真实位置，再以所述的两个纹理坐标从各自对应的二维单通道码表纹理中获取 码元矢量，将当前体素在两个码元矢量中的两个解码值以及该体素所在的一级 数据块的体数据值的算术平均值相加得到该体素的重构体数据值。逐体素完成 这一过程，HVQ解码完毕。

第二步是PSH的GPU解码，即获得当前体素的残差值。这一步的解码较简 单，使用当前体素初始位置对偏移表的长度进行模操作即可获得当前体素在偏 移表纹理中的纹理地址，称为偏移表纹理地址。再使用偏移表纹理地址从偏移 表纹理中获取偏移地址的值，将该值与当前体素位置相加得到新的体素位置， 然后使用新的体素位置对哈希表的长度进行模操作即可获得当前体素在哈希表 纹理中的纹理地址，称为哈希表纹理地址，该地址同样是当前体素在位置表纹 理中的纹理地址。然后先使用哈希表纹理地址从位置表纹理中进行取值，若取 出的值和当前体素初始位置一致，则再使用该地址从哈希表纹理中取得实际的 残差值；若不一致，则返回0值。最后将上述HVQ解码值和PSH解码值相加 作为最终的解码体数据值，整个解码过程结束。以最终的解码体数据值在GPU 中进一步进行标准的光线投射体绘制过程即可完成压缩域体绘制。

下面将给出本发明在一个真实体数据的应用结果以说明本发明技术的有效 性。实验平台为一台个人电脑，处理器为Intel Core(TM)i7-4770,主频3.4GHz， 内存16GB，显卡是NVIDIA GeForce GTX660，显存是2GB。从压缩性能角度， 对两个体数据进行了实验，性能参数对比结果如表1所示，HVQ的码表大小为 256，分块大小为4*4*4，索引的数据类型为无符号字符型。表1数据尺度格式 为x×y×z×t，x、y和z为数据空间尺度，t是时间步,压缩比是原始数据字节数 与压缩后数据字节数的比率。帧率是每个数据在每种方法下各运行十次，取算 术平均值，绘制窗口分辨率为512*512。

表1 压缩域体绘制性能参数对比结果

由实验结果可知，与已有的无损浮点体数据压缩方法(压缩比普遍不高，其 压缩率在2-3之间)相比，本发明方法压缩率更高；而与有损的基于矢量量化的 方法，如HVQ方法相比，本发明方法压缩比虽有所下降，但是在“近似无损” 的前提下，取得了可观的压缩比，尤其是面对稀疏体数据，如Hurricane Isabel  Cloud，压缩比更能说明本发明方法在面对此类数据时的巨大应用潜力。此外， 对实验数据，本发明的方法均保持不错的绘制帧率，达到实时绘制的要求。

图5、6、7是对Hurricane Isabel Cloud分别使用DVR、基于HVQ的CDVR 和本发明的方法的绘制结果。由图5-7可知，相对于HVQ方法，本发明方法解 码数据更接近原始数据，误差更小，失真明显减轻，绘制结果的视觉效果更加 光滑，与原始数据的DVR结果更加接近。同样地，对图6和图7的定量分析 的对比结果如表2所示，由此可见，相对于HVQ方法，本发明方法对该数据的 处理结果在质量上有了较大提升，红绿蓝三个颜色通道上的MSE均有明显下降， 而PSNR则均有提高。

表2 Hurricane Isabel Cloud压缩质量对比分析