172 2016，52（20） Computer Engineering and Applications 计算机工程与应用 ⦾ 图形图像处理 ⦾ 保持泊松噪声图像细节的快速变分去噪算法 杨 燕，金正猛，蒋晓连，刘 艳，张永燕 YANG Yan, JIN Zhengmeng, JIANG Xiaolian, LIU Yan, ZHANG Yongyan 南京邮电大学 理学院，南京 210046 School of Science, Nanjing University of Posts and Telecommunications, Nanjing 210046, China YANG Yan, JIN Zhengmeng, JIANG Xiaolian, et al. Fast variational algorithm based on detail preserving for Pois- son noise removal. Computer Engineering and Applications, 2016, 52（20）：172-176. Abstract：The removal problem of Poisson noise in the medical, astronomical images has been one of the hot topics until now. In this paper, it firstly analyzes the α -Le model of Poisson noise removal and develops a fast algorithm based on a box constraint to solve numerically the model by incorporating Alternating Direction Multiplier（ADMM） algorithm. Then the convergence of the fast algorithm is proved. Finally, numerical results are reported to show that the proposed algorithm, at the same time of denoising, not only preserves small detail characteristics in images, but also improves greatly the computational efficiency. Key words：image denoising; Poisson noise; Alternating Direction Method of Multipliers（ADMM）algorithm; details 摘 要：去除医学、天文图像中的泊松噪声一直是人们关注的热点问题之一。在充分分析泊松去噪 α -Le 模型的基 础上结合交替方向乘子（ADMM）算法，给出该模型一基于框式约束的快速求解算法，并证明了该算法的收敛性。数 值实验结果表明，该算法在去噪的同时，不仅能很好地保留图像中的边缘及小细节特征，还能大幅提高运算效率。 关键词：图像去噪；泊松噪声；交替方向乘子（ADMM）算法；细节 文献标志码：A 中图分类号：TN911.73 doi：10.3778/j.issn.1002-8331.1412-0330 1 引言 对图像进行处理，首要的任务就是对图像进行增强 信噪比的工作，即去除噪声和干扰，突出感兴趣的区域 或边缘，为后续处理奠定基础。根据噪声对图像信号的 影响，一般情况下，噪声大概可以被分为两大类：加性噪 声与乘性噪声。假设 f > 0 表示观察到的图像，u 表示 原始的清晰图像，n 是噪声，则加性噪声模型可表示为： （1） f = u + n 乘性噪声模型可简单地表示为： f = u × n （2） 加性噪声通过叠加的方式对信号产生干扰，与图像 信号无关。传统的加性噪声一般被认为满足高斯正态 分布。但是乘性噪声与信号关系密切，一般出现在相干 成像系统中。如在合成孔径雷达、声纳、超声波和激光 成像中最为常见。对于图像中泊松噪声的处理，传统的 加性噪声处理方法，如均值滤波方法、中值滤波方法 [1] 和 自 适 应 维 纳 滤 波 方 法 [2]等 处 理 后 效 果 均 不 是 很 好 。 2008 年，文献[3]提出基于最大似然估计的滤波方法去 除图像中的泊松噪声，取得了较好效果。但是该方法计 算太复杂，实用性不强。随后，Luisier 等在文献[4]中提 出基于 Haar 小波变换的快速泊松噪声消除算法，相比 于硬阈值和软阈值算法，该方法去除泊松噪声的效果较 好。其他的泊松噪声处理方法，可见参考文献[5-7]等。 本 文 假 设 图 像 中 的 泊 松 噪 声 满 足 加 性 噪 声 模 型 （式（1）），主要关注的是基于变分偏微分方程并保持泊 松图像细节的去噪方法。经典的变分去噪模型可追溯 基金项目：国家自然科学基金（No.11101218）；江苏省重点 STITP 项目（No.SZDG2014025）。 作者简介：杨燕（1989—），女，硕士研究生，主要研究领域为图像处理建模及快速优化算法等；金正猛（1982—），男，博士，副教授， 主要研究领域为非线性偏微分方程及其在图像处理中应用。 收稿日期：2014-12-24 修回日期：2015-06-19 文章编号：1002-8331（2016）20-0172-05 CNKI 网络优先出版：2015-07-03, http://www.cnki.net/kcms/detail/11.2127.TP.20150703.1558.001.html 

