图像边缘检测各种算子MATLAB实现以及实际应用

1、图像处理中的数学方法实验报告学生姓名：赵芳舟教师姓名：曾理学院：数学与统计学院专业：信息与计算科学学号：联系方式：梯度和拉普拉斯算子在图像边缘检测中的应用1、 数学方法边缘检测最通用的方法是检测灰度值的不连续性，这种不连续性用一阶和二阶导数来检测。1. （1）一阶导数：一阶导数即为梯度，对于平面上的图像来说，我们只需用到二维函数的梯度，即： f=gxgy=fxfy,该向量的幅值：f=magf=gx2+gy21/2=f/x2+f/y21/2,为简化计算，省略上式平方根，得到近似值fgx2+gy2；或通过取绝对值来近似，得到：fgx+gy。（2）二阶导数：二阶导数通常用拉普拉斯算子来计算，由二阶微

2、分构成：2fx,y=2f(x,y)x2+2f(x,y)y22. 边缘检测的基本思想：（1） 寻找灰度的一阶导数的幅度大于某个指定阈值的位置；（2） 寻找灰度的二阶导数有零交叉的位置。3. 几种方法简介（1） Sobel边缘检测器：以差分来代替一阶导数。Sobel边缘检测器使用一个33邻域的行和列之间的离散差来计算梯度，其中，每行或每列的中心像素用2来加权，以提供平滑效果。f=gx2+gy21/2=z7+2z8+z9-(z1+2z2+z3)2+z3+2z6+z9-z1+2z4+z721/2-1-21000121 -101-202-101 （2） Prewitt边缘检测器:使用下图所示模板来数字化

3、地近似一阶导数。与Sobel检测器相比，计算上简单一些，但产生的结果中噪声可能会稍微大一些。gx=z7+z8+z9-z1+z2+z3gy=z3+z6+z9-z1-z4-z7-1-1-1000111-101-101-101（3） Roberts边缘检测器：使用下图所示模板来数字化地将一阶导数近似为相邻像素之间的差，它与前述检测器相比功能有限（非对称，且不能检测多种45倍数的边缘）。gx=z9-z5gy=z8-z6-10010-110（4） Laplace边缘检测器：二维函数f(x,y)的拉普拉斯是一个二阶的微分定义：2fx,y=2f(x,y)x2+2f(x,y)y2模板算子可分为四邻域和八邻域，

4、如下：0101-41010（四邻域）1111-81111（八邻域）（5） LoG边缘检测器由于噪声点（灰度与周围点相差很大的像素点）对边缘检测有一定的影响，所以效果更好的是LoG算子，即Laplacian-Guass算子。引入高斯函数来平滑噪声：该函数的Laplace算子：它把Guass平滑滤波器和Laplace锐化滤波器结合起来，先平滑掉噪声，再进行边缘检测，所以效果比单用Laplace算子要更为平滑，效果更好。（6） Canny边缘检测器主要分为以下几个步骤：使用具有指定标准差的一个高斯滤波器来平滑图像，以减少噪声；在每个点处计算局部梯度和边缘方向；对步骤中确定的边缘点产生梯度中的脊线顶部

5、进行追踪，并将实际山不在脊线顶部的像素设置为零，从而在输出中给出一条细线（非最大值抑制），然后使用滞后阈值处理法对这些脊线像素进行阈值处理。最后进行边缘连接。2、 实验结果原始图像：1.Roberts算子2.Prewitt算子3.Sobel算子4.Laplace算子（4邻域）（8邻域）5. LoG四种结果的比较：（Laplace采用8邻域模板）图像来源：/show/.html（网行天下首页设计图库文化艺术插画集：牧羊犬灰度）图像规格：800677 格式：SVG 模式：RGB硬件条件：软件条件：运行环境：Matlab r2015b三、讨论1. 各方法

