1 引言
近年来，随着多媒体通信和3D采集技术的发展，3D点云数据广泛应用于增强现实、游戏、远程通信、沉浸式交流等新兴产业。3D点云(3D Point Cloud)是由大量拥有具体三维位置信息的点集组成，其中每个点具备一个或多个特征(如：颜色，法线等)。作为新型空间数据类型，3D点云适合用来表示3D模型和空间，且计算处理起来非常高效。尽管3D点云数据结构有着许多优势，但是其存在着数据量巨大的问题。比如，对于一般的点云来说，点的个数通常是百万数量级的。3D点云巨大的数据量成为3D点云实际应用的瓶颈问题。因此，在有限的带宽和存储空间下，如何设计高效的3D点云压缩算法就显的尤为重要[1-2]。
但是与传统的图像视频不同，3D点云的特征是不规则的。这种不规则性主要表现为3D点云的每个点并不是规则的落在标准的三维坐标空间上，而且对于3D点云视频序列，每一帧的点的个数也都是不同的，这些都给3D点云的压缩带来了极大难度。为了避免点云数据的不规则性，通常采用体素的形式(Voxel Representation)来表示不规则的3D点云[3-5]。具体过程如下：取一个固定的步长对3D点云的位置坐标进行量化，则点云的三维坐标便落在维度为2L×2L×2L的规则且轴对称的三维网格上，其中L 是深度。如果一个体素(Voxel)至少包含一个点，那么就称这个体素被占用。一个体素的几何位置ν∈
3×1，即对应于体素中心的三维坐标。一个未被占用的体素是透明的，不具备任何特征属性。体素化(Voxelized)的3D点云可以高效地排列为八叉树结构[6]。为了处理起来更有效率，我们通常把一个3D点云分成许多个大小为k×k×k的块，以每个块作为处理单元。同样的我们只对被占用的块进行处理。3D点云数据的压缩任务主要包括位置信息(或者称几何信息)和点云特征的压缩，点云特征主要是针对点云的颜色信息。位置信息的压缩算法研究比较成熟，目前基本都在沿用基于八叉树的位置压缩方法[6]，于是对于 3D点云的颜色信息压缩成为今年重点研究方向。Zhang等人[2]运用图形变换(Graph Transform，GT)对3D点云数据进行去相关。首先，图可由一个3D块内被占用的体素得到。然后，通过距离倒数模型可以得到图形拉普拉斯矩阵(Graph Laplacian Matrix)，其特征向量矩阵就是最终的变换矩阵。Queiroz and Chou[7]提出了一种区域自适应分级变换算法(Region-Adaptive Hierarchical Transform, RAHT)。该变换类似于Haar小波的自适应变化的分级子带变换，可以取得与GT差不多的压缩效果，但是计算复杂度却比GT减少很多。Cohen[5]等人将现有图像压缩算法(Shape Adaptive DCT, SA-DCT)[8]拓展到体素化的3D点云数据上，也取得了不错的压缩效果。Queiroz and Chou[9]进一步提出了类似于GT的基于高斯过程(Gaussian Process)的去相关算法，简称GPT，其性能优于上述所提算法。本文中，我们针对GT算法图形构造的不足提出了一种改进的图形变换算法使得每个块只产生唯一的子图。GT,RAHT和GPT这三个压缩算法中变换系数的幅值都偏大，针对这一不足，本文中用每个块所有点的均值作为直流系数，并且对去平均后的颜色值进行图形变换。熵编码部分，GT,RAHT和GPT都是将量化后所有交流系数进行编码，而在实验中我们发现量化后的交流系数有很大比例都为0值，针对这一点我们采用Run-Level编码方法对非零系数进行编码。对于直流系数，我们采用类似于H.264的预测编码方法对直流系数进行预测编码。实验结果表明与GT,RAHT和GPT相比，本文算法能取得不错的压缩效果。
2 改进的图形变换
2.1 图形变换
图形变换是一种线性变换，由潜在的图形结构决定，广泛应用于压缩网格的几何结构、深度图等方面[10-11]。图形变换简单介绍如下：
对于一个任意的图G=(V={1,…,n},ε)，其中V表示图中一系列的顶点，ε表示图中的边。我们在图的顶点上定义一个多元高斯随机向量X=(x1,…,xn)T，称之为高斯马尔科夫随机场(Gaussian Markov Random Field)。假定X的均值为μ，精度矩阵为Q，则X的概率密度为：
(1)
其中Qij≠0，Q是X协方差矩阵的逆矩阵。值得一提的是，精度矩阵Q的特征向量矩阵实际上就等价于Karhunen-Loeve Transform[12]。
2.2 改进的图形变换
图1 GT构造的图和图对应的相邻矩阵
Fig.1 Graph based on GT and the adjacency matrix
如图1所示，(a)是GT算法中构造的图,(b)是由构造的图形得到的相邻矩阵A，其中ε是点与点之间的距离。由A我们可以得到图形拉普拉斯矩阵Q，而它的特征向量矩阵Φ即为进行图形变换的变换矩阵。3D点云经过图形变换后将会产生一系列的变换系数。这些系数可分为两类：一类是幅值较大，包含较多信息的低频系数，称为直流系数；另一类是幅值趋于零，包含信息较少的高频系数，称为交流系数，具体的细节可以参考文献[2]。但是文献[2]构造图的方法可能会在每个块中产生多个子图。然而，子图的个数与变换后的直流系数的个数一致。也就是说，如果一个块中包含s个子图的话，变换后也将产生s个直流系数，这不利于编码。为此，我们提出一种改进方法，对每个块中的多个子图进行连接，使得每个块中只会存在一个图，从而将直流系数的个数固定为一。文献[13]中的实验结果表明，随着连接点的个数增加，图形变换的压缩效果会下降。因此，我们需要用较少的边将这些独立的子图连接起来。因此，受此启发，改进的图形变换如下：
对于一个包含N个被占用体素的块，用GT的方法将会得到一个N×N的相邻矩阵A。由于A是个对角对称矩阵，且对角线上的元素都为0，我们只需要分析矩阵A的左下半区即可。我们对A进行逐行处理(第一行除外)，对于第i行，计算该行前i个元素的和，如果和为0，则A(i,i-1)=1,A(i-1,i)=1。从第二行到最后一行，一共需要扫描N-1次。也就是说，如果一个点和之前扫描的都不相邻的话，那么当前的点将与上一个点连接起来，且令两点间的距离为1。图2为改进后的图。由图可见，当扫描到第5行的时候，第5个点与之前4个点都未连接，即第5个点与之前4个点是独立的，则将4,5两个点连接起来形成一个图。对于一个有s个独立子图的块，通过本文所提改进图形变换则只需添加s-1条边即可将所有子图连接为1个图，实现了用尽可能少的边来构造一个图。
图2 多子块连接后的图
Fig.2 The graph through multiple sub-block connection
3 量化和熵编码
3.1 量化
对于每一个块，改进后的N×N变换矩阵Φ可以由被占用的体素的位置信息得到。落在N个顶点的颜色信号可以记为三个N×1的向量(对应于YUV的三个通道)。我们对YUV三通道分别进行编码。假设y=[y1,…,yn]是Y通道的信号，对y进行变换可得变换系数：
c=Φ(y-mean(y))
(2)
给定量化步长Qa，量化后的系数为：
cq=round(c/Qa)
(3)
接下来会对量化后的系数进行熵编码，等待每个块的Y分量都编码完毕后，再对U,V分量进行相同的编码过程。
3.2 交流分量熵编码
由于交流系数的幅值较小，量化后0值将会占很大的比重。为了更直观的说明，我们对不同PSNR下的YUV三通道非零系数百分比的均值进行了统计。如图3所示，4个测试模型在PSNR较高的情况下，非零系数的比重依然没超过25%，比较适合用Run-Level编码方法来对量化后的交流系数进行编码。为此，我们将量化系数表示为2个一维的阵列，分别记为Run和Level。由于每个块的非零系数的个数不同，本文引入NumNonzeroCoef来记录每个块中的非零系数的个数。对于第i个块，只有当NumNonzeroCoef[i]不等于0的时候才会对(runi，leveli)进行编码。图4描述了交流系数熵编码的主要过程。对于量化后的系数ci，如果不含非零系数，则NumNonzeroCoef[i]为0，表示没有系数需要编码。如果ci中有Ni个非零系数，先把Ni存储在变量NumNonzeroCoef [i]中，然后存储runi和leveli到Run{i}和Level{i}中。当所有的块都处理完毕，进一步采用霍夫曼编码[14]对NumNonzeroCoef, Run和Level进行熵编码。
图3 4个测试模型PSNR对比非零系数百分比曲线
Fig.3 Curves of PSNR vs percentage of nonzero coefficients for four test models
图4 交流系数熵编码流程图
Fig.4 The flowchart of AC coefficients entropy coding
3.3 直流分量熵编码
考虑到点云的相邻块在空间上存在相关性，本文采用类似于H.264帧内编码的方法对直流系数进行预测编码[15]。在H.264中，当前像素可由其相邻像素预测得到，然后对残差进行编码。同样的，当前块的直流系数可以由其相邻的被占用的块进行预测得到，然后对残差进行量化后编码。假设一个点云被分为m×m×m个块，每个块的空间三维坐标为(xi,yi,zi)，其中1≤xi≤m，1≤yi≤m，1≤zi≤m。我们按照块坐标从小到大的顺序，并且从x,y,z三个维度依次对每个被占用的块进行预测编码。以x这个维度为例，具体的编码步骤如下：
1)对被占用的块的x坐标进行扫描，找出x坐标为1的块，直接量化，量化步长为Qd。
2)对于2≤x≤m的占用块，按照x从小到大的顺序逐级进行预测编码。对于当前块(xi,yi,zi)，我们对其上一级的7个相邻的块进行扫描。所谓相邻块，假设空间中有一个2×2×2的立方体，把当前块放置在该立方体的右下角，则其余剩下的7个块即为当前块的相邻块。如果这7个块都没有被占用，则对当前块直接进行量化，量化后的值多为当前块的预测残差；如果这7个相邻块存在被占用的块，则计算这些被占用块的均值，与当前块的直流值做差，对残差进行量化。
3)采用霍夫曼编码器对量化后的预测残差进行熵编码。
4 实验结果与分析
本文使用的实验数据是从动态点云数据库中提取的某一帧[16]。如图5所示，这四个上半身模型分别是(a)Ricardo, (b)Sarah, (c)Phil, (d)Andrew,每个模型的点数量分别是207077，286934，325077和301336。这些点云数据均经过体素化处理，由实时的高分辨率稀疏体素算法[17](Real-time High Resolution Sparse Voxelization algorithm)采集得到。采用深度L为9的体素化处理方法将产生一个512×512×512的体素空间，进一步对这个体素空间进行划分，产生64×64×64个8×8×8块。在实验中，Qd设置为2，Qa的范围是5到45，每个被占用体素(bpv)所需要的比特数来表示码率，Y分量的峰值信噪比(PSNR, 单位 dB)来表示失真程度。
图5 测试点云数据
Fig.5 The test point cloud data
本文将所提算法与GT, RAHT, GPT点云压缩算法进行比较。本文所提算法和GT算法的dmax值均为1。图6展示了4种算法RD曲线的实验结果。由图可见：(1)在Ricardo，Sarah和Andrew这个点云模型上，我们的算法始终比GT， RAHT和GPT的性能要好，尤其是在Ricardo和Sarah这两个模型上。主要是由于GT，RAHT，GPT这三个算法都是在所有系数都服从拉普拉斯分布的假设前提下对所有系数都进行了编码，然而这几个模型的非零系数比重比较小，采用这种编码方法很大程度上局限了压缩效率。而本文采用Run-Level的方法可以较好地节约比特消耗。(2)在Phil这个模型上，我们的算法比GT和RAHT的性能要好，基本取得与GPT相当的压缩性能。这主要原因是Phil这个模型的非零系数的占比较大，run占用的比特也迅速增加。针对这个问题，我们可以通过加大交流系数的量化步长来降低非零系数的比重，同时调节直流系数的量化步长来确保获得较好的率失真性能。(3)在所有的实验中，RAHT的压缩性能皆低于其他三个算法，这主要是因为RAHT是一种基于Harr变换的低复杂度算法，不会根据信号的统计特性自适应调整。
图6 Ricardo, Andrew, Phil, Sarah的RD曲线
Fig.6 RD curves for Ricardo, Andrew, Phil and Sarah
为了更直观对点云压缩结果进行对比，以测试点云Andrew为例，我们给出了在不同失真条件下重构后的点云主观效果图。如图7所示，左边是原始点云，没有经过压缩处理；中间图是压缩后重构的点云，其bpv为1.43，PSNR为36.0134 dB;右边图也是压缩后重构的点云，bpv为0.5626，PSNR为30.5543 dB。由图可见，本文所提算法能够在保证重建点云图质量的前提下实现高效的点云压缩。
图7 重建点云的主观比较图
Fig.7 Subjective comparison of reconstructed 3D point cloud
5 结论
本文提出改进图形变换以实现3D点云的高效压缩。本文主要贡献在于：(1)通过改进图形构造，将每个块中的所有子图连接为一个图，从源头减少直流系数个数；(2)引入块平均值降低了直流系数的幅值，并对去平均的颜色特征进行图形变换。(3)在熵编码部分，采用Run-Level和预测编码方法分别对交流和直流系数进行编码，提高了熵编码效率。实验结果表明，本文所提算法相比于现有3D点云压缩算法具有更高的压缩效率。
[1] Zhang C, Cai Q, Chou P A, et al. Viewport: A distributed, immersive teleconferencing system with infrared dot pattern[J]. IEEE MultiMedia, 2013, 20(1): 17-27.
[2] Zhang C, Florêncio D, Loop C. Point cloud attribute compression with graph transform[C]∥Proceedings of IEEE International Conference on Image Processing, 2014: 2066-2070.
[3] Mekuria R, Blom K, Cesar P. Design, implementation and evaluation of a point cloud codec for teleimmersive video[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2017, 27(4): 828- 842.
[4] Thanou D, Chou P A, Frossard P. Graph-based motion estimation and compensation for dynamic 3D point cloud compression[C]∥Proceedings of IEEE International Conference on Image Processing, 2015: 3235-3239.
[5] Cohen R, Tian D, Vetro A. Point cloud attribute compression using 3-d intra prediction and shapeadaptive transforms[C]∥Proceedings of Data Compression Conference (DCC), 2016: 141-150.
[6] Schnabel R, Klein R. Octree-based point-cloud compression[C]∥Proceedings of Eurographics Symposium on Point-Based Graphics, 2006: 111-121.
[7] Queiroz R L D, Chou P A. Compression of 3D point clouds using a region-adaptive hierarchical transform[J]. IEEE Transactions on Image Processing, 2016, 25(8): 3947-3956.
[8] Sikora T, Makai B. Shape-adaptive dct for generic coding of video[J]. IEEE Transactions on Circuits and Systems for Video Technology, 1995, 5(1): 59- 62.
[9] Queiroz R L D, Chou P A. Transform Coding for Point Clouds Using a Gaussian Process Model[J]. IEEE Transactions on Image Processing, 2017, 26(7): 3507-3517.
[10] Shuman D I, Narang S K, Frossard P, et al. The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[J]. Signal Processing Magazine, IEEE, 2013, 30(3): 83-98.
[11] Hu W, Cheung G, Li X, et al. Depth map compression using multi-resolution graph-based transform for depth-image-based rendering[C]∥Proceedings of IEEE International Conference on Image Processing, 2012: 2066-2070.1297-1300.
[12] Zhang C, Florencio D. Analyzing the optimality of predictive transform coding using graph-based models[J]. IEEE Signal Processing Letters, 2013, 20(1): 106-109.
[13] Hou J H, Chau L P, He Y, et al. Sparse representation for colors of 3D point cloud via virtual adaptive sampling[C]∥Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017: 2926-2930.
[14] Skretting K, Husøy J H, Aase S O. Improved huffman coding using recursive splitting[C]∥Proceedings of Norwegian Signal Processing, 1999.
[15] 胡伟军, 李克非, 成建波. H.264/AVC视频编码标准的技术特点及应用分析[J]. 信号处理, 2005, 21(1): 79- 85.
Hu Weijun, Li Kefei, Cheng Jianbo. Features and Applications of H.264/AVC[J]. Signal Processing, 2005, 21(1): 79- 85.(in Chinese)
[16] Loop C, Cai Q, Escolano S O, et al. Microsoft voxelized upper bodies-a voxelized point cloud dataset[C]∥ISO/IEC JTC1/SC29 Joint WG11/WG1 (MPEG/JPEG) Input Document m38673/M72012, 2016.
[17] Loop C, Zhang C, Zhang Z. Real-time high-resolution sparse voxelization with application to image-based modeling[C]∥Proceedings of the 5th HighPerformance Graphics Conference. ACM, 2013: 73-79.
