基于法向光滑和曲面拟合的高阶去噪模型求解-凯发娱乐官网

期刊菜单

基于法向光滑和曲面拟合的高阶去噪模型求解
higher-order denoising model solving based on normal smoothing and surface fitting

doi: , , html, ,
作者: 周筱婷, 杨奋林^*：吉首大学数学与统计学院，湖南吉首
关键词: ；；；；；；；

摘要: 高阶变分方程在有效抑制tv模型去噪中的阶梯效应方面表现出色，但其求解过程往往较为复杂。本文针对单位法向量梯度的l₁范数做正则项的高阶变分模型求解问题，提出利用法向光滑和曲面拟合两步将其转化为两个低阶变分模型来求解。对这两个低阶变分模型构造了有效的数值方法，实验结果表明本文提出的方法在峰值信噪比、结构相似性以及计算时间上均优于传统的tv模型。峰值信噪比与结构相似性比tv模型分别高0.3 db、0.02 db，在时间上也表现出了较大的差别。

abstract: higher-order variational equations are excellent in effectively suppressing the step effect in denoising tv models, but their solution process is often complicated. in this paper, for the problem of solving the higher-order variational model with the l1-parameter of the unit normal vector gradient as a regular term, it is proposed to use two steps of normal smoothing and surface fitting to transform it into two lower-order variational models to solve it. an effective numerical method is constructed for these two low-order variational models, and the experimental results show that the method proposed in this paper outperforms the traditional tv model in terms of peak signal-to-noise ratio, structural similarity, and computation time. the peak signal-to-noise ratio and structural similarity are 0.3 db and 0.02 db higher than the tv model, respectively, and also show a large difference in time.

文章引用：周筱婷, 杨奋林. 基于法向光滑和曲面拟合的高阶去噪模型求解[j]. 应用数学进展, 2024, 13(11): 4823-4830.

1. 引言

图像中存在噪声是不可避免的，这可能源于图像形成过程、记录或传输等环节。去噪的目标是从观察到的图像中恢复出原始图像，通常用表达式(1)表示，其中 $z$ 指的就是最终呈现在眼中的图像， $h$ 表示的是真实的图像， $η$ 表示的是噪声。

$z (x, y) = h (x, y) η (x, y) .$ (1)

在过去的二十年中，研究人员开发了多种方法来处理上述逆问题，其中变分法是应用最广泛的技术之一。vogel和oman [1]提出了一个经典的图像增强变分模型，该模型通过总变差(tv)方法在去除图像噪声的同时有效保留了物体的边缘和轮廓，取得了令人满意的效果。

最初的tv范数滤波器[2]是为了获得恢复的图像作为约束优化问题的解，以下是通常的表达形式

$\inf_{h} \int_{ω} | \nabla h | d x s .t . \int_{ω} {| h - z |}^{2} d x = σ^{2} .$ (2)

其中 $σ^{2}$ 表示噪声方差，受到[3]中提出四阶方法的启发，为克服阶梯效应按照[3]中的建议将(7)中的成本函数换为

$\int_{ω} (| h_{x x} | | h_{y y} |) d x 或 \int_{ω} \sqrt{{| h_{x x} |}^{2} {| h_{x y} |}^{2} {| h_{y x} |}^{2} {| h_{y y} |}^{2}} d x .$ (3)

然而，在某些情况下，该方法也表现出一些不利特性，例如图像对比度损失、边角模糊和阶梯效应。为了弥补tv模型的不足，为了克服这些缺点，研究者们尝试了多种改进策略，包括提出了许多其他非线性滤波器和引入其他正则化项，如l₁范数和总变差的加权版本，以增强对比度和细节。近年来，高阶偏微分方程组引起了特别的关注[2]-[4]。本文所使用的方法与高阶偏微分方程组的方法有一定关系，同时受到[5]中处理流场的启发，将求解过程分为两步：首先光滑法向然后曲面拟合得出去噪结果图像[6]。

在本文涉及了函数约束问题其中特别是映射函数，[7]将任意v函数投影到s圆上求解这个问题，这种方法不仅可以用于映射函数的约束问题，还被广泛应用于图像处理、信号恢复等领域。在处理向量值函数时，通过这样的约束投影，可以有效地保持数据的结构特征，同时满足特定的物理或几何约束条件。

$\inf_{h} \int_{ω} | \nabla \frac{\nabla h}{| \nabla h |} | d x, \int_{ω} {| h - z |}^{2} d x = σ^{2} .$ (4)

