图像去噪是数字图像处理的一个重要分支，主要目的是从观测到的带噪声图像中得到真实图像。在进行图像增强、图像分类和图像重构等算法之前，采用适当方法去除噪声是一个重要的处理过程^[1]。目前，如何在有效去除图像噪声的同时保持图像的边缘和细节成为研究的热点。尤其在强噪声环境下，例如计算机断层扫描(CT)，已广泛用于临床诊断和工业无损检测。然而，暴露在高剂量的辐射下可能会导致细胞病变，体检时高剂量的辐射可能会增加患者患癌症的几率，因而需要使用低剂量的辐射进行扫描，这会导致采集的原始数据中存在强噪声^[2]，进而降低CT图像的质量，破坏组织的边界，使得病灶难以辨别，增加了计算机辅助诊断(CAD)等后续任务^[3-5]的难度。在视频监控等应用中，例如无人停车场汽车过闸监控拍摄过程中存在较恶劣的光照环境，使获取的图像难以识别^[6]。

传统的图像去噪方法可分为计算调和分析和偏微分方程(PDE)2种^[7]。在计算调和分析领域，根据图像频谱的分布规律，从频率上将图像中的真实信息与噪声分开，例如小波变换对2维图像进行小波分解，进而去除噪声^[8]。这类方法认为噪声集中在频谱的高频部分，图像的真实信息分布在低频的一个有限区域内。然而，频谱的高频部分也包含了图像的真实信息，例如图像边缘和纹理细节，由于小波变化对2维图像几何结构的非稀疏表示，导致结果图像出现振铃现象^[9]。一些改进的方法采用多尺度几何分析作为计算调和工具，例如Contourlet^[10]、分形小波^[11]和双树复小波^[12]等变换代替小波变换，取得了更好的去噪效果.但这些方法对图像噪声在频谱上分布的基本假设不变，不能有效地将真实信息与噪声分开，图像边缘和纹理部分过度平滑，造成去噪后的图像清晰度降低，不能很好地去除噪声^[13]。因此，许多学者开始考虑空域上基于变分原理的图像去噪方法。PDE图像去噪模型通过最小化能量函数对图像信息进行约束，采用梯度、散度、曲率等参数进行求解，易于数值实现且具有良好的灵活性和普适性，在图像去噪和增强中得到广泛应用。Rudin等人^[14]提出的全变分(TV)模型以及Shen等人^[15]提出的改进算法通过最小化图像全变分得到分段平滑图像，能很好地保护图像的边缘信息。随后，Zhu等人^[16]提出了两个新的贝叶斯假设方法，分别利用均值曲率和高斯曲率进行先验正则化，较好地克服了TV模型产生分片常数效应的问题。为了获取较好的去噪效果同时保持图像边缘，文献[17-18]分别提取使用结构张量和移动帧方法，获取图像的局部细节，进而得到更好的去噪效果，但这些方法都是显式地计算图像曲率，效率较低。为了快速近似求解变分问题，Gong等人^[19]提出了一种离散滤波方法来计算正则化主导变分模型的降能图像。但曲率滤波算法在强噪声环境下，去噪效果较差，因其在计算图像曲率时，仅考虑局部最小三角切平面投影算子，获取的图像仍会存在部分斑点。而能量函数的正则项只考虑曲率正则，在噪声较大的情况下，曲率正则能量很难显著下降，难以得到较好的结果。

本文在曲率滤波算法的基础上，首先对强噪声图像的特征进行分析，假设强噪声图像中的噪声满足加性噪声模型，进而使用曲率滤波算法去噪并保持图像边缘与纹理细节。其次通过实验发现传统的曲率滤波算法在处理强噪声图像时，存在严重的去噪不完全现象，考虑增强去噪能力的投影算子，使用半窗三角切平面组合最小三角切平面投影算子，抑制结果图像的斑点。然后根据增强去噪能力的投影算子，修正能量函数的正则项，添加局部方差能量，进一步约束斑点的产生。最后与其他算法相比较，本文提出的强噪声图像去噪算法在峰值信噪比(PSNR)和结构相似性(SSIM)等客观度量上明显优于现有算法，具有更好的去噪效果和更强的边缘细节保持能力，特别适用于强噪声图像。

在每个变分模型中，都需要找到一个最小化的能量函数，其能量公式为

