基于张量特征分解的三维网格模型降噪

， ，

(浙江理工大学理学院，杭州 310018)

0 引 言

近年来，随着数字扫描仪的普遍化使用，人们可以获取到各种各样的三维模型数据。然而由于人为因素或者扫描仪的缺陷，可能导致获取的三维数据带有噪声，而噪声数据往往会影响后续相关数字几何处理。因此三维模型数据降噪是数字几何处理中的重要步骤。

关于三维模型降噪的方法已有很多[1-2]。第一类是基于滤波技术的三维网格降噪方法。例如Taubin[3]提出了一种表面信号低通滤波算法，该算法可以应用于具有拓扑结构的三维网格模型，与其他一些耗时的光顺降噪算法不同的是，该算法拥有较低的时间复杂度以及空间复杂度；Fleishman等[4]受图像双边滤波算法的启发，提出了基于三维网格模型的双边滤波降噪算法，该算法从三维网格中滤除了噪声，并成功地保留了三维模型原有的几何特征。

第二类是基于法向信息的三维网格降噪方法。如Yu等[5]提出了一种新颖有效的两阶段特征保持网格降噪方法，使用二次误差测度(Quadric error mactrics，QEM)[6]技术提高顶点更新时的特征保持能力，QEM技术能够实时获取需要更新的顶点并减小顶点更新时的误差，但是该方法不能精确处理和恢复浅层特征；Wei等[7]基于双向法向域的分段一致性，设计出一种双向法向过滤算法，解决了特征区域中噪声网格顶点法向的估值问题。

第三类是于曲面重建的三维网格降噪方法。如Wang等[8]提出了一种通用性好、鲁棒性强的网格降噪算法，该方法将双边滤波、特征识别、各向异性邻域搜索、曲面拟合等技术混合在一起；Huang等[9]提出了一种基于重采样的网格降噪方法，采用边缘感知方法对三维模型上的孤立点集进行降噪处理，但是该方法对于边界开放的三维模型，降噪效果并不理想；Wang等[10]首先采用全局的拉普拉斯降噪方法得到一个预估计的基平面，然后提出了一种压缩感知优化方案，用来保持以及恢复模型的形状特征，但是该方法是基于高斯白噪声的，对于高斯噪声、椒盐噪声、瑞利噪声等其他噪声模型，降噪效果并不理想。

虽然网格降噪方法已有很多，但是目前对三维模型特征保持的降噪研究仍然存在较多不足。本文为较好地保持模型特征，提出了一种基于张量的三维模型降噪方法，首先为三维模型上的每个顶点找到高斯曲率、形状指数、形状直径函数三个几何特征值；其次将这三个几何特征进行协同处理，构造三阶特征张量，并通过张量范数搜索给定区域的局部邻域相似匹配块；然后根据相似块内顶点法向量构造四阶张量，并通过张量分解，对法向进行修正；最后将顶点沿着修正后的法向方向进行相应偏移，实现降噪效果。本文提出的三维模型降噪方法可以应用于三维模型重建、模型分割、特征提取等领域中。

1 形状相似性分析

为了解决张量分解中子模型块的相似块匹配问题，本文给出了三维模型相似性分析方法，首先采用Golovinskiy等[11]提出的Randomizedcuts算法，将三维模型分割成若干个块，并根据每个分割块计算出GS描述符、SI描述符、SDF描述符；然后通过这三个三维模型的几何特征构建了本文的张量描述符；最后，通过张量的二范数进行形状相似性分析。

1.1 常用几何特征属性估算

现有的三维模型主要有网格和点云两种表示形式，三维模型的常用几何特性主要包括曲率、形状指数、形状直径函数等。本小节描述了这些几何特征在三维模型中的应用。

1.1.1 三维模型的曲率估算

三维模型表面上的曲率用来反映该点的弯曲程度。常用的曲率形式包括主曲率、平均曲率和高斯曲率等。主曲率描述的是曲面上所有经过给定点的曲线曲率的最值，主曲率可分为最大主曲率和最小主曲率，反映了当前顶点的平均弯曲程度。平均曲率通过求取当前点的最大主曲率和最小主曲率的平均值来计算。高斯曲率在数值上等于该点的两个主曲率的乘积，它可以描述曲面在一点处的弯曲程度。针对三维网格模型，本文使用Meye等[12]提出的Voronoi方法来估算模型的高斯曲率，得到Bunny模型的高斯曲率图如图1(b)所示，其中：深色表示高斯曲率值处于中等水平的点，浅色表示高斯曲率值较大的部分点。

