﻿ 基于交替方向乘子法的图像盲复原
 文章快速检索 高级检索
 山东大学学报(工学版)  2017, Vol. 47 Issue (4): 14-18  DOI: 10.6040/j.issn.1672-3961.0.2017.008 0

### 引用本文

LI Zhenwei, CUI Guozhong, GUO Congzhou, YU Changhao. Blind image restoration using alternating direction method of multipliers[J]. Journal of Shandong University (Engineering Science), 2017, 47(4): 14-18. DOI: 10.6040/j.issn.1672-3961.0.2017.008.

### 文章历史

1. 中国人民解放军信息工程大学理学院, 河南 郑州 450001;
2. 中国人民解放军信息工程大学指挥军官基础教育学院, 河南 郑州 450001

Blind image restoration using alternating direction method of multipliers
LI Zhenwei1, CUI Guozhong1, GUO Congzhou1, YU Changhao2
1. School of Science, The PLA Information Engineering University, Zhengzhou 450001, Henan, China;
2. School of Command Officer Basic Education, The PLA Information Engineering University, Zhengzhou 450001, Henan, China
Abstract: In order to overcome the low operating efficiency and poor reconstruction quality in the total variation blind image restoration model of the regularization theory, an iterative algorithm of blind restoration based on alternating direction method of multipliers algorithm was proposed. The restored image and the point spread function were estimated alternatively by alternating iteration to improve the running speed and reconstruction quality through a way without updating the penalty term. The normalization and threshold constraint condition of the point spread function, and the positive definite condition of the image were added while calculating. In the numerical experimentation, the blind restoration of the images with different fuzzy types were carried out, and it was compared with other existing blind image restoration methods.
Key words: blind image restoration    total variation regularization    alternating direction method of multipliers algorithm    point spread function    peak signal to noise ratio    structural similarity index
0 引言

 $Ax+\varepsilon =b,$ (1)

 $Ax=b,$ (2)

 $\text{min}\ f\left( Wx \right)~\quad \text{s}.\text{t}.~\ Ax=b,\text{ }x\in D\left( W \right),$ (3)

 $\underset{x\in D\left( W \right)}{\mathop{\text{min}}}\,\left\{ \frac{1}{2}\|Ax-b{{\|}^{2}}+\alpha f\left( Wx \right)+\beta f\left( WA \right) \right\},$ (4)

