【图像去噪】基于matlab全变分算法图像去噪【含Matlab源码 419期】

二、图像去噪及滤波简介

1 图像去噪

1.1 图像噪声定义

噪声是干扰图像视觉效果的重要因素，图像去噪是指减少图像中噪声的过程。噪声分类有三种：加性噪声，乘性噪声和量化噪声。我们用f(x,y）表示图像，g(x,y）表示图像信号，n(x,y）表示噪声。

图像去噪是指减少数字图像中噪声的过程。现实中的数字图像在数字化和传输过程中常受到成像设备与外部环境噪声干扰等影响，称为含噪图像或噪声图像。去噪是图像处理研究中的一个重点内容。在图像的获取、传输、发送、接收、复制、输出等过程中，往往都会产生噪声，其中的椒盐噪声是比较常见的一种噪声，它属于加性噪声。

1.2 图像噪声来源

（1）图像获取过程中

图像传感器CCD和CMOS采集图像过程中受传感器材料属性、工作环境、电子元器件和电路结构等影响，会引入各种噪声。

（2）图像信号传输过程中

传输介质和记录设备等的不完善，数字图像在其传输记录过程中往往会受到多种噪声的污染。

1.3 噪声分类

噪声按照不同的分类标准可以有不同的分类形式：

基于产生原因：内部噪声，外部噪声。

基于噪声与信号的关系：

加性噪声：加性噪声和图像信号强度是不相关的，这类带有噪声的图像g可看成为理想无噪声图像f与噪声n之和：

g = f + n；

乘性嗓声：乘性噪声和图像信号是相关的，往往随图像信号的变化而变化，载送每一个象素信息的载体的变化而产生的噪声受信息本身调制。在某些情况下，如信号变化很小，噪声也不大。为了分析处理方便，常常将乘性噪声近似认为是加性噪声，而且总是假定信号和噪声是互相统计独立。

g = f + f*n

按照基于统计后的概率密度函数：

是比较重要的，主要因为引入数学模型这就有助于运用数学手段去除噪声。在不同场景下噪声的施加方式都不同，由于在外界的某种条件下，噪声下图像-原图像（没有噪声时）的概率密度函数（统计结果）服从某种分布函数，那么就把它归类为相应的噪声。下面将具体说明基于统计后的概率密度函数的噪声分类及其消除方式。

1.4 图像去噪算法的分类

（1）空间域滤波

空域滤波是在原图像上直接进行数据运算，对像素的灰度值进行处理。常见的空间域图像去噪算法有邻域平均法、中值滤波、低通滤波等。

（2）变换域滤波

图像变换域去噪方法是对图像进行某种变换，将图像从空间域转换到变换域，再对变换域中的变换系数进行处理，再进行反变换将图像从变换域转换到空间域来达到去除图像嗓声的目的。将图像从空间域转换到变换域的变换方法很多，如傅立叶变换、沃尔什-哈达玛变换、余弦变换、K-L变换以及小波变换等。而傅立叶变换和小波变换则是常见的用于图像去噪的变换方法。

（3）偏微分方程

偏微分方程是近年来兴起的一种图像处理方法，主要针对低层图像处理并取得了很好的效果。偏微分方程具有各向异性的特点，应用在图像去噪中，可以在去除噪声的同时，很好的保持边缘。偏微分方程的应用主要可以分为两类:一种是基本的迭代格式，通过随时间变化的更新，使得图像向所要得到的效果逐渐逼近，这种算法的代表为Perona和Malik的方程，以及对其改进后的后续工作。该方法在确定扩散系数时有很大的选择空间，在前向扩散的同时具有后向扩散的功能，所以，具有平滑图像和将边缘尖锐化的能力。偏微分方程在低噪声密度的图像处理中取得了较好的效果，但是在处理高噪声密度图像时去噪效果不好，而且处理时间明显高出许多。

（4）变分法

另一种利用数学进行图像去噪方法是基于变分法的思想，确定图像的能量函数，通过对能量函数的最小化工作，使得图像达到平滑状态，现在得到广泛应用的全变分TV模型就是这一类。这类方法的关键是找到合适的能量方程，保证演化的稳定性，获得理想的结果。

形态学噪声滤除器将开与闭结合可用来滤除噪声，首先对有噪声图像进行开运算，可选择结构要素矩阵比噪声尺寸大，因而开运算的结果是将背景噪声去除；再对前一步得到的图像进行闭运算，将图像上的噪声去掉。据此可知，此方法适用的图像类型是图像中的对象尺寸都比较大，且没有微小细节，对这类图像除噪效果会较好。

2 基于TV的图像去噪

2.1 TV图像去噪模型

全变分(TV) 图像去噪模型是由Rudin、Osher and Fatemi[4] 提出的， 并已成为图像去噪以及图像复原中最为成功的方法之一。TV图像去噪模型的成功之处就在于利用了自然图像内在的正则性，易于从噪声图像的解中反映真实图像的几何正则性，比如边界的平滑性15。

令f为原始的清晰图像，fo为被噪声污染的图像，即

式中n为具有零均值、方差为o²的随机噪声。Q表示图像的定义域，像素点(x，y)EΩ。通常有噪声图像的全变分比无噪声图像的全变分明显大，最小化全变分(TV)可以消除噪声，因此基于全变分的图像降噪可以归结为如下最小化问题：

满足约束条件:

最小化式 (1.2) 可以等价于最小化下式:

式中，第1项为数据保真项，它主要起保留原图像特性和降低图像失真度的作用；第2项为正则化项，参数为入规整参数，对平衡去噪与平滑起重要作用，它依赖于噪声水平。其导出的欧拉-拉格朗日方程为：

2.2 TV去噪的数值实现

则求解方程 (1.5) 的离散迭代格式为:

Step 1读入带噪图像f0;

Step 2初始化参数:

三、部分源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=imread('toys.bmp'); % load image
I=double(I(20:120,10:105)); % cut a piece, convert to double
%%% Parameters
std_n=10; var_n=std_n^2;  % Gaussian noise standard deviation
reduced_pw = 1.5*var_n;   % power to reduce in first phase
sig_w = 5; ws=4*sig_w+1;  % window size
%%%%%%%%%%%%%%

%%% Add noise
In = randn(size(I))*std_n;
I0 = I + In;  % noisy input image
% show original and noisy images
figure(1); imshow(uint8(I)); title('Original')
figure(2); imshow(uint8(I0)); title('Noisy image')
snr_noisy=db(I,I0)

% run normal tv - strong denoising
tic;
J=I0; 
ep_J=0.1; % minimum mean change in image J
J_old=0;
lam=0; iter=10; dt=0.2; ep=1; 

while (mean(mean(abs(J - J_old))) > ep_J),  % iterate until convergence
   J_old=J;
   J=tv(J,iter,dt,ep,lam,I0); % scalar lam
   lam = calc_lam(J,I0,reduced_pw); % update lambda (fidelity term)
end % while
% figure(3); imshow(uint8(J)); title('residue TV')
% snr_residue= db(I,J)

Ir=I0-J;  % Ir scalar
Pr = mean(mean(Ir.^2)); % power of residue
LV = loc_var(Ir,ws,sig_w^2);  % local variance (local  power of the residue )
% P=mean(mean(LV));
Pxy=1*(var_n^2)./LV;  %%  Sxy    inverse proportional to the LV

%%% Varying Lambda
lamxy=zeros(size(I0));
J=I0; J_old=0; 
ep_J=0.001;
%eps=0.01;

while (mean(mean(abs(J - J_old))) > ep_J),  % iterate until convergence
   J_old=J;
   J=tv(J,iter,dt,ep,lamxy,I0); % adaptive lam
   %J=tv(J,iter,dt,ep_J,lamxy,I0);
   lamxy = calc_lamxy(J,I0,Pxy,sig_w); % update lambda (fidelity term)
end % while

figure(3); imshow(uint8(J)); title('Adaptive TV')
snr_adap= db(I,J)
toc
t1=toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run scalar TV denoising for comparision 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
J=I0; 
% params
ep_J = 0.01; % minimum mean change in image J
lam=0; J_old=0;
i=0;
while (mean(mean(abs(J - J_old))) > ep_J),  % iterate until convergence
   J_old = J;
   J=tv(J,iter,dt,ep,lam,I0);     % scalar lam
   lam = calc_lam(J,I0,var_n,ep); % update lambda (fidelity term) 
end % for i
% Ir=I0-J;  % Ir scalar
% Pr = mean(mean(Ir.^2)); % power of residue

function Ig=gauss(I,ks,sigma2)
%private function: gauss (by Guy Gilboa):
% Ig=gauss(I,ks,sigma2)
% ks - kernel size (odd number)
% sigma2 - variance of Gaussian

[Ny,Nx]=size(I);
hks=(ks-1)/2;  % half kernel size
if (Ny<ks)   % 1d convolutin
   x=(-hks:hks); %x:1x ks
   flt=exp(-(x.^2)/(2*sigma2));  % 1D gaussian
   flt=flt/sum(sum(flt));  % normalize
   % expand
   x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx));%x0,xn :1 x hks; 
   eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)];   % ???
   Ig=conv(eI,flt);
   Ig=Ig(:,ks+hks+1:Nx+ks+hks);  % truncate tails of convolution   
else
   %% 2-d convolution
   x=ones(ks,1)*(-hks:hks); y=x'; 
   flt=exp(-(x.^2+y.^2)/(2*sigma2));  % 2D gaussian
   flt=flt/sum(sum(flt));  % normalize
   % expand
   if (hks>1)
      xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')';% xL,xR :Ny x 1
   else
      xL=I(:,1); xR=I(:,Nx);
   end
   eI=[xL*ones(1,hks) I xR*ones(1,hks)];   % Ny x  Nx+2hks
   if (hks>1)
      xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:));  % xU,xD: 1 x Nx+2hks
   else
      xU=eI(1,:); xD=eI(Ny,:);   
   end

四、运行结果