空间域滤波器(实验报告)

数字图像处理作业

——空间域滤波器

摘要

在图像处理的过程中,消除图像的噪声干扰是一个非常重要的问题。本文利用matlab软件,采用空域滤波的方式,对图像进行平滑和锐化处理。平滑空间滤波器用于模糊处理和减小噪声，经常在图像的预处理中使用；锐化空间滤波器主要用于突出图像中的细节或者增强被模糊了的细节。本文使用的平滑滤波器有中值滤波器和高斯低通滤波器，其中，中值滤波器对去除椒盐噪声特别有效，高斯低通滤波器对去除高斯噪声效果比较好。使用的锐化滤波器有反锐化掩膜滤波、Sobel边缘检测、Laplacian边缘检测以及Canny算子边缘检测滤波器。不同的滤波方式，在特定的图像处理应用中有着不同的效果和各自的优势。

1、分别用高斯滤波器和中值滤波器去平滑测试图像test1和2，模板大小分别是3x3 ， 5x5 ，7x7；利用固定方差 sigma=1.5产生高斯滤波器. 附件有产生高斯滤波器的方法。

实验原理分析：

空域滤波是直接对图像的数据做空间变换达到滤波的目的。它是一种邻域运算，其机理就是在待处理的图像中逐点地移动模板，滤波器在该点地响应通过事先定义的滤波器系数与滤波模板扫过区域的相应像素值的关系来计算。如果输出像素是输入像素邻域像素的线性组合则称为线性滤波（例如最常见的均值滤波和高斯滤波），否则为非线性滤波（中值滤波、边缘保持滤波等）。

空域滤波器从处理效果上可以平滑空间滤波器和锐化空间滤波器：平滑空间滤波器用于模糊处理和减小噪声，经常在图像的预处理中使用；锐化空间滤波器主要用于突出图像中的细节或者增强被模糊了的细节。

模板在源图像中移动的过程中，当模板的一条边与图像轮廓重合后，模板中心继续向图像边缘靠近，那么模板的某一行或列就会处于图像平面之外，此时最简单的方法就是将模板中心点的移动范围限制在距离图像边缘不小于（n－1）/2个像素处，单处理后的图像比原始图像稍小。如果要处理整幅图像，可以在图像轮廓边缘时用全部包含于图像中的模板部分来滤波所有图像，或者在图像边缘以外再补上一行和一列灰度为零的像素点（或者将边缘复制补在图像之外）。

①中值滤波器的设计：

中值滤波器是一种非线性统计滤波器，它的响应基于图像滤波器包围的图像区域中像素的排序，然后由统计排序的中间值代替中心像素的值。它比小尺寸的线性平滑滤波器的模糊程度明显要低，对处理脉冲噪声（椒盐噪声）非常有效。中值滤波器的主要功能是使拥有不同灰度的点看起来更接近于它的邻近值，去除那些相对于其邻域像素更亮或更暗，并且其区域小于滤波器区域一半的孤立像素集。

在一维的情况下，中值滤波器是一个含有奇数个像素的窗口。在处理之后，位于窗口正中的像素的灰度值，用窗口内各像素灰度值的中值代替。例如若窗口长度为5，窗口中像素的灰度值为80、90、200、110、120，则中值为110，因为按小到大（或大到小）排序后，第三位的值是110。于是原理的窗口正中的灰度值200就由110取代。如果200是一个噪声的尖峰，则将被滤除。然而，如果它是一个信号，则滤波后就被消除，降低了分辨率。因此中值滤波在某些情况下抑制噪声，而在另一些情况下却会抑制信号。

将中值滤波推广到二维的情况。二维窗口的形式可以是正方形、近似圆形的或十字形等。本次作业使用正方形模板进行滤波，它的中心一般位于被处理点上。窗口的大小对滤波效果影响较大。

根据上述算法利用MATLAB软件编程，对源图像test1和test2进行滤波处理，结果如下图：

中值滤波后的test1.pgm（3x3

）

中值滤波后的test1.pgm（5x5

）

中值滤波后的test1.pgm（7x7）

中值滤波后的test2.tif（3x3

）

中值滤波后的test2.tif（5x5）

中值滤波后的test2.tif（7x7）

可见，窗口的大小对滤波效果影响较大。窗口越大，平滑效果越明显，图像细节越模糊。

②高斯滤波器的设计：

高斯滤波是一种根据高斯函数的形状来选择模板权值的线性平滑滤波方法。高斯平滑滤波器对去除服从正态分布的噪声是很有效果的。一维零均值高斯函数为 。其中，高斯分布参数决定了高斯滤波器的宽度。对图像来说，常用二维零均值离散高斯函数作平滑滤波器，函数表达式如下：

（1）

高斯函数具有5个重要性质：

（1）二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的。一般来说一幅图像的边缘方向是不知道的。因此，在滤波之前是无法确定一个方向比另一个方向上要更多的平滑的。旋转对称性意味着高斯滤波器在后续的图像处理中不会偏向任一方向。

（2）高斯函数是单值函数。高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点的权值是随着该点与中心点距离单调递减的。这一性质是很重要的，因为边缘是一种图像局部特征。如果平滑运算对离算子中心很远的像素点仍然有很大的作用，则平滑运算会使图像失真。

（3）高斯函数的傅立叶变换频谱是单瓣的。这一性质是高斯函数傅立叶变换等于高斯函数本身这一事实的直接推论。图像常被不希望的高频信号所污染，而所希望的图像特征，既含有低频分量，又含有高频分量。高斯函数傅立叶变换的单瓣意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需要的信号。

（4）高斯滤波器的宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的。σ越大，高斯滤波器的频带就越宽，平滑程度就越好。通过调节平滑程度参数σ，可在图像特征分量模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷。

（5）由于高斯函数的可分离性，大高斯滤波器可以有效实现。通过二维高斯函数的卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积的结果与方向垂直的相同一维高斯函数进行卷积。因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长。这些性质使得它在早期的图像处理中特别有用，表明高斯平滑滤波器无论在空间域还是在频率域都是十分有效的低通滤波器。