杨 燕，金正猛，蒋晓连，等：保持泊松噪声图像细节的快速变分去噪算法 2016，52（20） 173 u Ω  min |f - u|2dx 到 Rudin，Osher 和 Fatemi 于 1992 年提出的 ROF 模型（也 被称为 TV 模型）[8]： |Ñu|dx + λ 其中，u Î BV (Ω) ， 权重系数。该模型在去噪的同时，能很好地保护图像中 的边缘信息。2007 年，Le 等[9]在 TV 模型的基础上，采用 最大后验估计的方法推导出去除泊松噪声的变分模型 （Le 模型）： Ω |Ñu| 表示函数 u 的全变差，λ > 0 是 Ω （3） )  |Ñu|dx + λ ( Ω Ω u - f ln u dx min u > 0 随后，在文献[10]中作者结合自适应的 TV 正则项[8]， 提出保持小细节特征的泊松图像去噪模型（本文中称为 α -Le 模型）： （4）  Ω min u > 0 α(x)|Ñu|dx + λ Ω (u - f ln u)dx （5） 这里 α(x) 表示不同区域扩散速度控制函数。在灰度值 低的区域 α(x) ® 0 ，在灰度值高的区域 α(x) ® 1 。在文 献[10]中，作者证明了变分问题（5）BV 解的存在唯一性 和解的上下界估计，并采用梯度下降法来求解该模型的 数值解。注意到该方法收敛速度较慢且要求 u > 0 以保 证该模型的凸性。在本文中，利用该模型解的上下界估 计作为解的一框式约束（Box Constraint）限制，已达到 松弛 u > 0 这一硬性约束。在此基础上，结合交替方向乘 子 算 法（Alternating Direction Method of Multipliers， ADMM），本文给出基于框式约束的快速全变差图像泊 松去噪算法（CADMM 算法），并给出该算法的收敛性结 果。最后，通过数值实验验证该快速算法的可行性与有 效性。 2 模型介绍 这里记： E(u): =  在文献[10]中，作者构造如下的扩散速度控制函数： α(x) = æ ç è α(x)|Ñu|dx + λ 1 1 + k|G (u - f ln u)dx 1 + kM 2 *f (x)|2 （6） （7） kM 2 1 - ö ÷ ø Ω Ω σ 其中： M = sup x Î Ω (G σ *f )(x) (x) = (1/4πσ)exp{-|x|2/4σ 2} G σ σ > 0 和 k > 0 均为参数。从式（7）可以看出 α(x) 是单调 增函数，0 < α(x) < 1 ，由于在低灰度值的局域 α(x) 较小 (α(x) ® 0) ，所以在这些区域，一些小细节将尽可能地不 被平滑，从而保持一些细节信息；而在高灰度值的区域 α(x) 较大 (α(x) ® 1) ，噪声被平滑，达到去噪的效果。 在文献[10]中，作者已给出 α -Le 模型解的存在唯一 性定理，如下所示。 定理 1 假设式（2）中 f (x) ÎL¥(Ω) ，且 inf f (x) > 0 ， Ω 则变分问题： min {u Î BV (Ω)u > 0} E(u) 存在唯一解 u(x) 且对任意的 x Î Ω ，满足： Ω f (x) f (x)  u(x)  sup 0 < inf 如果记：a = inf （8） f (x) 。通过简单分 析定理 1 中的结论，发现 α -Le 模型等价于如下带框式 约束的变分模型： f (x) ，b = sup Ω Ω Ω min u Î BV (Ω)  Ω α(x)|Ñu|dx + λ Ω (u - f ln u)dx （9） 满足对任意的 x Î Ω ， a  u(x)  b （10） 由于当 u > 0 时，α -Le 模型是凸的。所以当式（9） 中的泛函满足式（10）时，该等价模型是严格凸的。因 此，可以利用基于凸优化的交替方向乘子（ADMM）算 法，给出模型（9）~（10）的快速数值求解算法，并证明该 算法的收敛性。 3 算法求解与收敛性分析 本文算法是在经典 ADMM 算法 [11]基础上，通过引 入辅助变量，将其转化为几个交替能量泛函极值子问 题，把原本复杂的运算转化为简单形式，提高了该类模 型的计算速度。下面结合 ADMM 算法，给出式（9）和 式（10）的一个快速数值求解方法（称为 CADMM 算法） 的步骤： 首先引入两个辅助变量 d ，z ，式（9）和式（10）可等 价转化为：  α(x)|d|dx + λ Ω (z - f ln z)dx min dzu s.t. d = Ñuz = ua  z  b （11） 将上述前两个等式约束问题转化成无约束问题，其 Ω 增广 Lagrange 函数为： ) = min dzu a  z  b L(dzuλ λ 1 2  Ω α(x)|d|dx + λ Ω (z - f ln z)dx + (d - Ñu) + λT 1 (z - u) + λ 2 γ 2 2 γ 1 2    d - Ñu 2 2 +  z - u 2 2 （12） 2 λ ) ÎRl 是 Lagrange 乘子向量，γ 1 这里的 λ = (λ > 0 1 是处罚参数。下面应用 ADMM[11-12]的迭代方法分别对 三个子问题 d，z，u 进行求解：开始 u = uk，λ 。 1 > 0 ，γ = λk 2 = λk 1 ，λ 2 2 （1）求解关于 d 的子问题： d k + 1 ¬ arg min α(x)|d|dx +  d Ω T(d - Ñuk) + λk 1  γ 1 2 d - Ñuk  2 2 （13） 根据文献[13-14]求解 Shrinkage 算子方法，得到： d k + 1(x) = Ñuk - | ||Ñuk - | λk 1 γ 1 λk 1 γ 1 | || | max | ì ||Ñuk - ï í | ï î | || | λk 1 γ 1 - α(x) γ 1 0 （14） ü ï ý ï þ 

174 2016，52（20） Computer Engineering and Applications 计算机工程与应用 （2）求解关于 z 的子问题 zk + 1 ¬ arg min λ (z - f ln z)dx + z Ω (z - uk) + λk 2  γ 2 2 z - uk  2 2 （15） 且 z 满足：a  z  b 。由于极小值点要在导数为零处取 得，关于 z 求导解得： k + 1 z 2 (x) = sk (x)/2 + (sk (x)/2)2 + λ γ 2 f (x) （16） 其中，sk = uk - k 2 λ + λ γ 2 。进一步， zk + 1(x) = P [ab] k + 1 2) (z 其中，P (z) = max{amin(P [ab] [ab] （3）求解关于 u 的子问题： (zb))} 。 uk + 1 ¬ min u T(d k + 1 - Ñu) + λk 1  γ 1 2  d k + 1 - Ñu 2 2 + (zk + 1 - u) + λk 2  γ 2 2  zk + 1 - u 2 2 F æ çç è div(d k + 1 + 利用快速 Fourier 变换及其性质求出 u ： ö ÷÷ ø ö ÷ ÷ ÷ ÷ ÷÷ ÷ ÷ ø 这里 F ，F -1 分别为 Fourier 变换及其逆变换。 uk + 1 = F -1 F(D) - zk + 1 - æ ç ç ç ç çç ç ç è λk 2 γ 1 λk 1 γ 1 γ 2 γ 1 γ 2 γ 1 ) - （6）k = k + 1 步骤 3 否则，迭代停止，求出 u 。 结合经典的 ADMM 算法的收敛性分析 [11-12，15]，下面 给出本文中 CADMM 算法的收敛性结果。由于其证明 过程与带框式约束的 ADMM 算法[16]的收敛性证明完全 类似，故这里省略其证明过程。 > 0 ，γ > 0 和任意的初始值 ) ，由 CADMM 算法产生的迭代序列 {(d kzkuk ) 的 定理 2 对于给定的 γ 1 λ0 (u0λ0 2 1 λk )} 收 敛 到 其 增 广 Lagrange 函 数 L(dzuλ λ λk 2 1 1 鞍点 (u*λ* λ* 2 1 ) 。 2 2 （17） 4 数值实验与结果分析 为了对所提出的 CADMM 算法进行验证，检验其去 除泊松噪声的效果，本章将对多幅测试图像进行仿真实 验，并同现有算法进行比较。这里，用峰值信噪比（Peak to Noise Ratio，PSNR）和 均 方 根 误 差（Mean Signal Square Error，MSE）评价去噪的效果。其定义形式分别为： PSNR = 10 lg （18） MSE = å Ω æ ç è max(u 0 MSE )2 0 (u - u )2 ö ÷ ø （19） MN 本文 CADMM 算法迭代终止条件为：  uk + 1 - uk uk < tol   算法中，tol 是用来控制迭代终止的量，文中取 tol = 10-3 。 由 式（14），式（16）和 式（18）得 出 ，CADMM 的 算 法 步 骤为： 步骤 1 k = 0 ，初始赋值 d = λ 1 = 0 ，z = λ = 0 。 2 步骤 2 当 > tol 时，重复以下子步骤：   uk + 1 - uk uk  （1）d k + 1(x) = Ñuk - | ||Ñuk - | λk 1 γ 1 λk 1 γ 1 | || | max | ì ||Ñuk - ï í | ï î | || | λk 1 γ 1 - α(x) γ 1 ü 0 ï ý ï þ （2）zk + 1(x) = P[ ] ab æ çç è sk + 1(x)/2 + (sk + 1(x)/2)2 + λ γ 其中，P[ ] ab (t) = max{amin(tb)} ，sk + 1 = uk - （3）uk + 1 = F -1 æ çç è div(d k + 1 + ) - λk 1 γ 1 γ 2 γ 1 zk + 1 - F(Δ) - γ 2 γ 1 （4）λk + 1 （5）λk + 1 1 2 = λk 1 = λk 2 (d k + 1 - Ñuk + 1) (zk + 1 - uk + 1) F æ ç ç ç ç çç ç ç è + γ 1 + γ 2 f (x) ö ÷÷ ø 。 2 λ + λk 2 γ 2 λk 2 γ 1 ö ÷÷ ø ö ÷ ÷ ÷ ÷ ÷÷ ÷ ÷ ø 0 其中，u 是去噪后的图像，u 是原始清晰图像，N 和 M 为图像尺寸大小，且式中 MSE 为归一化的均方根误差。 PSNR 和 MSE 是最普遍的评价去噪效果的客观测量方 法。 PSNR 越大，表示图像失真越少，相反，MSE 越小， 去噪效果越好，表明去噪后的图像与原图像越接近。表 1 列出了实验分别采用梯度下降法，Split Bregman 方法， ADMM 算法和本文的 CADMM 算法对应的 PSNR 值和 MSE 值。 表 1 梯度下降法、Split Bregman 算法、ADMM 算法和 CADMM 算法去噪后图像的 PSNR 和 MSE 值，以及各自的计算速度 图像 算法 去噪后 PSNR/dB 去噪后 的 MSE 梯度下降法 Split Bregman ADMM CADMM 梯度下降法 Split Bregman ADMM CADMM 梯度下降法 Split Bregman ADMM CADMM 18.995 7 19.416 9 19.506 6 19.570 1 18.724 0 19.499 4 19.391 7 21.500 9 17.890 8 18.128 3 18.120 1 19.102 9 0.012 6 0.011 4 0.011 2 0.011 0 0.013 4 0.011 2 0.011 5 0.007 1 0.016 3 0.015 4 0.015 4 0.012 3 1 2 3 迭代 次数 603 19 18 15 106 19 18 15 运行 时间/s 29.75 0.93 1.05 0.88 9.78 1.86 1.95 1.83 166 10.18 29 25 23 1.85 1.89 1.78 实验中采用的三幅图均是 0~255 灰度级的图像，分 别对图像 1、2 和 3 添加方差大小不等的泊松噪声，然后 

杨 燕，金正猛，蒋晓连，等：保持泊松噪声图像细节的快速变分去噪算法 2016，52（20） 175 采用本文所提的方法进行去噪。在本文数值实验部分， 图像的灰度值限制在[0，1]范围内。另外，本文的算法 还同文献[10]中梯度下降方法、现行的 Split Bregman 快 速算法[17]和 ADMM 算法进行比较。图中进行去噪处理 时所选参数及框式约束值如表 2。 表 2 三幅图选取的参数及框式约束值 图像 1 2 3 σ 值 0.1 0.1 0.1 k 值 3 0.002 0.01 λ 值 3 4 20 值 γ 1 值 γ 2 5 5 5 8 8 8 框式约束值 a = 0.070 6; b = 1 a = 0.134 3; b = 1 a = 0; b = 1 从表 1 数据可看出：对于人工生成的泊松图像 1、图 像 2 和图像 3，比较 PSNR 值，本文方法都高于梯度下降 法、SplitBregman 方法和 ADMM 算法；对于 MSE 值，本 文方法都低于梯度下降法、Split Bregman 方法和 ADMM 算法。从图 1、图 2 和图 3 可以看出，本文 CADMM 方法 的去噪保真效果更好，它在去除噪声的同时，很好地保 持了图像的边缘和小细节特征。另外。从表 1 中可以 （a）清晰图像 （b）泊松噪声图像 （c）梯度下降法 （d）Split Bregman （e）ADMM （f）CADMM 方法 方法 方法 图 1 4 种算法去噪效果的比较 （a）清晰图像 （b）泊松噪声图像 （c）梯度下降法 看出本文算法的收敛速度明显优于传统的梯度下降算 法，而且也优于 Split Bregman 算法和 ADMM 算法。这 是由于本文算法添加了框式约束使得在保证图像质量 的同时加快了能量泛函的收敛速度。 （a）清晰图像 （b）泊松噪声图像 （c）梯度下降法 （d）Split Bregman （e）ADMM （f）CADMM 方法 方法 方法 图 3 4 种算法去噪效果的比较 5 结束语 本文针对图像中的泊松噪声，从 α -Le 模型出发，给 出 了 一 种 基 于 框 式 约 束 的 全 变 差 泊 松 噪 声 去 除 算 法 （CADMM 算法）。实验结果表明：针对去除泊松噪声的 α -Le 模型，采用 CADMM 算法求解是可行的，且该方法 在有效去除泊松图像噪声的同时，不仅能很好地保持图 像的边缘及一些小细节特征，还能大幅提高该类模型的 计算速度。 参考文献： [1] Nodes T，Gallagher N.Median filters：some modifications and their properties[J].IEEE Transactions on Acoustics， Speech，and Signal Processing，1982，30（5）：739-746. [2] Zhang Huipin.Spatially adaptive-Wiener filtering for image transform[R].Hous- denoising using undecimated wavelet ton：Rice University，1999. [3] Lefkimmiatis S，Papandreous G，Maragos P.Photon-limited image denoising by inference on multiscale models[J].IEEE Transactions on Image Processing，2008，15：2332-2335. [4] Luisier F，Vonesch C，Blu T.Fast interscale wavelet denoising of Poisson-corrupted images[J].Signal Processing，2010， 90（2）：415-427. [5] Chinna R B，Madhavi L M.A combination of wavelet and image denoising technique[J].International Journal fractal of Electronics Engineering，2010，2（2）：259-264. [6] 孙玉宝，韦志辉，吴敏，等.稀疏性正则化的图像泊松去噪 （d）Split Bregman （e）ADMM （f）CADMM 算法[J].电子学报，2011，39（2）：285-290. 方法 方法 方法 [7] 白键，冯象初.一种基于积分微分方程的泊松噪声去除算 图 2 4 种算法去噪效果的比较 法[J].电子与信息学报，2013，35（2）：451-456. 

176 2016，52（20） Computer Engineering and Applications 计算机工程与应用 [8] Rudin L I，Osher S，Fatemi E.Nonlinear total variation based noise removal algorithms[J].Physica D，1992，60（1/4）： 259-268. [9] Le T，Chartrand R，Asaki T J.A variational approach to reconstructing images corrupted by Poisson noise[J].Journal of Mathematical Imaging and Vision，2007，27（3）：257-263. [10] Dong Gang，Guo Zhichang，Wu Boying.A convex adap- tive total variation model based on the gray level indi- cator for multiplicative noise removal[C]//Abstract and Applied Analysis，2013. [11] Gabay D，Mercier B.A dual algorithm for the solution of nonlinear variational problems via finite-element approxi- mations[J].Computers and Mathematics with Applica- tions，1976，2（1）：17-40. [12] Glowinski R，Marrocco A.Sur lapproximation par elements finis dordre un，et al resolution par penalisation-dualite dune classe de problemes de Dirichlet nonlineaires[J]. Revue Francaise D'Automatique，Informatique，Recherche Operationnelle，1975，9（2）：41-76. [13] Wang Y，Yang J，Yin W，et al.A new alternating mini- mization algorithm for total variation image reconstruc- tion[J].SIAM Journal of Mathematical Imaging and Vision， 2008，1（3）：948-951. [14] Yang J，Yin W，Zhang Y.A fast algorithm for edge- preserving variational multichannel image restoration[J]. SIAM Journal of Mathematical Imaging and Vision，2009， 2（2）：569-592. [15] Boyd S，Parikh N，Chu E.Distributed optimization and learning via the alternating direction method in Machine statistical of multipliers[J].Foundations and Trends Learning，2010，3（1）：1-122. [16] Chan R H，Tao M，Yuan X.Constrained total variation deblurring models and fast algorithms based on alter- nating direction method of multipliers[J].SIAM Journal on Imaging Sciences，2013，6（1）：680-697. [17] Goldstein T，Osher S.The split Bregman method for L1 regularized problems[J].SIIMS Journal on Imaging Sciences，2009，2（2）：323-343. （上接 102 页） [6] 余伟峰，钱夕元.基于 KNN 图的两阶段孤立点检测及应用 研究[J].计算机工程与应用，2008，44（2）：186-189. [7] 杨福萍，王洪国，董树霞，等.基于聚类划分的两阶段离群 点检测算法[J].计算法应用研究，2013，30（7）：1942-1945. [8] Ester M，Kriegel H P，Sander J，et al.A density-based algorithm for discovering clusters in large spatial databases with noise[C]//Proc 2nd Int Conf on Knowledge Discovery and Data Mining（KDD-96）.Portland：ACM Press，1996： 226-231. [9] 于亚飞，周爱武.一种改进的 DBSCAN 密度算法[J].计算机 技术与发展，2011，21（2）：30-33. [10] Daszykowski M，Walczak B，Massart D L.Looking for natural patterns in data[J].Chemometrics and Intelligent Laboratory Systems，2001，56（2）：83-92. [11] Hawkins D.Identification of outliers[M].London：Chapman and Hall，1980. [12] Knorr E M，Ng R T，Tucakov V.Distance-based outliers： algorithms and applications[J].VLDB Journal：Very Large Databases，2000：237-253. [13] Ramaswamy S，Rastogi R，Shim K.Efficient algorithms for mining outliers from large data sets[C]//Proceedings of the ACM SIGMOD Conference，2000：437-438. [14] Angiulli F，Pizzuti C.Fast outlier detection in high dimen- sional spaces[C]//Proceedings of the Sixth European Con- ference on the Principle of Data Mining and Knowledge Discovery，2002：15-16. [15] 张悦，刘杰，李航.一种基于概率的孤立点检测方法[J].计 算机工程，2013，39（3）：46-55. （上接 107 页） [15] M andic D P，Wu Z，Huang N E.Empirical mode decomposition-based time-frequency analysis of multivariate signals：the power of adaptive data analysis[J].IEEE Signal Processing Magazine，2013，30（6）：74-86. [16] Mackey G，Sehrish S，Wang J.Improving metadata man- agement for small files in HDFS[C]//Proceedings of IEEE International Conference on Cluster Computing and Work- shops，2009：1-4. [17] Shafer J，Rixner S，Cox A L.The hadoop distributed file- system：balancing portability and performance[C]//Pro- ceedings of 2010 IEEE International Symposium on Per- formance Analysis of Systems & Software（ISPASS）， 2010：122-133. [18] Tian C，Zhou H，He Y，et al.A dynamic mapreduce scheduler for heterogeneous workloads[C]//Proceedings of Eighth International Conference on Grid and Cooperative Com- puting，2009：218-224. [19] O'Neil E J，O'Neil P E，Weikum G.The LRU-K page re- placement algorithm for database disk buffering[J].ACM SIGMOD Record，1993，22：297-306. [20] Xin R S，Gonzalez J E，Franklin M J，et al.Graphx：a resilient distributed graph system on spark[C]//First Inter- national Workshop on Graph Data Management Experi- ences and Systems，2013.