计算机视觉实验二：图像滤波

1.滤波器实验

 1.1 滤波器实验目的与要求

（1）了解线性滤波器在数字图像处理中的应用。

（2）理解图像滤波器的处理过程和处理前后图像的变化关系。

1.2 滤波器原理

   线性低通滤波器是最常用的线性平滑滤波器。这种滤波器的所有系数都是正的。对3*3的模板来说，最简单的是取所有系数都为1。为保证输出图像仍在原来的灰度范围，在计算出算子后要将其除以9再进行赋值。这种方法称为邻域平均法。

  维纳滤波也是经典的线性降噪滤波器。维纳滤波是一种自适应滤波，它能根据图像的局部方差调整滤波器的输出。局部方差越大，滤波器的平滑作用越强。它的最终目标是使恢复图像f*(与原图像f(x,y)的均方误差     最小。

中值滤波器是最常用的非线性平滑滤波器。它是一种邻域运算，类似于卷积，但计算的不是加权求和，而是把邻域中的像素按灰度级进行排序，然后选择该组的中间值作为输出像素值。

具体步骤是：

（1）将模板在图像中漫游，并将模板中心与图像中某个像素的位置重合；

（2）读取模板下各对应像素的灰度值；

（3）将这些灰度值从小到大排成一列，找出这些值里排在中间的一个；

（4）将这个中间值赋给对应模板中心位置的像素。

   由此可以看出，中值滤波器的主要功能就是让与周围像素灰度值的差比较大的像素改取与周围的像素值接近的值，从而可以消除孤立的噪声点。中值滤波器不像均值滤波器那样，它在衰减噪声的同时不会使图像的边界模糊，这也是中值滤波受欢迎的主要原因。中值滤波器去噪声的效果依赖于两个要素：邻域的空间范围和中值计算中所涉及的像素数。一般来说，小于中值滤波器面积一半的亮或暗的物体基本上会被滤掉，而较大的物体则几乎会原封不动地保存下来。因此，中值滤波器的空间尺寸必须根据现有的问题来进行调整。较简单的模板是 的方形（这里N常是奇数），计算使用到的所有 各像素点。另外，我们也可以使用稀疏分布的模板来节省时间。对于不同的模板，中值也要依情况而定。计算中值所使用的像素的增加跟去噪时的效果之间的关系是非线性的。如果现在的问题需要用到大尺度的中值滤波器，那么用稀疏分布的中值模板或许能得到令人满意的效果。

1.3 实验内容和步骤

（1）对参考图像实现均值滤波器滤波操作，需要添加相关噪声。可参考程序，也可自行编写，程序参考代码如下：（实现线性平滑均值过滤器）

I=imread('lena.png');
J=imnoise(I,'salt',0.02);
figure ,subplot(2,2,1),imshow(I),title('原图像');
subplot(2,2,2),imshow(J),title('加入椒盐噪声的图像');
K1=filter2(fspecial('average',3),J)/255;
%的均值滤波处理结果
K2=filter2(fspecial('average',5),J)/255;
%的均值滤波处理结果
subplot(2,2,3),imshow(K1),
title('的均值滤波处理结果');
subplot(2,2,4),imshow(K2),
title('的均值滤波处理结果');

（2）使用高斯滤波进行降噪处理，可参考程序，也可自行编写，显示结果。

I=imread('lena.png’);
J=imnoise(I,'gaussian',0,0.005);
h=fspecial('gaussian');
K=filter2(h,J)/255;
figure, subplot(2,2,1),imshow(I),title('原图像');
subplot(2,2,2),imshow(J),title('加入高斯噪声的图像');
subplot(2,2,3),imshow(K),title('高斯低通滤波的结果');

（3）使用例子图片实现中值滤波器，可参考程序，也可自行编写，此处采用的是3*3的中值滤波器。

%非线性平滑滤波器
%加入椒盐噪声并采用模板实现中值过滤器
I=imread('peppers.png’);
J=imnoise(I,'salt',0.02);
figure ,subplot(2,2,1),imshow(I),title('原图像');
subplot(2,2,2),imshow(J),title('加入椒盐噪声后的图像');
K1=medfilt2(J,[3,3]);% 用模板实现中值过滤器
subplot(2,2,3),imshow(K1),title('中值滤波后的结果');

（4）读入一张灰度图像，比较不同窗口（模板）大小，至少3种（如5*5、9*9、15*15，25*25大小等），分别加入高斯噪声，椒盐噪声 和其它噪声如：指数分布噪声，泊松噪声，乘性噪声，比较不同模板大小，不同噪声下的使用均值滤波器，高斯滤波器和中值滤波的效果。

% 读入图像并转为灰度图 
img = imread('peppers.png'); 
img_gray = rgb2gray(img);
figure ,subplot(2,2,1),imshow(img_gray),title('原图像');
% 添加高斯噪声 
img_noise_gaussian = imnoise(img_gray, 'gaussian', 0, 0.01); 
% 添加椒盐噪声 
img_noise_saltpepper = imnoise(img_gray, 'salt & pepper', 0.02);
% 均值滤波器 
filter_mean = ones(5,5) / 25; 
% 5*5窗口 
img_filter_mean = imfilter(img_noise_gaussian, filter_mean);
subplot(2,2,2),imshow(img_filter_mean),title('均值滤波后的结果');

% 高斯滤波器 
filter_gaussian = fspecial('gaussian', [9 9], 0.5); 
% 9*9窗口，sigma = 0.5 
img_filter_gaussian = imfilter(img_noise_gaussian, filter_gaussian); 
subplot(2,2,3),imshow(img_filter_gaussian),title('高斯滤波后的结果');

% 中值滤波器 
img_filter_median = medfilt2(img_noise_saltpepper, [15 15]); 
subplot(2,2,4),imshow(img_filter_median),title('中值滤波后的结果');
% 15*15窗口

（5）尝试不使用matlab函数，实现高斯滤波器或者中值滤波器，并与 matlab自带函数对比分析是否一样。

%自编的均值滤波函数。x是需要滤波的图像,n是模板大小(即n×n)
function d=avefilt(x,n)   
a(1:n,1:n)=1;   %a即n×n模板,元素全是1
p=size(x);   %输入图像是p×q的,且p>n,q>n
x1=double(x);
x2=x1;
%A(a:b,c:d)表示A矩阵的第a到b行,第c到d列的所有元素
for i=1:p(1)-n+1
    for j=1:p(2)-n+1
        c=x1(i:i+(n-1),j:j+(n-1)).*a;  %取出x1中从(i,j)开始的n行n列元素与模板相乘
        s=sum(sum(c));                 %求c矩阵(即模板)中各元素之和
        x2(i+(n-1)/2,j+(n-1)/2)=s/(n*n); %将模板各元素的均值赋给模板中心位置的元素
    end
end
%未被赋值的元素取原值
d=uint8(x2);


%自编的中值滤波函数。x是需要滤波的图像,n是模板大小(即n×n)
function d=midfilt(x,n)   
p=size(x);   %输入图像是p×q的,且p>n,q>n
x1=double(x);
x2=x1;
for i=1:p(1)-n+1
    for j=1:p(2)-n+1
        c=x1(i:i+(n-1),j:j+(n-1));  %取出x1中从(i,j)开始的n行n列元素,即模板(n×n的)
        e=c(1,:);      %是c矩阵的第一行
        for u=2:n
            e=[e,c(u,:)];     %将c矩阵变为一个行矩阵   
        end
        mm=median(e);      %mm是中值
        x2(i+(n-1)/2,j+(n-1)/2)=mm;   %将模板各元素的中值赋给模板中心位置的元素
    end
end
%未被赋值的元素取原值
d=uint8(x2);



%自编的高斯滤波函数，S是需要滤波的图象，n是均值，k是方差
function d=gaussfilt(k,n,s)
Img = double(s);
n1=floor((n+1)/2);%计算图象中心
for i=1:n
    for j=1:n
      b(i,j) =exp(-((i-n1)^2+(j-n1)^2)/(4*k))/(4*pi*k);
    end
end
%生成高斯序列b。
Img1=conv2(Img,b,'same'); %用生成的高斯序列卷积运算，进行高斯滤波
d=uint8(Img1);

2. 图像的傅立叶变换

2.1 实验目的与要求

(1)了解图像变换的意义和手段;

(2)熟悉傅里叶变换的性质；

(3)热练掌握FFT方法的应用；

(4)通过实验了解二维频谱的分布特点；

(5)通过本实验掌握利用MATLAB编程实现数字图像的傅立叶变换。

2.2 实验原理及知识点

(1)应用傅立叶变换进行图像处理

傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。通过实验培养这项技能，将有助于解决大多数图像处理问题。对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。

(2)傅立叶（Fourier）变换的定义

对于二维信号，二维Fourier变换定义为：

二维离散傅立叶变换为：

图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目，有关傅立叶变换的快速算法的程序不难找到。实际上，现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。

(3)利用MATLAB软件实现数字图像傅立叶变换、DCT变换的程序。

2.3 实验仪器与软件

(1) PC计算机

(2) MatLab或opencv 

(3) 实验所需要的图片

2.4 实验报告要求

描述实验的基本步骤，用数据和图片给出各个步骤中取得的实验结果和源代码，并进行必要的讨论，必须包括原始图像及其计算/处理后的图像。

2.5 实验内容及步骤

(1)傅立叶变换。

a) 绘制一个二值图像矩阵,并将其傅立叶函数可视化。

f=zeros(30,30);
f(5:24,13:17)=1;
imshow(f,'notruesize')

F=fft2(f);
F2=log(abs(F));
figure,imshow(F2,[-1 5],'notruesize');colormap(jet); 

F=fft2(f,256,256); %零填充为256×256矩阵
figure,imshow(log(abs(F)),[-1 5],'notruesize');colormap(jet); 

F2=fftshift(F);    %将图像频谱中心由矩阵原点移至矩阵中心
figure,imshow(log(abs(F2)),[-1 5],'notruesize');colormap(jet);

b)利用傅立叶变换分析两幅图像的相关性，定位图像特征。读入图像‘text.png'，抽取其中的字母‘a’。

bw=imread('text.png');
a=bw(59:71,81:91);（注意:这里的坐标有点问题，坐标需要自行重新定位字母a的位置）
imshow(bw);

figure,imshow(a);

C=real(ifft2(fft2(bw).*fft2(rot90(a,2),256,256)));%求相关性
figure,imshow(C,[]);

thresh=max(C(:));
figure,imshow(C>thresh-10)

figure,imshow(C>thresh-15)

(2) 离散余弦变换(DCT)

a) 使用dct2对图像‘autumn.tif’进行DCT变换(2) 离散余弦变换(DCT)

