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

2024-03-21 17:40

本文主要是介绍Split Bregman迭代进行TV降噪介绍及MATLAB实现,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

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} uminxu+yu+2μuf22
为了能够使用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. umindx+dy+2μuf22s.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,umindx+dy+2μμf22+2λdxxu22+2λdyyu22
利用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,umindx+dy+2μμf22+2λdxxubxk22+2λdyyubyk22
为了解决上述问题的求解,则需要先解决
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μμf22+2λdxxubxk22+2λdyyubyk22
其中, ( μ 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(dxkbxk)+λyT(dykbyk)
为了实现最佳的效率,需要利用快速迭代的方法进行近似求解。利用高斯-赛德尔迭代(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+ui1,jk+ui,j+1k+ui,j1k+dx,i1,jkdx,i,jk+dy,i,j1kdy,i,jkbx,i1,jk+bx,i,jkby,i,j1k+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 ukuk122>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+1dxk+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+1dyk+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μuf22
同样地,将上述问题变为 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μμf22+2λdxxubx22+2λdyyuby22
其中 ∥ ( 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,jdx,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λdxxubx22+2λdyyuby22
在这里引入一个广义的收缩公式
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(sk1/λ,0)skxuk+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(sk1/λ,0)skyuk+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+bxk2+yuk+byk2
因此,最终利用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 ukuk122>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(sk1/λ,0)skxuk+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(sk1/λ,0)skyuk+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+1dxk+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+1dyk+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>toltx=dx-bx;ty=dy-by;up=u;%Update uu=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 dxtmpx=ux+bx;dx=sign(tmpx).*max(Z,abs(tmpx)-1/lambda);%Update dytmpy=uy+by;dy=sign(tmpy).*max(Z,abs(tmpy)-1/lambda);%Update bx and bybx=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 uu=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 dxdx=thresh.*tmpx;%Update dydy=thresh.*tmpy;%Update bx and bybx=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.

这篇关于Split Bregman迭代进行TV降噪介绍及MATLAB实现的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/833414

相关文章

python实现svg图片转换为png和gif

《python实现svg图片转换为png和gif》这篇文章主要为大家详细介绍了python如何实现将svg图片格式转换为png和gif,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录python实现svg图片转换为png和gifpython实现图片格式之间的相互转换延展:基于Py

Python利用ElementTree实现快速解析XML文件

《Python利用ElementTree实现快速解析XML文件》ElementTree是Python标准库的一部分,而且是Python标准库中用于解析和操作XML数据的模块,下面小编就来和大家详细讲讲... 目录一、XML文件解析到底有多重要二、ElementTree快速入门1. 加载XML的两种方式2.

Java的栈与队列实现代码解析

《Java的栈与队列实现代码解析》栈是常见的线性数据结构,栈的特点是以先进后出的形式,后进先出,先进后出,分为栈底和栈顶,栈应用于内存的分配,表达式求值,存储临时的数据和方法的调用等,本文给大家介绍J... 目录栈的概念(Stack)栈的实现代码队列(Queue)模拟实现队列(双链表实现)循环队列(循环数组

redis过期key的删除策略介绍

《redis过期key的删除策略介绍》:本文主要介绍redis过期key的删除策略,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录第一种策略:被动删除第二种策略:定期删除第三种策略:强制删除关于big key的清理UNLINK命令FLUSHALL/FLUSHDB命

C++如何通过Qt反射机制实现数据类序列化

《C++如何通过Qt反射机制实现数据类序列化》在C++工程中经常需要使用数据类,并对数据类进行存储、打印、调试等操作,所以本文就来聊聊C++如何通过Qt反射机制实现数据类序列化吧... 目录设计预期设计思路代码实现使用方法在 C++ 工程中经常需要使用数据类,并对数据类进行存储、打印、调试等操作。由于数据类

Python实现图片分割的多种方法总结

《Python实现图片分割的多种方法总结》图片分割是图像处理中的一个重要任务,它的目标是将图像划分为多个区域或者对象,本文为大家整理了一些常用的分割方法,大家可以根据需求自行选择... 目录1. 基于传统图像处理的分割方法(1) 使用固定阈值分割图片(2) 自适应阈值分割(3) 使用图像边缘检测分割(4)

Android实现在线预览office文档的示例详解

《Android实现在线预览office文档的示例详解》在移动端展示在线Office文档(如Word、Excel、PPT)是一项常见需求,这篇文章为大家重点介绍了两种方案的实现方法,希望对大家有一定的... 目录一、项目概述二、相关技术知识三、实现思路3.1 方案一:WebView + Office Onl

C# foreach 循环中获取索引的实现方式

《C#foreach循环中获取索引的实现方式》:本文主要介绍C#foreach循环中获取索引的实现方式,本文给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录一、手动维护索引变量二、LINQ Select + 元组解构三、扩展方法封装索引四、使用 for 循环替代

Spring Security+JWT如何实现前后端分离权限控制

《SpringSecurity+JWT如何实现前后端分离权限控制》本篇将手把手教你用SpringSecurity+JWT搭建一套完整的登录认证与权限控制体系,具有很好的参考价值,希望对大家... 目录Spring Security+JWT实现前后端分离权限控制实战一、为什么要用 JWT?二、JWT 基本结构

Java实现优雅日期处理的方案详解

《Java实现优雅日期处理的方案详解》在我们的日常工作中,需要经常处理各种格式,各种类似的的日期或者时间,下面我们就来看看如何使用java处理这样的日期问题吧,感兴趣的小伙伴可以跟随小编一起学习一下... 目录前言一、日期的坑1.1 日期格式化陷阱1.2 时区转换二、优雅方案的进阶之路2.1 线程安全重构2