CN108447038A - 一种基于非局部全变分算子的网格去噪方法 - Google Patents

Abstract

本发明公开了一种基于非局部全变分算子的网格去噪方法，主要包括获取待测对象的三角网格测量数据；分别为每条边在包含它的两个三角面上指定方向；计算每条边的非局部全变分算子；通过求解目标公式得到整个三角网格滤波后的面法向量；更新三角网格中所有的顶点坐标，更新后的坐标即为滤波后的顶点坐标。本发明将传统的二维图像去噪的全变分和非局部方法扩展到三角网格上，在网格去噪中考虑到非局部性和几何位置信息，通过定义非局部全变分算子，最终可以简单高效地达到保持特征的网格去噪的效果。

Description

一种基于非局部全变分算子的网格去噪方法

技术领域

本发明属于计算机图形学领域，涉及一种基于非局部全变分算子的网格去噪方法，特别适用于对三维数据采集过程中的噪声和扰动，特征模糊、丢失等进行处理。

背景技术

三角网格是计算机图形学、计算机视觉以及虚拟现实等领域中几何对象的基本表示模型。一般情况下，三角网格曲面是由扫描设备获取数据后，使用三角化方法对数据进行三角化得到的，其结构如图2所示。在获取网格曲面数据的过程中，由于采集设备、环境等因素的影响，所获取的网格曲面，不可避免地存在着误差。所以，大部分网格曲面均含有不同程度的噪声，这不仅无法客观反映数据中的真实信息，还对后续处理与计算造成影响。因此，三角网格曲面的去噪对获取高质量的网格数据是非常必要的。

全变分(total variation，TV)模型是上世纪九十年代出现的一种图像处理方法，该方法通过最小化全变分能量函数实现对图像保持特征的去噪。由于全变分能量函数能够刻画自然图像良好的保持特征的性质，全变分模型及其推广形式已经被广泛地应用在图像恢复中。在图像恢复中，已经证明全变分模型能够较好地保持图像的特征，比如尖锐的边、点等。

非局部模型是另一类有效的图像处理方法，该方法充分利用了图像内容的自相似性，融合了非局部块匹配思想，基于图像块之间的欧氏距离定义了像素点之间的相似度权重，然后通过像素点间的加权平均达到去噪的目的，考虑到局部邻域和非局部邻域之间的相似性，提升了图像去噪的效果。

上述两种方法适用于处理二维图像去噪问题，二维图像的噪声通常是图像像素块的灰度值产生了误差，而在三角网格中噪声是指通过设备采集后的物体的坐标偏离了真实的坐标，这两种噪声是有区别的。二维图像中的像素块是规则排列的，但是三角网格结构是不规则的，本身的数据存在很大差异，所以，图像去噪的方法不适合直接作用在三角网格上，需要重新对方法进行扩展。

发明内容

本发明针对解决二维图像去噪方法如何应用到三角网格去噪的问题，将全变分和非局部的图像去噪方法从像素块规则排列的二维图像重新扩展到不规则的三角网格去噪，最终达到网格去噪的效果，提出了一种基于非局部全变分算子的网格去噪方法,适用于解决封闭网格去噪问题。

三角网格是全部由三角形组成的多边形网格，用来模拟复杂物体的表面。三角网格数据的获取是通过数据采集设备对真实世界的物体扫描，按照三角划分方法得到。三角网格数据的常用表示形式为其中，和分别表示三角网格中的顶点、边和三角面的集合，v_i，e_i和τ_i分别表示第i个顶点、第i条边和第i个三角面，V表示顶点的个数，E表示边的个数，T表示三角面的个数。是第i个三角面τ_i的带噪面法向量，即测量值，是第i个三角面τ_i的待求解面法向量，即滤波值。

待测对象三角网格中所有的三角面都有对应的法向量，然而通过扫描设备获得的每个三角面的法向量都是带有噪声的，本发明的主要思想就是首先对进行滤波，得到滤波后的然后在第二个阶段根据滤波后的调整顶点的位置，最终得到的顶点位置是求解的目标值。

网格噪声的本质是由于采集后的三角网格的顶点坐标和真实的顶点坐标发生了偏离，去噪就是把顶点恢复到真实的位置上。

本发明采用了如下的技术方案和实现步骤：

步骤(1)：获取待测对象的三角网格测量数据其中，和分别表示三角网格中的顶点、边和三角面的集合，v_i，e_i和τ_i分别表示第i个顶点、第i条边和第i个三角面，V表示整个三角网格中顶点的个数，E表示整个三角网格中边的个数，T表示整个三角网格中三角面的个数；

