CN104361562A - 一种变分pde图像保边保结构复原方法 - Google Patents

Abstract

本发明公开了一种变分PDE图像保边保结构复原方法，在模型中引入结构特征检测函数来惩罚模型中的不同泛函空间，从而使模型具有自适应性，进而得到高质的复原图像。在模型算法实现的过程中，我们首先将问题转化为易求解的多个子问题，并采用周期差分的格式，从而数值模型具有可以利用快速傅里叶变化的特性，进而可以高效稳定的求解。本发明中的软件包侧重于将模型转化为可快速求解的子问题，由于子问题具有可分性，且采用循环差分格式可利用快速傅里叶变换，因此可以高效稳定地求解且具有可并行性。

Description

一种变分PDE图像保边保结构复原方法

技术领域

本发明属于图像复原技术技术领域，涉及一种变分PDE图像保边保结构复原方法。

背景技术

图像在形成、记录、处理和传输过程中，由于成像系统、记录设备和传输介质等不完善，经常在图像中引入噪声、模糊或者部分信息丢失等质量下降问题。例如，大气湍流的扰动效应，环境条件的变化和传感元器件自身的质量原因产生的噪声干扰，被摄物与成像设备间的相对运动造成的运动模糊等。图像复原是利用图像的某些先验知识来重建图像从而改善图像质量的技术。然而,随着对图像质量和图像分辨率要求的提高，图像复原算法的代价和硬件实现的复杂度显著地增加,且实际应用中对图像复原算法通常有实时性要求,这就对现有数字图像复原技术提出了更大的挑战。

传统的图像处理技术主要分为基于概率统计的方法、基于小波理论的方法和基于变分偏微分方程的方法。最近二三十年来,基于变分偏微分方程的图像复原模型由于具有良好的数学理论和物理意义,同时可以开发高效的并行算法和软件实现,从而在图像复原领域内得到了广泛的关注,同时也预示了其广阔的应用前景。一个经典的变分偏微分方程方法是由Rudin、Osher和Fatemi提出的全变分模型(ROF模型):

$\underset{u}{\min }{\∫}_{\mathrm{\Ω}}|\▿u|\mathrm{dx}+\frac{\mathrm{\λ}}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_2^2,$

其中K为模糊算子,f为退化图像,u为待复原图像。该模型由于含有全变分项,可以保持图像的边缘；同时是建立在能量泛函的基础上,尤其是最近几年可以利用分裂算法快速求解；因此在图像复原领域得到了广泛的应用和关注。然而ROF模型经常在图像的光滑渐变区域引入阶梯现象,可以利用高阶模型保持图像光滑渐变区域的特性来缓解阶梯现象:

$\underset{u}{\min }{\∫}_{\mathrm{\Ω}}\left|{\▿}^{2}u\right|\mathrm{dx}+\frac{\mathrm{\λ}}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_2^2,$

其中但是高阶模型在复原图像的过程中经常导致图像边界模糊。为了有效地保持图像的边缘和光滑渐变区域等结构特征,有效的方法是耦合上述两种模型的正则项。然而,现有的耦合模型并没有考虑图像的局部结构特征,并且主要针对图像复原问题的某一类问题建模(如:图像去噪或者图像去模糊等),缺乏统一的框架。

发明内容

本发明的目的在于克服上述技术存在的保持图像棱角和边缘及其光滑渐变区域等结构特征比较差、计算效率比较低和不利于处理大数据的的弱点，提供一种变分PDE图像保边保结构复原方法，该方法结合全变分型泛函空间和自适应的参数选取准则，采用将模型的目标函数转化为可高效求解多个子问题的策略，建立新的复原模型和开发出稳定高效的应用软件来复原高质量的图像，从而为遥感图像处理、医学图像处理、视频目标跟踪以及图像结构特征提取中图像分割、融合、补全、配准和识别等工作提供强有力的理论支持和技术支撑。