图1 Bunny网格的曲率颜色映射图

1.1.2 形状指数

Bradford等[13]提出的形状指数客观地反映了三维模型表面任意曲面的形状，可作为刻画三维模型表面凹凸程度的参数。针对三维网格模型，本文采用Bradford提出的形状指数估算Bunny模型表面的形状指数，得到Bunny模型的形状指数颜色映射图如图2所示，其中：浅色块表示模型中较凹的部分，深色块表示模型中较凸的部分。

图2 Bunny模型形状指数颜色映射图

1.1.3 网格形状直径函数

Shapira等[14]提出的形状直径函数(Shape diameter function，SDF)反映了三维模型的体积与模型边界之间的联系，可以作为描述三维模型局部特征的参数。针对网格模型，本文采用Shapira方法求取Bunny网格模型各个点的形状直径函数，得到Bunny模型形状直径函数颜色映射图如图3所示，其中：浅色表示形状直径函数值较小部分；深色表示形状直径函数值介于两点之间的曲面块。

图3 Bunny模型形状直径函数颜色映射图

1.2 基于Randomizedcuts的过分割

本文去噪算法是用来处理三维模型的子模型块，因此本节采用描述子模型块的构建方法。首先根据Golovinskiy等[11]提出的Randomizedcuts算法将整个网格模型划分为若干分割块，再通过相似性匹配算法分析这些分割块的相似性。Golovinskiy把网格模型构造成一个完整的无向图(其中：网格面片抽象成图节点；面片与面片间的边作为图节点间的边)，并为图节点的每条边构造两个权值代价(遍历代价和分割代价)；然后，用k-means聚类方法将图聚类n次，同时，统计n次分割中每条边成为分割边界的概率；最后，找出概率最大且相等的边作为最后的分割线。Bunny网格模型的过分割图如图4所示，分割块为200，黑线表示分割线。

图4 Bunny网格模型的过分割图

1.3 基于张量描述符的相似块匹配

张量是向量、矩阵在维度上的扩展，是描述客观存在的物理量，具有坐标不变性，并且在不同的参考坐标下有不同的分量。张量技术已经得到广泛运用，如：Shen等[15]将图像子块及其相似块构造成四阶张量，并利用张量分解去除图像中高频的噪声数据；Zhang等[16]利用张量对光谱进行分析。本文使用张量对网格模型的分割面片进行相似性匹配。

1.3.1 张量的构造

通过采用1.1节描述的几何特征的求解方式，本文提取了三维模型的特征描述符，即GS描述符、SI描述符和SDF描述符，并且利用三维网格模型的所有点的GS特征、SI特征以及SDF特征分别构成一个n×k维矩阵(n是模型顶点数，k是顶点以及顶点邻域个数总和)。首先，为分割块内的每一个顶点找到最近k-1个邻域点；接着，以当前顶点作为第一个元素，其他邻域点按照距离递增(相对当前顶点距离)构造一个k维列向量：

GSj=(gs1，gs2，…，gsk)T

(1)

其中：GSj表示分割块内第j个点及其周围k-1个点的GS描述符构成的列向量；gs1表示当前顶点的高斯曲率值；gs2，gs3，…，gsk是当前顶点的k-1个邻域点的高斯曲率值。同理可求得分割块内所有顶点的高斯曲率(GS)列向量、形状指数(SI)列向量、形状直径函数(SDF)列向量。然后，将得到的特征列向量(1)按照顶点个数扩展成特征矩阵

GS=(GS1,GS2,…,GSn)

(2)

其中：GSj表示每个子模型块的高斯曲率矩阵；GS1，GS2，…，GSn顺序是按照每个列向量第一个元素的高斯曲率值排序的。最后，将高斯曲率(GS)、形状指数(SI)、形状直径函数(SDF)的k×n维矩阵构造成一个k×n×3的三阶张量。

1.3.2 张量二范数分析形状相似性

由于不能保证模型的每个分割块顶点数一致，本文采用PCA进行降维处理。将顶点个数多的张量降维成顶点个数少的张量；然后再对分割块计算张量范数(如二范数)判定给定区域的相似块。Bunny、人马模型相似块匹配图如图5所示，其中：浅色块作为参考块，深色块代表参考块局部邻域的相似块。

