图像边缘是图像局部特征不连续性（灰度突变、纹理结构突变等）的反映，它标志着一个区域的终结和另一个区域的开始。边缘的特征是沿边缘走向的像素变化平缓，而垂直于边缘方向的像素变化剧烈。在图像边缘的较小领域中的像素集，它们的灰度值是不连续的，会有较大的跃变。因此可以通过导数来判别图形中像素的灰度值是否存在突变，从而检测出边缘。一般可以用一阶导数和二阶导数来实现。

     1.基于一阶导数的边缘检测

      1）检测原理

       设f(x)是一个离散的一元函数，则其一阶导数可定义为：，由于图像的灰度值函数f(x,y)是一个离散的二元函数，因此通常用梯度来表示图像的一阶导数。所谓梯度指两个由一阶导数组成的向量，其定义为：

      梯度向量的模值大小为：

      为了减小计算量，梯度向量模值可用绝对值近似表示为：

       梯度向量的方向角大小为：

         在边缘检测中，梯度常指梯度向量的模值。基于一阶导数的边缘检测的基本原理：令坐标为（x，y）

 的像素值的梯度值为G(x,y)，若满足G(x,y)>=T，T为阈值，则(x,y)为边缘点所在的位置。

         综上，像素位置(x,y)的梯度为：

         可以用下图中的掩模与点(x,y)的对应邻域内的像素来进行空间卷积得到，若要计算整幅图像所有像素的梯度，只要在图像上方移动掩模，依次进行卷积操作即可。

  代码：   

function h=EdgeDetect(filename)
%基于一阶的边缘检测
%使用的是一阶梯度算子
f=imread(filename);
figure,imshow(f);
[width,height]=size(f);
g=zeros(width,height);
h=zeros(width,height);
df=im2double(f);

wx=[-1 0;1 0];
wy=[-1 1;0 0];
for i=2:width-1
    for j=2:height -1
        gw=[df(i,j) df(i,j+1);df(i+1,j) df(i+1,j+1)];
        g(i,j)=sqrt((sum(sum(wx.*gw)))^2+(sum(sum(wy.*gw)))^2);
    end
end

T=0.25*max(g(:));
h=g>=T;
figure,imshow(h);

%edgedetect(‘ct.bmp‘)

 常用的基于一阶导数的边缘检测器有Roberts/Prewitt、Sobel、Robinson、Kirsch，它们的实现机理与前面介绍的类似，最主要的区别在于实现导数的方法和掩模的选取不同.

       Roberts检测器：

               导数公式：

               掩模及邻域：

        Prewitt检测器：

                导数公式：

               掩模及邻域：

        Sobel检测器：

                导数公式：

                掩模及邻域：

         Robinson检测器：

                 导数公式：

                 掩模及邻域：

          Kirsch检测器：

                  导数公式：

                   掩模及邻域：

         2.基于二阶导数的边缘检测

         1）检测原理

             设f(x)是一个离散的一元函数，其二阶导数可定义为：

                  对于二元的图像灰度函数f(x,y)，其二阶导数一般使用Laplacian算子进行计算。拉普拉斯算子定义为：

                    可得：

                    根据上式可得掩模：

                   通用掩摸如下（α表示锐化程度，取值范围为[0,1]）：

                 基于二阶导数的边缘检测的基本原理：在边缘两侧的像素点，其二阶导数符号相反，他们连线的中点就是边缘的中心。或者说在二阶导数零交叉                  的地方就是边缘中心。

              从检测原理上看，基于二阶导数的边缘检测包含两部分处理过程：一是计算二阶导数，二十寻找零交叉点。由于二阶导数对噪声相对敏感，因此在               计算前需要对图像进行平滑处理。在选择平滑滤波器时需要满足两个标准：第一，滤波器应该是平滑的且在频域中大致上是有限带宽的，以便减少               会导致函数变化的可能频率数；第二，空间定位的约束要求滤波器的响应需来自于图像中邻近的点。这两个标注是矛盾的，但是通过使用高斯分布               可以同时得到优化。

               高斯函数基本形式：

             由于高斯函数图形是关于中心对称的，为计算导数方便，可令r²=x²+y².

             则其一阶导数为：

             二阶导数为：

               得到导数计算公式：该公式通常称为高斯型的拉普拉斯算子(LoG),由于其图形形状的原               因，该算子也被称为墨西哥草帽(Mexicanhat)。

              综上，可以得到计算图像f(x,y)二阶导数的两种是实现方法：

                      1）

                      2）  

                上述两个公式是干货，第一种方法是先用高斯与图像做卷积（用一次掩模），再求卷积结果的二阶导数（再用一次卷积掩模）。第二种方法更简单              些利用了卷积与微分运算的可交换性，只需要一次卷积掩模。

                可由下式掩模(LoG掩模)实现：

                   综上可得图像二阶导数对应的LoG图像，下面为寻找零交叉点 两种方法：                                                

                         1）设LoG图像的所有正值区域为白色，负值区域为黑色。在白色与黑色之间的就是零交叉点。这种方法实现简单,但会出现空心粉效应。

                         2）使用一个2*2的窗口检测零交叉点，当两种极性的LoG图像数值同时出现在该窗口内时，就将边缘标签任意赋给一个角点，比如说左                              上点；当窗口内的数值都是正的或负的时，就不给边缘标签。这种方法灵敏度很高，会检测出非常弱的边缘的零交叉点，因此一                                      般需要通过一阶导数(下面代码中使用了Sobel算子)来消除这些点。

   代码：