对于给定的图像 $h$ ， $n = \nabla h / | \nabla h |$ 是水平曲线 $h$ 的单位法向量，如果直接最小化上述泛函的欧拉拉格朗日方程，由于所需的时间步长有限会难以数值求解。在本文中使用这个模型用于平滑法向量，将去噪分为两步求解，第一步法向光滑，第二步需要找到一个合适的曲面来拟合上步所得的法向量。

2. 法向光滑

对于噪声图像z，本文尝试平滑法向量 $n_{0} = \nabla z / | \nabla z |$ 。为了避免在低梯度下的问题分母为零的情况，在梯度幅度计算中添加一个小常数，即 $| \nabla h |$ 将被替代为 $\sqrt{{| \nabla h |}^{2} ε}$ ，所以求解过程将变成一个凸过程[8]。因为法向量可能是不连续的向量函数，所以使用tv范数进行平滑。

与[9]相似在模型中增加一个保真度项来平滑和捕获边缘，具体来讲使用的是双参数方法，一个参数控制流场平滑，另一个参数控制表面拟合，这样与初始状态有较大偏差的时候就会受到惩罚。第一步是求解以下最小化问题得到平滑流场(即：法向量)

$\min_{| n = 1 |} {\int_{ω} | \nabla n | d x \frac{λ}{2} \int_{ω} {| n - n_{0} |}^{2} d x} .$ (5)

其中 $λ > 0$ 是个权重参数， $n_{0}$ 指的是给定图像的法向量， $n$ 是需要求解图像的法向量。 $λ$ 可以通过动态变化选择，但是在本文中将 $λ$ 固定为一个常数，对于它的取值将在今后工作中展开。

$\min_{θ} {\int_{ω} | \nabla θ | d x \frac{λ}{2} {| n - n_{0} |}^{2} d x} .$ (6)

上述问题 $θ$ 的最优性条件为

$\nabla \cdot \frac{\nabla θ}{| \nabla θ |} - λ (n - n_{0}) \cdot \frac{\partial n}{\partial θ} = 0.$ (7)

为求解该问题引入一个人工时间变量 $t$ ，因此这相当于假设 $θ$ 是时间(以及空间)的函数

$θ_{t} = \nabla \cdot \frac{\nabla θ}{| \nabla θ |} - λ (n - n_{0}) \cdot \frac{\partial n}{\partial θ} .$ (8)

直接求解上述方程可能会得到多个解为避免此情况的出现，在此借鉴文献[6]中用极坐标化简的内容将求解(8)转化为

$θ = \arctan \frac{q}{p},$ (9)

$\nabla θ = \frac{p \nabla q - q \nabla p}{p^{2} q^{2}},$ (10)

$\frac{\partial n}{\partial θ} = (- \sin θ, \cos θ) = (- p, q) .$ (11)

有 $n_{0} = (p_{0}, q_{0})$ ，将上述三个式子代入(8)化简可以得到

$\begin{array}{l} \frac{p q_{t} - q p_{t}}{p^{2} q^{2}} = \nabla \cdot (\frac{p \nabla q - q \nabla p}{| p \nabla q - q \nabla p |}) - λ [(p - p_{0}) q - (q - q_{0}) p] \\ = \nabla \cdot (\frac{p \nabla q - q \nabla p}{| p \nabla q - q \nabla p |}) - λ [q p_{0} - p q_{0}] . \end{array}$ (12)

对于坐标的两个分量明显有 $p^{2} q^{2} = 1$ ，对它两边求导得到 $p p_{t} q q_{t} = 0$ ，依据此可以从(12)中将 $p_{t}, q_{t}$ 分离出来，即得到表达式

$p_{t} = - q d i v (\frac{p \nabla q - q \nabla p}{| p \nabla q - q \nabla p |}) λ q (q p_{0} - p q_{0}),$ (13)

$q_{t} = p d i v (\frac{p \nabla q - q \nabla p}{| p \nabla q - q \nabla p |}) - λ p (q p_{0} - p q_{0}) .$ (14)

接下来对(13) (14)离散化求解，为简化标记将 $p^{n}$ 表示 $p_{i, j}^{n}$ ， $q^{n}$ 表示 $q_{i, j}^{n}$ ，同时为方便计算定义 ${[d i v]}^{n} = {[d i v (\frac{p \nabla q - q \nabla p}{| p \nabla q - q \nabla p |})]}^{n}$ 及 ${[p q]}^{n} = {[q p_{0} - p q_{0}]}^{n}$ 。然后分别对(13) (14)有限差分离散化，可以写出半隐式表达为

