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

摘要

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

主权项

一种变分PDE图像保边保结构复原方法，其特征在于，包括以下步骤：(1)输入一副大小为N×N的退化图像f(x,y),令f(i,j)表示图像在第i行第j列的像素值,利用Matlab中的diff函数定义函数g(x,y)的梯度算子和海森算子${▿}^{2}g={({▿}_{\mathrm{xx}}g,{▿}_{\mathrm{xy}}g;{▿}_{\mathrm{yx}}g,{▿}_{\mathrm{yy}}g)}^T:$${▿}_{x}g=[\mathrm{diff}(g,\mathrm{1,2}),g(:,1)-g(:,N)];$${▿}_{y}g=[\mathrm{diff}(g,1,1),g(1,:)-g(m,:)];$${▿}_{\mathrm{xx}}g=[{▿}_{x}g(:,1)-{▿}_{x}g(:,N);\mathrm{diff}({▿}_{x}g,\mathrm{1,2})];$${▿}_{\mathrm{yy}}g=[{▿}_{y}g(1,:)-{▿}_{y}g(N,:);\mathrm{diff}({▿}_{y}g,1,1)];$${▿}_{\mathrm{xy}}g=[\mathrm{diff}({▿}_{x}g,1,1,):{▿}_{x}g(1,:)-{▿}_{x}g(N,:)].$选取初始值u⁰＝f,w₀＝0,v⁰＝0,z⁰＝0,α⁰＝0,β⁰＝0,λ⁰＝0.；(2)利求解离散化模型$\underset{u\in 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}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}\left|{▿u}_{i,j}\right|=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}\sqrt{{\left({▿}_x{u}_{i,j}\right)}^2+{\left({▿}_y{u}_{i,j}\right)}^2},$${|{▿}^2▿|}_{l^1}=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}\left|{{▿}^{2}u}_{i,j}\right|=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}\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}.$即引入约束条件$w=u,v=▿u,z={▿}^{2}u.$将其转化为下面极大极小的鞍点问题$\underset{u,w,v,z}{\min }\underset{\alpha ,\beta ,\lambda }{\max }L(u,w,v,z,\alpha ,\beta ,\lambda ),$其中$\begin{array}{c}L(u,w,v,z,\alpha ,\beta ,\lambda )=\frac{A_I}{2}{\left|\right|\mathrm{Ku}-f\left|\right|}_{l^2}^2{+⟨\alpha ,w-u⟩}_{l^2}+\frac{c_1}{2}{\left|\right|w-u\left|\right|}_{l^2}^2+{\delta }_C\left(w\right)+{\left|\mathrm{bv}\right|}_{l^1}\\ +{⟨\beta ,v-▿u⟩}_{l^2}+\frac{c_2}{2}{\left|\right|v-▿u\left|\right|}_{l^2}^2+{\left|(1-b)z\right|}_{l^1}+{⟨\lambda ,z-{▿}^{2}u⟩}_{l^2}+\frac{c_3}{2}{\left|\right|z-{▿}^{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{{\alpha }^k,{\beta }^k,{\lambda }^k}{\max }L(u,w^k,v^k,z^k,{\alpha }^k,{\beta }^k,{\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-{⟨{\alpha }^k,u⟩}_{l^2}+\frac{c_1}{2}{\left|\right|w^k-u\left|\right|}_{l^2}^2-{⟨{\beta }^k,▿u⟩}_{l^2}\\ +\frac{c_2}{2}{\left|\right|v^k-▿u\left|\right|}_{l^2}^2-{⟨{\lambda }^k,{▿}^{2}u⟩}_{l^2}+\frac{c_3}{2}{\left|\right|z^k-{▿}^{2}u\left|\right|}_{l^2}^2.\end{array}$求出其对应的显式解，即对应的Euler‑Lagrangian方程为$\begin{array}{c}(A_I{K}^{T}K-c_{1}I-c_2\Delta +c_3{\mathrm{div}}^2{▿}^2)u^{k+1}=A_I{K}^{T}f+{\alpha }^k+c_1{w}^k-{\mathrm{div}\beta }^k\\ -c_2{\mathrm{div}v}^k+{\mathrm{div}}^2{\lambda }^k-c_3{\mathrm{div}}^2{z}^k.\end{array}$其中div²为的伴随矩阵.由(1)中选取的为循环边界差分格式,同时利用散度定理中的$⟨p,▿q⟩=-⟨\mathrm{div}p,q⟩$和$⟨s,{▿}^{2}t⟩=⟨{\mathrm{div}}^{2}s,t⟩,$因此可以利用快速FFT变化求解u^k+1.即为：$u=F^{-1}\left(\frac{A_{I}F\left(K^T\right)F\left(f\right)+F\left({\alpha }^k\right)+c_{1}F\left(w^F\right)-F\left(\mathrm{div}\right)F\left({\beta }^k\right)-c_{2}F\left(\mathrm{div}\right)F\left(v^k\right)+F\left({\mathrm{div}}^2\right)F\left({\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(\Delta \right)+c_{3}F\left({\mathrm{div}}^2{▿}^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({\beta }^k\right)=\mathrm{conjoDx}.*\mathrm{fft}2\left({\beta }_x^k\right)+\mathrm{conjoDy}.*\mathrm{fft}2\left({\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{▿}^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({\lambda }^k\right)=\mathrm{conjoDxx}.*\mathrm{fft}2\left({\lambda }_{\mathrm{xx}}^k\right)+\mathrm{conjoDxy}.*\mathrm{fft}2\left({\lambda }_{\mathrm{xy}}^k\right)\\ +140\mathrm{ptconjoDyx}.*\mathrm{fft}2\left({\lambda }_{\mathrm{yx}}^k\right)+\mathrm{conjoDyy}.*\mathrm{fft}2\left({\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{{\alpha }^k,{\beta }^k,{\lambda }^k}{\max }L(u^{k+1},w,v^k,z^k,{\alpha }^k,{\beta }^k,{\lambda }^k)$舍去与w无关项,则对应的子问题为$w^{k+1}=\underset{w}{\arg \min }{⟨\alpha ,w⟩}_{l^2}+\frac{c_1}{2}{\left|\right|w-u^{k+1}\left|\right|}_{l^2}^2+{\delta }_C\left(w\right).$该问题可以用投影梯度法快速求解,因此有$w^{k+1}{=P}_C({\alpha }^k+c_1{u}^{k+1}):=\left\{\begin{array}{cc}{\alpha }^k+c_1{u}^{k+1},& \mathrm{if}({\alpha }^k+c_1{u}^{k+1})\in C\\ 0,& \mathrm{otherwise}.\end{array}\right.$求解子问题v^k+1.即考虑$\underset{u^{k+1},w^{k+1},v,z^k}{\min }\underset{{\alpha }^k,{\beta }^k,{\lambda }^k}{\max }L(u^{k+1},w^{k+1},v,z^k,{\alpha }^k,{\beta }^k,{\lambda }^k)$舍去与v无关项,则对应的子问题为$v^{k+1}=\underset{v}{\arg \min }\frac{c_2}{2}{\left|\right|v-({▿u}^{k+1}-\frac{{\beta }^k}{c_2})\left|\right|}_{l^2}^2+b{\left|v\right|}_{l^1}.$该问题对应经典的压缩阈值问题,因此其显式解为$v^{k+1}=\frac{{▿u}^{k+1}-\frac{{\beta }^k}{c_2}}{{|{▿u}^{k+1}-\frac{{\beta }^k}{c_2}|}_{l^1}}\max \{{|{▿u}^{k+1}-\frac{{\beta }^k}{c_2}|}_{l^1}-\frac{b}{c_2},0\}$求解子问题z^k+1.由于子问题z与v具有类似的形式,因此同样可写成其对应的显式解为$z^{k+1}=\frac{{{▿}^{2}u}^{k+1}-\frac{{\lambda }^k}{c_3}}{{|{{▿}^{2}u}^{k+1}-\frac{{\lambda }^k}{c_3}|}_{l^1}}\max \{{|{{▿}^{2}u}^{k+1}-\frac{{\lambda }^k}{c_3}|}_{l^1}-\frac{1-b}{c_3},0\}$由于拉格朗日乘子对应对应的问题是平凡的,因此其更新变量即为α^k+1＝α^k+c₁(w^k+1‑u^k+1)；${\beta }^{k+1}={\beta }^k{+c}_2(v^{k+1}-▿u^{k+1});$${\lambda }^{k+1}={\lambda }^k{+c}_3(z^{k+1}-{▿}^2{u}^{k+1}).$终止标准.由交替方向乘子法的收敛性定理可知,当${\left|\right|w^{k+1}-u^{k+1}\left|\right|}_2^2,{\left|\right|v^{k+1}-{▿u}^{k+1}\left|\right|}_2^2,{\left|\right|z^{k+1}-{▿}^2{u}^{k+1}\left|\right|}_2^2$收敛到零时,算法的解收敛到初始问题的解.因此在算法中选取终止标准定为当算法终止时,输入复原图像u:＝u^k+1。