其具体技术方案为：

一种变分PDE图像保边保结构复原方法，包括以下步骤：

(1)输入一副大小为N×N的退化图像f(x,y),令f(i,j)表示图像在第i行第j列的像素值,利用Matlab中的diff函数定义函数g(x,y)的梯度算子▽g＝(▽_xg,▽_yg)^T和海森算子

▽²g＝(▽_xxg,▽_xyg；▽_yxg,▽_yyg)^T；

▽_xg＝[diff(g,1,2),g(:,1)-g(:,N)]；

▽_yg＝[diff(g,1,1)；g(1,:)-g(m,:)]；

▽_xxg＝[▽_xg(:,1)-▽_xg(:,N)；diff(▽_xg,1,2)]；

▽_yyg＝[▽_yg(1,:)-▽_yg(N,:)；diff(▽_yg,1,1)]；

▽_xyg＝[diff(▽_xg,1,1)；▽_xg(1,:)-▽_xg(N,:)].

选取初始值u⁰＝f,w₀＝0,v⁰＝0,z⁰＝0,α⁰＝0,β⁰＝0,λ⁰＝0.

(2)求解发明的变分PDE模型模型

$\underset{u\∈C}{\min }\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2+b{|\▿|}_{l^1}+(1-b){\left|{\▿}^{2}u\right|}_{l^1},$

其中

${|\▿|}_{l^1}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|{\▿u}_{i,j}\right|=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{{\left({\▿}_x{u}_{i,j}\right)}^2+{\left({\▿}_y{u}_{i,j}\right)}^2},$

${|{\▿}^2\▿|}_{l^1}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left|{\▿}^2{u}_{i,j}\right|=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{{\left({\▿}_{\mathrm{xx}}{u}_{i,j}\right)}^2+{\left({\▿}_{\mathrm{xy}}{u}_{i,j}\right)}^2+{\left({\▿}_{\mathrm{yx}}{u}_{i,j}\right)}^2+{\left({\▿}_{\mathrm{yy}}{u}_{i,j}\right)}^2.}$

由于上述问题是含有非光滑的l¹范数的l¹-l²优化问题,从而在数值实现中给直接求解带来了巨大的挑战。最近几年,利用算子分裂的方法在求解l¹-l²问题领域得到了广泛的关注。其中一类有效的算子分裂方法是用增广拉格朗日方法,即引入约束条件

w＝u,v＝▽u,z＝▽²u.

将其转化为下面极大极小的鞍点问题

$\underset{u,w,v,z}{\min }\underset{\mathrm{\α},\mathrm{\β},\mathrm{\λ}}{\max }L(u,w,v,z,\mathrm{\α},\mathrm{\β},\mathrm{\λ}),$

其中