$\begin{array}{l} p^{n 1} = p^{n} \frac{δ t}{2} \times (q^{n} \frac{δ t}{2} (p^{n 1} p^{n}) ({[d i v]}^{n} - λ {[p q]}^{n}) q^{n}) \times (- {[d i v]}^{n} λ {[p q]}^{n}) \\ = p^{n} q^{n} δ t (- {[d i v]}^{n} λ {[p q]}^{n}) - {(\frac{δ t}{2})}^{2} \times {({[d i v]}^{n} - λ {[p q]}^{n})}^{2} (p^{n 1} p^{n}) . \end{array}$ (15)

$\begin{array}{l} q^{n 1} = q^{n} \frac{δ t}{2} \times (p^{n} \frac{δ t}{2} (q^{n 1} q^{n}) ({[d i v]}^{n} - λ {[p q]}^{n}) p^{n}) \times (- {[d i v]}^{n} λ {[p q]}^{n}) \\ = q^{n} p^{n} δ t (- {[d i v]}^{n} λ {[p q]}^{n}) - {(\frac{δ t}{2})}^{2} \times {({[d i v]}^{n} - λ {[p q]}^{n})}^{2} (q^{n 1} q^{n}) . \end{array}$ (16)

3. 曲面拟合

计算出流场后，需要找到与流场相匹配的图像，需要找到与流场相匹配的图像。为此，本文尝试求解以下最小化问题，以找到一个合适的曲面来拟合法向量[6]。

$\min_{\int_{ω} {| h - z |}^{2} d x = σ^{2}} \int_{ω} (| \nabla h | - \nabla h \cdot n) d x .$ (17)

其中 $n$ 是通过(5)求得的平滑法向量， $σ$ 是l₂范数中的噪声约束，通常使用统计的方法去估计噪声方差估计[10]。在噪声水平未知或不能用统计方法估计的情况下，向最小泛函添加保真度项并丢弃噪声水平约束。文中假设 $σ$ 已知给定为常数，且具有良好的实验结果。

当得到光滑的法向量后，接下来需要寻找与此相匹配的曲面。对于与约束问题(17)直接求解困难且复杂，所以采用拉格朗日乘子处理噪声约束 $\int_{ω} {| h - z |}^{2} d x = σ^{2}$ ，那么(17)的拉格朗日函数定义表示为

$l (h, μ) = \int_{ω} (| \nabla h | - \nabla h \cdot n) d x \frac{μ}{2} (\int_{ω} {| h - z |}^{2} d x - σ^{2}) .$ (18)

因为(17)是非凸的，直接求解非常困难，所以求解上述(17)的最小化则转换为需要找到鞍点求解它，故鞍点在区域内和边界的最优条件为

$\begin{array}{l} \nabla \cdot (\frac{\nabla h}{| \nabla h |} - n) - μ (h - z) = 0 (p, q) \in ω \\ (\frac{\nabla h}{| \nabla h |} - n) \cdot u = 0 (p, q) \in \partial ω \end{array}$ (19)

$\int_{ω} {| h - z |}^{2} d x = σ^{2} .$ (20)

$u$ 是边界 $\partial ω$ 的单位外法向量。直接求解可能会导致数值不稳定，为方便离散化和迭代求解上述方程引入不动点算法[11]，该算法具有便于并行计算能够提高计算效率等优点。因此可以将方程重写为以下形式首先 $μ$ 固定为常数，迭代公式表示如下

$\nabla \cdot (\frac{\nabla h^{n 1}}{| \nabla h^{n} |} - n) - μ (h^{n 1} - z) = 0.$ (21)

用算法的形式表示出求解过程，就可以写成如下形式

算法1 不动点算法求解

步骤1 初始化，给定初始值 $h^{0} = z$ ， $μ$ 为常数， $n = 0$ ， $ε_{n} = 10^{- 4}$ ；

步骤2 计算残差 $r (h) = μ a (\frac{\nabla u}{| \nabla u |}) h - z$ ，其中a为一个线性算子；

