matlab数字图像处理：图像复原及重建--常见的噪声模型以及滤波方法_图像复原,自适应局部噪声消除滤波器,matlab

---

前言

数字图像处理c++ opencv（VS2019 opencv4.53）持续更新

---

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

---

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

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

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

---

二、常见空间域噪声模型

1.高斯噪声

高斯噪声的概率密度函数表示为：
$$p(z)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-(z-z^{{}^{\prime }}{)}^2/(2\delta )}$$

其中$z$表示灰度值，$z^{{}^{\prime }}$表示灰度均值，$\delta$表示灰度的标准差。

---

高斯噪声的灰度分布表示曲线如下：

---

在图像中添加高斯噪声的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) 
1
2
3
4
5
6
7
8
9
10
11
12
13

---

结果如下：

---

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​(z−a)e−(z−a)2/b,0,​z>=az<a​
其中：
$$z^{{}^{\prime }}=a+\sqrt{\frac{\pi b}{4}}$$
$${\delta }^2=\frac{b(4-\pi )}{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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

---

结果：

---

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)={(b−1)!a2zb−1​e−az,0,​z>=az<a​
其中：
$$z^{{}^{\prime }}=-\frac{b}{a}$$
$${\delta }^2=\frac{b}{a^2}$$

---

伽马噪声的灰度分布表示曲线如下：

---

在图像中添加伽马噪声的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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

---

结果：

---

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)={ae−a2,0,​z>=0z<0​
其中：
$$z^{{}^{\prime }}=-\frac{1}{a}$$
$${\delta }^2=\frac{1}{a^2}$$

---

指数噪声的灰度分布表示曲线如下：

---

在图像中添加伽马噪声的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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

---

结果：

---

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)={ae−a2,0,​a=<z=<b其他​
其中：
$$z^{{}^{\prime }}=-\frac{1}{b-a}$$
$${\delta }^2=\frac{(b-a{)}^2}{12}$$

---

均匀分布噪声的灰度分布表示曲线如下：

---

在图像中添加均匀分布噪声的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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

---

结果：

---

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),1−p(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
2
3
4
5
6
7
8
9
10
11
12

---

结果：

---

三、图像中噪声判别

对于上述六种噪声，椒盐噪声与其他噪声图像差别较大，能够区分，其他噪声图像则不能够凭人眼区分。其他噪声区分方法为取原图像中灰度较为平滑的区域，然后根据灰度分布图来判断，如：

根据部分区域的灰度分布图结合前面各种噪声的灰度分布可以初步估计图中噪声的类型。

---

四、空间滤波去噪

---

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

算数均值滤波器即用周围灰度值的平均值代替滤波中心像素的灰度值：
$$g(x,y)=\frac{1}{mn}\sum _{(s,t)\in 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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

---

对添加高斯噪声的图像进行算数均值滤波，结果如下：

---

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

几何滤波器与算数均值滤波器相比，丢失的图像细节较少，公式如下：
$$g(x,y)=[\prod _{(s,t)\in S(x,y)}f(s,t){]}^{\frac{1}{mn}}$$

---

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

对高斯噪声滤波结果如下：

---

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

谐波均值滤波器对处理高斯噪声较好，滤波器公式如下：
$$g(x,y)=\frac{mn}{{\sum }_{(s,t)\in S(x,y)}\frac{1}{f(s,t)}}$$

---

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

对高斯噪声图像滤波结果如下：

---

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

逆谐波均值滤波器公式为：
$$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}$$

当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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

---

对椒盐噪声图像进行逆谐波均值滤波，结果如下：
如上图所示，处理结果并不是很好，特别是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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

---

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

---

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

---

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

---

以高斯噪声为例：

---

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

---

以高斯噪声为例：

---

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