$\begin{array}{c}L(u,w,v,z,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2+{<\mathrm{\&alpha;},w-u>}_{l^2}+\frac{c_1}{2}{\left|\right|w-u\left|\right|}_{l^2}^2+{\mathrm{\&delta;}}_C\left(w\right)+{\left|\mathrm{bv}\right|}_{l^1}\\ +{<\mathrm{\&beta;},v-\&dtri;u>}_{l^2}+\frac{c_2}{2}{\left|\right|v-\&dtri;u\left|\right|}_{l^2}^2+{\left|(1-b)z\right|}_{l^1}+{<\mathrm{\&lambda;},z-{\&dtri;}^{2}u>}_{l^2}+\frac{c_3}{2}{\left|\right|z-{\&dtri;}^{2}u\left|\right|}_{l^2}^2\end{array}$

这里α,β,λ分别为对应的拉格朗日乘子，c₁,c₂,c₃分别为对应正的罚参数，δ_C(w)为示性函数,即w∈C，则δ_C(w)＝0；否则，δ_C(w)＝+∞。

(3)极大极小的鞍点问题中耦合了七个变量,有效的方法是求解某一个变量的同时固定其它六个变量,即为去耦合方法,亦称之为交替方向乘子法。交替方法乘子法在求解凸优化问题时具有计算速度比较快和可并行实现的优点,并且在理论上能保持算法的收敛性和稳定性,从而适宜于求解大规模的图像复原问题。考虑交替方法乘子法求解的具体步骤,在算法中应用Gauss-Steidl迭代格式。

●求解子问题u^k+1。即考虑

$\underset{u,w^k,v^k,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u,w^k,v^k,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k)$

由于w^k,v^k,z^k,α^k,β^k,λ^k可看作常变量,因此舍弃与u无关的变量便有

$\begin{array}{c}{u}^{k+1}:=\underset{u}{\arg \min }\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2-{<{\mathrm{\&alpha;}}^k,u>}_{l^2}+\frac{c_1}{2}{\left|\right|w^k-u\left|\right|}_{l^2}^2-{<{\mathrm{\&beta;}}^k,\&dtri;u>}_{l^2}\\ +\frac{c_2}{2}{\left|\right|v^k-\&dtri;u\left|\right|}_{l^2}^2-{<{\mathrm{\&lambda;}}^k,{\&dtri;}^{2}u>}_{l^2}+\frac{c_3}{2}{\left|\right|z^k-{\&dtri;}^{2}u\left|\right|}_{l^2}^2.\end{array}$

上述优化问题为光滑优化问题,因此可求出其对应的显式解。即对应的Euler-Lagrangian方程为

(A_IK^TK-c₁I-c₂Δ+c₃div²▽²)u^k+1＝A_IK^Tf+α^k+c₁w^k-divβ^k

                       -c₂divv^k+div²λ^k-c₃div²z^k.

其中div²为▽²的伴随矩阵。由于上述等式中选取的为循环边界差分格式,同时利用散度定理中的<p,▽q>＝-<divp,q>和<s,▽²t>＝<div²s,t>,因此可以利用快速FFT变化求解u^k+1。即为：

$u=F^{-1}\left(\frac{A_{I}F\left(K^T\right)F\left(f\right)+F\left({\mathrm{\&alpha;}}^k\right)+c_{1}F\left(w^k\right)-F\left(\mathrm{div}\right)F\left({\mathrm{\&beta;}}^k\right)-c_{2}F\left(\mathrm{div}\right)F\left(v^k\right)+F\left({\mathrm{div}}^2\right)F\left({\mathrm{\&lambda;}}^k\right)-c_{3}F\left({\mathrm{div}}^2\right)F\left(z^k\right)}{(A_{I}F\left(K^{T}K\right)-c_{1}F\left(I\right)-c_{2}F\left(\mathrm{\&Delta;}\right)+c_{3}F\left({\mathrm{div}}^2{\&dtri;}^2\right))}\right)$

其中F表示快速傅里叶变化,F^-1表示快速傅里叶变化的逆变换。上述算子具体实现在程序中定义为:

●求解子问题w^k+1。即考虑

$\underset{u^{k+1},w,v^k,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u^{k+1},w,v^k,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k).$

舍去与w无关项,则对应的子问题为

$w^{k+1}=\underset{w}{\arg \min }{<\mathrm{\&alpha;},w>}_{l^2}+\frac{c_1}{2}{\left|\right|w-u^{k+1}\left|\right|}_{l^2}^2+{\mathrm{\&delta;}}_C\left(w\right).$

该问题可以用投影梯度法快速求解,因此有

$w^{k+1}=P_C({\mathrm{\&alpha;}}^k+c_1{u}^{k+1}):=\left\{\begin{array}{cc}{\mathrm{\&alpha;}}^k+c_1{u}^{k+1},& \mathrm{if}({\mathrm{\&alpha;}}^k+c_1{u}^{k+1})\&Element;C\\ 0,& \mathrm{otherwise}.\end{array}\right.$

●求解子问题v^k+1。即考虑

$\underset{u^{k+1},w^{k+1},v,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u^{k+1},w^{k+1},v,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k)$

舍去与v无关项,则对应的子问题为

$v^{k+1}=\underset{v}{\arg \min }\frac{c_2}{2}{\left|\right|v-({\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2})\left|\right|}_{l^2}^2+b{\left|v\right|}_{l^1}.$

该问题对应经典的压缩阈值问题,因此其显式解为

$v^{k+1}=\frac{{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}}{{|{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}|}_{l^1}}\max \{{|{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}|}_{l^1}-\frac{b}{c_2},0\}.$

●求解子问题z^k+1。由于子问题z与v具有类似的形式,因此同样可写成其对应的显式解为

$z^{k+1}=\frac{{\&dtri;}^2{u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}}{{|{\&dtri;}^2{u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}|}_{l^1}}\max \{{|{\&dtri;}^2{u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}|}_{l^1}-\frac{1-b}{c_3},0\}$

由于拉格朗日乘子对应对应的问题是平凡的,因此其更新变量即为

α^k+1＝α^k+c₁(w^k+1-u^k+1)；

β^k+1＝β^k+c₂(v^k+1-▽u^k+1)；

λ^k+1＝λ^k+c₃(z^k+1-▽²u^k+1).

●终止标准。由交替方向乘子法的收敛性定理可知,当 ${\left|\right|w^{k+1}-u^{k+1}\left|\right|}_2^2,{\left|\right|v^{k+1}-{\&dtri;u}^{k+1}\left|\right|}_2^2,{\left|\right|z^{k+1}-{\&dtri;}^2{u}^{k+1}\left|\right|}_2^2$ 收敛到零时,算法的解收敛到初始问题的解。因此在算法中选取终止标准定为当算法终止时,输入复原图像u:＝u^k+1。

与现有技术相比，本发明的有益效果为：

本发明中采有保持图像低阶全变分泛函空间耦合保持图像光滑渐变区域的高阶全变分泛函空间的方法建立新的变分PDE图像复原模型，从而建立的模型从理论上可以保持图像的棱角和光滑区域等结构特征；

本发明通过引入判断图像尺度自适应性的函数惩罚不同泛函空间的在模型中的权重，进而使得建立的模型能有效地刻画图像的局部信息，从而可以复原高质量的图像；

本发明中的软件包侧重于将模型转化为可快速求解的子问题，由于子问题具有可分性，且采用循环差分格式可利用快速傅里叶变换，因此可以高效求解且具有可并行性。

附图说明

图1是本发明基于耦合全变分泛函空间的图像保边保结构复原方法的流程示意图；

图2是图像去噪实验效果图，其中，图2(a)初始图像，图2(b)初始图像第255行像素值，图2(c)初始图像部分放大，图2(a1)噪声图像，图2(b1)噪声图像第255行像素值，图2(c1)噪声图像部分放大，图2(a2)ROF模型复原图像，图2(b2)复原图像第255行像素值，图2(c2)复原图像部分放大，图2(a3)高阶模型复原图像，图2(b3)复原图像第255行像素值，图2(c3)复原图像部分放大，图2(a4)建议模型复原图像，图2(b4)复原图像第255行像素值，图2(c4)复原图像部分放大；

图3是图像去模糊实验效果图，其中，图3(a)待修补图像，图3(b)ROF模型复原图像，图3(c)高阶模型复原图像，图3(d)建议模型复原图像；

图4是图像修补实验效果图，其中，图4(a)待修补图像，图4(b)ROF模型复原图像，图4(c)高阶模型复原图像，图4(d)建议模型复原图像。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合附图和具体实施例，进一步阐述本发明。

数值实验：

下述数值实验是应用Matlab7.6软件进行仿真实验，实验环境为Acer3820GT的个人笔记本电脑，Intel Core i5processing520M,2.40GHz,4GB内存。为了评估复原图像的质量，我们采用信噪比SNR来刻画：

$\mathrm{SNR}=10\log \left(\frac{{\left|\right|f-\stackrel{\&OverBar;}{f}\left|\right|}_{l^2}^2}{{\left|\right|\mathrm{\&eta;}-\stackrel{\&OverBar;}{\mathrm{\&eta;}}\left|\right|}_{l^2}^2}\right),$

其中表示f的均值,表示η的均值,η为复原图像与参考图像的差异。显然,当图像具有较高的SNR时预示着较好的复原图像。

实施例1图像去噪实验：

如图2所示，通过对比不难发现，发明的模型能明显的保持图像的结构特征，如图像的边缘区域和光滑渐变区域。其中ROF模型复原图像的信噪比为17.7677，高阶模型的信噪比为17.9998，发明的模型信噪比为18.3694。

实施例2图像去模糊实验：

如图3所示，发明的模型可以较好地复原退化图像，尤其在保持图像的细节特征和结构特征更为有效。ROF模型复原图像的信噪比为14.3778，高阶模型的信噪比为14.4778，发明的模型信噪比为14.5118。

实施例3图像去模糊实验：

如图4所示，通过局部放大图像可以看出，ROF模型中在辣椒茎部含有明显的字母信息,高阶模型导致复原图像的边界模糊，然而建议的模型能很好的克服这两类模型的缺陷。对应地，ROF模型复原图像的信噪比为15.4055，高阶模型的信噪比为16.4336，发明的模型信噪比为16.6034。

以上所述，仅为本发明最佳实施方式，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。

Claims (1)

1.一种变分PDE图像保边保结构复原方法，其特征在于，包括以下步骤：

(1)输入一副大小为N×N的退化图像f(x,y),令f(i,j)表示图像在第i行第j列的像素值,利用Matlab中的diff函数定义函数g(x,y)的梯度算子和海森算子

${\&dtri;}^{2}g={({\&dtri;}_{\mathrm{xx}}g,{\&dtri;}_{\mathrm{xy}}g;{\&dtri;}_{\mathrm{yx}}g,{\&dtri;}_{\mathrm{yy}}g)}^T:$

${\&dtri;}_{x}g=[\mathrm{diff}(g,\mathrm{1,2}),g(:,1)-g(:,N)];$

${\&dtri;}_{y}g=[\mathrm{diff}(g,1,1),g(1,:)-g(m,:)];$

${\&dtri;}_{\mathrm{xx}}g=[{\&dtri;}_{x}g(:,1)-{\&dtri;}_{x}g(:,N);\mathrm{diff}({\&dtri;}_{x}g,\mathrm{1,2})];$

${\&dtri;}_{\mathrm{yy}}g=[{\&dtri;}_{y}g(1,:)-{\&dtri;}_{y}g(N,:);\mathrm{diff}({\&dtri;}_{y}g,1,1)];$

${\&dtri;}_{\mathrm{xy}}g=[\mathrm{diff}({\&dtri;}_{x}g,1,1,):{\&dtri;}_{x}g(1,:)-{\&dtri;}_{x}g(N,:)].$

选取初始值u⁰＝f,w₀＝0,v⁰＝0,z⁰＝0,α⁰＝0,β⁰＝0,λ⁰＝0.；

(2)利求解离散化模型

$\underset{u\&Element;C}{\min }\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2+b{|\&dtri;|}_{l^1}+(1-b){\left|{\&dtri;}^{2}u\right|}_{l^1},$

其中

${|\&dtri;|}_{l^1}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|{\&dtri;u}_{i,j}\right|=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\sqrt{{\left({\&dtri;}_x{u}_{i,j}\right)}^2+{\left({\&dtri;}_y{u}_{i,j}\right)}^2},$

${|{\&dtri;}^2\&dtri;|}_{l^1}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|{{\&dtri;}^{2}u}_{i,j}\right|=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\sqrt{{\left({\&dtri;}_{\mathrm{xx}}{u}_{i,j}\right)}^2+{\left({\&dtri;}_{\mathrm{xy}}{u}_{i,j}\right)}^2+{\left({\&dtri;}_{\mathrm{yx}}{u}_{i,j}\right)}^2+{\left({\&dtri;}_{\mathrm{yy}}{u}_{i,j}\right)}^2}.$

即引入约束条件

$w=u,v=\&dtri;u,z={\&dtri;}^{2}u.$

将其转化为下面极大极小的鞍点问题

$\underset{u,w,v,z}{\min }\underset{\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;}}{\max }L(u,w,v,z,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;}),$