步骤3 计算迭代矩阵 $m = μ a d a^{'} i$ ,其中d是对角元为 $\frac{1}{‖ \nabla h ‖}$ 的度角矩阵，i是单位矩阵；

步骤4 计算迭代 $h^{n 1} = h^{n} - {(μ a d a^{'} i)}^{- 1} (μ a (\frac{\nabla h^{n}}{‖ \nabla h^{n} ‖}) h^{n} - z)$ ；

步骤5 $n = n 1$ ，判断 $r (h) < ε_{n}$ 是否成立，转至步骤2。

4. 数值实验

4.1. 评判指标与参数的选取

在本节中首先介绍了改进模型的一些结果，然后进行多组实验，对带有不同信息、不同尺寸图像的去噪结果与tv模型去噪结果作比较。所有的实验均在华硕core(tm) i5-12500h 2.50 ghz，内存16 gb，windows 11系统matlab r2018b中进行。利用指标对图像进行评定，在文中选取了具有代表的psnr (峰值信噪比)和ssim (结构相似性) [12]。

$psnr = 10 \log {\frac{255^{2}}{\frac{1}{m \times n} \sum_{i = 1}^{m} \sum_{j = 1}^{n} {[h (i, j) - z (i, j)]}^{2}}} .$ (22)

$ssim = \frac{(2 μ_{h} μ_{z} c_{1}) (2 σ_{h z} c_{2})}{(μ_{h}^{2} μ_{z}^{2} c_{1}) (σ_{h}^{2} σ_{z}^{2} c_{2})} .$ (23)

设对应于一幅 $m \times n$ 的灰度图像 $h (i, j)$ ，噪声图像 $z (i, j)$ 。 $μ_{h}, μ_{z}$ 分别表示图像 $z, h$ 的平均值， $σ_{z}, σ_{h}$ 分别表示图像 $z, h$ 的方差， $σ_{h z}$ 表示图像 $z, h$ 的像素协方差。 $c_{1} = {(k_{1} l)}^{2}$ ， $c_{2} = {(k_{2} l)}^{2}$ 是两个常数用于避免分母为零的情况，一般默认 $k_{1}, k_{2}$ 取0.01和0.03。其中psnr的单位为db， $ssim \in (0, 1)$ ，ssim值越接近1，两图像的相似程度就越高，表示恢复效果越好。

实验结果主要为二维图像去噪，并对每个示例列出其稳态下的原始噪声图像和干净图像以及参与图像。下面将介绍图像的实验结果，收敛的误差取值 $ε_{n} = 10^{- 4}$ ，为避免低梯度问题取值 $ε = 4$ 。文中所有实验均是在添加了 $σ = 10$ 的随机高斯噪声下进行的仿真实验。为了便于标注，由于该模型是从[6]中改进而来，所以将改进后的模型记为msx模型。

4.2. 实验结果与分析

在本文中测试了许多不同的图像。首先，在图1中使用了一幅“house (128 × 128)”的照片。

truth noisy (28.17 db) tv model (33.55 db) msx model (33.92 db)

figure 1. the graphs of denoising results of different models of “house”

图1. “house”不同模型去噪结果图

从结果可以看出，两种模型都可以有效去除噪声。tv模型在房子的边界处，尤其是烟囱位置，出现了明显的阶梯状失真。改进后的模型在边缘处仍能保持清晰的边界，说明该模型具有很好的边缘保持能力。同时，tv模型的结果显示整个图像偏暗，对比度有所降低，而改进后的模型在色彩保持方面表现较好。

然后在图2中采用的是经典的去噪图像“lena (256 × 256)”。该图像中包含大量的卡通元素(人物的脸部)、帽子和发丝等纹理细节。首先观察人物面部，两个模型都能很好地消除面部噪声。整体来看，改进后的模型使人物脸部更加清晰。在两个结果图中并未显示出明显的差异。为了进一步比较图像中的细节，第二行截取了部分区域进行放大对比。

truth noisy (28.12 db) tv model (32.70 db) msx model (32.82 db)

a region in the original image noisy image in the region the result of the tv model the result of msx model

figure 2. the denoising results of different models of “lena”

图2. “lena”不同模型去噪结果图

紧接着，对比两张经过去噪后放大的脸部图像发现，放大后的tv模型脸部出现了明显的块状现象，而改进后的模型脸部则相对比较平滑。在右侧的头发部分，tv模型的结果中能看出明显的阶梯状结构，右下部分出现了多处白色区域。相较之下，改进模型的去噪结果看起来更加自然还原。对于小尺度的纹理部分，人物的发丝无法与噪声区别开来，这是因为纹理部分具有与噪声接近的信号尺度。改进模型在处理光滑区域方面，相较于tv模型表现得更为适宜。接下来展示的是对非纹理图像“peppers (256 × 256)”的去噪结果，其中多处区域是光滑的。

通过观察图3可以看出改进模型修补的辣椒表面光滑，而tv模型修补的辣椒表面则显得明显块状。整体而言，经过去噪处理后的改进模型辣椒更接近真实辣椒。从局部去噪结果图中可以看出，改进模型的表面依然保持光滑，而tv模型的表面则出现了斑驳的小白块。

truth noisy (28.15 db) tv model (33.54 db) msx model (33.98 db)

a region in the original image noisy image in the region the result of the tv model the result of msx model

figure 3. the denoising results of different models of “peppers”

图3. “peppers”不同模型去噪结果图

为了量化比较结果，表1中列出了两种模型在上述实验中的psnr和ssim。

table 1. table of denoising results of different images under different models

表1. 不同图像在不同模型下的去噪结果表

图像模型	house			lena			peppers
图像模型	psnr	ssim	时间	psnr	ssim	时间	psnr	ssim	时间
tv model	33.55	0.90	1.12	32.70	0.89	1.58	33.54	0.90	1.49
msx model	33.92	0.92	0.87	32.82	0.91	2.01	33.98	0.93	1.79

通过数值表格可以发现，两个模型的去噪结果在psnr和ssim上相近，且改进模型的结果略高于tv模型的去噪结果。在简单的“house”图像中，虽然数值上相差不大，但结果图中可以看到tv模型在边界处出现了阶梯效应。在纹理细节丰富的卡通图像“lena”中，两个模型的psnr和ssim值非常接近，同时结合去噪结果图可以发现，两模型对于纹理细节的处理都不太理想。在光滑区域较多的“peppers”图像中，改进模型不仅在数值上优于tv模型，去噪图像也更接近真实图像，没有出现块状的阶梯效应。整体来看，改进模型的去噪结果优于经典的tv模型。

同时，通过对比时间能够发现，尽管本文求解的是两个约束问题，但在时间方面，不仅收敛速度极快，而且相比经典的tv模型所花费的时间更短，这极大地提高了收敛效率。总体而言，改进后的模型不但在去噪效果上优于tv模型，而且在去噪效率(时间维度)上也高于tv模型。

5. 总结

本文采用不动点迭代算法改进了[6]模型，文中将去噪过程分为两步，首先采用引入人工时间演化方法求解偏微分方程获得光滑法向量，接着使用不动点算法求解第二个偏微分方程寻找与光滑法向量匹配的图像。通过仿真实验说明改进模型不仅取得了良好的去噪结果可以很好地保持边缘边界保持对比度，其中改进模型的psnr值基本上比tv模型高0.3 db，尽管两模型的ssim值很接近但改进模型值仍然比tv模型高0.02，同样在时间上改进模型平均比tv模型快0.3秒。大大地降低了计算时间去噪效率提高，虽然与tv模型相似对处理纹理细节的图像没有很好的效果，但是弥补了tv模型在光滑区域出现的块状现象。接下来考虑如何提高对纹理细节的处理，提高模型对复杂图像的去噪效果。

notes

^*通讯作者。

参考文献

[1]	vogel, c.r. and oman, m.e. (1996) iterative methods for total variation denoising. siam journal on scientific computing, 17, 227-238.
[2]	rudin, l.i., osher, s. and fatemi, e. (1992) nonlinear total variation based noise removal algorithms. physica d: nonlinear phenomena, 60, 259-268.
[3]	lysaker, m., lundervold, a. and tai, x.-c. (2003) noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. ieee transactions on image processing, 12, 1579-1590.
[4]	you, y.-l. and kaveh, m. (2000) fourth-order partial differential equation for noise removal. ieee transactions on image processing, 9, 1723-1730.
[5]	tasdizen, t., whitaker, r., burchard, p. and osher, s. (2003) geometric surface processing via normal maps. acm transactions on graphics, 22, 1012-1033.
[6]	lysaker, m., osher, s. and tai, x.-c. (2004) noise removal using smoothed normals and surface fitting. ieee transactions on image processing, 13, 1345-1357.
[7]	vese, l.a. and osher, s.j. (2002) numerical methods for p-harmonic flows and applications to image processing. siam journal on numerical analysis, 40, 2085-2104.
[8]	zhu, w. and chan, t. (2012) image denoising using mean curvature of image surface. siam journal on imaging sciences, 5, 1-32.
[9]	rudin, l.i., osher, s. and fatemi, e. (1992) nonlinear total variation based noise removal algorithms. physica d: nonlinear phenomena, 60, 259-268.
[10]	meer, p., jolion, j.-m. and rosenfeld, a. (1990) a fast parallel algorithm for blind estimation of noise variance. ieee transactions on pattern analysis and machine intelligence, 12, 216-223.
[11]	acar, r. and vogel, c.r. (1994) analysis of bounded variation penalty methods for ill-posed problems. inverse problems, 10, 1217-1229.
[12]	邱豪骏. 基于变分偏微分方程的图像去噪算法研究[d]: [硕士学位论文]. 深圳: 深圳大学, 2020.

为你推荐

凯发娱乐官网的友情链接