数字图像处理学习笔记9:图像复原及重建1(常见噪声及滤波方法、噪声判别方法)


前言

图像增强处理是根据自己需要,按照主观要求去改善图像。图像复原则是一个客观的过程,将图像恢复到退化前的原图像。


一、图像退化/复原过程的模型

在这里插入图片描述


如上图所示,图像退化过程认为是对原图像 f ( x , y ) f(x,y) f(x,y)用一个退化函数进行处理后加上一个加性噪声,退化后的图像为 g ( x , y ) g(x,y) g(x,y)。图像复原就是对于退化后的图像 g ( x , y ) g(x,y) g(x,y),使用一个复原滤波器进行处理,使得到的图像接近于原图像 f ( x , y ) f(x,y) f(x,y)

空间域退化图像由下列公式表示
g ( x , y ) = h ( x , y ) ∗ f ( x , y ) + n ( x , y ) g(x,y)=h(x,y)*f(x,y)+n(x,y) g(x,y)=h(x,y)f(x,y)+n(x,y)
‘ ∗ ’ ‘*’ 表示卷积

频率域表示为
G ( x , y ) = G ( x , y ) F ( x , y ) + N ( x , y ) G(x,y)=G(x,y)F(x,y)+N(x,y) G(x,y)=G(x,y)F(x,y)+N(x,y)


二、常见空间域噪声模型

1.高斯噪声