其中

$\begin{array}{c}L(u,w,v,z,\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&lambda;})=\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2{+\&lang;\mathrm{\&alpha;},w-u\&rang;}_{l^2}+\frac{c_1}{2}{\left|\right|w-u\left|\right|}_{l^2}^2+{\mathrm{\&delta;}}_C\left(w\right)+{\left|\mathrm{bv}\right|}_{l^1}\\ +{\&lang;\mathrm{\&beta;},v-\&dtri;u\&rang;}_{l^2}+\frac{c_2}{2}{\left|\right|v-\&dtri;u\left|\right|}_{l^2}^2+{\left|(1-b)z\right|}_{l^1}+{\&lang;\mathrm{\&lambda;},z-{\&dtri;}^{2}u\&rang;}_{l^2}+\frac{c_3}{2}{\left|\right|z-{\&dtri;}^{2}u\left|\right|}_{l^2}^2\end{array}$

这里α,β,λ分别为对应的拉格朗日乘子，c₁,c₂,c₃分别为对应的正的罚参数，δ_C(w)为示性函数,即w∈C，则δ_C(w)＝0；否则，δ_C(w)＝+∞；

(3)极大极小的鞍点问题中耦合了七个变量,在算法中应用Gauss-Steidl迭代格式，

求解子问题u^k+1.即考虑

