图像边缘提取算法源程序.doc

图像边缘提取算法研究报告概述图像的边缘包含了图像最重要的信息。什么是边缘？一般是指图像灰度变化率最大的位置。从成因上看，一般图像边缘主要由四个方面的因素形成：(1)图像灰度在表面法向变化的不连续造成的边缘；(2)图像对像素在空间上不一致形成的边缘；(3)在光滑的表面上由于颜色的不一致形成的边缘：(4)物体的光影造成的边缘。图像边缘提取的作用有：(1)改良图像质量；(2)分离对象；(3)理解和重构视觉场景；(4)识别特征；(5)其他。图像边缘检测是图像处理与计算机视觉共同的基本课题，1960年以来，相继发展了一系列采用梯度算子和拉普拉斯算子的边缘检测技术；为了降低图像噪声对边缘检测算法的干扰，1980年以来，又建立了高斯低通滤波与拉普拉斯算子复合的过零点检测Marr-Hildreth理论；在另一个方向上，1980年代初期，Canny从信号处理的角度出发，使边缘检测算法更具有实用性。本报告主要介绍以上以上几个方面的内容，通过matlab程序实现以上几种算法，对比各种算法的性能。算法介绍及相应程序一、 基于微分算子的边缘检测检测图像边缘信息，可以把图像看做曲面，边缘就是图像的变化最剧烈的位置。这里所讲的边缘信息包含两个方面：一是边缘的具体位置，即像素的坐标；而是边缘的方向。微分算子有两个重要性质：定域性(或局部性)、敏感性(或无界性)。敏感性就是说，它对局部的函数值变化很敏感，但是因其对变化过于敏感又有了天然的缺陷不能抵抗噪声。局部性意思是指，每一点的导数只与函数在该点邻近的信息有关。主要有两大类基于微分算子的边缘检测技术：一阶微分算子边缘检测与二阶微分算子边缘检测。这些检测技术采用以下的基本步骤：(1) 将相应的微分算子简化为离散的差分格式，进而简化为模板(记为T)。(2) 利用模板对图像f(m,n)进行运算，获得模板作用后的结果Tf(m,n)。(3) 提出阈值h,在采用一阶微分算子情形记录下高于某个阈值h的位置坐标(而采用二阶微分算子情形，一般是对某个阈值确立 )(4) 对集合进行整理，同时调整阈值h。Roberts算子Roberts算子是一种利用局部差分算子寻找边缘的算子，两个模板分别为 则，= = 算法的步骤为：(1) 首先用两个模板分别对图像作用得到和；(2) 对，进行阈值判决，若大于阈值则相应的点位于便于边缘处。 对于阈值选取的说明：由于微分算子的检测性能受阈值的影响较大，为此，针对具体图像我们采用以下阈值的选取方法，对处理后的图像统计大于某一阈值的点，对这些数据求平均值，以下每个程序均采用此方法，不再做说明。具体程序如下：%-filename:Roberts.m-%-Useage:edge detecting by the Roberts operator-%-Writer： Subailong- function edgeRb=Roberts(oimage);xlen ylen = size(oimage); %read the sizeedgeX = zeros(xlen,ylen); %horizontal directionedgeY = zeros(xlen,ylen); %vertical directionedgeXY = zeros(xlen,ylen); %synthesize the two directions %-process the oringin image with the operator-for i = 1:xlen-1 for j = 1:ylen-1 edgeX(i,j) = oimage(i,j) - oimage(i+1,j+1); edgeY(i,j) = oimage(i+1,j) - oimage(i,j+1); endendedgeX = abs(edgeX);edgeY = abs(edgeY);edgeXY = sqrt(edgeX.*edgeX + edgeY.*edgeY); %-Thresh estimate-rsum=0;counter=0;for i = 2:xlen-1 for j = 2:ylen-1 if(edgeXY(i,j)15) rsum=rsum+edgeXY(i,j); counter=counter+1; end endendthreshold=rsum/counter; %threshold = 2*sum(sum(edgeXY(1:xlen-1,1:ylen-1)/(xlen-1)*(ylen-1); %-edge detecting-for i = 1:xlen-1 for j = 1:ylen-1 if (edgeXY(i,j)threshold) edgeimage(i,j) = 255; else edgeimage(i,j) = 0; end endendedgeRb = edgeimage;Sobel算子Sobel算子采用中心差分，但对中间水平线和垂直线上的四个邻近点赋予略高的权重。两个模板分别如下： 该算法的处理过程同Roberts算字，在这个程序中，我们采用两个方向分别进行阈值判断，程序如下： %-filename:Sobel.m-%-Useage:edge detecting by the Sobel operator-%-Writer： Subailong- function edgeS=Sobel(oimage);xlen ylen = size(oimage); %read the sizeedgeX = zeros(xlen,ylen); %horizontal directionedgeY = zeros(xlen,ylen); %vertical directionedgeXY = zeros(xlen,ylen); %synthesize the two directionsd= zeros(xlen,ylen); %direction of the conresponding pointsedgeimage= zeros(xlen,ylen); %the result imageH1 = -1 0 1;-2 0 2;-1 0 1; %horizontal operatorH2 = -1 -2 -1;0 0 0;1 2 1; %vertical operator %-process the oringin image with the operator-for i = 1:xlen for j = 1:ylen for m = 1:3 for n = 1:3 updateX = i-m+2; updateY = j-n+2; if (updateX=1) & (updateX=1) & (updateYedgeY(i,j) d(i,j)=1; %save the direction else d(i,j)=2; %save the direction end end end end endend%edgeX = abs(edgeX);%edgeY = abs(edgeY);%edgeXY = sqrt(edgeX.*edgeX + edgeY.*edgeY); %-Thresh estimate-rsum=0;counter=0;for i = 2:xlen-1 for j = 2:ylen-1 if(edgeX(i,j)20) rsum=rsum+edgeX(i,j); counter=counter+1; end endendthreshold=rsum/counter; %threshold = 4*sum(sum(edgeXY(2:xlen-1,2:ylen-1)/(xlen-2)*(ylen-2); %-edge detecting-for i = 2:xlen-1 for j = 2:ylen-1 if (d(i,j)=1&edgeX(i,j)edgeX(i,j+1)&edgeX(i,j)edgeX(i,j-1)&edgeX(i,j)threshold) edgeimage(i,j) = 1; elseif (d(i,j)=2&edgeY(i,j)edgeY(i+1,j)&(edgeY(i,j)edgeY(i-1,j)&edgeY(i,j)threshold) edgeimage(i,j) = 1; end endendedgeS=edgeimage;Prewitt算子 Prewitt算子也属于中心差分类型，但没有给最邻近点较高的权重，两个模板如下： 该算法的处理过程同Roberts算字，在这个程序中，我们采用两个方向分别进行阈值判断，程序如下： %-filename:Prewitt.m-%-Useage:edge detecting by the Prewitt operator-%-Writer： Subailong- function edgeP=Prewitt(oimage);xlen ylen = size(oimage); %read the sizeedgeX = zeros(xlen,ylen); %horizontal directionedgeY = zeros(xlen,ylen); %vertical directionedgeXY = zeros(xlen,ylen); %synthesize the two directionsedgeimage = zeros(xlen,ylen); %the result imaged=zeros(xlen,ylen); %direction of the conresponding pointsH1 = -1 0 1;-1 0 1;-1 0 1; %horizontal operatorH2 = 1 1 1;0 0 0;-1 -1 -1; %vertical operator %-process the oringin image with the operator-for i = 1:xlen for j = 1:ylen for m = 1:3 for n = 1:3 updateX = i-m+2; updateY = j-n+2; if (updateX=1) & (updateX=1) & (updateYedgeY(i,j) d(i,j)=1; %save the direction else d(i,j)=2; %save the direction end end end end endend%edgeX = abs(edgeX);%edgeY = abs(edgeY);%edgeXY = edgeX.*edgeX + edgeY.*edgeY;%threshold = 4*sum(sum(edgeXY(2:xlen-1,2:ylen-1)/(xlen-2)*(ylen-2);%-%-Thresh estimate-rsum=0;counter=0;for i = 2:xlen-1 for j = 2:ylen-1 if(edgeX(i,j)15) rsum=rsum+edgeX(i,j); counter=counter+1; end endendthreshold=rsum/counter; %-edge detecting-for i = 2:xlen-1 for j = 2:ylen-1 if (d(i,j)=1&edgeX(i,j)edgeX(i,j+1)&edgeX(i,j)edgeX(i,j-1)&edgeX(i,j)threshold) edgeimage(i,j) = 1; elseif (d(i,j)=2&edgeY(i,j)edgeY(i+1,j)&(edgeY(i,j)edgeY(i-1,j)&edgeY(i,j)threshold) edgeimage(i,j) = 1; end endend edgeP=edgeimage;Krisch算子Krisch算子采用以下8种形式的模板： ，图像中的每个位置都要经过8个模板的作用，最大值被选作输出，达到最大值的方向就是边缘的方向。实际上可采用以下快速算法：具体程序如下：%-filename:Kirsch.m-%-Useage:edge detecting by the Kirsch operator-%-Writer： Subailong- function edgeK=Kirsch(oimage);xlen ylen = size(oimage); %read the sizeedge = zeros(xlen,ylen); %the result imaged = zeros(xlen,ylen); %direction of the conresponding points %define three voriables to computer fastK = zeros(8,1);t = zeros(8,1);f = zeros(8,1); %-process the oringin image with the operator-for i = 2:xlen-1 for j = 2:ylen-1 temp = sum(sum(oimage(i-1:i+1,j-1:j+1)-oimage(i,j); K(1) = oimage(i-1,j-1)+oimage(i-1,j)+oimage(i-1,j+1); %North K(2) = oimage(i-1,j)+oimage(i-1,j+1)+oimage(i,j+1); %Northeast K(3) = oimage(i-1,j+1)+oimage(i,j+1)+oimage(i+1,j+1); %East K(4) = oimage(i,j+1)+oimage(i+1,j+1)+oimage(i+1,j); %Southeast K(5) = oimage(i+1,j+1)+oimage(i+1,j)+oimage(i+1,j-1); %South K(6) = oimage(i+1,j)+oimage(i+1,j-1)+oimage(i,j-1); %Southweast K(7) = oimage(i+1,j-1)+oimage(i,j-1)+oimage(i-1,j-1); %Weast K(8) = oimage(i,j-1)+oimage(i-1,j-1)+oimage(i-1,j); %NorthWest for m = 1:8 t(m) = temp-K(m); end f=abs(5*K-3*t); m c=max(f); edge(i,j)=m; %save the maximum d(i,j)=c; %save the direction of maximum endend %-Thresh estimate-rsum=0;counter=0;for i = 2:xlen-1 for j = 2:ylen-1 if(edge(i,j)200) rsum=rsum+edge(i,j); counter=counter+1; end endendthreshold=rsum/counter; %-edge detecting- for i = 2:xlen-1 for j = 2:ylen-1 if (d(i,j)=1)&(edge(i,j)edge(i-1,j)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=2)&(edge(i,j)edge(i-1,j+1)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=3)&(edge(i,j)edge(i,j+1)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=4)&(edge(i,j)edge(i+1,j+1)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=5)&(edge(i,j)edge(i+1,j)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=6)&(edge(i,j)edge(i+1,j-1)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=7)&(edge(i,j)edge(i,j-1)&(edge(i,j)threshold) edgeimage(i,j) = 255; elseif (d(i,j)=8)&(edge(i,j)edge(i-1,j-1)&(edge(i,j)threshold) edgeimage(i,j) = 255; end endendedgeK=edgeimage;Robinson算子Robinson算子也采用8形式的模板，其具体形式参加课本（课本中第一个算子有误），这里需要说明一下Robinson算子的特征，每个四个其相应元素的负号相反，就是说在计算中只需计算四个就行了。在本程序中为了比较数值（而不是绝对值，这样更好地体现方向性）的大小，仍计算8个。具体程序如下：%-filename:Robinson.m-%-Useage:edge detecting by the Robinson operator-%-Writer： by Subailong- function edgeR=Robinson(oimage);xlen ylen = size(oimage); %read the sizeedgeRobinson = zeros(xlen,ylen); edgeRobinsonimage = zeros(xlen,ylen); %the result imaged=zeros(xlen,ylen); %direction of the conresponding points %-process the oringin image with the operator-R = zeros(4,1);for i = 2:xlen-1 for j = 2:ylen-1 R(1) = oimage(i-1,j-1)+oimage(i-1,j+1)-oimage(i+1,j+1)-oimage(i+1,j-1). +oimage(i-1,j)-oimage(i+1,j)+oimage(i-1,j)-oimage(i+1,j); %North R(2) = oimage(i-1,j)+oimage(i,j+1)-oimage(i+1,j)-oimage(i,j-1). +oimage(i-1,j+1)-oimage(i+1,j-1)+oimage(i-1,j+1)-oimage(i+1,j-1); %Northeast R(3) = oimage(i-1,j+1)-oimage(i-1,j-1)+oimage(i+1,j+1)-oimage(i+1,j-1). +oimage(i,j+1)-oimage(i,j-1)+oimage(i,j+1)-oimage(i,j-1); %East R(4) = oimage(i,j+1)-oimage(i-1,j)+oimage(i+1,j)-oimage(i,j-1). +oimage(i+1,j+1)-oimage(i-1,j-1)+oimage(i+1,j+1)-oimage(i-1,j-1); %Southeast R(5) = oimage(i+1,j+1)-oimage(i-1,j-1)-oimage(i-1,j+1)+oimage(i+1,j-1). +oimage(i+1,j)-oimage(i-1,j)+oimage(i+1,j)-oimage(i-1,j); %South R(6) = oimage(i,j-1)-oimage(i-1,j)-oimage(i,j+1)+oimage(i+1,j). +oimage(i-1,j+1)-oimage(i+1,j-1)+oimage(i-1,j+1)-oimage(i+1,j-1); %Southweast R(7) = oimage(i-1,j-1)-oimage(i-1,j+1)-oimage(i+1,j+1)+oimage(i+1,j-1). +oimage(i,j-1)-oimage(i,j+1)+oimage(i,j-1)-oimage(i,j+1); %Weast R(8) = oimage(i,j-1)-oimage(i+1,j)+oimage(i-1,j)-oimage(i,j+1). +oimage(i-1,j-1)-oimage(i+1,j+1)+oimage(i-1,j-1)-oimage(i+1,j+1); %NorthWest % edgeRobinson(i,j) = max(abs(R); m c=max(R); %find the mamximum and the coresponding direction edgeRobinson(i,j)=m; %save the maximum d(i,j)=c; %save the direction endend %-edge detecting-for i = 2:xlen-1 for j = 2:ylen-1 if (d(i,j)=1)&(edgeRobinson(i,j)edgeRobinson(i-1,j)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=2)&(edgeRobinson(i,j)edgeRobinson(i-1,j+1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=3)&(edgeRobinson(i,j)edgeRobinson(i,j+1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=4)&(edgeRobinson(i,j)edgeRobinson(i+1,j+1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=5)&(edgeRobinson(i,j)edgeRobinson(i+1,j)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=6)&(edgeRobinson(i,j)edgeRobinson(i+1,j-1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=7)&(edgeRobinson(i,j)edgeRobinson(i,j-1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; elseif (d(i,j)=8)&(edgeRobinson(i,j)edgeRobinson(i-1,j-1)&(edgeRobinson(i,j)200) edgeRobinsonimage(i,j) = 255; end endend edgeR=edgeRobinsonimage;二阶微分算子采用一阶微分算子很难找到一个一致的阈值选择办法，保证检测出的图像有相对均匀的宽度，克服这个障碍的办法是改用二阶微分算子进行边缘检测定位。经常采用如下Laplace微分算子：并进而寻找的跨零点的位置(零点的局部正和负的取值都有)。当然实践中可以通过模板来实现，本程序采用如下模板：跨零点难得直接出现在某个像素上，往往需要进行插值等手段进行估计，文献“HUERTAS and MEDIONI, Detection of Intensity Changes with Subpixel Accuracy Using Laplacian-Gaussian Masks，1986”给出一种很好的估计方法，本程序采用定义的方法去检测跨零点。具体程序如下：%-filename:secondorder.m-%-Useage:edge detecting by the Laplace operator-%-Writer: by Xubin- function edgeSeO=secondorder(oimage); xlen ylen = size(oimage); %read the sizeedgeXY=zeros(xlen,ylen); edgeimage=zeros(xlen,ylen); H=1,1,1;1,-8,1;1,1,1; for i = 1:xlen for j = 1:ylen for m = 1:3 for n = 1:3 updateX = i-m+2; updateY = j-n+2; if (updateX=1) & (updateX=1) & (updateY180) sum=sum+edgeXY(i,j); counter=counter+1; end endendth=sum/counter; %-edge detecting by zero-crossing-for i = 2:xlen-1 for j = 2:ylen-1 if(edgeXY(i,j)0)&(edgeXY(i,j-1)-edgeXY(i,j)th) edgeimage(i,j)=1; end if(edgeXY(i,j+1)0)&(edgeXY(i,j+1)-edgeXY(i,j)th) edgeimage(i,j)=1; end if(edgeXY(i-1,j)0)&(edgeXY(i-1,j)-edgeXY(i,j)th) edgeimage(i,j)=1; end if(edgeXY(i+1,j)0)&(edgeXY(i+1,j)-edgeXY(i,j)th) edgeimage(i,j)=1; end elseif(abs(edgeXY(i,j)2) if(edgeXY(i,j-1)0)&(edgeXY(i,j+1)-edgeXY(i,j-1)2*th) edgeimage(i,j)=1; end if(edgeXY(i,j-1)0)&(edgeXY(i,j+1)2*th) edgeimage(i,j)=1; end if(edgeXY(i-1,j)0)&(edgeXY(i+1,j)-edgeXY(i-1,j)2*th) edgeimage(i,j)=1; end if(edgeXY(i-1,j)0)&(edgeXY(i+1,j)2*th) edgeimage(i,j)=1; end end endend edgeSeO=edgeimage;LOG：高斯低通与微分算子的复合边缘检测无论什么样的微分算子，直接用来进行边缘检测，会受到噪声很大的干扰。即使是二阶微分算子也不能克服噪声干扰。但是如果采用高斯低通滤波，所得的结果则比较好地保留了图像的边缘特征。Marr-Hildrech的LOG边缘检测算法：(1) 对图像先进行高斯低通滤波，再进行Laplace算子运算；(2) 保留一阶导数峰值的位置记录，然后从中寻找Laplace跨零点；(3) 采用插值方法对零点进行估计。实践当中，取高斯低通滤波器对图像，有容易计算出 具体程序如下：% Use the Log method to detect the edge of an image% the result is effected by the value of sigma_log% $ April.19th, 2010 $ %Writer:Wangxingcheng- % prepare data function edgeL=LOG(Gray_M);height width = size(Gray_M); % calculate the modulesigma_log = 1.02; delta_LOG = zeros(5,5); % initialize the convolution module for i = 1:5 for j=1:5 m_A = (i-3)2+(j-3)2; delta_LOG(i,j)=(exp(-m_A/(2*sigma_log2)*(m_A-2*sigma_log2);%/(2*pi*sigma_log4); endend % edge detection process result_log = conv2(delta_LOG, Gray_M); % 用模板和原始数据卷积result_log = result_log(3:height+2, 3:width+2); % 卷季后的尺寸会变大result_log = result_log + Gray_M; % 加上原来的值突出高频edge_log = zeros(height, width); rr = 2:height-1;cc = 2:width-1;thresh_log = mean2(abs(result_log(rr,cc); %设置阈值 rx,cx = find(result_log(rr,cc)0 & abs(result_log(rr,cc)-result_log(rr,cc+1)thresh_log);%- +的情况edge_log(rx+1)+cx*height)=1;rx,cx = find(result_log(rr,cc-1)0 & result_log(rr,cc)thresh_log);%+ -的情况edge_log(rx+1)+cx*height)=1;rx,cx = find(result_log(rr,cc)0 & abs(result_log(rr,cc)-result_log(rr+1,cc)thresh_log);%- +的情况edge_log(rx+1)+cx*height)=1;rx,cx = find(result_log(rr-1,cc)0 & result_log