步骤(2)：指定三角面三个顶点的读取顺序，指定每条边两个顶点的读取顺序，分别为每条边在包含它的两个三角面上指定方向，其中第j条边e_j在包含第j条边e_j的一个三角面上的方向用sgn(e_j,τ_i)表示，具体如下：

如图3所示，左侧三角面三个顶点的读取顺序是顺时针方向，为边AB指定的两个顶点的读取顺序是由A到B，则顶点A、B在左侧三角面中的读取顺序和在边中的读取顺序一致，因此边AB在左侧三角面上的方向为1。同理，边AB在右侧三角面上的方向为-1。

步骤(3)：计算每条边的非局部全变分算子,具体如下：

步骤(3.1)：求出每条边两侧的三角面的邻域的索引集合，其中第j条边两侧的三角面分别为τ_l,τ_r，三角面τ_l的领域的索引集合由三角网格上与三角面τ_l的质心距离不超过给定阈值d的三角面的索引组成，记为三角面τ_r的领域的索引集合由三角网格上与三角面τ_l的质心距离不超过给定阈值d的三角面的索引组成，记为具体表示如下：

其中，T表示整个三角网格中所有三角面的个数，d为阈值，d∈[2σ,3σ]，σ是相邻三角面质心之间的距离的平均值，分别表示三角面τ_k,τ_l,τ_r的质心。如图4所示，三角面τ_l的邻域的索引集合由虚线圆内的三角面的索引组成，代表和三角面τ_l几何位置最接近的三角面集合。

步骤(3.2)：计算每条边两侧的三角面的非局部面法向量，其中第j条边两侧的三角面τ_l,τ_r的非局部面法向量分别为和和的计算方法一样，的计算公式如下：

其中，是第i个三角面τ_i的滤波后的面法向量。

其中，分别表示三角面τ_i,τ_l的质心；

σ_s是质心的高斯核函数的方差，取值范围：(0，2)；

分别表示三角面τ_i,τ_l的输入的带噪法向量；

σ_r是带噪面法向量的高斯核函数的方差，取值范围：(0，2)；

和的计算方法一样，分别表示三角面τ_i,τ_l的邻域集合内三角面包含的顶点的高斯曲率通过直方图统计后得到的特征向量，是把[f₁,...,f_i,...,f_n]标准化后的向量，n表示直方图中划分的区间总数，f_i表示内三角面包含的顶点落入直方图中第i个区间的个数；

σ_f是特征向量的高斯核函数的方差，取值范围：(0，10)；

计算非局部面法向量实质是对三角面邻域中的三角面的法向量进行加权平均，权重是三角面之间的相似度。

步骤(3.3)：计算三角网格中所有边的非局部全变分梯度算子其中第j条边的非局部全变分梯度算子记为具体计算公式如下：

这里τ_l，τ_r是指共享边e_j的两个三角面；

步骤(4)：通过求解目标公式得到整个三角网格滤波后的面法向量目标公式具体如下:

其中，全变分项计算公式为:

其中，E表示整个三角网格中边的个数，||·||₂是2-范数，是指共享边e_j的两个三角面的带噪法向量，length(e_j)是指边e_j的长度，是步骤2中得到的第j条边的非局部全变分梯度算子的模长，

误差项具体公式如下：

其中，Nⁱⁿ是三角面带噪面法向量，T是整个三角网格中三角面的个数，是第i个三角面的滤波后的面法向量和第i个三角面的带噪法向量差的2-范数的平方，是采集后得到的三角面的法向量，area(τ_i)是第i个三角面的面积，α＞0，用来平衡两项。

全变分项实质是对步骤(3)中定义的非局部全变分梯度算子的一个加权平均。误差项是为了使得测量值Nⁱⁿ和滤波后的N不会偏离太多。滤波后的面法向量N还是和三角面垂直的法向量，为了更新顶点求解的。

步骤(5)：更新三角网格中所有的顶点坐标，更新后的坐标即为滤波后的顶点坐标，其中顶点i的更新后的坐标v′_i的计算公式如下：

其中，v_i是顶点i的测量坐标，F_v(i)表示顶点i的1-邻域，即有公共顶点i的三角形集合，|F_v(i)|表示构成1-邻域的三角面的个数，N_k表示第k个三角面滤波后的面法向量，c_k表示三角面k的质心。本步骤需要迭代计算，迭代次数:[5,20]。

