Split Bregman迭代进行TV降噪介绍及MATLAB实现

最新推荐文章于 2023-08-23 19:05:21 发布

YHCANDOU

最新推荐文章于 2023-08-23 19:05:21 发布

阅读量4.3k

点赞数 12

分类专栏：图像信号处理文章标签： matlab 算法

本文链接：https://blog.csdn.net/yhcwjh/article/details/114293548

版权

图像信号处理专栏收录该内容

6 篇文章 3 订阅

订阅专栏

TV（Total Variation）降噪

一般情况下，利用全变分降噪能够很好地去除噪声，但是存在着计算复杂度大等问题。这里利用Split Bregman迭代进行TV降噪，不仅实现简单，而且能够大大降低计算的复杂度。

考虑各向异性情况（Anisotropic case）

对于优化问题
$\underset{u}{\mathop{\min }}\,\left| {{\nabla }_{x}}u \right|+\left| {{\nabla }_{y}}u \right|+\frac{\mu }{2}\left\| u-f \right\|_{2}^{2}$
为了能够使用Split Bregman迭代，需要将上述问题转化为
$\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.$
则
$\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}$
利用Bregman迭代可以得到
$\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}$
为了解决上述问题的求解，则需要先解决
${{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}$
其中， $\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)$
为了实现最佳的效率，需要利用快速迭代的方法进行近似求解。利用高斯-赛德尔迭代（Gauss-Seidel）。
令 $u_{i,j}^{k+1}=G_{i,j}^{k}$
则，
$\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}$
因此，最终利用Split Bregman迭代进行各向异性TV降噪可以写作
初始化： ${{u}^{0}}=f,d_{x}^{0}=d_{y}^{0}=b_{x}^{0}=b_{y}^{0}=0$
While $\left\| {{u}^{k}}-{{u}^{k-1}} \right\|_{2}^{2}>threshold$
${{u}^{k+1}}={{G}^{k}}$
$d_{x}^{k+1}=shrink\left( {{\nabla }_{x}}{{u}^{k+1}}+b_{x}^{k},1/\lambda \right)$
$d_{y}^{k+1}=shrink\left( {{\nabla }_{y}}{{u}^{k+1}}+b_{y}^{k},1/\lambda \right)$
$b_{x}^{k+1}=b_{x}^{k}+\left( {{\nabla }_{x}}{{u}^{k+1}}-d_{x}^{k+1} \right)$
$b_{y}^{k+1}=b_{y}^{k}+\left( {{\nabla }_{y}}{{u}^{k+1}}-d_{y}^{k+1} \right)$
end

考虑各向同性情况（Isotropic case）

对于优化问题
$\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}$
同样地，将上述问题变为 ${{l}_{1}}$ 范数和 ${{l}_{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}$
其中 ${{\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}}}$
同理，在解决上述时需要提前解决
$\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}$
在这里引入一个广义的收缩公式
$d_{x}^{k+1}=\max \left( {{s}^{k}}-1/\lambda ,0 \right)\frac{{{\nabla }_{x}}{{u}^{k}}+b_{x}^{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}}}$
其中 ${{s}^{k}}=\sqrt{{{\left| {{\nabla }_{x}}{{u}^{k}}+b_{x}^{k} \right|}^{2}}+{{\left| {{\nabla }_{y}}{{u}^{k}}+b_{y}^{k} \right|}^{2}}}$
因此，最终利用Split Bregman迭代进行各向同性TV降噪可以写作
初始化： ${{u}^{0}}=f,d_{x}^{0}=d_{y}^{0}=b_{x}^{0}=b_{y}^{0}=0$
While $\left\| {{u}^{k}}-{{u}^{k-1}} \right\|_{2}^{2}>threshold$
${{u}^{k+1}}=G_{i,j}^{k}$
$d_{x}^{k+1}=\max \left( {{s}^{k}}-1/\lambda ,0 \right)\frac{{{\nabla }_{x}}{{u}^{k}}+b_{x}^{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}}}$
$b_{x}^{k+1}=b_{x}^{k}+\left( {{\nabla }_{x}}{{u}^{k+1}}-d_{x}^{k+1} \right)$
$b_{y}^{k+1}=b_{y}^{k}+\left( {{\nabla }_{y}}{{u}^{k+1}}-d_{y}^{k+1} \right)$
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