$ E\left( \mathit{\boldsymbol{U}} \right) = {E_{{\mathit{\Phi }_{\rm{d}}}}}\left( {\mathit{\boldsymbol{U}},\mathit{\boldsymbol{I}}} \right) + \lambda {E_{{\mathit{\Phi }_{\rm{r}}}}}\left( \mathit{\boldsymbol{U}} \right) $

(1)

式中，$E$表示总能量，包含数据拟合项$ {E_{{\mathit{\Phi }_{\rm{d}}}}}$，用于度量返回图像$\mathit{\boldsymbol{U}} $与输入图像$\mathit{\boldsymbol{I}}$的匹配程度，以及正则能量项${E_{{\mathit{\Phi }_{\rm{r}}}}} $，实现形式化地使用关于图像$\mathit{\boldsymbol{U}}$的先验知识，$λ$为正则能量的权重。

高斯曲率滤波假设原始图像构成的曲面是分块可展的，即高斯曲率处处为零，因而得到边缘保护作用的图像平滑效果。为了避免显式计算高斯曲率的复杂计算量，通过直接调整各点像素值使其位于邻域像素的切平面来满足分段可展的假设，进而隐式使用图像曲率信息，可以较好地保护图像的细节和边缘。采用最小距离调整原则，在所有邻域像素构成的切平面中，寻找与当前像素距离最近的面对像素值进行调整。对于($i, j$)处的像素，计算它与邻域内多个切平面的距离，得到它与多个切平面的最小距离，用于更新($i, j$)处的像素值。

令$\mathit{\boldsymbol{x}} = \left( {i, j} \right) \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $，表示图像空间坐标和2D图像域。用$v(i, j)$表示($i, j$)处的图像灰度值，$ \mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right)$表示当前重构图像，由式(1)可得$ {\mathit{\boldsymbol{U}}_0}\left( \mathit{\boldsymbol{x}} \right) = \mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{x}} \right)$。通过以上定义，图像$\mathit{\boldsymbol{I}}$的高斯曲率为

$ K\left( {\mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{x}} \right)} \right) = \frac{{{\mathit{\boldsymbol{I}}_{xx}}{\mathit{\boldsymbol{I}}_{yy}} - \mathit{\boldsymbol{I}}_{xy}^2}}{{{{\left( {1 + \mathit{\boldsymbol{I}}_x^2 + \mathit{\boldsymbol{I}}_y^2} \right)}^2}}} $

(2)

式中，下标表示不同意义的偏导数。在$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$上对$K$进行积分运算，可得$\mathit{\boldsymbol{I}}$的总高斯曲率。通过Gauss-Bonnet定理可知，任何表面的总高斯曲率与表面拓扑有关。因此，若使总的高斯曲率是拓扑不变量，只能使总绝对高斯曲率最小化，进而将式(1)更新为

$ \begin{array}{*{20}{c}} {E\left( \mathit{\boldsymbol{U}} \right) = {E_{{\mathit{\Phi }_{\rm{d}}}}}\left( {\mathit{\boldsymbol{U}},\mathit{\boldsymbol{I}}} \right) + \lambda E_{{\mathit{\Phi }_{\rm{r}}}}^{{\rm{GC}}}\left( \mathit{\boldsymbol{U}} \right) = }\\ {\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{x}} \right)} \right|{\rm{d}}\mathit{\boldsymbol{x}}} + \lambda \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {K\left( \mathit{\boldsymbol{U}} \right)} \right|{\rm{d}}\mathit{\boldsymbol{x}}} } \end{array} $

(3)

式中，上标GC表示使用高斯曲率先验作为能量函数的正则约束，每个局部的最小化计算$ {{\mathit{\boldsymbol{\hat S}}}_i}\left( \mathit{\boldsymbol{x}} \right), i = 1, \cdots , {N_s}$，通过

$ {\mathit{\boldsymbol{S}}_{\min }}\left( \mathit{\boldsymbol{x}} \right) = \arg \mathop {\min }\limits_{{{\mathit{\boldsymbol{\hat S}}}_i}\left( \mathit{\boldsymbol{x}} \right)} \left| {{{\mathit{\boldsymbol{\hat S}}}_i}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right)} \right| $

(4)

