【图像去噪】基于matlab小波域双重局部维娜滤板图像去噪【含Matlab源码 1642期】

一、图像去噪及滤波简介

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 基于邻域的小波阈值
(1) Bayes Shrink阈值
假设信号为{fj，i，j=1，…，N}，N是2的幂，它被一加性随机噪声污染，记为：

其中aj是独立同分布(iid independence identical distribution) 的均值为零的正态分布，并且独立于fy。去噪的目标是得到fy的估计值f^，使均方误差(MSE Mean Square Error) 最小。

令g={9mhj， f={f hj， e={a jhj，又令Y=Wg， X=Wf， V=We，其中W表示二进正交小波算子。于是可得：Y=X+V(3)
文献[1] 中，在贝叶斯估计理论框架下，假设小波系数服从广义高斯分布(GGD) ，即：

其中，

二、部分源代码


clear;close all;%---------Load the original image and show it------------------------------
im_original=imread('images\lena.png','png');% load the original image
im_original=double(im_original);            % convert it to double
figure(1);imshow(im_original,[]);title('The original image');% show it%--------Generate a noisy image, the noise is assumed to be additive,------
%--------white,Gaussian and independent of the original image--------------
sigma=25; % Standard deviation of the noise
im_noisy=im_original+sigma*randn(size(im_original));% Generate noisy image%--------Show the noisy image----------------------------------------------
figure(2);imshow(im_noisy,[]);title('The noisy image');%--------Set the parameters of the denoising algorithm---------------------
% Attention: the parameters is set for reproducing the result of our IEEE
%            Signal Processing Letters paper, do not change them unless you
%            are an advanced user with a deep understanding of the theory.
wav_base1='db4'; % the wavelet base used in the first LWFDW
wav_base2='sym8';% the wavelet base used in the second LWFDW
level=5;         % the wavelet decomposition level
pad_mode=0;% 0:stands for the period extention mode% 1:stands for the symmetric extention modesymbol=0; % 0: stands for denoising using Decimated wavelet transform% 1: stands for denoising using Undecimated wavelet transform% If symbol=0, uncomment the following 4 code lines for reproducing the
% results of denoising image using maximal decimated wavelet transform
r1=[5,4,4,3,3];a1=2;  % the parameters of the directional window used in
%                       the first LWFDW under Decimated wavelet transform
r2=[3,2,2,1,1];a2=1.5;% the parameters of the directional window used in
%                       the second LWFDW% If symbol=1, uncomment the following 4 code lines for reproducing the
% results of denoising image using Undecimated wavelet transform
% r1=[7,8,8,9,9];a1=2;  % the parameters of the directional window used in
%                       % first LWFDW under Undecimated wavelet transform
% r2=[3,4,4,5,5];a2=2;% the parameters of the directional window used in
%                       % the second LWFDW%--------Call the denoising function to denoise the noisy image------------
tic;
im_denoised=denoise_DLWFDW(im_noisy,wav_base1,wav_base2,level,r1,a1...,r2,a2,pad_mode,symbol);
toc;%--------Compute the PSNR and show the result------------------------------
PSNR=10*log10(255^2/mse(im_original-im_denoised))
figure(3);imshow(im_denoised,[]);title('The denoised image using DLWFDW');function im_denoised=DLWFDW_DWT(im_noisy,wav_base1,wav_base2,level,...r1,a1,r2,a2);
% DLWFDW using the maximal decimated wavelet transform
% written by Peng-lang Shui and Zuo-feng Zhou
% last revision:16/8/2007%-------maximal decimated wavelet decomposition using the DWT2-------------
t=im_noisy;index=1;
% estimate the standard deviation of the noise using the MAD estimator
[tt,ht,vt,dt]=dwt2(im_noisy,wav_base1,'mode','per');
sigma=median(median(abs(dt)))/.6745;% wavelet decomposition and the first local Wiener filtering
% with directional windows
for i=1:level[t,h,v,d]=dwt2(t,wav_base1,'mode','per');% for horizontal subbandW=elliptic_window(r1(i),1/a1,1);% generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,h.^2-sigma^2)/M);% the signal varianceh=(var./(var+sigma^2)).*h;      % local Wiener filtering% for vertical subbandW=elliptic_window(r1(i),a1,1);% generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,v.^2-sigma^2)/M);% the signal variancev=(var./(var+sigma^2)).*v;      % local Wiener filtering% for diagonal subbandW=elliptic_window(r1(i),a1,2);% generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,d.^2-sigma^2)/M);% the signal varianced=(var./(var+sigma^2)).*d;      % local Wiener filteringcoeff{index}=h;coeff{index+1}=v;coeff{index+2}=d;index=index+3;
end% wavelet reconstruction to get the pilot image
for i=1:levelt=idwt2(t,coeff{(level-i+1)*3-2},coeff{(level-i+1)*3-1},...coeff{(level-i+1)*3},wav_base1,'mode','per');
end
im_pilot=t;% get the pilot image% The second local Wiener filtering with directional windows
t_p=im_pilot;t_n=im_noisy;index=1;
for i=1:level[t_p,h_p,v_p,d_p]=dwt2(t_p,wav_base2,'mode','per');[t_n,h_n,v_n,d_n]=dwt2(t_n,wav_base2,'mode','per');% for horizontal subbandW=elliptic_window(r2(i),1/a2,1);% generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,h_p.^2)/M); % the signal varianceh_n=(var./(var+sigma^2)).*h_n;  % local Wiener filtering% for vertical subbandW=elliptic_window(r2(i),a2,1);  % generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,v_p.^2)/M); % the signal variancev_n=(var./(var+sigma^2)).*v_n;  % local Wiener filtering% for diagonal subbandW=elliptic_window(r2(i),a2,2);  % generate the directional windowM=sum(sum(W));                  % the size of the windowvar=max(0,filter2(W,d_p.^2)/M); % the signal varianced_n=(var./(var+sigma^2)).*d_n;      % local Wiener filteringcoeff{index}=h_n;coeff{index+1}=v_n;coeff{index+2}=d_n;index=index+3;
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.