$\underset{u,w^k,v^k,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u,w^k,v^k,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k),$

由于w^k,v^k,z^k,α^k,β^k,λ^k可看做常变量,因此舍弃与u无关的变量便有

$\begin{array}{c}{u}^{k+1}:=\underset{u}{\arg \min }\frac{A_I}{2}\left|\right|{\mathrm{Ku}-\left|\right|}_{l^2}^2-{\&lang;{\mathrm{\&alpha;}}^k,u\&rang;}_{l^2}+\frac{c_1}{2}{\left|\right|w^k-u\left|\right|}_{l^2}^2-{\&lang;{\mathrm{\&beta;}}^k,\&dtri;u\&rang;}_{l^2}\\ +\frac{c_2}{2}{\left|\right|v^k-\&dtri;u\left|\right|}_{l^2}^2-{\&lang;{\mathrm{\&lambda;}}^k,{\&dtri;}^{2}u\&rang;}_{l^2}+\frac{c_3}{2}{\left|\right|z^k-{\&dtri;}^{2}u\left|\right|}_{l^2}^2.\end{array}$

求出其对应的显式解，即对应的Euler-Lagrangian方程为

$\begin{array}{c}(A_I{K}^{T}K-c_{1}I-c_2\mathrm{\&Delta;}+c_3{\mathrm{div}}^2{\&dtri;}^2)u^{k+1}=A_I{K}^{T}f+{\mathrm{\&alpha;}}^k+c_1{w}^k-{\mathrm{div}\mathrm{\&beta;}}^k\\ -c_2{\mathrm{div}v}^k+{\mathrm{div}}^2{\mathrm{\&lambda;}}^k-c_3{\mathrm{div}}^2{z}^k.\end{array}$

