MATLAB实现遥感图像预处理——图像滤波

前言

这些代码均是使用最基础的方法，通过一步一步迭代过程来理解算法的原理及实现过程，并不采用于实用工程，读者以此作为学习参考即可。

观测图像

观测图像=真实图像+噪声：

常用的空域滤波方法

常用的频域滤波方法

频域滤波的步骤：

（1）中心变换，用(−1^(𝑥+𝑦)乘以输入图像
（2）计算图像的 DFT，即𝐹(𝑢, 𝑣)
（3）用滤波器函数𝐻(𝑢, 𝑣)乘以𝐹(𝑢, 𝑣)
（4）计算 3 中结果的反 DFT
（5）取 4 中结果的实部
（6）用(−1)^(𝑥+𝑦)乘以 (5) 中结果

频域滤波函数

运行结果

1、对灰度图像分别加入均值为 0，标准差为 0.1、0.2、0.5 的高斯噪声进行对比

2、对灰度图像分别加入噪声密度为0.1、0.2、0.5的椒盐噪声进行对比

结果分析：结合Figure2与Figure3，可以得出高斯噪声的标准差越大或者椒盐噪声的噪声密度越高，图像中的噪声成分越多，图像就越模糊。由于高斯噪声的标准差与椒盐噪声的噪声密度都是其噪声能量的体现方式，所以可以得出，噪声能量越大，图像中的噪声成分含量就越多，图像越模糊。

3、对均值为0，标准差为0.1的高斯噪声图像进行空域滤波，分别用大小为33、55、7*7的滑窗对图像进行处理，并对结果进行分析

结果分析：从图中三幅图来看，均值滤波能够有效的滤除噪声，而其滑窗的大小决定了处理后图像的平滑程度，滑窗越大，处理过后的图像就越平滑，其边缘部分就变得十分模糊。

4、对均值为0，标准差为0.1的高斯噪声图像进行频域滤波，分别用阈值大小为50、80、100、150的低通滤波器对图像进行处理，并对结果进行分析

结果分析：图像频域中高频部分对应图像中变化较快的部分，而低频对应图片中变化较慢的部分，从图中四幅图来看，频域滤波也能够有效的滤除噪声，而其阈值的大小影响处理后图像的平滑程度，阈值越大，处理过后的图像中变化较快的部分就越少，图像里边缘部分看起来就越清晰。

源代码

clc;
clearvars
close all;

Im1 = imread('ygtx.jpg');
Im2 = rgb2gray(Im1);
figure(1)
subplot(2,1,1),imshow(Im1),title('原图像');
subplot(2,1,2),imshow(Im2),title('原图像灰度图');
%高斯噪声
[width,height,z]=size(Im2);
av1=0;
std1=0.1;
u11=rand(width,height);
u21=rand(width,height);
x1=std1*sqrt(-2*log(u11)).*cos(2*pi*u21)+av1;
gs1=double(Im2)/255+x1;
gs1=uint8(255*gs1);
figure(2)
subplot(2,2,1),imshow(Im2);title('没有加入高斯噪声的图像');
subplot(2,2,2),imshow(gs1);title('加入均值为0，标准差为0.1的高斯噪声图像');

av2=0;
std2=0.2;
u12=rand(width,height);
u22=rand(width,height);
x2=std2*sqrt(-2*log(u12)).*cos(2*pi*u22)+av2;
gs2=double(Im2)/255+x2;
gs2=uint8(255*gs2);
subplot(2,2,3),imshow(gs2);title('加入均值为0，标准差为0.2的高斯噪声图像');

av3=0;
std3=0.5;
u13=rand(width,height);
u23=rand(width,height);
x3=std3*sqrt(-2*log(u13)).*cos(2*pi*u23)+av3;
gs3=double(Im2)/255+x3;
gs3=uint8(255*gs3);
subplot(2,2,4),imshow(gs3);title('加入均值为0，标准差为0.5的高斯噪声图像');
%椒盐噪声
figure(3)
subplot(2,2,1),imshow(Im2);title('没有加入椒盐噪声的图像');

jy1=Im2;
k11=0.1;
k21=0.3;
a11=rand(width,height)<k11;
a21=rand(width,height)<k21;
jy1(a11&a21)=0;
jy1(a11& ~a21)=255;
subplot(2,2,2),imshow(jy1),title('加入椒0.03盐0.07噪声图像');

jy2=imnoise(Im2,'salt & pepper',0.2);
subplot(2,2,3),imshow(jy2),title('加入噪声密度：0.2的椒盐噪声图像');
jy3=imnoise(Im2,'salt & pepper',0.5);
subplot(2,2,4),imshow(jy3),title('加入噪声密度：0.5的椒盐噪声图像');

%空域滤波
ky=double(gs1);
[x,y,z1]=size(ky);
%3*3窗口均值滤波
klb3=zeros(x,y);
for i= 2:x-1
            for j= 2:y-1
                s3=ky(i-1:i+1,j-1:j+1);
                klb3(i,j)=mean(s3(:));
            end
end
figure(4)
subplot(2,2,1),imshow(gs1);title('均值为0，标准差为0.1的高斯噪声图像');
subplot(2,2,2),imshow(klb3,[]),title('3*3窗口滤波结果');
%5*5窗口均值滤波
klb5=zeros(x,y);
for i= 3:x-2
             for j= 3:y-2
                  s5=ky(i-2:i+2,j-2:j+2);
                  klb5(i,j)=mean(s5(:));
             end
end
subplot(2,2,3),imshow(klb5,[]),title('5*5窗口滤波结果');
%7*7窗口均值滤波
klb7=zeros(x,y);
for i= 4:x-3
            for j= 4:y-3
                s7=ky(i-3:i+3,j-3:j+3);
                klb7(i,j)=mean(s7(:));
            end
end
subplot(2,2,4),imshow(klb7,[]),title('7*7窗口滤波结果');

%频域滤波
py=double(gs1);
figure(5)
subplot(3,2,1),imshow(gs1);title('均值为0，标准差为0.1的高斯噪声图像');
py_f=fftshift(fft2(py));
py_F=log10(abs(py_f)+1);
subplot(3,2,2),imshow(py_F,[]);title('傅里叶变换频域图像');
[m n]=size(py_f);
m_mid=fix(m/2);
n_mid=fix(n/2);
%阈值为50的低通滤波
d0=50;
plb1=zeros(m,n);
for i=1:m
            for j=1:n
                d=sqrt((i-m_mid)^2+(j-n_mid)^2);
                if d<=d0
                    h(i,j)=1;
                else
                    h(i,j)=0;
                end
                plb1(i,j)=h(i,j)*py_f(i,j);
             end
end
plb1=ifftshift(plb1);
plb1=uint8(real(ifft2(plb1)));
subplot(3,2,3),imshow(plb1);title('阈值50低通滤波图像');
%阈值为80的低通滤波
d0=80;
plb2=zeros(m,n);
for i=1:m
             for j=1:n
                 d=sqrt((i-m_mid)^2+(j-n_mid)^2);
                 if d<=d0
                    h(i,j)=1;
                 else
                    h(i,j)=0;
                 end
                 plb2(i,j)=h(i,j)*py_f(i,j);
             end
end
plb2=ifftshift(plb2);
plb2=uint8(real(ifft2(plb2)));
subplot(3,2,4),imshow(plb2);title('阈值80低通滤波图像');
%阈值为100的低通滤波
d0=100;
plb3=zeros(m,n);
for i=1:m
            for j=1:n
                d=sqrt((i-m_mid)^2+(j-n_mid)^2);
                if d<=d0
                    h(i,j)=1;
                else
                    h(i,j)=0;
                end
                plb3(i,j)=h(i,j)*py_f(i,j);
             end
end
plb3=ifftshift(plb3);
plb3=uint8(real(ifft2(plb3)));
subplot(3,2,5),imshow(plb3);title('阈值100低通滤波图像');
%阈值为150的低通滤波
d0=150;
plb4=zeros(m,n);
for i=1:m
            for j=1:n
                d=sqrt((i-m_mid)^2+(j-n_mid)^2);
                if d<=d0
                    h(i,j)=1;
                else
                    h(i,j)=0;
                end
                plb4(i,j)=h(i,j)*py_f(i,j);
            end
end
plb4=ifftshift(plb4);
plb4=uint8(real(ifft2(plb4)));
subplot(3,2,6),imshow(plb4);title('阈值150低通滤波图像');

---

最开始写的代码比较随意，图片也是随意找的，没有注释，给大家磕了。。