最终，使三角网格的顶点恢复到真实的顶点坐标位置就达到了网格去噪的效果。

有益效果

本发明将传统的二维图像去噪的全变分和非局部方法扩展到三角网格上，在网格去噪中考虑到非局部性和几何位置信息，通过定义非局部全变分算子，最终可以简单高效地达到保持特征的网格去噪的效果。

附图说明

图1为本发明整体流程示意图；

图2为三角网格曲面数据结构；

图3为边在相邻两个三角面上方向的示意图；

图4为三角面在邻域的索引集合示意图；

图5a为原始的fandisk模型；

图5b为通过Kinect扫描采集到的fandisk模型；

图5c为双边滤波对fandisk模型网格去噪的效果图；

图5d为非迭代fandisk模型网格去噪的效果图；

图5e为L₁-压缩感知对fandisk模型网格去噪的效果图；

图5f本发明对扫描的fandisk模型网格去噪的效果图。

具体实施方式

本实验中所用实验数据是用Kinect摄影机对模型fandisk进行扫描，fandisk模型如图5a所示，该模型有6475个顶点，19419条边，12946个三角面。测量数据包含：顶点的三维坐标，边的两个顶点的索引号，三角面的三个顶点的索引号。

在MATLAB程序实现网格去噪的具体实施步骤如下：

步骤(1)：求出每条边在共享它的两个三角面的方向：

指定三角面三个顶点的索引号的读取顺序为顺时针方向，随机指定每条边两个顶点的索引号读取顺序并固定该顺序，分别为每条边在包含它的两个三角面上指定方向，其中第j条边e_j在包含第j条边e_j的一个三角面上的方向用sgn(e_j,τ_i)表示，具体如下：

每条边被两个三角面共享，所以有sgn(e_j,τ_l)和sgn(e_j,τ_r)，τ_l，τ_r是指每条边的两个三角面。

步骤(2)：计算每条边的非局部全变分梯度算子：

步骤(2.1)：找到每条边的两个三角面，记为τ_l,τ_r，然后，求τ_l,τ_r的 公式如下：

d＝0.24。

步骤(2.2)：计算每条边两侧的三角面的非局部面法向量：

对于上述步骤得到的每条边两侧三角面的邻域集合和计算相似度：

取σ_s＝0.35,σ_r＝0.3,σ_f＝5。和分别表示三角面τ_i,τ_l的邻域集合内三角面包含的顶点的高斯曲率通过直方图统计后得到的特征向量。比如，求首先计算出中顶点的高斯曲率，然后将直方图划分的区间数指定为100，f_i表示内三角面包含的顶点落入直方图中第i个区间的个数，是把[f₁,...,f_i,...,f₁₀₀]标准化后的向量。

计算的非局部面法向量，本质是对中三角面的法向量进行加权平均计算，权重是上述求得的相似度，所以有；

是非局部面法向量。同理，可计算出

步骤(2.3)：结合上述步骤得到的每条边的sgn(e_j,τ_i)和计算三角网格中所有边的非局部全变分梯度算子其中，第j条边的非局部全变分梯度算子记为具体计算公式如下：

这里τ_l，τ_r是指共享边e_j的两个三角面；

步骤(3)：求解目标公式：

步骤(3.1):计算全变分项:计算出每条边的边长，结合上述步骤(2.3)的对19419条边进行累加计算，得到：

其中，length(e_j)表示边长，表示非局部全变分向量的模长。

步骤(3.2):误差项:计算每个三角面的面积，把12946个面的面法向量的测量值和滤波后的值的差值的2-范数的平方进行加权平均，三角面的面积作为权重。

其中，area(τ_i)是指三角面τ_i的面积。

步骤(3.3)：求解目标公式：

取α＝1000，使用增广拉格朗日方法和变量分裂法求解滤波后的面法向量，然后对该法向量进行标准化得到

步骤(4)：根据上述步骤得到的滤波后的面法向量调整测量的顶点坐标的值:对于6475个点都进行以下计算：

其中，v_i是顶点i的测量坐标，F_v(i)表示顶点i的1-邻域，即有公共顶点i的三角形集合，|F_v(i)|表示构成1-邻域的三角面的个数，N_k表示第k个三角面滤波后的面法向量，c_k表示三角面k的质心。该步骤迭代计算35次，得到最终更新后的顶点坐标。

最终，使三角网格的顶点恢复到真实的顶点坐标位置就达到了网格去噪的效果。