根据上述分析，利用MATLAB软件设计高斯滤波器，对源图像test1和test2进行滤波处理，结果如下图：

高斯滤波平滑后的test1.pgm（3x3）

50

100

150

200

250

[1**********]250 高斯滤波平滑后的test1.pgm（5x5）

50

100

150

200

250

[1**********]250

高斯滤波平滑后的test1.pgm（7x7）

50

100

150

200

250

[1**********]250

高斯滤波平滑后的test2.tif（3x3）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

高斯滤波平滑后的test2.tif（5x5）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

高斯滤波平滑后的test2.tif（7x7）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

可见，对于高斯滤波器，模板的大小对滤波效果影响不大。高斯滤波虽然能够在一定程度上去掉噪声，但也使得图象变得模糊不清，效果并不能令人满意。

2、利用高通滤波器滤波测试图像test3,4：包括unsharp masking, Sobel edge detector, and Laplace edge detection；Canny algorithm.

实验原理分析：

锐化滤波能减弱或消除图像中的低频率分量，但不影响高频率分量。因为低频分量对应图像中灰度值缓慢变化的区域，因而与图像的整体特性，如整体对比度和平均灰度值等有关。锐化滤波将这些分量滤去可使图像反差增加，边缘明显。在实际应用中，锐化滤波可用于增强被模糊的细节或者低对比度图像的目标边缘。

图像锐化的主要目的有两个：一是增强图像边缘，使模糊的图像变得更加清晰，颜色变得鲜明突出，图像的质量有所改善，产生更适合人眼观察和识别的图像；二是希望经过锐化处理后，目标物体的边缘鲜明，以便于提取目标的边缘、对图像进行分割、目标区域识别、区域形状提取等，为进一步的图像理解与分析奠定基础。

由于锐化使噪声受到比信号还要强的增强，所以要求锐化处理的图像有较高的信噪比；否则，锐化后图像的信噪比更低。

①反锐化掩膜图像增强（unsharp masking）

图像的反锐化掩蔽算法可以表示为：

(1) 其中fs(x,y)表示经过反锐化掩蔽得到的锐化图像，是f(x,y)的模糊形式。反锐化掩蔽进一步的普遍形式称为高提升滤波。在图像中任何一点(x,y)处，高提升滤波后的图像fhb可定义如下：

(2) 其中A≥1，与前式一样，是f的模糊形式，此式也可以写成：

(3)

结合式(1)，可以得到：

这一表达式可计算高提升滤波图像。

如果选择拉普拉斯变换，式(4)变成： (4)

(5) 高提升滤波处理可以通过任何一个图1所示的掩模得以实现。当A=1时，高提升滤波处理就是标准的拉普拉斯变换。随着A超过l不断增大，锐化处理的效果越来越不明显。最终，当A足够大时，高提升图像将近似等于经常数调制的图像。

图1 高频提升滤波技术可以用其中一种掩膜来实现（A1） 本文采用的反锐化掩膜滤波器中A=1。根据以上分析，利用MATLAB软件设计反锐化掩膜滤波器对test3和test4进行滤波，结果如下图：

原图像——test3c

orrupt.pgm

反锐化掩膜后的test3corrupt.pgm

原图像

——test4 copy.bmp反锐化掩膜后的test4 copy.bmp

可见，经过反锐化掩膜滤波后，图像的边缘得到了增强，细节更加明显，但

同时也带来了高频噪声的影响。

②索贝尔边缘检测（Sobel edge detector）

索贝尔算子（Sobel operater）主要用作边缘检测，在技术上，它是一离散性差分算子，用来运算图像亮度函数的灰度之近似值。在图像的任何一点使用此算子，将会产生对应的灰度矢量或其法矢量。

该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像灰度值，其公式如下：

1Gx2

1

000

12*A1

1

Gy0

1

21



00A*

21

图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的

大小。

G

然后可用以下公式计算梯度方向。

arctan(

GyGx

)

如果角度Θ等于零，即代表图像在该处拥有纵向边缘，左方较右方暗。

根据以上算法分析，利用MATLAB软件设计Sobel边缘检测滤波器，对源

图像test3和test4进行滤波，结果如下图示：

源图像——test3c

orrupt.pgmSobel 边缘检测的test3c

orrupt.pgm

源图像——test3c

orrupt.pgm

Sobel 边缘检测——test4 copy.bmp

③Laplacian边缘检测（Laplacian edge detector）

拉普拉斯算子是最简单的各向同性微分算子，具有旋转不变性。一个二维图像函数

的拉普拉斯变换是各向同性的二阶导数，定义为：

（6）

为了更适合于数字图像处理，将该方程表示为离散形式：

（7）

另外，拉普拉斯算子还可以表示成模板的形式，如图2所示。图2（a）表示离散拉普拉斯算子的模板，图2（b）表示其扩展模板，图2（c）和（d）则分别表示其他两种拉普拉斯的实现模板。从模板形式容易看出，如果在图像中一个较暗的区域中出现了一个亮点，那么用拉普拉斯运算就会使这个亮点变得更亮。因为图像中的边缘就是那些灰度发生跳变的区域，所以拉普拉斯锐化模板在边缘检测中很有用。

一般增强技术对于陡峭的边缘和缓慢变化的边缘很难确定其边缘线

的位置。但该算子却可用二次微分正峰和负峰之间的过零点来确定，对孤立点或端点更为敏感，因此特别适用于以突出图像中的孤立点、孤立线或线端点为目的的场合。同梯度算子一样，拉普拉斯算子也会增强图像中的噪声，有时用拉普拉

（a（b

图2 拉普拉斯的4种模板

拉普拉斯锐化的基本方法可以由下式表示：

