数字图像处理第六次作业——图像噪声和恢复

---

一、基本概念及原理

1. 高斯噪声：

高斯噪声是指它的概率密度函数服从高斯分布的一类噪声。如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。高斯白噪声的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。

产生原因：

1. 图像传感器在拍摄时市场不够明亮、亮度不够均匀；

2. 电路各元器件自身噪声和相互影响；

3. 图像传感器长期工作，温度过高。

概率密度函数：
$$p(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-(z-\bar{z}{)}^2/2{\sigma }^2}$$

2. 椒盐噪声：

椒盐噪声又称脉冲噪声，它随机改变一些像素值，是由图像传感器，传输信道，解码处理等产生的黑白相间的亮暗点噪声。

产生原因：椒盐噪声往往由图像切割引起。

概率密度函数：
$$P(z)=\{\begin{array}{rlr}& P_a& z=a\\ & P_b& z=b\\ & 1-P_a-P_b& else\end{array}$$

椒盐噪声具体可以分为白盐噪声和胡椒噪声。白盐噪声又称白噪声，是在图像中添加一些随机的白色像素点(255)；胡椒噪声是在图像中添加一些随机的黑色像素点(0)；而通常说的椒盐噪声则是在图像中既有白色像素点，又有黑色像素点。

3. 维纳滤波器：

维纳滤波也称为最小均方误差滤波，它能处理被退化函数退化和噪声污染的图像。

该滤波方法建立在图像和噪声都是随机变量的基础之上，目标时找到未污染图像 $f(x,y)$ 的一个估计 $\hat{f}$ ，使它们之间的均方误差最小，即$e^2=E(f-\hat{f}{)}^2$，其中 E { ⋅ } E\{·\} E{⋅} 是参数期望值。

在假设噪声和图像不相关，其中一个或另一个有零均值，且估计中的灰度级是退化图像中灰度级的线性函数的条件下，均方误差函数的最小值在频率域由如下表达式给出：

$$\begin{array}{rl}\hat{F}(u,v)& =[\frac{H^{\ast }(u,v)S_f(u,v)}{S_f(u,v)|H(u,v){|}^2+S_{\eta }(u,v)}]G(u,v)\\ & =[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}]G(u,v)\\ & =[\frac{1}{H(u,v)}\frac{H(u,v{)}^2]}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}]G(u,v)\end{array}$$

其中，$H(u,v)$为退化函数，$H^{\ast }(u,v)$为$H(u,c)$的复共轭。

从上面的公式可以发现，如果没有噪声，即$S_{\eta }(u,v)=0$，此时维纳滤波变为直接逆滤波。
如果有噪声，那么$S_{\eta }(u,v)$如何估计将成问题，同时$S_f(u,v)$的估计也成问题。

在实际应用中假设退化函数已知，如果噪声为高斯白噪声，则$S_{\eta }(u,v)$为常数，但$S_f(u,v)$通常难以估计。一种近似的解决办法是用一个系数K代替$S_{\eta }(u,v)/S_f(u,v)$，因此上面的公式变为如下式所示的形式：

$$\hat{F}(u,v)=[\frac{1}{H(u,v)}\frac{|H(u,v){|}^2}{|H(u,v){|}^2+K}]G(u,v)$$

在实际应用中，根据处理的效果选取合适的K值。

---

二、实现过程和结果分析

第一题、第二题

1. 在测试图像上产生高斯噪声lena图-需能指定均值和方差；并用多种滤波器恢复图像，分析各自优缺点；
2. 在测试图像lena图加入椒盐噪声（椒和盐噪声密度均是0.1）；用学过的滤波器恢复图像；在使用反谐波分析Q大于0和小于0的作用；

1. 高斯噪声和椒盐噪声：

分析： 从对lena添加高斯噪声和椒盐噪声可以基本看出两种噪声对图像的影响。其中，高斯噪声的效果与均值、方差两个参数密切相关，方差越大时越模糊，颗粒感更大、更重，标准差越大时，图像越亮、越偏白。椒盐噪声的两个参数分别表示胡椒噪声和白盐噪声的密度。

MATLAB代码：

%加高斯噪声：
function Img_out=GaussianNoise(Img,av,std)
[M,N]=size(Img);
u1=rand(M,N);   u2=rand(M,N);
x=std*sqrt(-2*log(u1)).*cos(2*pi*u2)+av;
Img_out=uint8(255*(double(Img)/255+x));
end
%加椒盐噪声：
function Img_out=SaltPepperNoise(Img,a,b)
[M,N]=size(Img);
x=rand(M,N);
Img_out=Img;
Img_out(find(x<=a))=0;
Img_out(find(x>a&x<(a+b)))=255;
end