得到，即在有限集合中选择变化最小的元素，更新$ \mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right) = {\mathit{\boldsymbol{S}}_{\min }}\left( \mathit{\boldsymbol{x}} \right)$。

由于图像空间是离散的，式(3)可使用求和代替积分操作，进而得到

$ E\left( \mathit{\boldsymbol{U}} \right) = \sum\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} {\left| {\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{I}}\left( \mathit{\boldsymbol{x}} \right)} \right|} + \lambda \sum\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} {\left| {\mathit{\boldsymbol{K}}\left( {\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right)} \right)} \right|} $

(5)

根据噪声对图像信号的干扰模式，噪声分为加性噪声和乘性噪声两大类。假设$\mathit{\boldsymbol{I}} $表示观察到的图像，$\mathit{\boldsymbol{U}}$表示原始的清晰图像，$n_{1}，n_{2}$分别表示乘性噪声和加性噪声，则噪声模型可表示为

$ \mathit{\boldsymbol{I}} = {n_1} \times \mathit{\boldsymbol{U}} + {n_2} $

(6)

当$n_{1}=0, n_{2}≠0$时表示加性噪声，当$n_{1}≠0$时表示乘性噪声。从式(6)可以看出，加性噪声与图像原始信息无关，而是通过叠加的方式对图像信号产生破坏，常见的椒盐噪声、高斯白噪声和泊松噪声等均属于加性噪声^[20]。乘性噪声与信号强度相关, 在声纳、超声波和激光成像中较为常见。

本文考虑的强噪声图像主要有以下特点：数据中噪声强度大，分布难以预计，且可能存在复合噪声的现象。这些噪声与图像原始信息无关，属于加性噪声。例如来自计算机断层扫描的原始数据、电子元器件内部电子随机热运动产生的电子噪声、高频器件产生的电磁干扰噪声，以及在低剂量辐射强度下产生的较大噪声。由于噪声的随机性呈类泊松分布，但噪声并不是精确的泊松分布，因此很难假设图像噪声的原始分布。加之它的均值和方差相互关联，但并不严格相等，都给图像去噪带来了很大挑战。此外，在视频监控应用中，无人停车场汽车过闸的监控拍摄过程存在较弱的自然光源或其他车辆的强光意外照射等较恶劣光照环境的影响，给图像带来严重的噪声污染。再者雨天也可能给图像带来高斯噪声和椒盐噪声。因而需要一种能有效去除噪声并保留图像细节的方法。

高斯曲率滤波在处理弱噪声图像时可以取得较好的边缘保持效果，但在去除强噪声图像时优势并不明显。原因在于处理强噪声图像时，仅考虑了8邻域内12种(8个)最小三角切平面得到的高斯曲率，且正则项仅考虑了曲率能量的迭代终止条件，存在严重的去噪不完全现象。如图 1所示，高斯曲率滤波较好地保持了图像(例如帽子部分)的边缘信息，但是在纹理、帽羽、头发等细节区域存在去噪能力不足的现象。为了得到保留良好边缘和纹理细节的强噪声图像去噪方法，本文提出通过扩展的增强去噪算子和松弛迭代终止条件以达到对强噪声图像去噪目的。

高斯曲率具有良好的几何特性，用于许多常见的滤波器。其处处为零的曲面可以在形变和扭曲的前提下等距映射到一个平面上，也就是所谓的分段可展，最小化高斯曲率就是使图像尽可能地分段可展。但是使用高斯曲率作为正则能量项依然存在2个问题: 1)基于扩散流的最小化高斯曲率算法收敛缓慢；2)如果显式地迭代计算图像的高斯曲率，需要图像满足2阶可微才能进行计算操作。Gong等人^[19]提出的松弛方法放宽了上述限制，通过隐式计算减小高斯曲率能量时不需要图像的2阶可微性假设，在计算速度上明显优于传统的偏分方程求解方式，展示了在图像去噪处理中的高效性。

为了弥补高斯曲率滤波去噪能力的不足，本文通过扩展的增强去噪投影算子，即半窗三角切平面和最小三角切平面的组合，以及局部方差修正能量函数的正则项，增强了曲率滤波的去噪能力。

定义:若任意可展表面$\mathit{\boldsymbol{S}}$可以用局部切平面$\mathit{\boldsymbol{TS}}$近似表示，则