为了说明本发明的效果，做了以下对比试验：

在三角网格模型上进行实验，并与双边滤波方法、非迭代方法、L₁-压缩感知方法进行比较。实验采用峰值信噪比(PSNR)进行评价，PSNR值越大，则去噪效果较好。

从图5a-图5f中我们可以看出，传统的网格去噪方法对网格去噪后，棱边或者尖锐特征处理不好，或者过于光滑，而本发明在进行网格去噪的同时较好地保持了尖锐特征。表1是不同网格用不同去噪方法得到的PSNR的对比，非局部全变方法不但完成了对三维模型的去噪，而且PSNR有明显的提升。

表1.不同去噪方法得到的PSNR

网格模型	双边滤波	非迭代	L1-压缩感知	非局部全变分
					fandisk	47.0403	52.4706	47.7518	56.9178

Claims (1)

1.一种基于非局部全变分算子的网格去噪方法，其特征在于包括以下步骤：

步骤(1)：获取待测对象的三角网格测量数据其中，和分别表示三角网格中的顶点、边和三角面的集合，v_i，e_i和τ_i分别表示第i个顶点、第i条边和第i个三角面，V表示整个三角网格中顶点的个数，E表示整个三角网格中边的个数，T表示整个三角网格中三角面的个数；

步骤(2)：指定三角面三个顶点的读取顺序，随机指定每条边两个顶点的读取顺序，分别为每条边在包含它的两个三角面上指定方向，其中第j条边e_j在包含第j条边e_j的一个三角面上的方向用sgn(e_j,τ_i)表示，具体如下：

步骤(3)：计算每条边的非局部全变分算子,具体如下：

步骤(3.1)：求出每条边两侧的三角面的邻域的索引集合，其中第j条边两侧的三角面分别为τ_l,τ_r，三角面τ_l的领域的索引集合由三角网格上与三角面τ_l的质心距离不超过给定阈值d的三角面的索引组成，记为三角面τ_r的领域的索引集合由三角网格上与三角面τ_l的质心距离不超过给定阈值d的三角面的索引组成，记为

步骤(3.2)：计算每条边两侧的三角面的非局部面法向量，其中第j条边两侧的三角面τ_l,τ_r的非局部面法向量分别为和 和的计算方法一样，的计算公式如下：

其中，是第i个三角面τ_i的滤波后的面法向量，

其中，分别表示三角面τ_i,τ_l的质心，

σ_s是质心的高斯核函数的方差，取值范围：(0，2)，

分别表示三角面τ_i,τ_l的带噪法向量，

σ_r是带噪面法向量的高斯核函数的方差，取值范围：(0，2)，

和的计算方法一样，分别表示三角面τ_i,τ_l的邻域集合内三角面包含的顶点的高斯曲率通过直方图统计后得到的特征向量，是把[f₁,...,f_i,...,f_n]标准化后的向量，n表示直方图中划分的区间总数，f_i表示内三角面包含的顶点落入直方图中第i个区间的个数；

σ_f是特征向量的高斯核函数的方差，取值范围：(0，2)；

步骤(3.3)：计算三角网格中所有边的非局部全变分梯度算子其中第j条边的非局部全变分梯度算子记为具体计算公式如下：

这里τ_l，τ_r是指共享边e_j的两个三角面；

步骤(4)：通过求解目标公式得到整个三角网格滤波后的面法向量目标公式具体如下:

其中，全变分项R_nltv(▽N^nltv)计算公式为:

其中，E表示整个三角网格中边的个数，||·||₂是2-范数，是指共享边e_j的两个三角面的带噪法向量，length(e_j)是指边e_j的长度，|▽(e_j)|是步骤2中得到的第j条边的非局部全变分梯度算子的模长，

误差项具体公式如下：

其中，Nⁱⁿ是三角面的输入的带噪的面法向量，T是整个三角网格中三角面的个数，是第i个三角面的滤波后的面法向量和第i个三角面的带噪法向量差的2-范数，area(τ_i)是第i个三角面的面积；

α＞0；

步骤(5)：更新三角网格中所有的顶点坐标，更新后的坐标即为滤波后的顶点坐标，其中顶点i的更新后的坐标v′_i的计算公式如下：

其中，v_i是顶点i的测量坐标，F_v(i)表示顶点i的1-邻域，即有公共顶点i的三角形集合，|F_v(i)|表示构成1-邻域的三角面的个数，N_k表示第k个三角面滤波后的面法向量，c_k表示三角面k的质心。