2. 用高斯滤波和中值滤波进行恢复：

分析： 通过中值滤波和高斯滤波对加了两种噪声的图像进行恢复的效果可以看出，两种滤波方式的恢复效果基本相当，高斯滤波略优一些，得到的图像稍微光滑一些。另外，中值滤波得到的图像亮度要稍高于高斯滤波得到的图像。

MATLAB代码：

% 高斯滤波：
function Img_out=GaussianFilter(Img,masksize,sigma)
for i=1:masksize
    for j=1:masksize
        x=i-ceil(masksize/2);
        y=j-ceil(masksize/2);
        h(i,j)=exp(-(x^2+y^2)/(2*sigma^2))/(2*pi*sigma^2);
    end
end
Img_out=uint8(conv2(Img,h,'same'));
end
% 中值滤波：
function Img_out=MedianFilter(Img,masksize)
exsize=floor(masksize/2);   %各方向扩展大小
Imgex=padarray(Img,[exsize,exsize],'replicate','both'); %扩展图片
[m,n]=size(Img);
Img_out=Img;    %将Img_out准备为和Img相同的size
for i=1:m
    for j=1:n
        neighbor=Imgex(i:i+masksize-1,j:j+masksize-1);  %截取邻域
        Img_out(i,j)=median(neighbor(:));   %中值滤波
    end
end
end

3. 反谐波分析及其中Q的作用：

　　分析： 通过上述处理可以看出，Q为正数时，对胡椒噪声有较好的恢复效果，而对白盐噪声则恢复效果很差，几乎将整幅图像处理成白色。当Q为辐射时，对白盐噪声有很好的恢复效果，而对胡椒噪声的恢复效果很差，几乎将整幅图像处理成黑色。

MATLAB代码：

%反谐波均值滤波：
function Img_out=Contraharmonic(Img,Q)
[M,N]=size(Img);
ImgSize=3;   ImgSize=(ImgSize-1)/2;
Img_out=Img;
for x=1+ImgSize:1:M-ImgSize
    for y=1+ImgSize:1:M-ImgSize
        is=Img(x-ImgSize:1:x+ImgSize,y-ImgSize:1:y+ImgSize);
        Img_out(x,y)=sum(double(is(:)).^(Q+1))/sum(double(is(:)).^(Q));
    end
end
end

---

第三题

推导维纳滤波器并实现下边要求；
a) 实现模糊滤波器如方程Eq. (5.6-11).
b) 模糊lena图像：45度方向，T=1；
c) 再模糊的lena图像中增加高斯噪声，均值= 0 ，方差=10 pixels 以产生模糊图像；
d) 分别利用方程 Eq. (5.8-6)和(5.9-4)，恢复图像；并分析算法的优缺点.

分析： 通过上述处理可以看出，用课本中所给的模糊滤波器函数处理得到的图像的模糊效果与运动模糊并不相同，它有着明显的方框感。而运动模糊得到的图像更加柔和一些，是将图像向指定方向均匀抹开的感觉。用课本上的维纳滤波器和约束最小二乘滤波对运动模糊加噪声后的图像进行恢复，可以看出，维纳滤波器恢复得到的颗粒感较大，但在一些细节处更清除，而约束最小二乘法得到的图像则较为平滑，但细节处也更加模糊。

MATLAB代码：

%模糊：
function Img_out=Blur(Img_in,a,b,T)
Img_in=double(imread('lena.bmp'));
Img_fft_shift=fftshift(fft2(Img_in));
[M,N]=size(Img_fft_shift);
for i=1:M
    for j=1:N
        H(i,j)=(T/(pi*(i*a+j*b)))*sin(pi*(i*a+j*b))*exp(-sqrt(-1)*pi*(i*a+j*b));
        G(i,j)=H(i,j)*Img_fft_shift(i,j);
    end
end
Img_out=ifft2(ifftshift(G));
Img_out=256.*Img_out./max(max(Img_out));
Img_out=uint8(real(Img_out));
end

 

---

参考资料

[1] 冈萨雷斯.数字图像处理(第三版)北京:电子工业出版社,2011

---

 
 
 

数字图像处理作业链接：
数字图像处理第一次作业——Bmp格式与基本变换
数字图像处理第二次作业——图像仿射变换
数字图像处理第三次作业——基于直方图的图像空域操作
数字图像处理第四次作业——空域滤波与边缘检测
数字图像处理第五次作业——频域滤波器
数字图像处理第六次作业——图像噪声和恢复
数字图像处理第七次作业——图像直线检测