$ \begin{array}{*{20}{c}} {\forall \mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{S}},\forall \varepsilon > 0,\exists {\mathit{\boldsymbol{x}}_0} \in \mathit{\boldsymbol{S}}}\\ {{\rm{s}}.{\rm{t}}.\;\;\;0 < \left| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_0}} \right| < \varepsilon \;且\;\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{TS}}\left( {{\mathit{\boldsymbol{x}}_0}} \right)} \end{array} $

(7)

即在表面$\mathit{\boldsymbol{S}}$上，$\mathit{\boldsymbol{x}}_{0}$附近任意一点$\mathit{\boldsymbol{x}}$的值可由$\mathit{\boldsymbol{x}}_{0}$切平面上的点近似得到。

根据微分几何的定义，2D表面的高斯曲率为2个主曲率的乘积，记为$k_{1}k_{2}$。式(5)要求总的绝对高斯曲率最小化，即高斯曲率为零，得到尽可能分段可展的表面。则可以转化为求解$k_{1}k_{2}=0$，进一步推导出使得$|k_{1}|、|k_{2}|$中最小值为零。因此通过隐式计算最小化高斯曲率能量，可以较好地实现高斯曲率约束。

对输入图像进行邻域分解时，因为需要考虑邻域构造的曲面，可以用8邻域作为最小切平面的构造区域。为了弥补曲率滤波去噪能力不足的缺陷，本文使用半窗切平面组合和最小三角切平面作为基础投影算子，将领域分成4种半窗三角切平面和4种最小三角切平面，如图 2所示。这样分解有以下优点：1)使用半窗切平面增强了邻域像素的联系，抑制了结果图像中的斑点，弥补了传统高斯曲率滤波算法去噪能力不足的缺陷。2)考虑最小三角切平面，避免了结果图像过度平滑、丢失图像的边缘和纹理细节。

通过枚举中心像素到所有切平面的距离，找到中心点到切平面的最小距离$d_{\min} $。本文考虑局部邻域中8个不同的距离$d_{i}$，其中4种表示半窗三角切平面，如图 2(a)所示；另4种表示最小三角切平面，如图 2(b)所示。为了简化计算，使用中心像素的灰度值与切平面元素集合灰度值均值的差值作为中心像素到切平面的距离，最小距离$d_{\min} $表示$\mathit{\boldsymbol{x}}$从当前图像$\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $到目标图像$\mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right) $的距离，$ \mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right)$=$\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $+$d_{\min} $表示对$\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $的更新，从而减小了局部的高斯曲率能量，进而减小了$ \min \{ |{k_1}|, |{k_2}|\} $。

在迭代图像空间中每个邻域时，通过滑动窗口的方式计算数据拟合项和正则项能量，实现更新每个邻域的运算。通过多次实验，发现原始的高斯曲率滤波算法的能量函数正则项仅考虑了高斯曲率能量的约束。然而对于强噪声图像，数据的原始信息被严重污染，较大程度上破环了图像的几何结构。若使用增强去噪能力的投影算子在原始的能量函数进行迭代，能量函数在第5次迭代内即开始增加，如图 3(a)所示，影响模型的去噪能力，使得算法无法收敛到较好的结果。

通过对传统曲率滤波算法结果图像的特征分析，发现在结果图像中存在较多斑点，且大多出现在图像的平坦区域。受双边滤波器^[21]的启发，在考虑局部几何空间的同时，考虑邻域内相邻像素的灰度值差距，增强算法对噪声的抑制能力，但曲率滤波是通过空间距离关系和像素值差距的非线性组合，过于复杂和耗时，所以本文考虑通过高斯曲率能量和局部像素值差距的线性组合来适应强噪声图像的特征，提出一个新的正则能量约束如下

$ E_{{\mathit{\Phi }_{\rm{r}}}}^{{\rm{GC}}}\left( \mathit{\boldsymbol{U}} \right) = \sum\limits_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left( {\left| {K\left( {\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right)} \right)} \right| + \left| {\sigma \left( \mathit{\boldsymbol{x}} \right)} \right|} \right)} $

(8)