这种简单的锐化方法既可以产生拉普拉斯锐化处理的效果，同时又能保留背景信息，将原始图像叠加到拉普拉斯变换的处理结果中去，可以使图像中的各灰度值得到保留，使灰度突变处的对比度得到增强，最终结果是在保留图像背景的前提下，突现出图像中小的细节信息。

根据以上算法分析，利用MATLAB软件设计Laplacian边缘检测滤波器，

对源图像test3和test4进行滤波，结果如下图示：

源图像—test3corrupt.pgm

Laplacian边缘检测-test3corrupt.pgm

源图像

—test4 copy.bmpLaplacian边缘检测（未滤波）

—test4 copy.bmp

由上图可知，运用Laplacian算子对test4 copy.bmp进行边缘检测后，图像的边缘信息没有被很好的检测出来。这是由于拉普拉斯梯算子会增强图像中的噪声，因此本文再用拉普拉斯算子对test4 copy.bmp进行边缘检测时，先将图像进行高斯平滑滤波处理。处理结果如下图：

源图像

—test4 copy.bmp

Laplacian边缘检测(高斯滤波后)—test4 copy.bmp

④Canny算子边缘检测（Canny algorithm）

在图像边缘检测中，抑制噪声和边缘精确定位是无法同时满足的。边缘检测算法通过平滑滤波去除图像噪声的同时，也增加了边缘定位的不确定性；反之，提高边缘检测算子对边缘敏感性的同时，也提高了对噪声的敏感性。Canny算子力图在抗噪声干扰和精确定位边缘之间寻求最佳折中方案。用Canny算子检测图像边缘的步骤如下：

step1:用高斯滤波器平滑图象；

step2:用一阶偏导的有限差分来计算滤波后图像梯度的幅值和方向； step3:对梯度幅值进行非极大值抑制，其过程为找出图像梯度中的局部极大值点，把其他非局部极大值点置零以得到细化的边缘。

step4:用双阈值算法检测和连接边缘。使用两个阈值T1和T2（T1>T2），T1用来找到每条线段，T2用来在这些线段的两个方向上延伸寻找边缘的断裂处，并连接这些边缘。 利用上述原理和MATLAB软件设计滤波器，对test3和test4进行边缘检测，结果如下图所示：

源图像-test3c

orrupt.pgmCanny边缘检测-test3c

orrupt.pgm

Canny边缘检测-test4 copy.bmp

可见，经过Canny算子较好的检测出了图像的所有边缘信息，具有较强的去噪能力。

附录

一、参考文献

[1] 冈萨雷斯著.数字图像处理(第三版).北京：电子工业出版社，2010

[2] 杨杰 李庆著.数字图像处理及MATLAB实现——学习与实验指导.北京：电子

工业出版社，2010

[3] 苏金明 王永利著.MATLAB图形图像. 北京：电子工业出版社，2005

[4] 朱习军 隋思涟等著.MATLAB在信号与图像处理中的应用. 北京：电子工业

出版社，2009

[5] 张秀兰著.基于MATLAB的数字图像的边缘检测.吉林化工学院学报，2010 [6] 杨先花 黎粤华著.基于MATLAB图像边缘检测算法效果对比.机电产品开发

与创新，2010

二、源代码：

第一题

(1) 平滑滤波器

1. 中值滤波器（以3x3的模板大小为例）

①I=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test1.pgm','pgm'); n=3; a=ones(n,n); p=size(I);

x1=double(I);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)); e=c(1,:); for u=2:n

e=[e,c(u,:)]; end

mm=median(e);

x2(i+(n-1)/2,j+(n-1)/2)=mm; end end

I2=uint8(x2); imshow(I2)

title('中值滤波后的test1.pgm（3x3）')

②I=imread('E:\大三下\图像处理英文课件\作业\第四次作业\test2.tif','tif'); n=3;

a=ones(n,n); p=size(I);

x1=double(I);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)); e=c(1,:); for u=2:n

e=[e,c(u,:)]; end

mm=median(e);

x2(i+(n-1)/2,j+(n-1)/2)=mm; end end

I2=uint8(x2); imshow(I2)

title('中值滤波后的test2.tif（3x3）')

2、高斯滤波器（以3x3的模板大小为例）

① n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0;

[I,map]=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test1.pgm','pgm'); r=[cos(theta) -sin(theta); sin(theta) cos(theta)];

for i = 1 : n2 for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2

*sigma2^2))/(sigma2*sqrt(2*pi));

end

end

h = h / sqrt(sum(sum(h.*h))); f1=conv2(I,h,'same'); figure(1);

imagesc(I);title('test1.pgm'); colormap(gray); figure(2);

imagesc(f1);title(' '高斯滤波平滑后的test1.pgm(3x3)'); colormap(gray);

② n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0;

[I,map]=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test2.tif','tif');

r=[cos(theta) -sin(theta); sin(theta) cos(theta)]; for i = 1 : n2 for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2 end end

*sigma2^2))/(sigma2*sqrt(2*pi));

h = h / sqrt(sum(sum(h.*h))); f1=conv2(I,h,'same'); figure(1);

imagesc(I);title('test2.tif'); colormap(gray);

figure(2);

imagesc(f1);title(' '高斯滤波平滑后的test2.tif(3x3)'); colormap(gray);

（2）锐化滤波器：

①反锐化掩膜（以处理test3_corrupt.pgm为例）

pic = imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); I=pic;

picSize = size(pic);

pic = cast(pic,'int32');

t = zeros(picSize(1) + 2 , picSize(2) + 2); t = cast(t,'int32');

t(2:picSize(1) + 1 , 2:picSize(2) + 1) = pic; t(: , 1) = t(: , 2);

t(: , picSize(2) + 2) = t(: , picSize(2) + 1); t(1 , :) = t(2 , :);

t(picSize(1) + 2 , :) = t(picSize(1) + 1 , :); A = 1.0;

for i=2:1:picSize(1)+1

for j=2:1:picSize(2)+1