其中div²为的伴随矩阵.由(1)中选取的为循环边界差分格式,同时利用散度定理中的 $\&lang;p,\&dtri;q\&rang;=-\&lang;\mathrm{div}p,q\&rang;$ 和 $\&lang;s,{\&dtri;}^{2}t\&rang;=\&lang;{\mathrm{div}}^{2}s,t\&rang;,$ 因此可以利用快速FFT变化求解u^k+1.即为：

$u=F^{-1}\left(\frac{A_{I}F\left(K^T\right)F\left(f\right)+F\left({\mathrm{\&alpha;}}^k\right)+c_{1}F\left(w^F\right)-F\left(\mathrm{div}\right)F\left({\mathrm{\&beta;}}^k\right)-c_{2}F\left(\mathrm{div}\right)F\left(v^k\right)+F\left({\mathrm{div}}^2\right)F\left({\mathrm{\&lambda;}}^k\right)-c_{3}F\left({\mathrm{div}}^2\right)F\left(z^k\right)}{(A_{I}F\left(K^{T}K\right)-c_{1}F\left(I\right)-c_{2}F\left(\mathrm{\&Delta;}\right)+c_{3}F\left({\mathrm{div}}^2{\&dtri;}^2\right))}\right)$ 其中F表示快速傅里叶变化,F^-1表示快速傅里叶变化的逆变换.上述算子具体实现在程序中定义为:

