基于非线性全变分的去噪算法
摘要: 本文提出一种用于图像去噪的数值最优化约束的算法。通过约束噪声的统计,图像的全变分被最优化。约束使用拉格朗日乘数。该方案通过使用梯度投影法得到。这相当于在一个被约束决定的流形上求解一个与时间有关的偏微分方程。当t→∞t\rightarrow\inftyt→∞,该方案转换成一个稳定的状态,即一个无噪声的图像。该数值算法简单且相应地快。该结果对于高噪声的图像时高水平的。该方法是无损的,会为图像产
摘要: 本文提出一种用于图像去噪的数值最优化约束的算法。通过约束噪声的统计,图像的全变分参数被最优化。约束使用拉格朗日乘数。该方法通过使用梯度投影法得到。这相当于在一个被约束决定的流形上求解一个与时间有关的偏微分方程。当t→∞t\rightarrow\inftyt→∞,该方法转换成一个稳定的状态,即一个无噪声的图像。该数值算法简单且相应地快。该结果对于高噪声的图像时高水平的。该方法是无损的,会为图像产生锐利的边缘。该技术可以理解为首先将图像向平滑变化,然后将图像投影到约束集。
一、介绍
噪声这种随机的失真让图像的处理变得困难。评估噪声中的真实信号的最广泛使用的方法基于最小二乘法。其原理为统计中最小二乘估计是所有可能的图像中最佳的。该方法依赖L2范数,然而已经证实图像的合适范数是全变分【total variation,TV】范数,而不是L2范数。TV范数本质上是L1范数的延申,因此L1估计方法更加接近于图像估计(复原)的条件。边缘全变分函数空间在离散精准估计的方法中有着重要的作用。
本方法通过最小化估计结果的TV范数为图像去噪。一个约束由噪声统计决定、依赖时间的非线性偏导数可以得到一个约束的最优化算法。
以往的方法在图像处理过程之前抑制或去除噪声,这也是本文的方法。然而,TV或L1的思想用来与其他噪声敏感的图像处理任务的去噪结合设计算法。
二、基于非线性偏微分方程的去噪算法
考虑观察得到的强度函数u0(x,y)u_0(x, y)u0(x,y)定义了噪声图像x,y∈Ωx, y \in Ωx,y∈Ω的每个像素值。使用u(x,y)u(x, y)u(x,y)定义期望的清晰图像,则u0(x,y)=u(x,y)+n(x,y)(2.1)u_0(x, y) = u(x, y) + n(x, y) \tag{2.1}u0(x,y)=u(x,y)+n(x,y)(2.1)其中n(x,y)n(x, y)n(x,y)是加性噪声。
当u0u_0u0重建为uuu的过程中,大多传统的变分方法包含L2项以解决线性方程。在二维连续空间中,这种受约束的最优化问题为min∫Ω(uxx+uyy)2(2.2a)min \int_Ω (u_{xx} + u_{yy})^2\tag{2.2a}min∫Ω(uxx+uyy)2(2.2a)约束包括平均值∫Ωu=∫Ωu0(2.2b)\int_Ω u = \int_Ω u_0\tag{2.2b}∫Ωu=∫Ωu0(2.2b)与标准差∫Ω(u−u0)2=σ2(2.2c)\int_Ω (u - u_0)^2 = σ^2 \tag{2.2c}∫Ω(u−u0)2=σ2(2.2c)线性系统使用现代线性代数易解。然而,结果不能令人满意,但优于MEM。
L1范数一般不被使用,因为变分的表达式如∫Ω∣u∣dx\int_Ω |u|dx∫Ω∣u∣dx产生了单一的分布作为系数,这不能单纯的使用代数解决。但是梯度的L1范数更为合适,这就是基本的边缘全变分【Bounded total variation,BV】函数空间。方便起见,我们去除了空间中虚假的震荡,保留了尖锐的信号。
本文的受约束的最优化问题为min∫Ωux2+uy2dxdy(2.3a)min \int_Ω \sqrt{u_x^2 + u_y^2}dxdy\tag{2.3a}min∫Ωux2+uy2dxdy(2.3a)约束包含u0u_0u0,形式同上∫Ωudxdy=∫Ωu0dxdy(2.3b)\int_Ω udxdy = \int_Ω u_0dxdy\tag{2.3b}∫Ωudxdy=∫Ωu0dxdy(2.3b)该约束表明了高斯噪声n(x,y)n(x, y)n(x,y)的均值为0,以及∫012(u−u0)2dxdy=σ2(2.3c)\int_0 \frac{1}{2}(u - u_0)^2 dxdy = σ^2\tag{2.3c}∫021(u−u0)2dxdy=σ2(2.3c)其使用一个先验信息,即高斯噪声n(x,y)n(x, y)n(x,y)的标准差是σσσ。
因此,存在一个线性约束与一个非线性约束。该方法同时考虑了约束的数量与形式。考虑欧拉-拉格朗日方程0=∂∂xuxux2+uy2+∂∂yuyux2+uy2−λ1−λ2(u−u0),x,y∈Ω(2.4a)0 = \frac{∂}{∂x}\frac{u_x}{\sqrt{u_x^2 + u_y^2}} + \frac{∂}{∂y}\frac{u_y}{\sqrt{u_x^2 + u_y^2}} - λ_1 - λ_2(u - u_0), x, y \in Ω \tag{2.4a}0=∂x∂ux2+uy2ux+∂y∂ux2+uy2uy−λ1−λ2(u−u0),x,y∈Ω(2.4a)∂u/∂n=0,u∈∂Ω(2.4b) ∂u/∂n = 0, u \in ∂Ω\tag{2.4b}∂u/∂n=0,u∈∂Ω(2.4b) 该方法使用ttt为参数的抛物方程,或等效的梯度下降法,即解决ut=∂∂xuxux2+uy2+∂∂yuyux2+uy2−λ(u−u0),t>0,x,y∈Ω(2.5a)u_t = \frac{∂}{∂x}\frac{u_x}{\sqrt{u_x^2 + u_y^2}} + \frac{∂}{∂y}\frac{u_y}{\sqrt{u_x^2 + u_y^2}} - λ_(u - u_0), t>0, x, y \in Ω \tag{2.5a}ut=∂x∂ux2+uy2ux+∂y∂ux2+uy2uy−λ(u−u0),t>0,x,y∈Ω(2.5a) u(x,y,0)(2.5b)u(x, y, 0) \tag{2.5b}u(x,y,0)(2.5b) ∂u/∂n=0,u∈∂Ω(2.5c) ∂u/∂n = 0, u \in ∂Ω \tag{2.5c}∂u/∂n=0,u∈∂Ω(2.5c)注意到第一约束被去除了,因为当u(x,y,0)=u0(x,y)u(x, y, 0) = u_0(x, y)u(x,y,0)=u0(x,y),其会被该方程自动执行。随着ttt的增加,图像将接近于无噪声。
考虑计算λ(t)λ(t)λ(t),通过不断的对式(2.5a)(2.5a)(2.5a)乘(u−u0)(u - u_0)(u−u0)并在整个ΩΩΩ上积分,其会达到一个稳定值,式(2.5a)(2.5a)(2.5a)的左侧将会消掉,得到λ=−12σ2∫[ux2+uy2−((u0)xuxux2+uy2+(u0)yuyux2+uy2)]dxdy(2.6)λ = -\frac{1}{2σ^2}\int[\sqrt{u_x^2 + u_y^2} - (\frac{(u_0)_xu_x}{\sqrt{u_x^2 + u_y^2}} + \frac{(u_0)_yu_y}{\sqrt{u_x^2 + u_y^2}})]dxdy \tag{2.6}λ=−2σ21∫[ux2+uy2−(ux2+uy2(u0)xux+ux2+uy2(u0)yuy)]dxdy(2.6)这给出了λ(t)λ(t)λ(t)的动态值,并当t→∞t\rightarrow \inftyt→∞时收敛。
该方法引入二维空间如下,定义xi=ih,yj=jh,i,j=0,1,...,N,Nh=1tn=nΔt,n=0,1,...uijn=u(xi,yj,tn)uij0=u0(ih,jh)+σφ(ih,jh)(2.7)x_i = ih, y_j = jh, i, j = 0, 1, ..., N, Nh = 1 \\ t_n = nΔt, n = 0, 1, ... \\ u_{ij}^n = u(x_i, y_j, t_n) \\ u_{ij}^0 = u_0(ih, jh) + σφ(ih, jh) \tag{2.7}xi=ih,yj=jh,i,j=0,1,...,N,Nh=1tn=nΔt,n=0,1,...uijn=u(xi,yj,tn)uij0=u0(ih,jh)+σφ(ih,jh)(2.7)选择初始数据,使约束被初始的满足,即φφφ的均值为0,而L2范数为1。
根据(2.5)(2.5)(2.5),(2.6)(2.6)(2.6)式,数值近似为uijn+1=uijn+Δth[Δ−x(Δ+xuijn(Δ+xuijn)2+(m(Δ+yuijn,Δ−yuijn))2)1/2)+Δ−y(Δ+yuijn(Δ+yuijn)2+(m(Δ+xuijn,Δ−xuijn))2)1/2)]−Δtλn(uijn−u0(ih,jh)),i,j=1,...,N(2.8a)\begin{aligned} u_{ij}^{n+1} =& u_{ij}^n + \frac{Δt}{h}[Δ^x_{-}(\frac{Δ^x_{+}u_{ij}^n}{(Δ^x_{+}u_{ij}^n)^2+(m(Δ^y_{+}u_{ij}^n, Δ^y_{-}u_{ij}^n))^2)^{1/2}}) \\&+ Δ^y_{-}(\frac{Δ^y_{+}u_{ij}^n}{(Δ^y_{+}u_{ij}^n)^2+(m(Δ^x_{+}u_{ij}^n, Δ^x_{-}u_{ij}^n))^2)^{1/2}})] \\& - Δtλ^n(u_{ij}^n - u_0(ih, jh)), i, j = 1, ..., N\end{aligned} \tag{2.8a}uijn+1=uijn+hΔt[Δ−x((Δ+xuijn)2+(m(Δ+yuijn,Δ−yuijn))2)1/2Δ+xuijn)+Δ−y((Δ+yuijn)2+(m(Δ+xuijn,Δ−xuijn))2)1/2Δ+yuijn)]−Δtλn(uijn−u0(ih,jh)),i,j=1,...,N(2.8a)及边缘条件u0jn=u1jn,uNjn=uN−1,jn,ui0n=ui1n,uiNn=ui,N−1n(2.8b)u_{0j}^n = u_{1j}^n, u_{Nj}^n = u_{N-1, j}^n, u_{i0}^n = u_{i1}^n, u_{iN}^n = u_{i, N-1}^n\tag{2.8b}u0jn=u1jn,uNjn=uN−1,jn,ui0n=ui1n,uiNn=ui,N−1n(2.8b)其中Δ∓x=∓(ui∓1,j−uij),similiarly for Δ∓ym(a,b)=(sgn(a)+sgn(b))min(∣a∣,∣b∣)/2,sgn(x)=x/∣x∣,x≠0λn=−h2σ2[∑i,j(Δ+xuijn)2+(Δ+yuijn)2−(Δ+xuij0)(Δ+xuijn)(Δ+xuijn)2+(Δ+yuijn)2−(Δ+yuij0)(Δ+yuijn)(Δ+xuijn)2+(Δ+yuijn)2)](2.9)\begin{aligned}Δ^x_{\mp} =& \mp(u_{i\mp1, j} - u_{ij}), similiarly\ for\ Δ^y_{\mp} \\ m(a, b) =& (sgn(a) + sgn(b))min(|a|, |b|)/2, sgn(x) = x/|x|,x \ne 0 \\ λ^n =& -\frac{h}{2σ^2}[\sum_{i, j}(\sqrt{Δ^x_{+}u_{ij}^n)^2 + (Δ^y_{+}u_{ij}^n)^2}\\ &- \frac{(Δ^x_{+}u_{ij}^0)(Δ^x_{+}u_{ij}^n)}{\sqrt{(Δ^x_{+}u_{ij}^n)^2 + (Δ^y_{+}u_{ij}^n)^2}} - \frac{(Δ^y_{+}u_{ij}^0)(Δ^y_{+}u_{ij}^n)}{\sqrt{(Δ^x_{+}u_{ij}^n)^2 + (Δ^y_{+}u_{ij}^n)^2}})]\end{aligned} \tag{2.9}Δ∓x=m(a,b)=λn=∓(ui∓1,j−uij),similiarly for Δ∓y(sgn(a)+sgn(b))min(∣a∣,∣b∣)/2,sgn(x)=x/∣x∣,x=0−2σ2h[i,j∑(Δ+xuijn)2+(Δ+yuijn)2−(Δ+xuijn)2+(Δ+yuijn)2(Δ+xuij0)(Δ+xuijn)−(Δ+xuijn)2+(Δ+yuijn)2(Δ+yuij0)(Δ+yuijn))](2.9)为了稳定,步长被限制为Δt/h2≤cΔt/h^2 \le cΔt/h2≤c
三、结果
图1:(a)分辨率图;(b)含有噪声的分辨率图,SNR=1.0;(c)维纳滤波复原;(d)TV复原
更多推荐


所有评论(0)