pic(i-1,j-1) = t(i,j)*A - (t(i,j)*(-8) + t(i,j-1) + t(i-1,j-1) + t(i-1,j)

+ t(i,j+1) + t(i+1,j) + t(i+1,j+1) + t(i+1,j-1) + t(i-1,j+1));

end

end

pic = cast(pic,'uint8');

imwrite(pic,['UnsharpMasking_A_', num2Str(10*A),'.png'],'png'); figure; subplot(1,2,1) imshow(I);

title('原图像——test4 copy.bmp'); subplot(1,2,2) imshow(pic);

title('反锐化掩膜后的test4 copy.bmp');

②Sobel 边缘检测（以处理test4_corrupt.pgm为例）

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp');

[N,M]=size(I) I=double(I);

h1=[-1,0,1;-2,0,2;-1,0,1]; h2=[-1,-2,-1;0,0,0;1,2,1]; Gx=conv2(I,h1,'same'); Gy=conv2(I,h2,'same'); F=abs(Gx)+abs(Gy); for i=1:N for j=1:M

I(i,j)=F(i,j); end end

I=uint8(I); figure; imshow(I);

title('Sobel 边缘检测——test4 copy.bmp');

③Laplace边缘检测

a）Laplace边缘检测——test3_corrupt.pgm

I = imread('E:\大三下\图像处理英文课件\作业\第四次作业\test3_corrupt.pgm','pgm'); [N,M]=size(I); figure;

subplot(1,2,1); imshow(I);

title('ԴͼÏñ¡ªtest3_corrupt.pgm'); t=ones(N,M); t=I;

for i=2:N-1

for j=2:M-1

I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j); end

end

subplot(1,2,2);

imshow(I);

title('Laplacian边缘检测-test3_corrupt.pgm');

b）Laplace边缘检测——test4 copy.bmp

未进行高斯滤波：

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); [N,M]=size(I); figure;

subplot(1,2,1);

imshow(I);

title('源图像—test4 copy.bmp');

t=ones(N,M);

t=I;

for i=2:N-1

for j=2:M-1

I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j);

end

end

subplot(1,2,2);

imshow(I);

title('Laplacian边缘检测(未滤波)-test4 copy.bmp');

加入高斯滤波：

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp');

[N,M]=size(I);

figure;

subplot(1,2,1);

imshow(I);

title('源图像—test4 copy.bmp ');

n1=3;sigma1=0.5;n2=3;sigma2=0.5;theta=0;

r=[cos(theta) -sin(theta); sin(theta) cos(theta)];

for i = 1 : n2

for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2*s end

end

h = h / sqrt(sum(sum(h.*h)));

f1=conv2(I,h,'same');

t=ones(N,M);

t=f1;

for i=2:N-1

for j=2:M-1

f1(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j);

end

end

subplot(1,2,2);

imshow(f1);

title(' Laplacian边缘检测(高斯滤波后)-test4 copy.bmp ');

igma2^2))/(sigma2*sqrt(2*pi));

④ Canny算子边缘检测（以处理test4 copy.bmp为例）

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); a = im2single(I);

[m,n] = size(a);

e = false(m,n);

GaussianDieOff = 0.0001;

PercentOfPixelsNotEdges = 0.7;

ThresholdRatio = 0.4;

sigma = 1;

thresh = [];

pw = 1:30;

ssq = sigma^2;

width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');

if isempty(width)

width = 1;

end

t = (-width:width);

gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);

[x,y]=meshgrid(-width:width,-width:width);

dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);

aSmooth=imfilter(a,gau,'conv','replicate');