function Zero_Crossing(filename)
%基于二阶导数的边缘检测
%使用Log掩膜 Laplacian of a Gaussian
f=imread(filename);
figure,imshow(f);
[width,height]=size(f);
g1=zeros(width,height);
g2=g1;
df=im2double(f);
w=[0 0 -1 0 0;0 -1 -2 -1 0;-1 -2 16 -2 -1;0 -1 -2 -1 0;0 0 -1 0 0];
for i=3:width-2
    for j=3:height-2
        gw=[df(i-2,j-2) df(i-2,j-1) df(i-2,j) df(i-2,j+1) df(i-2,j+2);df(i-1,j-2) df(i-1,j-1) df(i-1,j) df(i-1,j+1) df(i-1,j+2);df(i,j-2) df(i,j-1) df(i,j) df(i,j+1) df(i,j+2);df(i+1,j-2) df(i+1,j-1) df(i+1,j) df(i+1,j+1) df(i+1,j+2);df(i+2,j-2) df(i+2,j-1) df(i+2,j) df(i+2,j+1) df(i+2,j+2)];
        g2(i,j)=sum(sum(w.*gw));
    end
end

%使用Sobel算子检测器
wx=[-1 -2 -1;0 0 0;1 2 1];              
wy=[-1 0 1;-2 0 2;-1 0 1];
for i=2:width-1
    for j=2:height-1
                gw=[df(i-1,j-1) df(i-1,j) df(i-1,j+1);df(i,j-1) df(i,j) df(i,j+1);df(i+1,j-1) df(i+1,j) df(i+1,j+1)];
                g1(i,j)=sqrt((sum(sum(wx.*gw)))^2+(sum(sum(wy.*gw)))^2)/2;
    end
end

for i=1:width
    for j=1:height
        if g2(i,j)>0
            if(g2(i+1,j)<0 || g2(i,j+1)<0 || g2(i+1,j+1)<0)
                g2(i,j)=1;
            else
                g2(i,j)=0;
            end
        elseif g2(i,j)<0
            if(g2(i+1,j)>0 || g2(i,j+1)>0 || g2(i+1,j+1)>0)
                g2(i,j)=1;
            else
                g2(i,j)=0;
            end
        else
            g2(i,j)=0;
        end
    end
end

g1=g1>=0.25*max(g1(:));
g=g1.*g2;
figure,imshow(g);

边缘检测