图5 Bunny、人马模型相似块匹配

2 基于相似块和张量奇异值分解的三维模型降噪

2.1 高阶奇异值分解在图像上的运用

首先，为大小为m×m的图像子块寻找k-1个相似块；然后，将图像子块与k-1个相似块构成一个m×m×3×k的四阶张量T∈Rm×m×3×k。并对张量T进行分解，张量T的高阶奇异值分解为：

T=S×1U(1)×2U(2)×3U(3)×4U(4)

(3)

(4)

(5)

本文将基于张量方法的图像降噪方法应用到三维模型降噪。但是直接应用到三维模型降噪上会存在以下问题:a) 三维模型的顶点顺序并不是像图像那样规则排列，对于张量的构造以及相似块的查找带来一定的困难；b) 图像降噪的处理数据是R，G，B颜色信息，降噪的本质就是对一个像素点的R，G，B信息进行微调，所以三维模型需要找到一个位置偏移的量，运用降噪方法对这个位置偏移量进行微调。

2.2 三维子模型块匹配

本文的三维子模型块是以顶点以及顶点的k-1个邻域点组成的，所以对于三维模型块的相似性匹配，可看作顶点的k邻域相似性查找问题。首先，把每个顶点和顶点的k-1个邻域点的GS描述符、SI描述符、SDF描述符构造成一个k维列向量；然后，将得到的特征列向量按照顶点个数扩展成特征矩阵；最后，将GS、SI、SDF特征的k×n维矩阵构造成一个k×n×3的三阶张量，并两两比较子模型块张量二范数的值，值越小表示子模型块越相似。

2.3 三维模型张量构造

基于张量的图像降噪，是利用张量分解将图像的高频部分去除，保留图像中的低频部分，张量分解是对图像的R，G，B值进行处理。而三维模型中却不存在类似R，G，B的值，但是可以通过修正噪声点的法向方向，完成降噪过程。首先，以子模型块内点的法向量以及m-1个相似块为基础，构造一个k×k×3×m(第一个k表示矩阵中每个列向量第一个元素有k-1个最近点，第二个k表示一个子模型块中有k个顶点，3表示每一个顶点的法向量存在x，y，z轴三个分量，m表示子模型块及其相似块的个数)的四阶张量；然后，对张量进行分解，修正顶点的法向量；最后，在每个顶点以及顶点的2邻域构造一个最小二乘平面，求得顶点到平面的距离d，并将顶点沿着修正后的法向量偏移d距离，即得到降噪后的顶点坐标。

Shen等[15]利用图像非局部相似特性，在整个图像范围内，为噪声图像子块匹配若干个相似块，通过对该图像子块及其相似块进行协同滤波，可得到较好的降噪效果。本文为三维模型上的所有顶点寻找最近k-1个顶点，并以顶点以及顶点的最近k-1个点作为一个三维子模型块。三维子模型块具有无序性，还需要进一步处理。首先，把所有顶点的最近k-1个点按照距离递增排序，并将每个顶点以及顶点最近k-1个点的法向量构造成k维列向量：

Bix=(n1x，n2x，…，nkx)T

(6)

其中：Bix由子模型块内某一顶点及其最近k-1个点法向量的x轴分量构成。同理可求出分割块内其他顶点的Bix、Biy、Biz三轴列向量。将子模型块内所有顶点的Bix，按照Bix中第一个元素对应顶点的高斯曲率进行排序，可得到：

Mjx=(B1x，B2x，…，Bkx)

(7)

其中：B1x，B2x，…，Bkx是式(6)求得的列向量；Mjx内元素顺序是按照每个Bix第一个元素对应顶点高斯曲率由小到大的顺序；Mjx代表第j个分割块内所有顶点法向量沿x轴分量构成的法向量矩阵。同理可求得Mjy，Mjz。然后，将参考块的Mjx、Mjy、Mjz(j为参考块序号)以及p个相似块的Mjx、Mjy、Mjz构造成k×k×3×p的四阶张量，并对张量T进行张量分解，得到修正后的顶点法向量；最后，用最小二乘法拟合出每个顶点的最近2邻域所构成的平面，将顶点距离这个平面的距离记为d，则

(8)