aSmooth=imfilter(aSmooth,gau','conv','replicate');

ax = imfilter(aSmooth, dgau2D, 'conv','replicate');

ay = imfilter(aSmooth, dgau2D', 'conv','replicate');

mag = sqrt((ax.*ax) + (ay.*ay));

magmax = max(mag(:));

if magmax>0

mag = mag / magmax;

end

if isempty(thresh)

counts=imhist(mag, 64);

highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,... 1,'first') / 64;

lowThresh = ThresholdRatio*highThresh;

thresh = [lowThresh highThresh];

elseif length(thresh)==1

highThresh = thresh;

if thresh>=1

eid = sprintf('Images:%s:thresholdMustBeLessThanOne', mfilename); msg = 'The threshold must be less than 1.';

error(eid,'%s',msg);

end

lowThresh = ThresholdRatio*thresh;

thresh = [lowThresh highThresh];

elseif length(thresh)==2

lowThresh = thresh(1);

highThresh = thresh(2);

if (lowThresh >= highThresh) || (highThresh >= 1)

eid = sprintf('Images:%s:thresholdOutOfRange', mfilename); msg = 'Thresh must be [low high], where low

end

end

idxStrong = [];

for dir = 1:4

[m,n] = size(mag);

idx = find((ay-ay) | (ay>=0 & ax

if ~isempty(idx)

v = mod(idx,m);

extIdx = find(v==1 | v==0 | idx(n-1)*m));

idx(extIdx) = [];

end

ixv = ax(idx);

iyv = ay(idx);

gradmag = mag(idx);

d = abs(iyv./ixv);

gradmag1 = mag(idx+m).*(1-d) + mag(idx+m-1).*d;

gradmag2 = mag(idx-m).*(1-d) + mag(idx-m+1).*d;

idxLocMax = idx(gradmag>=gradmag1 & gradmag>=gradmag2); idxWeak = idxLocMax(mag(idxLocMax) > lowThresh);

e(idxWeak)=1;

idxStrong = [idxStrong; idxWeak(mag(idxWeak) > highThresh)]; end

if ~isempty(idxStrong)

rstrong = rem(idxStrong-1, m)+1;

cstrong = floor((idxStrong-1)/m)+1;

e = bwselect(e, cstrong, rstrong, 8);

e = bwmorph(e, 'thin', 1);

end

imshow(e);

title('Canny边缘检测-test4 copy.bmp');

数字图像处理作业

——空间域滤波器

摘要

在图像处理的过程中,消除图像的噪声干扰是一个非常重要的问题。本文利用matlab软件,采用空域滤波的方式,对图像进行平滑和锐化处理。平滑空间滤波器用于模糊处理和减小噪声，经常在图像的预处理中使用；锐化空间滤波器主要用于突出图像中的细节或者增强被模糊了的细节。本文使用的平滑滤波器有中值滤波器和高斯低通滤波器，其中，中值滤波器对去除椒盐噪声特别有效，高斯低通滤波器对去除高斯噪声效果比较好。使用的锐化滤波器有反锐化掩膜滤波、Sobel边缘检测、Laplacian边缘检测以及Canny算子边缘检测滤波器。不同的滤波方式，在特定的图像处理应用中有着不同的效果和各自的优势。

1、分别用高斯滤波器和中值滤波器去平滑测试图像test1和2，模板大小分别是3x3 ， 5x5 ，7x7；利用固定方差 sigma=1.5产生高斯滤波器. 附件有产生高斯滤波器的方法。

实验原理分析：

空域滤波是直接对图像的数据做空间变换达到滤波的目的。它是一种邻域运算，其机理就是在待处理的图像中逐点地移动模板，滤波器在该点地响应通过事先定义的滤波器系数与滤波模板扫过区域的相应像素值的关系来计算。如果输出像素是输入像素邻域像素的线性组合则称为线性滤波（例如最常见的均值滤波和高斯滤波），否则为非线性滤波（中值滤波、边缘保持滤波等）。

空域滤波器从处理效果上可以平滑空间滤波器和锐化空间滤波器：平滑空间滤波器用于模糊处理和减小噪声，经常在图像的预处理中使用；锐化空间滤波器主要用于突出图像中的细节或者增强被模糊了的细节。

模板在源图像中移动的过程中，当模板的一条边与图像轮廓重合后，模板中心继续向图像边缘靠近，那么模板的某一行或列就会处于图像平面之外，此时最简单的方法就是将模板中心点的移动范围限制在距离图像边缘不小于（n－1）/2个像素处，单处理后的图像比原始图像稍小。如果要处理整幅图像，可以在图像轮廓边缘时用全部包含于图像中的模板部分来滤波所有图像，或者在图像边缘以外再补上一行和一列灰度为零的像素点（或者将边缘复制补在图像之外）。

①中值滤波器的设计：

中值滤波器是一种非线性统计滤波器，它的响应基于图像滤波器包围的图像区域中像素的排序，然后由统计排序的中间值代替中心像素的值。它比小尺寸的线性平滑滤波器的模糊程度明显要低，对处理脉冲噪声（椒盐噪声）非常有效。中值滤波器的主要功能是使拥有不同灰度的点看起来更接近于它的邻近值，去除那些相对于其邻域像素更亮或更暗，并且其区域小于滤波器区域一半的孤立像素集。

在一维的情况下，中值滤波器是一个含有奇数个像素的窗口。在处理之后，位于窗口正中的像素的灰度值，用窗口内各像素灰度值的中值代替。例如若窗口长度为5，窗口中像素的灰度值为80、90、200、110、120，则中值为110，因为按小到大（或大到小）排序后，第三位的值是110。于是原理的窗口正中的灰度值200就由110取代。如果200是一个噪声的尖峰，则将被滤除。然而，如果它是一个信号，则滤波后就被消除，降低了分辨率。因此中值滤波在某些情况下抑制噪声，而在另一些情况下却会抑制信号。

将中值滤波推广到二维的情况。二维窗口的形式可以是正方形、近似圆形的或十字形等。本次作业使用正方形模板进行滤波，它的中心一般位于被处理点上。窗口的大小对滤波效果影响较大。

根据上述算法利用MATLAB软件编程，对源图像test1和test2进行滤波处理，结果如下图：

中值滤波后的test1.pgm（3x3

）

中值滤波后的test1.pgm（5x5

）

中值滤波后的test1.pgm（7x7）

中值滤波后的test2.tif（3x3

）

中值滤波后的test2.tif（5x5）

中值滤波后的test2.tif（7x7）

可见，窗口的大小对滤波效果影响较大。窗口越大，平滑效果越明显，图像细节越模糊。

②高斯滤波器的设计：

高斯滤波是一种根据高斯函数的形状来选择模板权值的线性平滑滤波方法。高斯平滑滤波器对去除服从正态分布的噪声是很有效果的。一维零均值高斯函数为 。其中，高斯分布参数决定了高斯滤波器的宽度。对图像来说，常用二维零均值离散高斯函数作平滑滤波器，函数表达式如下：

（1）

高斯函数具有5个重要性质：

（1）二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的。一般来说一幅图像的边缘方向是不知道的。因此，在滤波之前是无法确定一个方向比另一个方向上要更多的平滑的。旋转对称性意味着高斯滤波器在后续的图像处理中不会偏向任一方向。

（2）高斯函数是单值函数。高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点的权值是随着该点与中心点距离单调递减的。这一性质是很重要的，因为边缘是一种图像局部特征。如果平滑运算对离算子中心很远的像素点仍然有很大的作用，则平滑运算会使图像失真。

（3）高斯函数的傅立叶变换频谱是单瓣的。这一性质是高斯函数傅立叶变换等于高斯函数本身这一事实的直接推论。图像常被不希望的高频信号所污染，而所希望的图像特征，既含有低频分量，又含有高频分量。高斯函数傅立叶变换的单瓣意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需要的信号。

（4）高斯滤波器的宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的。σ越大，高斯滤波器的频带就越宽，平滑程度就越好。通过调节平滑程度参数σ，可在图像特征分量模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷。

（5）由于高斯函数的可分离性，大高斯滤波器可以有效实现。通过二维高斯函数的卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积的结果与方向垂直的相同一维高斯函数进行卷积。因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长。这些性质使得它在早期的图像处理中特别有用，表明高斯平滑滤波器无论在空间域还是在频率域都是十分有效的低通滤波器。

根据上述分析，利用MATLAB软件设计高斯滤波器，对源图像test1和test2进行滤波处理，结果如下图：

高斯滤波平滑后的test1.pgm（3x3）

50

100

150

200

250

[1**********]250 高斯滤波平滑后的test1.pgm（5x5）

50

100

150

200

250

[1**********]250

高斯滤波平滑后的test1.pgm（7x7）

50

100

150

200

250

[1**********]250

高斯滤波平滑后的test2.tif（3x3）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

高斯滤波平滑后的test2.tif（5x5）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

高斯滤波平滑后的test2.tif（7x7）

50

100

150

200

250

300

350

400

450

500

[***********][1**********]

可见，对于高斯滤波器，模板的大小对滤波效果影响不大。高斯滤波虽然能够在一定程度上去掉噪声，但也使得图象变得模糊不清，效果并不能令人满意。

2、利用高通滤波器滤波测试图像test3,4：包括unsharp masking, Sobel edge detector, and Laplace edge detection；Canny algorithm.

实验原理分析：

锐化滤波能减弱或消除图像中的低频率分量，但不影响高频率分量。因为低频分量对应图像中灰度值缓慢变化的区域，因而与图像的整体特性，如整体对比度和平均灰度值等有关。锐化滤波将这些分量滤去可使图像反差增加，边缘明显。在实际应用中，锐化滤波可用于增强被模糊的细节或者低对比度图像的目标边缘。

图像锐化的主要目的有两个：一是增强图像边缘，使模糊的图像变得更加清晰，颜色变得鲜明突出，图像的质量有所改善，产生更适合人眼观察和识别的图像；二是希望经过锐化处理后，目标物体的边缘鲜明，以便于提取目标的边缘、对图像进行分割、目标区域识别、区域形状提取等，为进一步的图像理解与分析奠定基础。

由于锐化使噪声受到比信号还要强的增强，所以要求锐化处理的图像有较高的信噪比；否则，锐化后图像的信噪比更低。

①反锐化掩膜图像增强（unsharp masking）

图像的反锐化掩蔽算法可以表示为：

(1) 其中fs(x,y)表示经过反锐化掩蔽得到的锐化图像，是f(x,y)的模糊形式。反锐化掩蔽进一步的普遍形式称为高提升滤波。在图像中任何一点(x,y)处，高提升滤波后的图像fhb可定义如下：

(2) 其中A≥1，与前式一样，是f的模糊形式，此式也可以写成：

(3)

结合式(1)，可以得到：

这一表达式可计算高提升滤波图像。

如果选择拉普拉斯变换，式(4)变成： (4)

(5) 高提升滤波处理可以通过任何一个图1所示的掩模得以实现。当A=1时，高提升滤波处理就是标准的拉普拉斯变换。随着A超过l不断增大，锐化处理的效果越来越不明显。最终，当A足够大时，高提升图像将近似等于经常数调制的图像。

图1 高频提升滤波技术可以用其中一种掩膜来实现（A1） 本文采用的反锐化掩膜滤波器中A=1。根据以上分析，利用MATLAB软件设计反锐化掩膜滤波器对test3和test4进行滤波，结果如下图：

原图像——test3c

orrupt.pgm

反锐化掩膜后的test3corrupt.pgm

原图像

——test4 copy.bmp反锐化掩膜后的test4 copy.bmp

可见，经过反锐化掩膜滤波后，图像的边缘得到了增强，细节更加明显，但

同时也带来了高频噪声的影响。

②索贝尔边缘检测（Sobel edge detector）

索贝尔算子（Sobel operater）主要用作边缘检测，在技术上，它是一离散性差分算子，用来运算图像亮度函数的灰度之近似值。在图像的任何一点使用此算子，将会产生对应的灰度矢量或其法矢量。

该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像灰度值，其公式如下：

1Gx2

1

000

12*A1

1

Gy0

1

21



00A*

21

图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的

大小。

G

然后可用以下公式计算梯度方向。

arctan(

GyGx

)

如果角度Θ等于零，即代表图像在该处拥有纵向边缘，左方较右方暗。

根据以上算法分析，利用MATLAB软件设计Sobel边缘检测滤波器，对源

图像test3和test4进行滤波，结果如下图示：

源图像——test3c

orrupt.pgmSobel 边缘检测的test3c

orrupt.pgm

源图像——test3c

orrupt.pgm

Sobel 边缘检测——test4 copy.bmp

③Laplacian边缘检测（Laplacian edge detector）

拉普拉斯算子是最简单的各向同性微分算子，具有旋转不变性。一个二维图像函数

的拉普拉斯变换是各向同性的二阶导数，定义为：

（6）

为了更适合于数字图像处理，将该方程表示为离散形式：

（7）

另外，拉普拉斯算子还可以表示成模板的形式，如图2所示。图2（a）表示离散拉普拉斯算子的模板，图2（b）表示其扩展模板，图2（c）和（d）则分别表示其他两种拉普拉斯的实现模板。从模板形式容易看出，如果在图像中一个较暗的区域中出现了一个亮点，那么用拉普拉斯运算就会使这个亮点变得更亮。因为图像中的边缘就是那些灰度发生跳变的区域，所以拉普拉斯锐化模板在边缘检测中很有用。

一般增强技术对于陡峭的边缘和缓慢变化的边缘很难确定其边缘线

的位置。但该算子却可用二次微分正峰和负峰之间的过零点来确定，对孤立点或端点更为敏感，因此特别适用于以突出图像中的孤立点、孤立线或线端点为目的的场合。同梯度算子一样，拉普拉斯算子也会增强图像中的噪声，有时用拉普拉

（a（b

图2 拉普拉斯的4种模板

拉普拉斯锐化的基本方法可以由下式表示：

这种简单的锐化方法既可以产生拉普拉斯锐化处理的效果，同时又能保留背景信息，将原始图像叠加到拉普拉斯变换的处理结果中去，可以使图像中的各灰度值得到保留，使灰度突变处的对比度得到增强，最终结果是在保留图像背景的前提下，突现出图像中小的细节信息。

根据以上算法分析，利用MATLAB软件设计Laplacian边缘检测滤波器，

对源图像test3和test4进行滤波，结果如下图示：

源图像—test3corrupt.pgm

Laplacian边缘检测-test3corrupt.pgm

源图像

—test4 copy.bmpLaplacian边缘检测（未滤波）

—test4 copy.bmp

由上图可知，运用Laplacian算子对test4 copy.bmp进行边缘检测后，图像的边缘信息没有被很好的检测出来。这是由于拉普拉斯梯算子会增强图像中的噪声，因此本文再用拉普拉斯算子对test4 copy.bmp进行边缘检测时，先将图像进行高斯平滑滤波处理。处理结果如下图：

源图像

—test4 copy.bmp

Laplacian边缘检测(高斯滤波后)—test4 copy.bmp

④Canny算子边缘检测（Canny algorithm）

在图像边缘检测中，抑制噪声和边缘精确定位是无法同时满足的。边缘检测算法通过平滑滤波去除图像噪声的同时，也增加了边缘定位的不确定性；反之，提高边缘检测算子对边缘敏感性的同时，也提高了对噪声的敏感性。Canny算子力图在抗噪声干扰和精确定位边缘之间寻求最佳折中方案。用Canny算子检测图像边缘的步骤如下：

step1:用高斯滤波器平滑图象；

step2:用一阶偏导的有限差分来计算滤波后图像梯度的幅值和方向； step3:对梯度幅值进行非极大值抑制，其过程为找出图像梯度中的局部极大值点，把其他非局部极大值点置零以得到细化的边缘。

step4:用双阈值算法检测和连接边缘。使用两个阈值T1和T2（T1>T2），T1用来找到每条线段，T2用来在这些线段的两个方向上延伸寻找边缘的断裂处，并连接这些边缘。 利用上述原理和MATLAB软件设计滤波器，对test3和test4进行边缘检测，结果如下图所示：

源图像-test3c

orrupt.pgmCanny边缘检测-test3c

orrupt.pgm

Canny边缘检测-test4 copy.bmp

可见，经过Canny算子较好的检测出了图像的所有边缘信息，具有较强的去噪能力。

附录

一、参考文献

[1] 冈萨雷斯著.数字图像处理(第三版).北京：电子工业出版社，2010

[2] 杨杰 李庆著.数字图像处理及MATLAB实现——学习与实验指导.北京：电子

工业出版社，2010

[3] 苏金明 王永利著.MATLAB图形图像. 北京：电子工业出版社，2005

[4] 朱习军 隋思涟等著.MATLAB在信号与图像处理中的应用. 北京：电子工业

出版社，2009

[5] 张秀兰著.基于MATLAB的数字图像的边缘检测.吉林化工学院学报，2010 [6] 杨先花 黎粤华著.基于MATLAB图像边缘检测算法效果对比.机电产品开发

与创新，2010

二、源代码：

第一题

(1) 平滑滤波器

1. 中值滤波器（以3x3的模板大小为例）

①I=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test1.pgm','pgm'); n=3; a=ones(n,n); p=size(I);

x1=double(I);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)); e=c(1,:); for u=2:n