式中，$σ(\mathit{\boldsymbol{x}})$表示在$\mathit{\boldsymbol{x}}$的8邻域像素几何灰度值的方差。因为斑点噪声表现为与周围邻域像素的灰度值差异很大，增大了局部方差，所以通过考虑局部方差可以较好地约束斑点的产生。而且同时考虑曲率能量和局部方差能量，可以较好地保护图像边缘等细节信息。从图 3(b)可以看出，通过使用增加局部方差能量的正则能量项，能量函数在多次迭代后可以收敛至较稳定的结果，证明了本文算法的有效性。

本文算法的具体步骤如下：

1) 对输入的噪声图像$\mathit{\boldsymbol{U}}$，构造组合半窗三角切平面和最小三角切平面的增强去噪投影算子，根据投影算子计算中心像素$\mathit{\boldsymbol{x}}$到8邻域内所有切平面的距离，再从中选出最小距离$d_{\min} $，最后使用得到的最小距离更新$\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $，通过滑动窗口实现噪声图像$\mathit{\boldsymbol{U}}$的更新($ update$($\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $))，具体如下：

输入：$\mathit{\boldsymbol{U}}(i, j)$

输出：$ \mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right)$=$\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right) $+$d_{\min} $

$ \begin{array}{l} {d_1} = \left( {\mathit{\boldsymbol{U}}\left( {i - 1,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,} \right.} \right.\\ \left. {\left. {j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j - 1} \right)} \right)/5 - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_2} = \left( {\mathit{\boldsymbol{U}}\left( {i - 1,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,} \right.} \right.\\ \left. {\left. {j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j + 1} \right)} \right)/5 - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_3} = \left( {\mathit{\boldsymbol{U}}\left( {i + 1,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,} \right.} \right.\\ \left. {\left. {j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j + 1} \right)} \right)/5 - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_4} = \left( {\mathit{\boldsymbol{U}}\left( {i + 1,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,} \right.} \right.\\ \left. {\left. {j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j - 1} \right)} \right)/5 - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_5} = \left( {\mathit{\boldsymbol{U}}\left( {i - 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j} \right)} \right)/3\\ - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_6} = \left( {\mathit{\boldsymbol{U}}\left( {i,j - 1} \right) + \mathit{\boldsymbol{U}}\left( {i + 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i,j + 1} \right)} \right)/3\\ - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_7} = \left( {\mathit{\boldsymbol{U}}\left( {i + 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j} \right)} \right)/3\\ - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ \begin{array}{l} {d_8} = \left( {\mathit{\boldsymbol{U}}\left( {i,j + 1} \right) + \mathit{\boldsymbol{U}}\left( {i - 1,j} \right) + \mathit{\boldsymbol{U}}\left( {i,j - 1} \right)} \right)/3\\ - \mathit{\boldsymbol{U}}\left( {i,j} \right) \end{array} $

$ {d_{\min }} = \min \left\{ {\left| {{d_1}} \right|,\left| {{d_2}} \right|, \cdots ,\left| {{d_8}} \right|} \right\} $

2) 获取更新图像后，使用式(8)计算原始图像和更新图像的能量。

3) 在迭代滤波过程中，比较去噪前后的能量变化，若去噪后图像能量大于去噪之前，则停止迭代。强噪声高斯曲率滤波流程如下：

输入：$\mathit{\boldsymbol{U}}(i, j)$

输出：$ \mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right)$

do{

$E\left( \mathit{\boldsymbol{U}} \right) = {E_{{\mathit{\Phi }_{\rm{d}}}}}\left( {\mathit{\boldsymbol{U}}, \mathit{\boldsymbol{I}}} \right) + \lambda {E_{{\mathit{\Phi }_{\rm{r}}}}}\left( \mathit{\boldsymbol{U}} \right) $

$\mathit{\boldsymbol{\hat U}}\left( \mathit{\boldsymbol{x}} \right) = update\left( {\mathit{\boldsymbol{U}}\left( \mathit{\boldsymbol{x}} \right)} \right) $

$tempE = {E_{{\mathit{\Phi }_{\rm{d}}}}}\left( {\mathit{\boldsymbol{\hat U}}, \mathit{\boldsymbol{I}}} \right) + \lambda {E_{{\mathit{\Phi }_{\rm{r}}}}}(\mathit{\boldsymbol{\hat U}}) $

$\} {\rm{while}}\left( {tempE - E\left( \mathit{\boldsymbol{U}} \right) < 0} \right) $