其中：xi表示噪声数据；yi表示降噪后的数据；ni表示该顶点的修正后的法向量。如果顶点的法向方向与平面法向同向，则n取偶数；如果反向，则n取奇数。

图6为fandisk模型降噪示例图像，如图中所示，fandisk噪声模型以网格的形式显示，其中：深色块内点为任意一顶点Vi，Vi周围深色区域是以Vi为中心的邻域；浅色块内的深色点为相似块的参考点，浅色块为深色块的相似块。具体步骤如下：

步骤1：为图6(a)的参考块查找相似块，得到图6(b)中的浅色块即为相似块。

步骤2：通过将参考块和相似块的法向量构造成4阶张量，并对张量进行分解，得到修正后的参考块法向量；

步骤3：令参考块内的点沿着各自修正后的法向量进行偏移，得到图6(c)降噪后的局部效果图。其中：图6(d)为整个fandisk模型的全局降噪效果图。

图6 fandisk模型降噪示例图像

3 实验分析与讨论

本文实验中计算机软硬件配置为：处理器是Intel(R) Core(TM) i7-4510U CPU@2.00 GHz 2.60 GHz，内存为6 G，操作系统是Windows7系统，计算软件是Visual Studio C++。本文实验的三维模型数据来源于三维公园(http://threepark.net/geometryplusplus/models)，在这个三维模型数据库中下载对应obj文件。

根据本文降噪方法可知，相似块个数取值会对实验效果产生影响。相似块少，降噪效果不明显；相似块个数过多，会导致彼此不相似的块也参与了参考块的法向修正过程，从而影响整个降噪效果。本文增加了对shark1、lamp1、giraffe、dancingChildren等模型相似块分别取5、10、15、20、25、30、35、40、45的实验结果，并对每个降噪模型的顶点位置作方差统计，如表1所示。通过表1发现当模型顶点数目较小(如6000)时，子块取20时本文的降噪效果较好；当模型顶点数目较大(如15000)时，子块取30时降噪效果较好。综合本文方法的降噪效果和运行时间，选择相似块个数为20作为本文降噪方法的相似块个数。

为了研究相似块对模型降噪的影响，本文统计了各个模型对于不同相似块的时间开销，如表2所示。

表1 τ=0.25时降噪前后模型的方差值

注：shark模型顶点个数5000个，面片个数9992个；lamp模型顶点个数8000个，面片个数17726个；giraffe模型顶点个数15000个，面片个数35302个；dancingChildren顶点个数20000个，面片个数42000个。

表2 τ=0.25时不同相似块降噪时间开销 s

注：shark模型顶点个数5000个，面片个数9992个；lamp模型顶点个数8000个，面片个数17726个；giraffe模型顶点个数15000个，面片个数35302个；children顶点个数20000个，面片个数42000个。

为了进一步说明本文相似块个数会对降噪效果产生影响，本文列出了shark1、lamp1、giraffe、children等模型分别取10个、20个相似块的实验结果，如图7所示。

然后，将本文降噪方法与Wei等[17]、Wang等[10]和Zheng等[18]方法进行对比，结果如图8所示。从图8中可以看出，Wei等[17]方法影响了边界的光顺结果，并留有少许噪点；Wang等[10]方法丢失了尖锐点的特征；zheng等[18]方法导致模型边界收缩，丢失了边界特征点。而本文方法降噪后的模型边界较为光滑，并且能够还原尖锐点的特征。

为了验证本文方法能够在处理噪点的同时做好特征点的保持工作，本文在planck模型上与其他三种降噪算法进行了对比，如图9所示，从原模型可以看出，planck模型鼻尖呈圆状，而Wei等[17]方法降噪后显示，planck模型鼻尖有一丝扁平状；zheng等[18]方法降噪后，模型鼻尖上移，整个鼻子收缩变小；而本文方法降噪后的模型能够很好的保持特征点。

本文分别统计了本文方法、Wei方法、Wang方法、Zheng方法在fandisk模型以及planck模型上的降噪运行时间，如表3所示。在表3中，可发现本文方法拥有较快的运行处理速度。

为了定量分析本文方法、Wei方法、Wang方法、Zheng等方法的降噪效果，本文增加了各个方法对原模型与降噪后模型的顶点位置作方差对比，计算结果如表4所示。

图7 本文方法降噪效果

图8 fandisk模型降噪效果比较

图9 planck模型降噪效果比较

注：fandisk模型顶点个数6475个，面片个数12946个；planck模型顶点个数5000个，面片个数9884个。

表4 本文方法降噪方差与其他降噪算法方差对比(τ=0.25，m=20)

注：fandisk模型顶点个数6475个，面片个数12946个；planck模型顶点个数5000个，面片个数9884个。

最后，为了说明本文降噪方法的鲁棒性，还分别对其他一些模型进行了降噪实验，结果如10图所示。实验表明，本文降噪算法可以较好的保持三维模型的特征。

图10 本文降噪模型效果

4 结 语

本文根据张量分解的滤波特性提出了一种三维模型降噪算法。主要通过模型几何属性获取子模型块相似块，对子模型块及其相似块的法向量进行张量处理；同时，通过顶点法向量调整，修正当前点的位置。基于张量的网格降噪过程中，现有方法主要通过修正法向方向进行降噪，法向偏移没有进行优化处理。如何设置合理的法向偏移量得到更有效的降噪结果将是今后的主要工作。另外本文现有方法只适用于三维网格模型，如何将该方法应用到三维点云模型的降噪也是今后的主要工作。

参考文献：

[1] 李楠楠，曹俊杰，李波，等.基于面法向规范化的重加权全局双边滤波[J].计算机辅助设计与图形学学报，2014，26(3):370-377．

[2] 江亮亮，杨付正，任光亮.用于网格降噪的自适应双边滤波器[J].华南理工大学学报(自然科学版)，2015，43(11):54-60．

[3] Taubin G. A signal processing approach to fair surface design[C]//Conference on Computer Graphics and Interactive Techniques. ACM，1995:351-358.

[4] Fleishman S， Drori I， Cohen-or D. Bilateral mesh denoising[J]. Acm Transactions on Graphics，2003，22(3):950-953．

[5] Yu J， Wei M， Qin J， et al. Feature-preserving mesh denoising via normal guided quadric error metrics[J]. Optics & Lasers in Engineering，2014，62(6):57-68．

[6] 韦虎，张丽艳，刘胜兰，等.基于支撑域的网格简化算法[J].中国图象图形学报，2011，16(5):892-897．

[7] Wei M， Shen W， Qin J， et al. Feature-preserving optimization for noisy mesh using joint bilateral filter and constrained Laplacian smoothing[J]. Optics & Lasers in Engineering，2013，51(11):1223-1234．

[8] Wang J， Zhang X， Yu Z. A cascaded approach for feature-preserving surface mesh denoising[J]. Computer-Aided Design，2012，44(7):597-610．

[9] Huang H， Wu S， Gong M， et al. Edge-aware point set resampling[J]. Acm Transactions on Graphics，2013，32(1):1-12．

[10] Wang R， Yang Z， Liu L， et al. Decoupling noise and features via weighted1-analysis compressed sensing[J]. Acm Transactions on Graphics，2014，33(2):1-12．

[11] Golovinskiy A， Funkhouser T. Randomized cuts for 3D mesh analysis[J]. Acm Transactions on Graphics，2008，27(5):145-156．

[12] Meyer M， Desbrun M， Schröder P， et al. Discrete differential-geometry operators for triangulated 2-manifolds[J]. Mathematics & Visualization，2002，6(8-9):35-57．

[13] Bradford J， Westhead D. Improved prediction of protein-protein binding sites using a support vector machines approach[J]. Bioinformatics，2005，21(8):1487-1494．

[14] Shapira L，Shamir A，Cohen-or D. Consistent mesh partitioning and skeletonisation using the shape diameter function[J]. The Visual Computer，2008，24(4):249-258．

[15] Shen Y， Han B， Braverman E. Adaptive frame-based color image denoising[J]. Applied and Computational Harmonic Analysis，2016，41(1):54-74．

[16] Zhang Q， Wang H， Plemmons R， et al. Tensor methods for hyperspectral data analysis: a space object material identification study[J]. Journal of the Optical Society of America A Optics Image Science & Vision，2008，25(12):3001-3012．

[17] Wei M， Yu J， Pang W， et al. Bi-normal filtering for mesh denoising[J]. IEEE Transactions on Visualization & Computer Graphics，2015，21(1):43-55．

[18] Zheng Y， Fu H， Au K， et al. Bilateral normal filtering for mesh denoising[J]. IEEE Transactions on Visualization & Computer Graphics，2011，17(10):1521-1530．