e=[e,c(u,:)]; end

mm=median(e);

x2(i+(n-1)/2,j+(n-1)/2)=mm; end end

I2=uint8(x2); imshow(I2)

title('中值滤波后的test1.pgm（3x3）')

②I=imread('E:\大三下\图像处理英文课件\作业\第四次作业\test2.tif','tif'); n=3;

a=ones(n,n); p=size(I);

x1=double(I);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)); e=c(1,:); for u=2:n

e=[e,c(u,:)]; end

mm=median(e);

x2(i+(n-1)/2,j+(n-1)/2)=mm; end end

I2=uint8(x2); imshow(I2)

title('中值滤波后的test2.tif（3x3）')

2、高斯滤波器（以3x3的模板大小为例）

① n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0;

[I,map]=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test1.pgm','pgm'); r=[cos(theta) -sin(theta); sin(theta) cos(theta)];

for i = 1 : n2 for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2

*sigma2^2))/(sigma2*sqrt(2*pi));

end

end

h = h / sqrt(sum(sum(h.*h))); f1=conv2(I,h,'same'); figure(1);

imagesc(I);title('test1.pgm'); colormap(gray); figure(2);

imagesc(f1);title(' '高斯滤波平滑后的test1.pgm(3x3)'); colormap(gray);

② n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0;

[I,map]=imread(' E:\大三下\图像处理英文课件\作业\第四次作业\test2.tif','tif');

r=[cos(theta) -sin(theta); sin(theta) cos(theta)]; for i = 1 : n2 for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2 end end

*sigma2^2))/(sigma2*sqrt(2*pi));