高斯噪声的概率密度函数表示为:
p ( z ) = 1 2 π δ e − ( z − z ′ ) 2 / ( 2 δ ) p(z)=\frac{1}{\sqrt{2π}δ}e^{-(z-z^{'})^2/(2δ)} p(z)=2π δ1e(zz)2/(2δ)

其中 z z z表示灰度值, z ′ z^{'} z表示灰度均值, δ δ δ表示灰度的标准差。


高斯噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加高斯噪声的MATLAB代码如下:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2);                   
a=0;                        %均值
b=0.002;                    %方差
N_Gau=a+sqrt(b)*randn(M,N); %高斯噪声
im3=im2+N_Gau;
im3=im2uint8(im3);
figure
imshow(im3) 

结果如下:
在这里插入图片描述


2.瑞利噪声

瑞利噪声的概率密度函数表示为:
p ( z ) = { 2 b ( z − a ) e − ( z − a ) 2 / b , z>=a 0 , z<a p(z) = \begin{cases} \frac{2}{b}(z-a)e^{-(z-a)^2/b}, & \text{z>=a} \\ 0, & \text{z<a} \end{cases} p(z)={b2(za)e(za)2/b,0,z>=az<a
其中:
z ′ = a + π b 4 z^{'}=a+\sqrt{\frac{πb}{4}} z=a+4πb
δ 2 = b ( 4 − π ) 4 δ^2=\frac{b(4-π)}{4} δ2=4b(4π)


瑞利噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加瑞利噪声的MATLAB代码为:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2); 
a=0;
b=0.06;
B=0.5;
N_Ray1=a+b*raylrnd(B,M,N); %瑞利噪声
im3=im2+N_Ray1;
im3=im2uint8(im3);
figure
imshow(im3)

结果:

在这里插入图片描述


3.伽马噪声

伽马噪声的概率密度函数表示为:
p ( z ) = { a 2 z b − 1 ( b − 1 ) ! e − a z , z>=a 0 , z<a p(z) = \begin{cases} \frac{a^2z^{b-1}}{(b-1)!}e^{-az}, & \text{z>=a} \\ 0, & \text{z<a} \end{cases} p(z)={(b1)!a2zb1eaz,0,z>=az<a
其中:
z ′ = − b a z^{'}=-\frac{b}{a} z=ab
δ 2 = b a 2 δ^2=\frac{b}{a^2} δ2=a2b


伽马噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加伽马噪声的MATLAB代码:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2); 
a=0;
b=0.03;
A=1;
B=1;
N_Gam=a+b*gamrnd(A,B,[M,N]);  %伽马噪声
im3=im2+N_Gam;
im3=im2uint8(im3);
figure
imshow(im3)

结果:
在这里插入图片描述


4.指数噪声

指数噪声的概率密度函数表示为:
p ( z ) = { a e − a 2 , z>=0 0 , z<0 p(z) = \begin{cases} ae^{-a^2}, & \text{z>=0} \\ 0, & \text{z<0} \end{cases} p(z)={aea2,0,z>=0z<0
其中:
z ′ = − 1 a z^{'}=-\frac{1}{a} z=a1
δ 2 = 1 a 2 δ^2=\frac{1}{a^2} δ2=a21


指数噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加伽马噪声的MATLAB代码:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2); 
a=0;
b=0.03;
A=1;
B=1;
N_Gam=a+b*gamrnd(A,B,[M,N]);  %伽马噪声
im3=im2+N_Gam;
im3=im2uint8(im3);
figure
imshow(im3)

结果:
在这里插入图片描述


5.均匀分布噪声

均匀分布噪声的概率密度函数表示为:
p ( z ) = { a e − a 2 , a=<z=<b 0 , 其他 p(z) = \begin{cases} ae^{-a^2}, & \text{a=<z=<b} \\ 0, & \text{其他} \end{cases} p(z)={aea2,0,a=<z=<b其他
其中:
z ′ = − 1 b − a z^{'}=-\frac{1}{b-a} z=ba1
δ 2 = ( b − a ) 2 12 δ^2=\frac{(b-a)^2}{12} δ2=12(ba)2


均匀分布噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加均匀分布噪声的MATLAB代码:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2); 
a=0;
b=0.05;
A=0;
B=2;
N_unif=a+b*unifrnd(A,B,[M,N]);  %均匀分布噪声
im3=im2+N_unif;
im3=im2uint8(im3);
figure
imshow(im3)

结果:
在这里插入图片描述


6.脉冲(椒盐)噪声

椒盐噪声的概率密度函数表示为:
p ( z ) = { p ( a ) , z=a p ( b ) , z=b 1 − p ( a ) − p ( b ) , 其他 p(z) = \begin{cases} p(a), & \text{z=a} \\ p(b), & \text{z=b}\\ 1-p(a)-p(b), & \text{其他} \end{cases} p(z)=p(a),p(b),1p(a)p(b),z=az=b其他


椒盐噪声的灰度分布表示曲线如下:
在这里插入图片描述


在图像中添加椒盐噪声的MATLAB代码:

im1=imread('1.jpg');
im2=rgb2gray(im1);
figure;
imshow(im2);
[M,N] = size(im2);                
im2 = im2double(im2); 
a=0.05;
im3 = imnoise(im2,'salt & pepper',a); %椒盐噪声
im3=im2uint8(im3);
figure
imshow(im3)


结果:
在这里插入图片描述


三、图像中噪声判别

对于上述六种噪声,椒盐噪声与其他噪声图像差别较大,能够区分,其他噪声图像则不能够凭人眼区分。其他噪声区分方法为取原图像中灰度较为平滑的区域,然后根据灰度分布图来判断,如:
在这里插入图片描述
根据部分区域的灰度分布图结合前面各种噪声的灰度分布可以初步估计图中噪声的类型。


四、空间滤波去噪

在这里插入图片描述


1.算数均值滤波器及MATLAB代码

算数均值滤波器即用周围灰度值的平均值代替滤波中心像素的灰度值:
g ( x , y ) = 1 m n ∑ ( s , t ) ∈ S ( x , y ) f ( s , t ) g(x,y)=\frac{1}{mn}\sum_{(s,t)\in S(x,y)}f(s,t) g(x,y)=mn1(s,t)S(x,y)f(s,t)
其中m,n为滤波模板大小,f为需要处理的图形,g为滤波后的图像


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=im2(i-r:i+r,j-r:j+r);
        im3(i,j)=sum(sum(h))/(d^2);       %算数均值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

对添加高斯噪声的图像进行算数均值滤波,结果如下:
在这里插入图片描述


2.几何均值滤波器及MATLAB代码

几何滤波器与算数均值滤波器相比,丢失的图像细节较少,公式如下:
g ( x , y ) = [ ∏ ( s , t ) ∈ S ( x , y ) f ( s , t ) ] 1 m n g(x,y)=[\prod_{(s,t)\in S(x,y)}f(s,t)]^{\frac{1}{mn}} g(x,y)=[(s,t)S(x,y)f(s,t)]mn1


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        h(h==0)=1;                %防止出现×0的情况
        im3(i,j)=prod(prod(h))^(1/d^2);   %几何均值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

对高斯噪声滤波结果如下:
在这里插入图片描述


3.谐波均值滤波器及MATLAB代码

谐波均值滤波器对处理高斯噪声较好,滤波器公式如下:
g ( x , y ) = m n ∑ ( s , t ) ∈ S ( x , y ) 1 f ( s , t ) g(x,y)=\frac{mn}{\sum_{(s,t)\in S(x,y)}\frac{1}{f(s,t)}} g(x,y)=(s,t)S(x,y)f(s,t)1mn


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        h(h==0)=1;                %防止出现×0的情况
        im3(i,j)=d^2/(sum(sum(1./h)));   %谐波均值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

对高斯噪声图像滤波结果如下:
在这里插入图片描述


4.逆谐波均值滤波器及MATLAB代码

逆谐波均值滤波器公式为:
g ( x , y ) = ∑ ( s , t ) ∈ S ( x , y ) f ( s , t ) Q + 1 ∑ ( s , t ) ∈ S ( x , y ) f ( s , t ) Q g(x,y)=\frac{\sum_{(s,t)\in S(x,y)}f(s,t)^{Q+1}}{\sum_{(s,t)\in S(x,y)}f(s,t)^Q} g(x,y)=(s,t)S(x,y)f(s,t)Q(s,t)S(x,y)f(s,t)Q+1

当Q=0时,变为算数均值滤波器。当Q=-1时,变为谐波均值滤波器。当Q为正时,可以去除胡椒噪声。当Q为负时,可以去除盐粒噪声。


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
Q=3;                                      %逆谐波滤波器阶数
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        h(h==0)=1;                %防止出现×0的情况
        im3(i,j)=(sum(sum(h.^(Q+1))))/(sum(sum(h.^Q)));   %逆谐波均值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

对椒盐噪声图像进行逆谐波均值滤波,结果如下:
在这里插入图片描述如上图所示,处理结果并不是很好,特别是Q为负时,暂不清楚原因。


5.中值滤波器及MATLAB代码

中值滤波器即将邻域像素灰度的中值作为中心点的灰度值


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        im3(i,j)=median(median(h));        %中值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

使用中值滤波器处理椒盐噪声,结果如下:

在这里插入图片描述


6.最大值、最小值滤波器及MATLAB代码

即以滤波窗口中排序最后一位或第一位代替窗口中心灰度值


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        im3(i,j)=max(max(h));        %最大值值滤波
        % im3(i,j)=min(min(h));      %最小值值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

7.中点滤波器及MATLAB代码

及以最大值与最小值中点的值为中心像素灰度值


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        im3(i,j)=1/2*(max(max(h))+min(min(h))) ; %中点滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

以高斯噪声为例:
在这里插入图片描述


8.修正的阿尔法均值滤波器及MATLAB代码

去掉部分灰度最高的点以及灰度值最低的点后,以剩下的像素灰度平均值作为结果。


MATLAB程序如下(令滤波模板为n*n大小,n为奇数):

im=imread('1.jpg');
im1=rgb2gray(im);
figure
imshow(im1)

d=3;                                      %滤波模板大小,d=3,5,7......
r=floor(d/2);                             %图像需要填充的大小
im2 = padarray(im1,[r,r],'symmetric');    %对原图像进行扩充,使得图像边界像素能够得到处理
[m,n]=size(im2);   
im3=zeros(size(m,n));
d2=3;                                      %去除最大最小值得个数
for i=r+1:m-r
    for j=r+1:n-r
        h=double(im2(i-r:i+r,j-r:j+r));
        h1=reshape(h,[1,d^2]);         %将h矩阵中的数重新排列为一维矩阵
        h2=sort(h1,'ascend');          %对矩阵h1进行升序排列,即从小到大排列
        h3=h2(1,d2+1:d^2-d2);          %去除部分最大最小值后的矩阵
        im3(i,j)=sum(sum(h3))/(d^2-2*d2) ;       %修正阿尔法均值滤波
    end
end
im3=im3(r+1:m-r,r+1:n-r);
im3=uint8(im3);
figure
imshow(im3)

以高斯噪声为例:
在这里插入图片描述


9.自适应局部降低噪声滤波器及MATLAB代码

  • 13
    点赞
  • 59
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

‭刘燚

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值