1 交替方向乘子法

 $\text{min}~f\left( y \right)\quad ~\text{s}.\text{t}.~\ Ax=b,\text{ }Wx=y,\text{ }x\in D\left( W \right)。$ (5)

 \begin{align} &{{L}_{{{\rho }_{1}},\text{ }{{\rho }_{2}}}}\left( x,\text{ }y,\text{ }\lambda ,\text{ }\mu \right)=f\left( y \right)+\left\langle \lambda ,\text{ }Ax-b \right\rangle + \\ &\left\langle \mu ,\text{ }Wx-y \right\rangle +\frac{{{\rho }_{1}}}{2}\|Ax-b{{\|}^{2}}+\frac{{{\rho }_{2}}}{2}\|Wx-y{{\|}^{2}}, \\ \end{align} (6)

 \left\{ \begin{align} &{{x}_{k+1}}\in \text{arg}\ \text{min}\ {{L}_{{{\rho }_{1}},\text{ }{{\rho }_{2}}}}\left( x,\text{ }{{y}_{k}};{{\lambda }_{k}},{{\mu }_{k}} \right), \\ &{{y}_{k+1}}\in \text{arg}\ \text{min}\ {{L}_{{{\rho }_{1}},\text{ }{{\rho }_{2}}}}\left( {{x}_{k+1}},\text{ }y;{{\mu }_{k}} \right), \\ &{{\lambda }_{k+1}}={{\lambda }_{k}}+{{\rho }_{1}}\left( A{{x}_{k+1}}-b \right), \\ &{{\mu }_{k+1}}={{\mu }_{k}}+{{\rho }_{2}}\left( W{{x}_{k+1}}-{{y}_{k+1}} \right) 。\\ \end{align} \right. (7)

(1) 输入y0, λ0, μ0, 正常数ρ1ρ2, Γ＞1, 噪声等级δ＞0。(δ代表的是测量值与实际数据的偏差, 例如:‖bδ-b‖≤δ)。

(2) 令y0δ=y0, λ0δ=λ0, μ0δ=μ0

(3) 对k=0, 1, 2…, 更新x, y和拉格朗日乘子:

 \begin{align} &\quad \quad \quad {{x}^{k+1}}=\underset{x\in D\left( W \right)}{\mathop{\text{arg}\ \text{min}}}\,\left\{ \left\langle \lambda _{k}^{\delta },\text{ }Ax \right\rangle +\left\langle \mu _{k}^{\delta },\text{ }Wx \right\rangle \right.+ \\ &\left. \quad \quad \quad \frac{{{\rho }_{1}}}{2}\|Ax-{{b}^{\delta }}{{\|}^{2}}+\frac{{{\rho }_{2}}}{2}\|Wx-y_{k}^{\delta }~{{\|}^{2}} \right\}, \\ &y_{k+1}^{\delta }=\underset{y\in Y}{\mathop{\text{arg }\!\!~\!\!\text{ min}}}\,\left\{ f\left( y \right)-\left\langle \mu _{k}^{\delta },y \right\rangle +\frac{{{\rho }_{2}}}{2}\|Ww_{k+1}^{\delta }-y{{\|}^{2}} \right\}, \\ &\quad \quad \quad \quad \lambda _{k+1}^{\delta }=\lambda _{k}^{\delta }+{{\rho }_{1}}\left( Ax_{k+1}^{\delta }-{{b}^{\delta }} \right), \\ &\quad \quad \quad \quad \mu _{k+1}^{\delta }=\mu _{k}^{\delta }+{{\rho }_{2}}\left( Wx_{k+1}^{\delta }-y_{k+1}^{\delta } \right) 。\\ \end{align}

(4) 终止条件:

 \begin{align} &\rho _{1}^{2}~\|Ax_{k}^{\delta }-{{b}^{\delta }}{{\|}^{2}}+\rho _{2}^{2}~\|Wx_{k}^{\delta }-y_{k}^{\delta }~{{\|}^{2}}\le \\ &\quad \quad \quad \quad \text{max}\left( \rho _{1}^{2},\text{ }\rho _{2}^{2} \right){{\tau }^{2}}{{\delta }^{2}} 。\\ \end{align} (8)
2 交替方向乘子盲复原算法

(1) 输入初始y0, λ0, μ0, 常数ρ1, ρ2, ρ3, ρ4均为正, 噪声等级Γ＞1, 噪声等级δ＞0。

(2) 令y0δ=y0, λ0δ=λ0, μ0δ=μ0

(3) 对k=0, 1, 2, …, 更新x, y和拉格朗日乘子:

 \begin{align} &{{A}^{k+1}}=\text{arg }\!\!~\!\!\text{ min}\left\{ \left\langle \lambda _{1k}^{\delta },\text{ }{{A}_{k}}x \right\rangle +\left\langle \mu _{1k}^{\delta },\text{ }W{{A}_{k}} \right\rangle \right.+ \\ &\left. \ \frac{{{\rho }_{1}}}{2}\|{{A}_{k}}{{x}_{k}}-{{b}^{\delta }}{{\|}^{2}}+\frac{{{\rho }_{2}}}{2}\|W{{A}_{k}}-y_{1k}^{\delta }~{{\|}^{2}} \right\}, \\ &{{x}^{k+1}}=\text{arg min}\left\{ \left\langle \lambda _{2k}^{\delta },\text{ }{{A}_{k+1}}x \right\rangle +\left\langle \mu _{2k}^{\delta },\text{ }Wx \right\rangle \right.+ \\ &\left. \ \frac{{{\rho }_{3}}}{2}\|{{A}_{k+1}}{{x}_{k}}-{{b}^{\delta }}{{\|}^{2}}+\frac{{{\rho }_{4}}}{2}\|W{{x}_{k}}-y_{2k}^{\delta }{{\|}^{2}} \right\}, \\ &\quad y_{1k+1}^{\delta }=\text{arg}\,\text{min}\left\{ f\left( {{y}_{1k}} \right)-\left\langle \mu _{1k}^{\delta },\text{ }{{y}_{1k}} \right\rangle \right.+ \\ &\left. \quad \quad \quad \quad \frac{{{\rho }_{2}}}{2}\|Wx_{k+1}^{\delta }-{{y}_{1k}}{{\|}^{2}} \right\}, \\ &\quad y_{2k+1}^{\delta }=\text{arg}\,\text{min}\left\{ f\left( {{y}_{2k}} \right)-\left\langle \mu _{2k}^{\delta },\text{ }{{y}_{2k}} \right\rangle \right.+ \\ &\quad \quad \quad \quad \left. \frac{{{\rho }_{4}}}{2}\|Wx_{k+1}^{\delta }-{{y}_{2k}}{{\|}^{2}} \right\}, \\ &\quad \quad \quad \quad \lambda _{1k+1}^{\delta }=\lambda _{1k}^{\delta }+{{\rho }_{1}}\left( {{A}_{k+1}}x_{k+1}^{\delta }-{{b}^{\delta }} \right), \\ &\quad \quad \quad \quad \lambda _{2k+1}^{\delta }=\lambda _{2k}^{\delta }+{{\rho }_{3}}\left( {{A}_{k+1}}x_{k+1}^{\delta }-{{b}^{\delta }} \right), \\ &\quad \quad \quad \quad \mu _{1k+1}^{\delta }=\mu _{1k}^{\delta }+{{\rho }_{2}}\left( WA_{k+1}^{\delta }-y_{k+1}^{\delta } \right), \\ &\quad \quad \quad \quad \mu _{2k+1}^{\delta }=\mu _{2k}^{\delta }+{{\rho }_{4}}\left( Wx_{k+1}^{\delta }-y_{2k+1}^{\delta } \right) 。\\ \end{align}

(4) 终止条件:

 \begin{align} &\rho _{1}^{2}\|{{A}_{k}}x_{k}^{\delta }-{{b}^{\delta }}{{\|}^{2}}+\rho _{2}^{2}~\|Wx_{k}^{\delta }-y_{1k}^{\delta }~{{\|}^{2}}\le \\ &\quad \quad \quad \quad \text{max}\left( \rho _{1}^{2},\text{ }\rho _{2}^{2} \right){{\tau }^{2}}{{\delta }^{2}}, \\ &\rho _{3}^{2}\|{{A}_{k}}x_{k}^{\delta }-{{b}^{\delta }}{{\|}^{2}}+\rho _{4}^{2}\|WA_{k}^{\delta }-y_{2k}^{\delta }{{\|}^{2}}\le \\ &\quad \quad \quad \quad \text{max}\left( \rho _{3}^{2},\text{ }\rho _{4}^{2} \right){{\tau }^{2}}{{\delta }^{2}} 。\\ \end{align} (9)

 $\begin{array}{l} {\left( {{\nabla _1}x} \right)_{i,{\rm{ }}j}} = \left\{ \begin{array}{l} {x_{i + 1,{\rm{ }}j}} - {x_{i,{\rm{ }}j}},{\rm{ }}i = 1, \cdots M - 1;{\rm{ }}j = 1, \cdots ,{\rm{ }}N,\\ {x_{1,{\rm{ }}j}} - {x_{M,{\rm{ }}j}},{\rm{ }}i = M;{\rm{ }}j = 1, \cdots ,{\rm{ }}N。\end{array} \right.\\ {\left( {{\nabla _2}x} \right)_{i,{\rm{ }}j}} = \left\{ \begin{array}{l} {x_{i,{\rm{ }}j + 1}} - {x_{i,{\rm{ }}j}},{\rm{ }}i = 1, \cdots ,{\rm{ }}M;{\rm{ }}j = 1, \cdots ,{\rm{ }}N - 1,\\ {x_{i,{\rm{ }}1}} - {x_{i,{\rm{ }}N}},{\rm{ }}i = 1, \cdots ,{\rm{ }}M;{\rm{ }}j = N。\end{array} \right. \end{array}$

 \begin{align} &f\left( u,v \right)=\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\sqrt{u_{i,\text{ }j}^{2}+v_{i,\text{ }j}^{2}}}}+ \\ &\quad \quad \frac{v}{2}\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{(u_{i,\text{ }j}^{2}+v_{i,\text{ }j}^{2})}} 。\\ \end{align}

 \begin{align} &A_{k+1}^{\delta }={{\left( {{\rho }_{1}}x_{k}^{*}{{x}_{k}}+{{\rho }_{2}}{{W}^{*}}w \right)}^{-1}}~\left( x_{k}^{*}\left( {{\rho }_{1}}{{b}^{\delta }}-\lambda _{1k}^{\delta } \right) \right.+ \\ &\quad \quad \quad \quad \quad \left. {{W}^{*}}\left( {{\rho }_{2}}y_{1k}^{\delta }-\mu _{1k}^{\delta } \right) \right) 。\\ \end{align} (10)

 $\begin{array}{l} {A^i}\left( {x,{\rm{ }}y} \right) = \left\{ \begin{array}{l} {A^i}\left( {x,{\rm{ }}y} \right),{\rm{ }}{A^i}\left( {x,{\rm{ }}y} \right) \ge \delta {\rm{max}}\left( {{A^i}} \right),\\ 0,{\rm{else,}} \end{array} \right.\\ \quad \quad \quad \quad \quad \int {{A^i}(x,{\rm{ }}y){\rm{d}}x{\rm{d}}y = 1} 。\end{array}$ (11)

 {{x}^{i}}\left( x,\text{ }y \right)=\left\{ \begin{align} &{{x}^{i}}\left( x,\text{ }y \right),\text{ }{{x}^{i}}\left( x,\text{ }y \right)\ge 0, \\ &0,\text{else} 。\\ \end{align} \right. (12)

3 试验结果及评价

 图 1 试验测试图“Lenna” Figure 1 Experimental test chart "Lenna"

 $\text{PSNR}=10\times \text{lo}{{\text{g}}_{10}}\left( \frac{({{2}^{n}}-1)}{\text{MSE}} \right),$ (13)

 $\text{SSIM}\left( x,\text{ }y \right)=\frac{(2{{\mu }_{x}}{{\mu }_{y}}+{{c}_{1}})(2{{\sigma }_{xy}}+{{c}_{2}})}{(\mu _{x}^{2}+\mu _{y}^{2}+{{c}_{1}})(\sigma _{x}^{2}+\sigma _{y}^{2}+{{c}_{2}})},$ (14)

3.1 不同模糊类型的图像复原

 图 2 不同模糊类型盲复原效果图 Figure 2 Effect of different fuzzy types of blind restoration

3.2 与其他图像盲复原算法的比较

 图 3 不同方法盲复原结果 Figure 3 Results of different methods of blind restoration

 图 4 不同算法估计的PSF Figure 4 PSF of different algorithms
3.3 其他图片的复原结果

 图 5 其他图像盲复原结果 Figure 5 Results of blind restoration of other images
4 结论

 [1] WROBEL P, CZYZYCKI M. Direct deconvolution approach for depth profiling ofelement concentrations in multi-layered materials by confocal micro-beam X-ray fluorescence spectrometry[J]. Talanta, 2013, 113: 62-67 DOI:10.1016/j.talanta.2013.03.087 [2] SHIN H C, PRAGER R, GOMERSALL H, et al. Estimation of average speed of sound using deconvolution of medical ultrasound data[J]. Ultrasound in Medicine & Biology, 2010, 36(4): 623-636 [3] BIOUCAS-DIAS J M, PLAZA A, CAMPS-VALLS G, et al. Hyperspectral remote sensing data analysis and future challenges[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(2): 6-36 DOI:10.1109/MGRS.2013.2244672 [4] 贺妍斐. 基于稀疏表示与自适应倒易晶胞的遥感图像复原方法研究[D]. 南京: 南京信息工程大学, 2015H. HE Yanfei. Research on image restoration method based on sparse representation and adaptive reciprocal cell[D].Nanjing: Nanjing University of Information Science and Technology, 2015. http://cdmd.cnki.com.cn/Article/CDMD-10300-1015639685.htm [5] ZHAO X Y, SUN D, TOH K C. A newton-CG augmented lagrangian method for semidefinite programming[J]. SIAM Journal on Optimization, 2010, 20(4): 1737-1765 DOI:10.1137/080718206 [6] 刘国良, 邵云龙. 基于伸缩梯度投影的天文图像复原改进算法[J]. 传感器与微系统, 2016, 35(4): 134-136 LIU Guoliang, SHAO Yunlong. Improved astronomical image restoration algorithm based on scaled gradient projection[J]. Transducer and Microsystem Technologies, 2016, 35(4): 134-136 [7] CHEN Q, MONTESINOS P, SUN Q S, et al. Adaptive total variation denoising based on difference curvature[J]. Image and Vision Computing, 2010, 28(3): 298-306 DOI:10.1016/j.imavis.2009.04.012 [8] 李旭超, 宋博. 原始-对偶模型的牛顿迭代原理与图像恢复[J]. 电子学报, 2015(10): 1984-1993 LI Xuchao, SONG Bo. Newton iterative principle of primal-dual model and image restoration[J]. Acta Electronica Sinica, 2015(10): 1984-1993 DOI:10.3969/j.issn.0372-2112.2015.10.016 [9] WEN Y W, CHAN R H, YIP A M. A primal-dual method for total-variation-based wavelet domain inpainting[J]. IEEE Transactions on Image Processing, 2012, 21(1): 106-114 DOI:10.1109/TIP.2011.2159983 [10] TAO M, YANG J. Alternating direction algorithms for total variation deconvolution in image reconstruction[EB/OL](2009-11-17).[2016-07-15]. http://www.optimization-online.org/DB_HTML/2009/11/2463.html. [11] LI W, LI Q, GONG W, et al. Total variation blind deconvolution employing split bregman iteration[J]. Journal of Visual Communication & Image Representation, 2012, 23(3): 409-417 [12] 李权利. 全变差正则化盲图像复原技术研究[D]. 重庆: 重庆大学, 2012L. LI Quanli.Total variation blind image restoration[D].Chongqing: Chongqing University, 2012. http://d.wanfangdata.com.cn/Thesis/Y2151362 [13] SHI Y, CHANG Q. Efficient algorithm for isotropic and anisotropic total variation deblurring and denoising[J]. Journal of Applied Mathematics, 2013, 2013: 1-14 [14] JIAO Y, JIN Q, LU X, et al. Alternating direction method of multipliers for linear inverse problems[J]. Siam Journal on Numerical Analysis, 2016, 54(4): 2114-2137 DOI:10.1137/15M1029308 [15] XU Y, HUANG T Z, LIU J, et al. Split bregman iteration algorithm for image deblurring using fourth-order total bounded variation regularization model[J]. Journal of Applied Mathematics, 2013, 2013(3): 417-433