conjoI＝conj(psf2otf([1,0],[N,N]))；

conjoDx＝conj(psf2otf([1,-1],[N,N]))；

conjoDy＝conj(psf2otf([1；-1],[N,N]))；

conjoDxx＝conj(psf2otf([1,-2,1],[N,N]))；

conjoDyy＝conj(psf2otf([1；-2；1],[N,N]))；

conjoDxy＝conj(psf2otf([1,-1；-1,1],[N,N]))；

otfH＝psf2otf(K,[N,N])；

F(K^T)F(K)＝abs(otfH).²；

F(I)＝abs(conjoI).²；

F(Δ)＝abs(conjoDx).²+abs(conjoDy).²；

F(K^T)F(f)＝conj(otfH).*fft2(f)；

F(α^k)＝fft2(α^k)；

$F\left(\mathrm{div}\right)F\left({\mathrm{\&beta;}}^k\right)=\mathrm{conjoDx}.*\mathrm{fft}2\left({\mathrm{\&beta;}}_x^k\right)+\mathrm{conjoDy}.*\mathrm{fft}2\left({\mathrm{\&beta;}}_y^k\right);$

$F\left(\mathrm{div}\right)F\left(v^k\right)=\mathrm{conjoDx}.*\mathrm{fft}2\left(v_x^k\right)+\mathrm{conjoDy}.*\mathrm{fft}2\left(v_y^k\right);$

$F\left({\mathrm{div}}^2{\&dtri;}^2\right)=\mathrm{abs}{\left(\mathrm{conjoDxx}\right).}^2+\mathrm{abs}{\left(\mathrm{conjoDyy}\right).}^2+2*\mathrm{abs}{\left(\mathrm{conjoDxy}\right).}^2\text{;}$

$\begin{array}{c}F\left({\mathrm{div}}^2\right)F\left({\mathrm{\&lambda;}}^k\right)=\mathrm{conjoDxx}.*\mathrm{fft}2\left({\mathrm{\&lambda;}}_{\mathrm{xx}}^k\right)+\mathrm{conjoDxy}.*\mathrm{fft}2\left({\mathrm{\&lambda;}}_{\mathrm{xy}}^k\right)\\ +140\mathrm{ptconjoDyx}.*\mathrm{fft}2\left({\mathrm{\&lambda;}}_{\mathrm{yx}}^k\right)+\mathrm{conjoDyy}.*\mathrm{fft}2\left({\mathrm{\&lambda;}}_{\mathrm{yy}}^k\right));\end{array}$

