TV（Total Variation）降噪
一般情况下，利用全变分降噪能够很好地去除噪声，但是存在着计算复杂度大等问题。这里利用Split Bregman迭代进行TV降噪，不仅实现简单，而且能够大大降低计算的复杂度。
考虑各向异性情况（Anisotropic case）
对于优化问题
min u ∣ ∇ x u ∣ + ∣ ∇ y u ∣ + μ 2 ∥ u − f ∥ 2 2 \underset{u}{\mathop{\min }}\,\left| {
{\nabla }_{x}}u \right|+\left| {
{\nabla }_{y}}u \right|+\frac{\mu }{2}\left\| u-f \right\|_{2}^{2} umin∣∇xu∣+∣∇yu∣+2μ∥u−f∥22
为了能够使用Split Bregman迭代，需要将上述问题转化为
{ min u ∣ d x ∣ + ∣ d y ∣ + μ 2 ∥ u − f ∥ 2 2 s . t . d x = ∇ x u , d y = ∇ y u \left\{ \begin{aligned} & \underset{u}{\mathop{\min }}\,\left| {
{d}_{x}} \right|+\left| {
{d}_{y}} \right|+\frac{\mu }{2}\left\| u-f \right\|_{2}^{2} \\ & s.t.\ \ {
{d}_{x}}=\ {
{\nabla }_{x}}u,\ \ {
{d}_{y}}={
{\nabla }_{y}}u \\ \end{aligned} \right. ⎩⎨⎧umin∣dx∣+∣dy∣+2μ∥u−f∥22s.t. dx= ∇xu, dy=∇yu
则
min d x , d y , u ∣ d x ∣ + ∣ d y ∣ + μ 2 ∥ μ − f ∥ 2 2 + λ 2 ∥ d x − ∇ x u ∥ 2 2 + λ 2 ∥ d y − ∇ y u ∥ 2 2 \underset{
{
{d}_{x}},{
{d}_{y}},u}{\mathop{\min }}\,\left| {
{d}_{x}} \right|+\left| {
{d}_{y}} \right|+\frac{\mu }{2}\left\| \mu -f \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{x}}-{
{\nabla }_{x}}u \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{y}}-{
{\nabla }_{y}}u \right\|_{2}^{2} dx,dy,umin∣dx∣+∣dy∣+2μ∥μ−f∥22+2λ∥dx−∇xu∥22+2λ∥dy−∇yu∥22
利用Bregman迭代可以得到
min d x , d y , u ∣ d x ∣ + ∣ d y ∣ + μ 2 ∥ μ − f ∥ 2 2 + λ 2 ∥ d x − ∇ x u − b x k ∥ 2 2 + λ 2 ∥ d y − ∇ y u − b y k ∥ 2 2 \underset{
{
{d}_{x}},{
{d}_{y}},u}{\mathop{\min }}\,\left| {
{d}_{x}} \right|+\left| {
{d}_{y}} \right|+\frac{\mu }{2}\left\| \mu -f \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{x}}-{
{\nabla }_{x}}u-b_{x}^{k} \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{y}}-{
{\nabla }_{y}}u-b_{y}^{k} \right\|_{2}^{2} dx,dy,umin∣dx∣+∣dy∣+2μ∥μ−f∥22+2λ∥∥dx−∇xu−bxk∥∥22+2λ∥∥dy−∇yu−byk∥∥22
为了解决上述问题的求解，则需要先解决
u k + 1 = min u μ 2 ∥ μ − f ∥ 2 2 + λ 2 ∥ d x − ∇ x u − b x k ∥ 2 2 + λ 2 ∥ d y − ∇ y u − b y k ∥ 2 2 {
{u}^{k+1}}=\underset{u}{\mathop{\min }}\,\frac{\mu }{2}\left\| \mu -f \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{x}}-{
{\nabla }_{x}}u-b_{x}^{k} \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{y}}-{
{\nabla }_{y}}u-b_{y}^{k} \right\|_{2}^{2} uk+1=umin2μ∥μ−f∥22+2λ∥∥dx−∇xu−bxk∥∥22+2λ∥∥dy−∇yu−byk∥∥22
其中， ( μ I − λ Δ ) u k + 1 = μ f + λ ∇ x T ( d x k − b x k ) + λ ∇ y T ( d y k − b y k ) \left( \mu I-\lambda \Delta \right){
{u}^{k+1}}=\mu f+\lambda \nabla _{x}^{T}\left( d_{x}^{k}-b_{x}^{k} \right)+\lambda \nabla _{y}^{T}\left( d_{y}^{k}-b_{y}^{k} \right) (μI−λΔ)uk+1=μf+λ∇xT(dxk−bxk)+λ∇yT(dyk−byk)
为了实现**的效率，需要利用快速迭代的方法进行近似求解。利用高斯-赛德尔迭代（Gauss-Seidel）。
令 u i , j k + 1 = G i , j k u_{i,j}^{k+1}=G_{i,j}^{k} ui,jk+1=Gi,jk
则，
G i , j k = λ μ + 4 λ ( u i + 1 , j k + u i − 1 , j k + u i , j + 1 k + u i , j − 1 k + d x , i − 1 , j k − d x , i , j k + d y , i , j − 1 k − d y , i , j k − b x , i − 1 , j k + b x , i , j k − b y , i , j − 1 k + b y , i , j k ) + λ μ + 4 λ f i , j \begin{aligned} & G_{i,j}^{k}\text{=}\frac{\lambda }{\mu +4\lambda }(u_{i+1,j}^{k}+u_{i-1,j}^{k}+u_{i,j+1}^{k}+u_{i,j-1}^{k} \\ & +d_{x,i-1,j}^{k}-d_{x,i,j}^{k}+d_{y,i,j-1}^{k}-d_{y,i,j}^{k}-b_{x,i-1,j}^{k}+b_{x,i,j}^{k}-b_{y,i,j-1}^{k}+b_{y,i,j}^{k})+\frac{\lambda }{\mu +4\lambda }{
{f}_{i,j}} \\ \end{aligned} Gi,jk=μ+4λλ(ui+1,jk+ui−1,jk+ui,j+1k+ui,j−1k+dx,i−1,jk−dx,i,jk+dy,i,j−1k−dy,i,jk−bx,i−1,jk+bx,i,jk−by,i,j−1k+by,i,jk)+μ+4λλfi,j
因此，最终利用Split Bregman迭代进行各向异性TV降噪可以写作
初始化： u 0 = f , d x 0 = d y 0 = b x 0 = b y 0 = 0 {
{u}^{0}}=f,d_{x}^{0}=d_{y}^{0}=b_{x}^{0}=b_{y}^{0}=0 u0=f,dx0=dy0=bx0=by0=0
While ∥ u k − u k − 1 ∥ 2 2 > t h r e s h o l d \left\| {
{u}^{k}}-{
{u}^{k-1}} \right\|_{2}^{2}>threshold ∥∥uk−uk−1∥∥22>threshold
u k + 1 = G k {
{u}^{k+1}}={
{G}^{k}} uk+1=Gk
d x k + 1 = s h r i n k ( ∇ x u k + 1 + b x k , 1 / λ ) d_{x}^{k+1}=shrink\left( {
{\nabla }_{x}}{
{u}^{k+1}}+b_{x}^{k},1/\lambda \right) dxk+1=shrink(∇xuk+1+bxk,1/λ)
d y k + 1 = s h r i n k ( ∇ y u k + 1 + b y k , 1 / λ ) d_{y}^{k+1}=shrink\left( {
{\nabla }_{y}}{
{u}^{k+1}}+b_{y}^{k},1/\lambda \right) dyk+1=shrink(∇yuk+1+byk,1/λ)
b x k + 1 = b x k + ( ∇ x u k + 1 − d x k + 1 ) b_{x}^{k+1}=b_{x}^{k}+\left( {
{\nabla }_{x}}{
{u}^{k+1}}-d_{x}^{k+1} \right) bxk+1=bxk+(∇xuk+1−dxk+1)
b y k + 1 = b y k + ( ∇ y u k + 1 − d y k + 1 ) b_{y}^{k+1}=b_{y}^{k}+\left( {
{\nabla }_{y}}{
{u}^{k+1}}-d_{y}^{k+1} \right) byk+1=byk+(∇yuk+1−dyk+1)
end
考虑各向同性情况（Isotropic case）
对于优化问题
min u ( ∇ x u ) i 2 + ( ∇ y u ) i 2 + μ 2 ∥ u − f ∥ 2 2 \underset{u}{\mathop{\min }}\,\sqrt{\left( {
{\nabla }_{x}}u \right)_{i}^{2}+\left( {
{\nabla }_{y}}u \right)_{i}^{2}}+\frac{\mu }{2}\left\| u-f \right\|_{2}^{2} umin(∇xu)i2+(∇yu)i2+2μ∥u−f∥22
同样地，将上述问题变为 l 1 {
{l}_{1}} l1范数和 l 2 {
{l}_{2}} l2范数两部分的组合
min d x , d y , u ∥ ( d x , d y ) ∥ 2 + μ 2 ∥ μ − f ∥ 2 2 + λ 2 ∥ d x − ∇ x u − b x ∥ 2 2 + λ 2 ∥ d y − ∇ y u − b y ∥ 2 2 \underset{
{
{d}_{x}},{
{d}_{y}},u}{\mathop{\min }}\,{
{\left\| \left( {
{d}_{x}},{
{d}_{y}} \right) \right\|}_{2}}+\frac{\mu }{2}\left\| \mu -f \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{x}}-{
{\nabla }_{x}}u-{
{b}_{x}} \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{y}}-{
{\nabla }_{y}}u-{
{b}_{y}} \right\|_{2}^{2} dx,dy,umin∥(dx,dy)∥2+2μ∥μ−f∥22+2λ∥dx−∇xu−bx∥22+2λ∥dy−∇yu−by∥22
其中 ∥ ( d x , d y ) ∥ 2 = ∑ i , j d x , i , j 2 + d y , i , j 2 {
{\left\| \left( {
{d}_{x}},{
{d}_{y}} \right) \right\|}_{2}}\text{=}\sum\limits_{i,j}^{
{}}{\sqrt{d_{x,i,j}^{2}+d_{y,i,j}^{2}}} ∥(dx,dy)∥2=i,j∑dx,i,j2+dy,i,j2
同理，在解决上述时需要提前解决
( d x k + 1 , d y k + 1 ) = min d x , d y ∥ ( d x , d y ) ∥ 2 + λ 2 ∥ d x − ∇ x u − b x ∥ 2 2 + λ 2 ∥ d y − ∇ y u − b y ∥ 2 2 \left( d_{x}^{k+1},d_{y}^{k+1} \right)=\underset{
{
{d}_{x}},{
{d}_{y}}}{\mathop{\min }}\,{
{\left\| \left( {
{d}_{x}},{
{d}_{y}} \right) \right\|}_{2}}+\frac{\lambda }{2}\left\| {
{d}_{x}}-{
{\nabla }_{x}}u-{
{b}_{x}} \right\|_{2}^{2}+\frac{\lambda }{2}\left\| {
{d}_{y}}-{
{\nabla }_{y}}u-{
{b}_{y}} \right\|_{2}^{2} (dxk+1,dyk+1)=dx,dymin∥(dx,dy)∥2+2λ∥dx−∇xu−bx∥22+2λ∥dy−∇yu−by∥22
在这里引入一个广义的收缩公式
d x k + 1 = max ( s k − 1 / λ , 0 ) ∇ x u k + b x k s k d_{x}^{k+1}=\max \left( {
{s}^{k}}-1/\lambda ,0 \right)\frac{
{
{\nabla }_{x}}{
{u}^{k}}+b_{x}^{k}}{
{
{s}^{k}}} dxk+1=max(sk−1/λ,0)sk∇xuk+bxk
d y k + 1 = max ( s k − 1 / λ , 0 ) ∇ y u k + b y k s k d_{y}^{k+1}=\max \left( {
{s}^{k}}-1/\lambda ,0 \right)\frac{
{
{\nabla }_{y}}{
{u}^{k}}+b_{y}^{k}}{
{
{s}^{k}}} dyk+1=max(sk−1/λ,0)sk∇yuk+byk
其中 s k = ∣ ∇ x u k + b x k ∣ 2 + ∣ ∇ y u k + b y k ∣ 2 {
{s}^{k}}=\sqrt{
{
{\left| {
{\nabla }_{x}}{
{u}^{k}}+b_{x}^{k} \right|}^{2}}+{
{\left| {
{\nabla }_{y}}{
{u}^{k}}+b_{y}^{k} \right|}^{2}}} sk=∣∇xuk+bxk∣2+∣∣∇yuk+byk∣∣2
因此，最终利用Split Bregman迭代进行各向同性TV降噪可以写作
初始化： u 0 = f , d x 0 = d y 0 = b x 0 = b y 0 = 0 {
{u}^{0}}=f,d_{x}^{0}=d_{y}^{0}=b_{x}^{0}=b_{y}^{0}=0 u0=f,dx0=dy0=bx0=by0=0
While ∥ u k − u k − 1 ∥ 2 2 > t h r e s h o l d \left\| {
{u}^{k}}-{
{u}^{k-1}} \right\|_{2}^{2}>threshold ∥∥uk−uk−1∥∥22>threshold
u k + 1 = G i , j k {
{u}^{k+1}}=G_{i,j}^{k} uk+1=Gi,jk
d x k + 1 = max ( s k − 1 / λ , 0 ) ∇ x u k + b x k s k d_{x}^{k+1}=\max \left( {
{s}^{k}}-1/\lambda ,0 \right)\frac{
{
{\nabla }_{x}}{
{u}^{k}}+b_{x}^{k}}{
{
{s}^{k}}} dxk+1=max(sk−1/λ,0)sk∇xuk+bxk
d y k + 1 = max ( s k − 1 / λ , 0 ) ∇ y u k + b y k s k d_{y}^{k+1}=\max \left( {
{s}^{k}}-1/\lambda ,0 \right)\frac{
{
{\nabla }_{y}}{
{u}^{k}}+b_{y}^{k}}{
{
{s}^{k}}} dyk+1=max(sk−1/λ,0)sk∇yuk+byk
b x k + 1 = b x k + ( ∇ x u k + 1 − d x k + 1 ) b_{x}^{k+1}=b_{x}^{k}+\left( {
{\nabla }_{x}}{
{u}^{k+1}}-d_{x}^{k+1} \right) bxk+1=bxk+(∇xuk+1−dxk+1)
b y k + 1 = b y k + ( ∇ y u k + 1 − d y k + 1 ) b_{y}^{k+1}=b_{y}^{k}+\left( {
{\nabla }_{y}}{
{u}^{k+1}}-d_{y}^{k+1} \right) byk+1=byk+(∇yuk+1−dyk+1)
end
仿真参数
| 参数名称 | 参数值 |
|---|---|
| 图像大小 | 512 |
| μ \mu μ | 10 |
仿真结果
从图中可以看出，经过两种方法处理后的图像均能够很好地降噪，更进一步地评价两种方法的优劣可以采用PSNR，MSE等指标，这里就不再叙述。
仿真代码如下：
主函数
clear; close all; clc; N = 512; %图像大小 n = N^2; ori_image = double(imread('Lena.png')); %读入原始图像 noise_image = ori_image(:) + 0.09*max(ori_image(:))*randn(n,1); %加噪后的图像 mu = 10; %正则化参数 noise_image=reshape(noise_image,N,N); ATV_image = ATV_ROF(noise_image,mu); ITV_image = ITV_ROF(noise_image,mu); figure; colormap gray; subplot(221); imagesc(ori_image); axis image; title('Original Image'); subplot(222); imagesc(reshape(noise_image,N,N)); axis image; title('Noisy Image'); subplot(223); imagesc(reshape(ATV_image,N,N)); axis image; title('Anisotropic TV denoising'); subplot(224); imagesc(reshape(ITV_image,N,N)); axis image; title('Isotropic TV denoising'); 讯享网
ATV_ROF.m
讯享网function u=ATV_ROF(f,mu) %============================================ % Anisotropic ROF denoising % This function performs the minimization of % u=arg min |D_x u|+|D_y u|+0.5*mu*||u-f||_2^2 % by the Split Bregman Iteration % f = noisy image % mu = regularization parameter % Reference：Goldstein T , Osher S . The Split Bregman Method for % L1-Regularized Problems[J]. %SIAM Journal on Imaging Sciences, 2009, 2(2):323-343. %============================================ [M,N]=size(f); f=double(f); dx=zeros(M,N); dy=zeros(M,N); bx=zeros(M,N); by=zeros(M,N); u=f; Z=zeros(M,N); Mask=zeros(M,N); Mask(1,1) = 1; FMask=fft2(Mask); lambda=100; %Fourier Laplacian mask initialization D = zeros(M,N); D([end,1,2],[end,1,2]) = [0,1,0;1,-4,1;0,1,0]; FD=fft2(D); %Fourier constant initialization FW=((mu/lambda)*abs(FMask).^2-real(FD)).^-1; FF=(mu/lambda)*conj(FMask).*fft2(f); err=1; tol=1e-3; while err>tol tx=dx-bx; ty=dy-by; up=u; %Update u u=real(ifft2(FW.*(FF-fft2(tx-tx(:,[1,1:N-1])+ty-ty([1,1:M-1],:))))); ux=u-u(:,[1,1:N-1]); uy=u-u([1,1:M-1],:); %Update dx tmpx=ux+bx; dx=sign(tmpx).*max(Z,abs(tmpx)-1/lambda); %Update dy tmpy=uy+by; dy=sign(tmpy).*max(Z,abs(tmpy)-1/lambda); %Update bx and by bx=tmpx-dx; by=tmpy-dy; err=sum(sum((up-u).^2)); end
ITV_ROF.m
function u=ITV_ROF(f,mu) % Isotropic ROF denoising % This function performs the minimization of % u=arg min sqrt(|D_x u|^2+|D_y u|^2)+0.5*mu*||u-f||_2^2 % by the Split Bregman Iteration % f = noisy image % mu = regularization parameter % Reference：Goldstein T , Osher S . The Split Bregman Method for % L1-Regularized Problems[J]. %SIAM Journal on Imaging Sciences, 2009, 2(2):323-343. %============================================ [M,N]=size(f); f=double(f); dx=zeros(M,N); dy=zeros(M,N); bx=zeros(M,N); by=zeros(M,N); s=zeros(M,N); u=f; Z=zeros(M,N); lambda=100; Mask=zeros(M,N); Mask(1,1) = 1; FMask=fft2(Mask); %Fourier Laplacian mask initialization D = zeros(M,N); D([end,1,2],[end,1,2]) = [0,1,0;1,-4,1;0,1,0]; FD=fft2(D); %Fourier constant initialization FW=((mu/lambda)*abs(FMask).^2-real(FD)).^-1; FF=(mu/lambda)*conj(FMask).*fft2(f); err=1; tol=1e-3; while err>tol %while K<Niter, tx=dx-bx; ty=dy-by; up=u; %Update u u=real(ifft2(FW.*(FF-fft2(tx-tx(:,[1,1:N-1])+ty-ty([1,1:M-1],:))))); ux=u-u(:,[1,1:N-1]); uy=u-u([1,1:M-1],:); tmpx=ux+bx; tmpy=uy+by; s=sqrt(tmpx.^2+tmpy.^2); thresh=max(Z,s-1/lambda)./max(1e-12,s); %Update dx dx=thresh.*tmpx; %Update dy dy=thresh.*tmpy; %Update bx and by bx=tmpx-dx; by=tmpy-dy; err=sum(sum((up-u).^2)); end 参考文献
[1] Goldstein T, Osher S. The split Bregman method for L1-regularized problems[J]. SIAM journal on imaging sciences, 2009, 2(2): 323-343.
[2]Gilles J, Osher S. Bregman implementation of Meyer’s G-norm for cartoon+ textures decomposition[J]. UCLA Cam Report, 2011: 11-73.
[3] Bush J. Bregman algorithms[J]. Senior Thesis. University of California, Santa Barbara, 2011.
[4] Yin W, Osher S, Goldfarb D, et al. Bregman iterative algorithms for ℓ 1 \ell_1 ℓ1-minimization with applications to compressed sensing[J]. SIAM Journal on Imaging sciences, 2008, 1(1): 143-168.
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容，请联系我们，一经查实，本站将立刻删除。
如需转载请保留出处：https://51itzy.com/kjqy/20276.html