h = h / sqrt(sum(sum(h.*h))); f1=conv2(I,h,'same'); figure(1);

imagesc(I);title('test2.tif'); colormap(gray);

figure(2);

imagesc(f1);title(' '高斯滤波平滑后的test2.tif(3x3)'); colormap(gray);

（2）锐化滤波器：

①反锐化掩膜（以处理test3_corrupt.pgm为例）

pic = imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); I=pic;

picSize = size(pic);

pic = cast(pic,'int32');

t = zeros(picSize(1) + 2 , picSize(2) + 2); t = cast(t,'int32');

t(2:picSize(1) + 1 , 2:picSize(2) + 1) = pic; t(: , 1) = t(: , 2);

t(: , picSize(2) + 2) = t(: , picSize(2) + 1); t(1 , :) = t(2 , :);

t(picSize(1) + 2 , :) = t(picSize(1) + 1 , :); A = 1.0;

for i=2:1:picSize(1)+1

for j=2:1:picSize(2)+1

pic(i-1,j-1) = t(i,j)*A - (t(i,j)*(-8) + t(i,j-1) + t(i-1,j-1) + t(i-1,j)

+ t(i,j+1) + t(i+1,j) + t(i+1,j+1) + t(i+1,j-1) + t(i-1,j+1));

end