$\begin{array}{c}F\left({\mathrm{div}}^2\right)F\left(z^k\right)=\mathrm{conjoDxx}.*\mathrm{fft}2\left(z_{\mathrm{xx}}^k\right)+\mathrm{conjoDxy}.*\mathrm{fft}2\left(z_{\mathrm{xy}}^k\right)\\ +140\mathrm{ptconjoDyx}.*\mathrm{fft}2\left(z_{\mathrm{yx}}^k\right)+\mathrm{conjoDyy}.*\mathrm{fft}2\left(z_{\mathrm{yy}}^k\right)).\end{array}$

求解子问题w^k+1.即考虑

$\underset{u^{k+1},w,v^k,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u^{k+1},w,v^k,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k)$

舍去与w无关项,则对应的子问题为

$w^{k+1}=\underset{w}{\arg \min }{\&lang;\mathrm{\&alpha;},w\&rang;}_{l^2}+\frac{c_1}{2}{\left|\right|w-u^{k+1}\left|\right|}_{l^2}^2+{\mathrm{\&delta;}}_C\left(w\right).$

该问题可以用投影梯度法快速求解,因此有

$w^{k+1}{=P}_C({\mathrm{\&alpha;}}^k+c_1{u}^{k+1}):=\left\{\begin{array}{cc}{\mathrm{\&alpha;}}^k+c_1{u}^{k+1},& \mathrm{if}({\mathrm{\&alpha;}}^k+c_1{u}^{k+1})\&Element;C\\ 0,& \mathrm{otherwise}.\end{array}\right.$

求解子问题v^k+1.即考虑

$\underset{u^{k+1},w^{k+1},v,z^k}{\min }\underset{{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k}{\max }L(u^{k+1},w^{k+1},v,z^k,{\mathrm{\&alpha;}}^k,{\mathrm{\&beta;}}^k,{\mathrm{\&lambda;}}^k)$

舍去与v无关项,则对应的子问题为

$v^{k+1}=\underset{v}{\arg \min }\frac{c_2}{2}{\left|\right|v-({\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2})\left|\right|}_{l^2}^2+b{\left|v\right|}_{l^1}.$

该问题对应经典的压缩阈值问题,因此其显式解为

$v^{k+1}=\frac{{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}}{{|{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}|}_{l^1}}\max \{{|{\&dtri;u}^{k+1}-\frac{{\mathrm{\&beta;}}^k}{c_2}|}_{l^1}-\frac{b}{c_2},0\}$

求解子问题z^k+1.由于子问题z与v具有类似的形式,因此同样可写成其对应的显式解为

$z^{k+1}=\frac{{{\&dtri;}^{2}u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}}{{|{{\&dtri;}^{2}u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}|}_{l^1}}\max \{{|{{\&dtri;}^{2}u}^{k+1}-\frac{{\mathrm{\&lambda;}}^k}{c_3}|}_{l^1}-\frac{1-b}{c_3},0\}$

由于拉格朗日乘子对应对应的问题是平凡的,因此其更新变量即为

α^k+1＝α^k+c₁(w^k+1-u^k+1)；

${\mathrm{\&beta;}}^{k+1}={\mathrm{\&beta;}}^k{+c}_2(v^{k+1}-\&dtri;u^{k+1});$

${\mathrm{\&lambda;}}^{k+1}={\mathrm{\&lambda;}}^k{+c}_3(z^{k+1}-{\&dtri;}^2{u}^{k+1}).$

终止标准.由交替方向乘子法的收敛性定理可知,当

${\left|\right|w^{k+1}-u^{k+1}\left|\right|}_2^2,{\left|\right|v^{k+1}-{\&dtri;u}^{k+1}\left|\right|}_2^2,{\left|\right|z^{k+1}-{\&dtri;}^2{u}^{k+1}\left|\right|}_2^2$ 收敛到零时,算法的解收敛到初始问题的解.因此在算法中选取终止标准定为当算法终止时,输入复原图像u:＝u^k+1。