试验结果表明，本文算法在执行10次内即可得到较好的去噪效果。

为了验证本文算法对强噪声图像的去噪效果，将本文算法与传统基于偏微分方程的ROF模型和高斯曲率滤波方法进行对比研究。实验选用经典的256灰度等级和512×512像素的Lena和airplane图像、包含纯色的圆和三角形以及渐变过渡的矩形基础图形图片、真实的CT肺部切片图像^[22]和视频监控中汽车正面带车牌的图像。

各算法的参数设置如下：ROF模型设置最大迭代次数为1 000，梯度下降步长为0.25，正则项权重为1。高斯曲率滤波算法和本文强噪声曲率滤波算法均设置最大迭代次数为10，正则项权重为1。根据1.2节描述的强噪声图像特点，许多真实场景下的强噪声可以近似为泊松分布和混合高斯分布，本文通过在图像中人为叠加泊松噪声、高斯白噪声和椒盐噪声的多种混合噪声来研究和比较图像去噪效果，具体参数设置在3.2节进行描述。

图 4与图 5是肺气肿患者和肺纤维化患者的肺部CT切片图像的去噪结果对比。图 4是密度为0.1的椒盐噪声叠加方差为50的高斯噪声的混合噪声，图 5是密度为0.1的椒盐噪声和混合方差为50的高斯噪声叠加泊松噪声的混合噪声。

由实验结果可见，ROF模型和高斯曲率滤波去噪后，图像在白色区域还存在一些噪声，而且肺组织部分有一些过平滑，很大程度上干扰了病灶的识别和提取。从残差图像可以看出，本文方法在肺组织和空白区域存在一些过平滑现象，但是在肺组织内部和白色区域很好地去掉了噪声，且没有过度平滑现象，可以很好地帮助特征的发现与提取，如图 4左肺下部与图 5左肺上部的病变区域。相比而言，本文算法取得了最优的去噪结果，光滑区域过渡自然，细节信息保持较完整，如图 4(e)和图 5(e)所示。

图 6是3种算法对5幅噪声图像去噪结果的主观视觉对比，其噪声为混合密度为0.1的椒盐噪声和方差为50的高斯噪声叠加泊松分布的噪声。可以看出，采用本文提出的去噪算法得到的图像，不仅能够有效地去除噪声，例如airplane图像中平坦区域和CT-lung图像中肺部区域的噪声去除，而且原图像中的许多边缘和纹理细节特征也得以很好地保留，例如Lena图像中帽羽、面部五官轮廓信息的保持和汽车图像的车牌区域内车牌号信息的保留都具有较高的视觉质量。与原始图像比较可以看出，利用传统的去噪算法得到的图像的噪声依旧很多，例如ROF对car图像去噪后，车牌部分仍存在较多噪声很难识别车牌信息，而且存在破坏图像边缘和纹理细节的现象。而采用本文去噪算法得到的图像在较好地去除噪声的同时，保留了图像的边缘和纹理细节。

为了定量衡量去噪效果，通过定义峰值信噪比(PSNR)和结构相似性(SSIM) 2项指标对去噪效果进行统计分析。对于一幅$m×n$像素的图像，定义去噪后图像的PSNR为

$ P = 10\lg {\left[ {\frac{{\max \mathit{\boldsymbol{O}}\left( {i,j} \right)}}{R}} \right]^2} $

(9)

式中，$\max \mathit{\boldsymbol{O}}\left( {i, j} \right), i \in \left\{ {1, 2, \cdots , n} \right\}, j \in \left\{ {1, 2, \cdots , m} \right\} $表示原始图像中最大的灰度值，$R$为均方误差百分比(RMSE)，计算公式为

$ R = {\rm{sqrt}}\left( {\sqrt {\frac{1}{{n \times m}}} \sum\limits_{i = 0}^{n - 1} {\sum\limits_{j = 0}^{m - 1} {{{\left[ {\mathit{\boldsymbol{\hat U}}\left( {i,j} \right) - \mathit{\boldsymbol{O}}\left( {i,j} \right)} \right]}^2}} } } \right) $

(10)

式中，$\mathit{\boldsymbol{\hat U}}\left( {i, j} \right), \mathit{\boldsymbol{O}}\left( {i, j} \right) $分别是去噪后和原始图像在位置($i, j$)上的像素点的灰度值。