end

pic = cast(pic,'uint8');

imwrite(pic,['UnsharpMasking_A_', num2Str(10*A),'.png'],'png'); figure; subplot(1,2,1) imshow(I);

title('原图像——test4 copy.bmp'); subplot(1,2,2) imshow(pic);

title('反锐化掩膜后的test4 copy.bmp');

②Sobel 边缘检测（以处理test4_corrupt.pgm为例）

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp');

[N,M]=size(I) I=double(I);

h1=[-1,0,1;-2,0,2;-1,0,1]; h2=[-1,-2,-1;0,0,0;1,2,1]; Gx=conv2(I,h1,'same'); Gy=conv2(I,h2,'same'); F=abs(Gx)+abs(Gy); for i=1:N for j=1:M

I(i,j)=F(i,j); end end

I=uint8(I); figure; imshow(I);

title('Sobel 边缘检测——test4 copy.bmp');

③Laplace边缘检测

a）Laplace边缘检测——test3_corrupt.pgm

I = imread('E:\大三下\图像处理英文课件\作业\第四次作业\test3_corrupt.pgm','pgm'); [N,M]=size(I); figure;

subplot(1,2,1); imshow(I);

title('ԴͼÏñ¡ªtest3_corrupt.pgm'); t=ones(N,M); t=I;

for i=2:N-1

for j=2:M-1

I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j); end

end

subplot(1,2,2);

imshow(I);

title('Laplacian边缘检测-test3_corrupt.pgm');

b）Laplace边缘检测——test4 copy.bmp

未进行高斯滤波：

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); [N,M]=size(I); figure;

subplot(1,2,1);

imshow(I);

title('源图像—test4 copy.bmp');

t=ones(N,M);

t=I;

for i=2:N-1

for j=2:M-1

I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j);

end

end

subplot(1,2,2);

imshow(I);

title('Laplacian边缘检测(未滤波)-test4 copy.bmp');

加入高斯滤波：

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp');

[N,M]=size(I);

figure;

subplot(1,2,1);

imshow(I);

title('源图像—test4 copy.bmp ');

n1=3;sigma1=0.5;n2=3;sigma2=0.5;theta=0;

r=[cos(theta) -sin(theta); sin(theta) cos(theta)];

for i = 1 : n2

for j = 1 : n1

u = r*[j-(n1+1)/2 i-(n2+1)/2]';

h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2*s end

end

h = h / sqrt(sum(sum(h.*h)));

f1=conv2(I,h,'same');

t=ones(N,M);

t=f1;

for i=2:N-1

for j=2:M-1

f1(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j);

end

end

subplot(1,2,2);

imshow(f1);

title(' Laplacian边缘检测(高斯滤波后)-test4 copy.bmp ');

igma2^2))/(sigma2*sqrt(2*pi));

④ Canny算子边缘检测（以处理test4 copy.bmp为例）

I= imread('E:\大三下\图像处理英文课件\作业\第四次作业\test4 copy.bmp','bmp'); a = im2single(I);

[m,n] = size(a);

e = false(m,n);

GaussianDieOff = 0.0001;

PercentOfPixelsNotEdges = 0.7;

ThresholdRatio = 0.4;

sigma = 1;

thresh = [];

pw = 1:30;

ssq = sigma^2;

width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');

if isempty(width)

width = 1;

end

t = (-width:width);

gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);

[x,y]=meshgrid(-width:width,-width:width);

dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);

aSmooth=imfilter(a,gau,'conv','replicate');

aSmooth=imfilter(aSmooth,gau','conv','replicate');

ax = imfilter(aSmooth, dgau2D, 'conv','replicate');

ay = imfilter(aSmooth, dgau2D', 'conv','replicate');

mag = sqrt((ax.*ax) + (ay.*ay));

magmax = max(mag(:));

if magmax>0

mag = mag / magmax;

end

if isempty(thresh)

counts=imhist(mag, 64);

highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,... 1,'first') / 64;

lowThresh = ThresholdRatio*highThresh;

thresh = [lowThresh highThresh];

elseif length(thresh)==1

highThresh = thresh;

if thresh>=1

eid = sprintf('Images:%s:thresholdMustBeLessThanOne', mfilename); msg = 'The threshold must be less than 1.';

error(eid,'%s',msg);

end

lowThresh = ThresholdRatio*thresh;

thresh = [lowThresh highThresh];

elseif length(thresh)==2

lowThresh = thresh(1);

highThresh = thresh(2);

if (lowThresh >= highThresh) || (highThresh >= 1)

eid = sprintf('Images:%s:thresholdOutOfRange', mfilename); msg = 'Thresh must be [low high], where low

end

end

idxStrong = [];

for dir = 1:4

[m,n] = size(mag);

idx = find((ay-ay) | (ay>=0 & ax

if ~isempty(idx)

v = mod(idx,m);

extIdx = find(v==1 | v==0 | idx(n-1)*m));

idx(extIdx) = [];

end

ixv = ax(idx);

iyv = ay(idx);

gradmag = mag(idx);

d = abs(iyv./ixv);

gradmag1 = mag(idx+m).*(1-d) + mag(idx+m-1).*d;

gradmag2 = mag(idx-m).*(1-d) + mag(idx-m+1).*d;

idxLocMax = idx(gradmag>=gradmag1 & gradmag>=gradmag2); idxWeak = idxLocMax(mag(idxLocMax) > lowThresh);

e(idxWeak)=1;

idxStrong = [idxStrong; idxWeak(mag(idxWeak) > highThresh)]; end

if ~isempty(idxStrong)

rstrong = rem(idxStrong-1, m)+1;

cstrong = floor((idxStrong-1)/m)+1;

e = bwselect(e, cstrong, rstrong, 8);

e = bwmorph(e, 'thin', 1);

end

imshow(e);

title('Canny边缘检测-test4 copy.bmp');