6、优缺点及比较：（1） Roberts算子：Roberts算子是22算子，对具有陡峭的低噪声图像响应最好，并且检测垂直边缘的效果好于斜向边缘，定位精度高。然而，它对噪声敏感，无法抑制噪声的影响。因此，它适用于边缘明显且噪声较少的图像分割。（2） Prewitt算子：Prewitt算子将两个点的各自一定领域内的灰度值求和，并根据两个灰度值和的差来计算x,y的偏导数。它是平均滤波，对噪声有抑制作用，对于灰度渐变的低噪声图像有较好的检测效果，但是像素平均相当于对图像的低通滤波，所以它对边缘的定位不如Roberts算子。对于混合多复杂噪声的图像，效果不太理想。（3） Sobel算子：Sobel算子是滤波

7、算子的形式，用于提取边缘，可以利用快速卷积函数， 简单有效，因此应用广泛。美中不足的是，Sobel算子并没有将图像的主体与背景严格地区分开来，即Sobel算子没有严格地模拟人的视觉生理特征，所以提取的图像轮廓有时并不能令人满意。（4） Laplace算子：Laplace算子是一种各向同性算子，在只关心边缘的位置而不考虑其周围的象素灰度差值时比较合适。Laplace算子对孤立象素的响应要比对边缘或线的响应要更强烈，因此只适用于无噪声图象。存在噪声情况下，使用Laplacian算子检测边缘之前需要先进行低通滤波。（5） LoG算子：由于二阶导数是线性运算，利用LoG卷积一幅图像与首先使用高斯型平滑

8、函数卷积改图像，然后计算所得结果的拉普拉斯是一样的。所以在LoG公式中使用高斯函数的目的就是对图像进行平滑处理，使用Laplacian算子的目的是提供一幅用零交叉确定边缘位置的图像；图像的平滑处理减少了噪声的影响并且它的主要作用还是抵消由Laplacian算子的二阶导数引起的逐渐增加的噪声影响。（6） Canny算子：是一个具有滤波，增强和检测的多阶段的优化算子，在进行处理前，Canny算子先利用高斯平滑滤波器来平滑图像以消除噪声，Canny分割算法采用一阶偏导的有限差分来计算梯度的幅值和方向。在处理过程中，Canny算法还将经过一个非极大值抑制的过程。最后Canny算法将采用两个阈值来连接边

9、缘。高定位精度、低误判率、抑制虚假边缘，适用于高噪声图像。四、实际应用1、概述：图像边缘检测广泛应用于车牌识别，人脸识别，地震带检测，生物医疗以及产品外观检测等方面。为了体现其实际应用，在此针对车牌识别来进行实验，实验流程如下：确定区域文字分割模板配对均值滤波边缘检测二值化灰度处理3、 实验结果原图像：灰度处理：将彩色图像转化为灰度图像，缩小图片占用的空间，减少处理时间边缘检测：车牌区域的边缘信息一般比其他部分丰富且集中，只要 背景不太复杂，没有过多干扰，边缘特征可将车牌区域与其他部分区别开来，从而消除图像中的无关信息。腐蚀和膨胀：腐蚀可以分割独立的图像元素，膨胀用于连接相邻的元素，更加利于图

10、像的分割最后将车牌号分割成七个单个字符，建立模板库，与其进行配对，最终得出结果。五、程序介绍（1） Roberts边缘检测 clear;sourcePic=imread(C:Users34899Desktoptimg.jpg); %读取原图像grayPic=mat2gray(sourcePic); %实现图像矩阵的归一化操作m,n=size(grayPic);newGrayPic=grayPic;%为保留图像的边缘一个像素robertsNum=0; %经roberts算子计算得到的每个像素的值robertThreshold=0.2; %设定阈值for j=1:m-1 %进行边界提取 for k