SSIM作为一种图像质量的评判标准，可以着重反映图像边缘结构的去噪效果。SSIM指标越接近1，说明去噪图像与原图像结构越接近，去噪效果越好。计算公式为

$ S\left( {\mathit{\boldsymbol{\hat U}},\mathit{\boldsymbol{O}}} \right) = \frac{{\left( {2{\mu _{\hat U}}{\mu _O} + {c_1}} \right)\left( {2{\sigma _{\hat UO}} + {c_2}} \right)}}{{\left( {\mu _{\hat U}^2 + \mu _O^2 + {c_1}} \right)\left( {\sigma _{\hat U}^2 + \sigma _O^2 + {c_2}} \right)}} $

(11)

式中，$ {\mu _{\hat U}}, {\mu _O}$分别是去噪后图像$ {\mathit{\boldsymbol{\hat U}}}$和原始图像$\mathit{\boldsymbol{O}}$的灰度平均值；$ {\sigma _{\hat U}}, {\sigma _O}$分别是去噪后图像$ {\mathit{\boldsymbol{\hat U}}}$和原始图像$\mathit{\boldsymbol{O}}$的灰度值方差；$ {\sigma _{\hat UO}}$是去噪后图像$ {\mathit{\boldsymbol{\hat U}}}$和原始图像$\mathit{\boldsymbol{O}}$的灰度值协方差；$ {c_1} = {({k_1}L)^2}, {c_2} = {({k_2}L)^2}$是维持稳定的常数，$L $是像素值的动态范围，$k_{1}=0.01，k_{2}=0.03$。

表 1列出了多幅图像在不同噪声组合下使用ROF模型、高斯曲率滤波(GCF)模型和本文去噪算法处理一幅噪声图像的平均用时、PSNR和SSIM的比较结果。其中S & P、GN和PN分别表示椒盐噪声、高斯白噪声和泊松噪声。一般方差大于30的高斯白噪声属于强噪声^[23]，本文主要考虑以下3种混合的强噪声：1)混合密度为0.02的椒盐噪声和方差为30的高斯白噪声；2)混合密度为0.1的椒盐噪声和方差为50的高斯白噪声；3)混合密度为0.1的椒盐噪声和方差为50的高斯噪声叠加泊松分布的混合噪声。

从表 1可以看出，本文算法的去噪效果明显优于其他算法，例如在椒盐噪声密度为0.02、高斯噪声方差为30时，传统的ROF模型虽然能够得到较好的去噪效果，但需要26 s的运行时间，无法实现实时的去噪目的。同时曲率滤波虽然能够在较短的时间内得到去噪图像，但去噪能力不足，无法很好地去除噪声。而本文算法可以在较短时间内很好地去除图像的噪声并得到很好的去噪结果，尤其在强噪声环境下，例如在椒盐噪声密度为0.1、高斯噪声方差为50，并叠加泊松噪声中，传统的ROF模型和曲率滤波都无法工作的情况下，本文强噪声图像去噪算法，仍可以得到较好的去噪效果并保持图像结构信息。

总的来说，无论是从主观视觉效果还是从PSNR和SSIM的客观评价方面来看，本文算法与传统的去噪算法相比，去噪效果有明显改进，既能够很好地消除噪声，又能够较好地保持图像边缘和纹理细节。

本文提出了一种基于曲率滤波的强噪声图像去噪算法。首先通过扩展的增强去噪投影算子，即半窗三角切平面和最小三角切平面的组合，计算每个中心像素到达邻域切平面的最小距离，并使用最小距离更新当前中心像素。使用曲率滤波中隐式计算曲率的方法，辅助算法快速收敛。最后修正正则项能量函数，考虑邻域内像素灰度值方差，从而达到增强去噪能力和保护图像边缘与细节的目的。实验结果表明，本文去噪方法明显优于实验结果中列出的现有算法。在去除噪声的同时，能有效保持图像的边缘和纹理特征，很好地保留图像的精细结构，取得了良好的去噪效果。在强噪声图像去噪和重构等领域具有重要的应用价值。未来的工作将致力于深入研究估计图像噪声强度并自适应调整邻域内切平面投影算子操作，并进一步研究一种并行加速的去噪方法，实现实时的去噪算法。