图像去噪基于matlab全变分算法（TV）图像去噪含Matlab源码 625期

Posted 2021-06-24 紫极神光

一、简介

全变分（Total variation），也称为全变差，是图象复原中常用的一个名词。本文简要介绍全变分的概念以及在图象去噪中的应用。

1 一维信号的全变分和去噪

1.1 一维连续函数的全变分
一维连续实函数f(x)f(x)在区间[a,b]⊂R[a,b]⊂R上的全变分定义为参数曲线x→f(x),x∈[a,b]x→f(x),x∈[a,b]的弧长。其表达式为
Vba(f)=∫ba|f′(x)|dx
Vab(f)=∫ab|f′(x)|dx
说白了，所谓的“变分”就是|f(x+Δx)−f(x)||f(x+Δx)−f(x)|，对于连续函数Δx→0Δx→0。而全变分是对函数定义的区间而言的，就是将“变分”在区间上累加起来。
一维离散信号的全变分
从上面连续实函数的全变分，我们可以很容易想到它的离散形式。给出信号序列{yi},i=1,..,n{yi},i=1,..,n，它的全变分定义为
V(y)=∑i=1n|yi+1−yi|
V(y)=∑i=1n|yi+1−yi|
用一句话来概括，全变分是前后项之差的绝对值之和。

1.2 一维信号去噪

当我们得到观察信号xixi，希望xixi变得平滑，也就是对xixi去噪。一种很直观的想法就是让信号的全变分变小。全变分对应的物理意义就是输入信号的平滑度。设得到的恢复信号为yiyi，它应该满足两个条件：
yiyi与观察信号xixi的差距不大。这个差距的常用数学表达式就是
E(x,y)=12∑i(xi−yi)2
E(x,y)=12∑i(xi−yi)2
yiyi的全变分不大。
将物理约束转化为数学模型，求解yy等价于求解下面这个优化问题：

minyE(x,y)+λV(y)
minyE(x,y)+λV(y)
其中参数λλ是正常数，用于调节两个约束的作用大小。注意到E(x,y)E(x,y)和V(y)V(y)都是凸函数，这是一个无约束凸优化问题，有很多经典方法可以求解。

2 二维离散信号（图象）的全变分和去噪

图象是典型的二维离散信号，Rudin在1992年将其全变分定义为

V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2−−−−−−−−−−−−−−−−−−−−−−−√
V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2
这个函数是各项同性的，但是不可微，也并不是凸函数。非凸函数的优化求解难度、速度和稳定性都无法与凸函数相比。因此二维全变分有另一种常用定义
V(y)=∑i,j|yi+1,j−yi,j|2−−−−−−−−−−−√+|yi,j+1−yi,j|2−−−−−−−−−−−√=∑i,j|yi+1,j−yi,j|+|yi,j+1−yi,j|
V(y)=∑i,j|yi+1,j−yi,j|2+|yi,j+1−yi,j|2=∑i,j|yi+1,j−yi,j|+|yi,j+1−yi,j|
这个函数是凸函数了。

仿照一维信号的去噪，基于全变分的图象去噪可以看成求解优化问题

minyE(x,y)+λV(y)
minyE(x,y)+λV(y)
其中E(x,y)E(x,y)作为数据误差项定义为

E(x,y)=12∑i,j(xi,j−yi,j)2
E(x,y)=12∑i,j(xi,j−yi,j)2
当VV有凸函数形式时，问题变为无约束凸优化问题，从而容易求解。

二、源代码

function image=bldconv(g)


a1=0.1;
a2=0.01;

PQ=paddedsize(size(g));
G=fft2(g,PQ(1),PQ(2));

[y x]=size(G);
htemp=ones(3);
h0=freqz2(htemp,PQ(1),PQ(2));
R=Rcreat(y,x);

H=h0;
for k=1:10 %计算IMG和psf
iMG=(conj(H).*G)./((conj(H).*H)+a1*R);
H=conj(iMG).*G./((conj(iMG).*iMG)+a2*R);
end

IMG=mat2gray(real(ifft2(iMG)));
image=IMG(1:size(g,1),1:size(g,2));
imwrite(image,\'testfile.tif\');
imshow(image,[]);
hold on;

end

%===========================
%计算R矩阵的函数Rcreat
%===========================

function R=Rcreat(y,x)
%R矩阵生成子函数
%by Realasking
%为bdeconv.m编制

i=1:y;
j=1:x;

RI=zeros(y,x);
RJ=zeros(y,x);
R=zeros(y,x);

for k=1:y %向量化代码生成R的矩阵
RI(k,i)=-2*cos(2*pi*i./y);
end

for k=1:x
RJ(j,k)=-2*cos(2*pi*j\'./x);
end
Img=imread(\'Baboon1.bmp\'); %读取图片
PSF=fspecial(\'motion\',3); %创建PSF
gb=imfilter(Img,PSF,\'circular\'); %创建退化图像
Img=imnoise(gb,\'gaussian\',0,0.01); %加噪声
figure,imshow(Img)
Img=double(Img);
Img0=Img;
PQ=paddedsize(size(Img));
IMG=fft2(Img,PQ(1),PQ(2));
IMG0=fft2(Img0,PQ(1),PQ(2));
[nrow,ncol]=size(Img);  % 获取图像尺寸大小
lamda1=0.02;
lamda2=0.02;
dt=0.28;   % 0.25-0.35为最佳
G=gauss(Img,7,3);
Ix = 0.5*(G(:,[2:ncol,ncol])-G(:,[1,1:ncol-1])); % x方向梯度
Iy = 0.5*(G([2:nrow,nrow],:)-G([1,1:nrow-1],:)); % y方向梯度
gradG = Ix.^2+Iy.^2;  % 梯度大小
P=1+1./(1+gradG);  %自适应滤波器
deltax=Img0; %zeros(nrow,ncol);  %产生与图像大小相同的矩阵
deltay=Img0; %zeros(nrow,ncol);
htemp=ones(3);
h0=freqz2(htemp,PQ(1),PQ(2));
R=Rcreat(PQ(1),PQ(2));
H=h0;



for M=1:10 % 设置迭代次数

for i=2:(nrow-1)
for j=2:(ncol-1)

deltax(i,j)=(Img(i+1,j)-Img(i,j))./(((Img(i+1,j)-Img(i,j)).^2+(Img(i,j+1)-Img(i,j-1)).^2/4+1).^(1-0.5.*P(i,j)))-(Img(i,j)-Img(i-1,j))./(((Img(i,j)-Img(i-1,j)).^2+(Img(i-1,j+1)-Img(i-1,j-1)).^2/4+1).^(1-0.5.*P(i-1,j)));
deltay(i,j)=(Img(i,j+1)-Img(i,j))./(((Img(i+1,j)-Img(i-1,j)).^2/4+(Img(i,j+1)-Img(i,j)).^2+1).^(1-0.5.*P(i,j)))-(Img(i,j)-Img(i,j-1))./(((Img(i+1,j-1)-Img(i-1,j-1)).^2/4+(Img(i,j)-Img(i,j-1)).^2+1).^(1-0.5.*P(i,j-1)));
end
end

div=deltax+deltay;
DIV=fft2(div,PQ(1),PQ(2));
IMG=IMG+dt*(-conj(H).*H.*IMG+conj(H).*IMG0+lamda1.*DIV);
H=conj(IMG).*IMG0./(conj(IMG).*IMG+lamda2.*R);

Img1=real(ifft2(IMG));
Img=Img1(1:size(Img0,1),1:size(Img0,2));
end

三、运行结果

四、备注

版本：2014a

完整代码或代写加1564658423