11、=1:n-1 robertsNum = abs(grayPic(j,k)-grayPic(j+1,k+1) + abs(grayPic(j+1,k)-grayPic(j,k+1); if(robertsNum robertThreshold) newGrayPic(j,k)=255; else newGrayPic(j,k)=0; end endendsubplot(1,2,1);imshow(sourcePic);title(原图); subplot(1,2,2);imshow(newGrayPic);title(Robert算子处理后图像); （2） Prewitt边缘检测 clear;s

12、ourcePic=imread(C:Users34899Desktoptimg.jpg); %读取原图像grayPic=mat2gray(sourcePic);%实现图像矩阵的归一化操作m,n=size(grayPic);newGrayPic=grayPic;%为保留图像的边缘一个像素PrewittNum=0;%经Prewitt算子计算得到的每个像素的值PrewittThreshold=0.5;%设定阈值for j=2:m-1 %进行边界提取for k=2:n-1 PrewittNum=abs(grayPic(j-1,k+1)-grayPic(j+1,k+1)+grayPic(j-1,k)-g

13、rayPic(j+1,k)+grayPic(j-1,k-1)-grayPic(j+1,k-1)+abs(grayPic(j-1,k+1)+grayPic(j,k+1)+grayPic(j+1,k+1)-grayPic(j-1,k-1)-grayPic(j,k-1)-grayPic(j+1,k-1); if(PrewittNum PrewittThreshold) newGrayPic(j,k)=255; else newGrayPic(j,k)=0; end endendsubplot(1,2,1);imshow(sourcePic);title(原图); subplot(1,2,2);ims

14、how(newGrayPic);title(Prewitt算子处理后图像); （3） Sobel边缘检测 clear all; close all; imag = imread(C:Users34899Desktoptimg.jpg); %读取关键帧 imag = rgb2gray(imag); %转化为灰度图 subplot(1,2,1);imshow(imag);title(原图); high,width = size(imag); % 获得图像的高度和宽度 F2 = double(imag); U = double(imag); uSobel = imag; for i = 2:high

15、 - 1 %sobel边缘检测 for j = 2:width - 1 Gx = (U(i+1,j-1) + 2*U(i+1,j) + F2(i+1,j+1) - (U(i-1,j-1) + 2*U(i-1,j) + F2(i-1,j+1); Gy = (U(i-1,j+1) + 2*U(i,j+1) + F2(i+1,j+1) - (U(i-1,j-1) + 2*U(i,j-1) + F2(i+1,j-1); uSobel(i,j) = sqrt(Gx2 + Gy2); end end subplot(1,2,2);imshow(im2uint8(uSobel):title(Sobel边缘检

16、测后); （4） Laplace边缘检测主函数：f=imread(C:Users34899Desktoptimg.jpg)fGray=rgb2gray(f);figure()imshow(fGray),title(灰度图像);fGray=double(fGray);T=60;LapModType=8;%设置laplace模板方式fLapEdge=LaplaceEdge(fGray,LapModType,T);fGrayLapEdge=uint8(fLapEdge);figure()imshow(fLapEdge),title(laplace边缘图像);四邻域算子和八邻域算子function f

17、LapEdge=LaplaceEdge(fGray,MoldType,Thresh)%-参数介绍-%输入参数：% fGray;输入的灰度图像% MoldType:模板类型，包括四邻域和八邻域% Thresh:边缘检测阈值%输出参数：% fEdge:边缘像素点，存储的是二值化图像r,c=size(fGray);fLapEdge=zeros(r,c); %四邻域拉普拉斯边缘检测算子if 4=MoldType for i=2:r-1 for j=2:c-1 Temp=-4*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1);

18、if TempThresh fLapEdge(i,j)=255; else fLapEdge(i,j)=0; end end endend%八邻域拉普拉斯边缘检测算子if 8=MoldType for i=2:r-1 for j=2:c-1 Temp=-8*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1)+fGray(i-1,j-1)+fGray(i+1,j+1)+fGray(i+1,j-1)+fGray(i-1,j+1); if TempThresh fLapEdge(i,j)=255; else fLapEdge(

19、i,j)=0; end end endEnd（5） 拉普拉斯高斯检测器clc; clear all; close all;hsize = 5 5;sigma = 0.5;h = fspecial(log, hsize, sigma);I = imread(C:Users34899Desktoptimg.jpg);bw = imfilter(I, h, replicate);figure;subplot(1, 2, 1); imshow(I, ); title(原图像, FontWeight, Bold);subplot(1, 2, 2); imshow(bw, ); title(边缘图像, F

20、ontWeight, Bold);（6） Canny算子tic % 记录CPU使用时间，开始计时I = imread(C:Users34899Desktoptimg.jpg); % 读入图像I=rgb2gray(I); % 转化为灰色图像% figure,imshow(I);title(原图) % 显示原图像% a = im2single(I); % 将图像矩阵转换为单精度类型% a = im2uint8(a); % 将图像矩阵转换为8位无符号整数类型% I = uint16(I); % 8位无符号整数转化为16位无符号整数 %* matlab的高斯平滑 *row,line = size(I)

21、; % 行列temp = zeros(row+8,line+8); % 申请空间imagedata = zeros(row,line); % 申请空间% 高斯核计算,核大小为9*9GaussianDieOff = .0001;sigma = 1;pw = 1:30; % possible widthsssq = sigma2;width = find(exp(-(pw.*pw)/(2*ssq)GaussianDieOff,1,last);if isempty(width) width = 1; % the user entered a really small sigmaendt = (-wi

22、dth:width);gau = exp(-(t.*t)/(2*ssq)/(2*pi*ssq); kernel = gau * gau; % 高斯核,大小为9*9kernel = im2single(kernel); % 图像边缘复制for q = 5:row+4 % 首先载入源图像信息 for p = 5:line+4 temp(q,p) = I(q-4,p-4); endendfor q = 1:row+8 % 复制列5:line+4 for p = 5:line+4 if(qrow+4) temp(q,p) = temp(row+4,p); end endendfor q = 1:row

23、+8 % 复制行1:row+8 for p = 1:line+8 if(pline+4) temp(q,p) = temp(q,line+4); end endendfor q = 1:row % 高斯卷积 for p = 1:line for n = 1:9 for m = 1:9 imagedata(q,p) = imagedata(q,p) + kernel(n,m) * temp(q+n-1,p+m-1); end end endend%* matlab的高斯平滑说明结束 * %* matlab求梯度 *dx = zeros(size(I); % 申请空间,存放x方向的幅值dy = z

24、eros(size(I); % 申请空间,存放y方向的幅值x,y=meshgrid(-width:width,-width:width);dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq)/(pi*ssq); % 核% 图像边缘复制for q = 5:row+4 % 首先载入源图像信息 for p = 5:line+4 temp(q,p) = imagedata(q-4,p-4); endendfor q = 1:row+8 % 复制列5:line+4 for p = 5:line+4 if(qrow+4) temp(q,p) = temp(row+4,p); end en

25、dendfor q = 1:row+8 % 复制行1:row+8 for p = 1:line+8 if(pline+4) temp(q,p) = temp(q,line+4); end endendfor q = 1:row % x方向卷积 for p = 1:line for n = 1:9 for m = 1:9 dx(q,p) = dx(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1); end end endenddgau2D = dgau2D;for q = 1:row % y方向卷积 for p = 1:line for n = 1:9 for m =

26、 1:9 dy(q,p) = dy(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1); end end endend%* matlab求梯度结束 * %* 自适应阈值计算 *value = sqrt(dx.2 + dy.2); % 幅值value = double(value);value = value/max(max(value);% 归一化 temp_value = value;counts=imhist(value, 64);high_threshold = find(cumsum(counts) 0.7*row*line,1,first) / 64;low

27、_threshold = 0.4 * high_threshold; % 低阈值clear counts%* 自适应阈值计算结束 * %* 方向与阈值判别 *num = 0; % 当前堆栈个数flag = zeros(80000,2); % 堆栈temp_flag = zeros(80000,2); % 临时堆栈imagedata = zeros(row,line); % 初始化% 方向：% 0为045% 1为4590% 2为90135% 3为135180direction = zeros(size(I); % 申请空间,存放小于高阈值而大于低阈值for q = 2:row-1 for p =

28、 2:line-1 if(dy(q,p)-dy(q,p) | (dy(q,p)=0 & dx(q,p)0 & -dy(q,p)=dx(q,p) | (dx(q,p)0 & -dy(q,p)=dx(q,p) d = abs(dx(q,p)/dy(q,p); % 4590方向 gradmag = temp_value(q,p); gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p+1)*d; gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p-1)*d; elseif(dx(q,p)dy

29、(q,p) | (dx(q,p)=0 & dx(q,p)dy(q,p) d = abs(dx(q,p)/dy(q,p); % 90135方向 gradmag = temp_value(q,p); gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p-1)*d; gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p+1)*d; elseif(dy(q,p)0 & dx(q,p)0 & dx(q,p)=dy(q,p) d = abs(dy(q,p)/dx(q,p); % 135180方向 g

30、radmag = temp_value(q,p); gradmag1 = temp_value(q,p-1)*(1-d) + temp_value(q-1,p-1)*d; gradmag2 = temp_value(q,p+1)*(1-d) + temp_value(q+1,p+1)*d; end if(gradmag=gradmag1 & gradmag=gradmag2) if(gradmag = high_threshold) % 高阈值判别 value(q,p) = 255; elseif(gradmag = low_threshold)% 低阈值判别 value(q,p) = 125

31、; else value(q,p) = 0; % 小于低阈值 end else value(q,p) = 0; % 非极大值抑制 end end endfor q = 2:row-1 % 高阈值判别，检查高阈值邻域8个方向范围内是否存在低阈值 for p = 2:line-1 if(value(q,p) = 255) imagedata(q,p) = 255; if(value(q-1,p-1) = 125) value(q-1,p-1) = 255; imagedata(q-1,p-1) = 255; if(q-1 1) & (p-1 1) num = num + 1; flag(num,1

32、) = q-1; flag(num,2) = p-1; end end if(value(q-1,p) = 125) value(q-1,p) = 255; imagedata(q-1,p) = 255; if(q-1 1) num = num + 1; flag(num,1) = q-1; flag(num,2) = p; end end if(value(q-1,p+1) = 125) value(q-1,p+1) = 255; imagedata(q-1,p+1) = 255; if(q-1 1) & (p+1 1) num = num + 1; flag(num,1) = q; fla

33、g(num,2) = p-1; end end if(value(q,p+1) = 125) value(q,p+1) = 255; imagedata(q,p+1) = 255; if(p+1 line) num = num + 1; flag(num,1) = q; flag(num,2) = p+1; end end if(value(q+1,p-1) = 125) value(q+1,p-1) = 255; imagedata(q+1,p-1) = 255; if(q+1 1) num = num + 1; flag(num,1) = q+1; flag(num,2) = p-1; e

34、nd end if(value(q+1,p) = 125) value(q+1,p) = 255; imagedata(q+1,p) = 255; if(q+1 row) num = num + 1; flag(num,1) = q+1; flag(num,2) = p; end end if(value(q+1,p+1) = 125) value(q+1,p+1) = 255; imagedata(q+1,p+1) = 255; if(q+1 row) & (p+1 1) & (p-1 1) num = num + 1; temp_flag(num,1) = q-1; temp_flag(n

35、um,2) = p-1; end end if(value(q-1,p) = 125) value(q-1,p) = 255; imagedata(q-1,p) = 255; if(q-1 1) num = num + 1; temp_flag(num,1) = q-1; temp_flag(num,2) = p; end end if(value(q-1,p+1) = 125) value(q-1,p+1) = 255; imagedata(q-1,p+1) = 255; if(q-1 1) & (p+1 1) num = num + 1; temp_flag(num,1) = q; temp_flag(num,2) = p-1; end end if(value(q,p+1) = 125) value(q,p+1) = 255; imagedata(q,p+1) = 255; if(p+1 line) num = num + 1; temp_flag(num,1) = q; temp_flag(num,2) = p+1;