RGB=imread('autumn.tif');
imshow(RGB)

I=rgb2gray(RGB); %转换为灰度图像
figure,imshow(I)

J=dct2(I);
figure,imshow(log(abs(J)),[]),colormap(jet(64));colorbar;

b) 将上述DCT变换结果中绝对值小于10的系数舍弃，使用idct2重构图像并与原图像比较。

RGB=imread('autumn.tif');
I=rgb2gray(RGB); 		%转换为灰度图像
J=dct2(I);
figure,imshow(I)

K=idct2(J);
figure,imshow(K,[0 255])

J(abs(J)<10)=0; 			%舍弃系数
K2=idct2(J);
figure,imshow(K2,[0 255]);

c) 利用DCT变换进行图像压缩。

I=imread('cameraman.tif');
I=im2double(I);
T=dctmtx(8);
B=blkproc(I,[8,8],'P1*x*P2',T,T');
mask=[1  1  1  1  0  0  0  0
      1  1  1  0  0  0  0  0
      1  1  0  0  0  0  0  0
      1  0  0  0  0  0  0  0
      0  0  0  0  0  0  0  0
      0  0  0  0  0  0  0  0
      0  0  0  0  0  0  0  0
      0  0  0  0  0  0  0  0];
B2=blkproc(B,[8 8],'P1.*x',mask);
I2=blkproc(B2,[8 8],'P1*x*P2',T',T);
imshow(I)
figure,imshow(I2);

（3）对(1)和(2)中经过FFT和DCT的图像进行反变换,试图保留不同维度的低频信息,将去掉的高频成分设置为0,显示保留不同维度的低频信息反变换后的图像,分析保留不同维度下反变换图像的质量区别.

F3 = iff2(F);
figure,imshow(log(abs(F3)),[-1 5],'InitialMagnification','fit');colormap(jet);

由于进行傅里叶变换的是矩阵，将空域转换成了频域，而傅里叶反变换则将空域 转换回频域，重新显现出原来的图形。

②对(2)中的图像进行反变换

K=idct2(J);
figure,imshow(K,[0 255])

(4)讨论不同的图像内容与FFT、DCT频谱之间的对应关系。

进行傅里叶变换（FFT2）后，频谱的低频部分主要集中在四周，中间则为高频部分； 而DCT则从左上角到右下角从低频到高频，系数逐渐减小。相同点为图像能量均集中在 低频部分，低频系数体现出图像中目标的轮廓和灰度分布特性，高频系数体现了目标形状细节。需要注意的是，频谱图上各点和原图像不存在一一对应的关系，即使是在不移频的情况下也没有。

3 频域增强实验

3.1 实验原理

  卷积理论是频域技术的基础。设函数f(x,y)与线性位不变算子h(x,y)的卷积结果是g(x,y)，即

g(x,y)=h(x,y)*f(x,y)                         

那么，根据卷积定理在频域有

G(u,v)=H(u,v)F(u,v)                        

其中，G(u,v)、H(u,v)、F(u,v)分别是 g(x,y)、h(x,y)、f(x,y)的傅立叶变换。频域增强的主要步骤是：

1. 计算所需增强图的傅立叶变换；
2. 将其与一个（根据需要设计的）转移函数相乘；
3. 在将结果进行傅立叶反变换以得到增强的图。

  常用的频域增强方法有低通滤波和高通滤波。下面分别介绍这两种滤波方法。

（1）低通滤波

图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分。因此能降低高频成分幅度的滤波器就能减弱噪声的影响。

Butterworth低通滤波器是一种物理上可以实现的低通滤波器。n阶截断频率为d0的Butterworth低通滤波器的转移函数为：

%  实现Butterworth低通滤波器
I1=imread('lena.png');
I2=imnoise(I1,'salt');
f=double(I2);
g=fft2(f);     %采用傅立叶变换
g=fftshift(g); %数据矩阵平移
[N1,N2]=size(g);
n=2;
d0=50;
d1=5;
n1=fix(N1/2);
n2=fix(N2/2);
for i=1:N1
    for j=1:N2
        d=sqrt((i-n1)^2+(j-n2)^2);
        %计算Butterworth低通变换函数
        h=1/(1+0.414*(d/d0)^(2*n));
        result(i,j)=h*g(i,j);
    end
end
result=ifftshift(result);
X2=ifft2(result);
X3=uint8(real(X2));

（2）高通滤波

   高通滤波也称高频滤波，它的频率值在0频率处为单位1，随着频率的增长，传递函数的值逐渐增加；当频率增加到一定值后，传递函数的值通常又回到0值或者降低到某个大于1的值。在前一种情况下，高频增强滤波器实际上是一种带通滤波器，只不过规定0频率处的增益为单位1。

在实际应用中，为了减少图像种面积大而且缓慢变化的成分的对比度，有时让0频率处的增益小雨单位1更合适。如果传递函数通过原点，则可称为laplacian滤波器。

% 实现Butterworth高通滤波器参考代码
%计算Butterworth高通变换函数
if d==0
    h=0;
else
    h=1/(1+(d1/d)^(2*n));
end
    result(i,j)=h*g(i,j);
result=ifftshift(result);
X4=ifft2(result);
X5=uint8(real(X2));
figure(1),subplot(2,2,1),imshow(I1),title('原图像');
subplot(2,2,2),imshow(I2),title('加噪图像');
subplot(2,2,3),imshow(X3),
title('Butterworth低通滤波器去噪图像');
subplot(2,2,4),imshow(X5),
title('Butterworth高通滤波器去噪图像');

3.2 实验内容和要求

分别用两张以上不同的图片， 实现Butterworth低通滤波和Butterworth高通滤波， 分别加入三种不同噪声如高斯噪声，椒盐噪声,以及指数分布噪声，泊松噪声，乘性噪声等，分析Butterworth高通和低通的滤波效果。分析比较Butterworth低通滤波器和上面几种空域平滑滤波器的差别，结合实验结果对比图，进行分析。

% 生成原始图像  
imageA = phantom();  
imageB = imread('cameraman.tif');  
  
% 添加高斯噪声  
imageA_noise = imnoise(imageA,'gaussian',0,0.01);  
imageB_noise = imnoise(imageB,'gaussian',0,0.01);  
imageA_noise1 = imnoise(imageA,'salt & pepper',0.02);  
imageB_noise1 = imnoise(imageB,'salt & pepper',0.02);  
imageA_noise2 = imnoise(imageA,'poisson');  
imageB_noise2 = imnoise(imageB,'poisson');  

  
% 添加其他类型的噪声，比如椒盐噪声、指数分布噪声、泊松噪声等，只需更改'gaussian'为'salt & pepper'、'exponential'、'poisson'等。  
  
% 显示噪声图像  
subplot(3,2,1), imshow(imageA_noise), title('高斯噪声')  
subplot(3,2,2), imshow(imageB_noise), title('高斯噪声')
subplot(3,2,3), imshow(imageA_noise1), title('椒盐噪声')  
subplot(3,2,4), imshow(imageB_noise1), title('椒盐噪声')
subplot(3,2,5), imshow(imageA_noise2), title('泊松噪声')  
subplot(3,2,6), imshow(imageB_noise2), title('泊松噪声')

实验总结：

实验中对两张不同的图片分别添加了高斯噪声、椒盐噪声、泊松噪声，并分析比较了这两种滤波器的效果。

实验结果表明，Butterworth低通滤波器在去除图像噪声方面具有良好的效果，尤其对于高频噪声，如椒盐噪声和乘性噪声等，去除效果更为明显。这是由于低通滤波器可以减弱高频成分幅度，进而降低噪声影响。然而，它也可能导致图像边缘的模糊化。

相比之下，Butterworth高通滤波器则能够增强图像的边缘特征，提高图像的清晰度。这是因为高通滤波器可以让高频成分通过，而这些成分通常对应图像的边缘和细节。但是，对于噪声较多或者质量较差的图像，高通滤波可能会增强噪声，影响图像质量。

与空域平滑滤波器相比，如均值滤波器和中值滤波器等，Butterworth低通和高通滤波器在处理图像噪声和增强图像细节方面具有一定的优势。它们能够在保护图像边缘和细节的同时，更有效地去除噪声。

Butterworth低通和高通滤波器在图像处理中具有重要的应用价值。在实际应用中，应根据具体需求和图像特点选择合适的滤波器类型和参数设置，以达到最佳的图像处理效果。同时，实验结果也验证了频域技术在图像处理中的有效性和优越性。