低层次特征提取（一）————边缘检测

    低层次特征是不需要任何形状/空间关系的信息就可以从图像中自动提取的基本特征。所有低层次方法都可以应用于高层次特征提取，从而在图像中找到形状。第一种低层次特征，称之”edge detection”。它的目的要是要制作一个线图。一阶检测算子相于一阶微分法，二阶边缘检测算子相当于高一阶微分处理。

    边缘检测：在视觉计算理论框架中，抽取二维图像上的边缘、角点、纹理等基本特征，是整个系统框架中的第一步。这些特征所组成的图称为基元图。在不同“尺度”意义下的边缘点，在一定条件下包含了原图像的全部信息。

定义：

    •目前，具有对边缘的描述性定义，即两个具有不同灰度的均匀图像区域的边界，即边界反映局部的灰度变化。

    •局部边缘是图像中局部灰度级以简单（即单调）的方式作极快变换的小区域。这种局部变化可用一定窗口运算的边缘检测算子来检测。

边缘的描述：

1)边缘法线方向——在某点灰度变化最剧烈的方向，与边缘方向垂直；

2)边缘方向——与边缘法线方向垂直，是目标边界的切线方向；

3) 
边缘强度——沿边缘法线方向图像局部的变化强度的量度。

    边缘检测的基本思想是通过检测每个像素和其邻域的状态，以决定该像素是否位于一个物体的边界上。如果一个像素位于一个物体的边界上，则其邻域像素的灰度值的变化就比较大。假如可以应用某种算法检测出这种变化并进行量化表示，那么就可以确定物体的边界。

边缘检测算法有如下四个步骤：

滤波：边缘检测算法主要是基于图像强度的一阶和二阶导数，但导数的计算对噪声很敏感，因此必须使用滤波器来改善与噪声有关的边缘检测器的性能．需要指出，大多数滤波器在降低噪声的同时也导致了边缘强度的损失，因此，增强边缘和降低噪声之间需要折衷．

增强：增强边缘的基础是确定图像各点邻域强度的变化值．增强算法可以将邻域（或局部）强度值有显著变化的点突显出来．边缘增强一般是通过计算梯度幅值来完成的．

检测：在图像中有许多点的梯度幅值比较大，而这些点在特定的应用领域中并不都是边缘，所以应该用某种方法来确定哪些点是边缘点．最简单的边缘检测判据是梯度幅值阈值判据．

定位：如果某一应用场合要求确定边缘位置，则边缘的位置可在子像素分辨率上来估计，边缘的方位也可以被估计出来．

   在边缘检测算法中，前三个步骤用得十分普遍。这是因为大多数场合下，仅仅需要边缘检测器指出边缘出现在图像某一像素点的附近，而没有必要指出边缘的精确位置或方向．边缘检测误差通常是指边缘误分类误差，即把假边缘判别成边缘而保留，而把真边缘判别成假边缘而去掉．边缘估计误差是用概率统计模型来描述边缘的位置和方向误差的．我们将边缘检测误差和边缘估计误差区分开，是因为它们的计算方法完全不同，其误差模型也完全不同．

边缘检测的三个共性准则：

•好的检测结果，或者说对边缘的误测率尽可能低，就是在图像边缘出现的地方检测结果中不应该没有；另一方面不要出现虚假的边缘；

•对边缘的定位要准确，也就是我们标记出的边缘位置要和图像上真正边缘的中心位置充分接近；

•对同一边缘要有尽可能低的响应次数，也就是检测响应最好是单像素的。

几种常用的边缘检测算子主要有Roberts边缘检测算子，Sobel算子、Prewitt算子、Krisch边缘算子，高斯-拉普拉斯算子。

边缘大致可以分为两种，一种是阶跃状边缘，边缘两边像素的灰度值明显不同；另一种为屋顶状边缘，边缘处于灰度值由小到大再到小的变化转折点处。
边缘检测的主要工具是边缘检测模板。我们以一个一维模板为例来考察边缘检测模板是如何作用的。假设有一个模板 和一幅图象

可以看出，图象中左边暗，右边亮，中间存在着一条明显的边缘，是一个典型阶跃状边缘。使用模板 进行模板操作后，结果如下

可以看出，边缘检测后的图象在原图象暗亮边缘处的灰度值高很多。观察时，就能发现一条很明显的亮边，其他区域都很暗，这样就起到了边缘检测的作用。
模板的作用是将右邻点的灰度值减去左邻点的灰度值作为该点的灰度值。在灰度相近的区域内，这么做的结果使得该点的灰度值接近于0；而在边缘附近，灰度值有明显的跳变，这么做的结果使得该点的灰度值很大，这样就出现了上面的结果。这种模板就是一种边缘检测器，它在数学上的涵义是一种基于梯度的滤波器，习惯上又称边缘算子。我们知道，梯度是有方向的，和边缘的方向总是垂直的。模板是水平方向的，而上面那幅图象的边缘恰好是垂直方向的，使用模板就可以将它检测出来。如果图象的边缘是水平方向的，我们可以用梯度是垂直
方向的模板 检测它的边缘。如果图象的边缘是45。方向的，我们可以用模板 检测它的边缘。

常用的边缘检测模板有Laplacian算子、Roberts算子、Sobel算子、log（Laplacian-Gauss）算子、Kirsch算子和Prewitt算子等。 

下面是算子的基本模板：

1、Roberts算子

如果我们沿如下图方向角度求其交叉方向的偏导数，则得到Roberts于1963年提出的交叉算子边缘检测方法。该方法最大优点是计算量小，速度快。但该方法由于是采用偶数模板，如下图所示，所求的(x,y)点处梯度幅度值，其实是图中交叉点处的值，从而导致在图像(x,y)点所求的梯度幅度值偏移了半个像素（见下图）。

 

Roberts算子的推导：

 

2、Sobel算子

Sobel算子也有两个，一个是检测水平边缘的模板 ，另一个是检测水平边缘的模板。与Prewitt算子相比，Sobel算子对于象素位置的影响作了加权，因此效果更好。

sobel算子的另一种形式是各向同性Sobel算子，也有两个模板组成，一个是检测水平边缘的 ，另一个是检测垂直边缘的。各向同性Sobel算子和普通Sobel算子相比，位置加权系数更为准确，在检测不同方向的边缘时梯度的幅度一致。

Sobel算子的推导：

3、Prewitt算子

Prewitt算子由两部分组成，检测水平边缘的模板

和检测垂直边缘的模板

prewitt算子一个方向求微分，一个方向求平均，所以对噪声相对不敏感。

 4、laplacian算子

Laplacian算子定义为

它的差分形式为

表示成模板的形式就是 。Laplacian算子另外一种形式是，也经常使用。Laplace算子是一种各向同性算子，在只关心边缘的位置而不考虑其周围的象素灰度差值时比较合适。Laplace算子对孤立象素的响应要比对边缘或线的响应要更强烈，因此只适用于无噪声图象。存在噪声情况下，使用Laplacian算子检测边缘之前需要先进行低通滤波。

5、canny算子

5.1. Canny边缘检测基本原理

     (1)图象边缘检测必须满足两个条件：一能有效地抑制噪声；二必须尽量精确确定边缘的位置。

     (2)根据对信噪比与定位乘积进行测度，得到最优化逼近算子。这就是Canny边缘检测算子。

     (3)类似与Marr（LoG）边缘检测方法，也属于先平滑后求导数的方法。

5.2 Canny边缘检测算法：

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

     step2:用一阶偏导的有限差分来计算梯度的幅值和方向；

     step3:对梯度幅值进行非极大值抑制；

     step4:用双阈值算法检测和连接边缘。

step1:高斯平滑函数

Step2：一阶微分卷积模板

step3:非极大值抑制仅仅得到全局的梯度并不足以确定边缘，因此为确定边缘，必须保留局部梯度最大的点，而抑制非极大值。（non-maxima suppression,NMS）

解决方法：利用梯度的方向。

 

图1非极大值抑制

四个扇区的标号为0到3，对应3*3邻域的四种可能组合。在每一点上，邻域的中心象素M与沿着梯度线的两个象素相比？？？？。如果M的梯度值不比沿梯度线的两个相邻象素梯度值大，则令M=0。

即： 

 Step4:阈值化

    减少假边缘段数量的典型方法是对N[i，j]使用一个阈值。将低于阈值的所有值赋零值。但问题是如何选取阈值？

    解决方法：双阈值算法。双阈值算法对非极大值抑制图象作用两个阈值τ1和τ2，且2τ1≈τ2，从而可以得到两个阈值边缘图象N1［i,j］和N2［i，j］。由于N2［i，j］使用高阈值得到，因而含有很少的假边缘，但有间断(不闭合)。双阈值法要在N2［i，j］中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在N1［i,j］的8邻点位置寻找可以连接到轮廓上的边缘，这样，算法不断地在N1［i,j］中收集边缘，直到将N2［i,j］连接起来为止。

5.3 canny算法程序实现

Canny算法程序中将上述的4个步骤再加以细分，分成以下7步：

l 生成高斯滤波系数；

l 用生成的高斯滤波系数对原图像进行平滑；

l 求滤波后图像的梯度；

l 进行非最大抑制；

l 统计图像的直方图，对阈值进行判定；

l 利用函数寻找边界起点；

l 根据第6步执行的结果，从一个像素点开始搜索，搜索以该像素点为边界起点的一条边界的一条边界的所有边界点；

各个算子在Matlab中使用格式和边缘检测方法：

1. Roberts算子

Roberts交叉算子为梯度幅值计算提供了一种简单的近似方法：

Roberts算子是该点连续梯度的近似值，而不是所预期的在点处的近似值。

在MATLAB中可以由edge函数实现。其语法格式如下：

BW=edge(I, ‘roberts’)

BW=edge(I, ‘roberts’,thresh)

[BW,thresh]=edge(I,’roberts’,…)

BW=edge(I,’roberts’)自动选择阈值用Robert算子进行边缘检测。

BW=edge(I,’roberts’,thresh)根据所指定的敏感度阈值thresh用Roberts算子进行边缘检测，它忽略了所有小于阈值的边缘。当thresh为空时，自动选择阈值。

[BW,thresh]=edge(I,’roberts’,…)返回阈值。

edge函数对灰度图像I进行边缘检测，返回与I同样大小的二值图像BW，其中1表示I的边缘，0表示非边缘。I是uint8型、uint16型或double型的，BW是uint8型的。

clc
close all
clear all
%%%生成高斯平滑滤波模板%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
hg=zeros(3,3); %设定高斯平滑滤波模板的大小为3*3
delta=0.5;
for x=1:1:3
for y=1:1:3
u=x-2;
v=y-2;
hg(x,y)=exp(-(u^2+v^2)/(2*pi*delta^2));
end
end
h=hg/sum(hg(:));
%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%读入图像%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%

f = imread('1111.tif'); % 读入图像文件
f=rgb2gray(im2double(f));
imshow(f)
title('原始图像');
[m,n]=size(f);
ftemp=zeros(m,n);
rowhigh=m-1;
colhigh=n-1;
%%%高斯滤波%%%
for x=2:1:rowhigh-1
for y=2:1:colhigh-1
mod=[f(x-1,y-1) f(x-1,y) f(x-1,y+1); f(x,y-1) f(x,y) f(x,y+1);f(x+1,y-1) f(x+1,y) f(x+1,y+1)];
A=h.*mod;
ftemp(x,y)=sum(A(:));
end
end
f=ftemp
figure,imshow(f)
title('通过高斯滤波器后的图像');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%利用roberts算子进行边缘检测%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sx=[-1 -2 -1;0 0 0;1 2 1];
sy=[-1 0 1;-2 0 2;-1 0 1];%%%%%你可以替换成其他算子，这里是罗伯特算子

% sx=[-1 -2 -1;0 0 0;1 2 1];
% sy=[-1 0 1;-2 0 2;-1 0 1];这个是Sobel算子，类似的，你可以替换成canny算子等等
for x=2:1:rowhigh-1
for y=2:1:colhigh-1
mod=[f(x-1,y-1) f(x-1,y) f(x-1,y+1); f(x,y-1) f(x,y) f(x,y+1);f(x+1,y-1) f(x+1,y) f(x+1,y+1)];
fsx=sx.*mod;
fsy=sy.*mod;
ftemp(x,y)=sqrt((sum(fsx(:)))^2+(sum(fsy(:)))^2);
end
end
fr=im2uint8(ftemp);
figure,imshow(fr)
title('用roberts算子边缘检测的原始图像');

%%%域值分割%%%
TH1=60; %设定阈值
for x=2:1:rowhigh-1
for y=2:1:colhigh-1
if (fr(x,y)>=TH1)&((fr(x,y-1) <= fr(x,y)) & (fr(x,y) > fr(x,y+1)) )
fr(x,y)=200;
elseif(fr(x,y)>=TH1)&( (fr(x-1,y) <=fr(x,y)) & (fr(x,y) >fr(x+1,y)))
fr(x,y)=200;
else fr(x,y)=50;
end
end
end
figure,imshow(fr)
title('用roberts算子边缘检测并细化后的图像'); 

2. Sobel算子

Sobel算子是边缘检测器中最常用的算子之一。采用3×3邻域可以避免在像素之间内插点上计算梯度。考虑如图10-5所示的点周围点的排列。

edge函数实现的语法格式如下：

BW=edge(I, ‘sobel’)

BW=edge(I, ‘sobel’,thresh)

BW=edge (I, ‘sobel’, thresh,direction)

[BW, thresh]=edge (I, ‘sobel’…)

BW=edge(I, ‘sobel’)自动选择阈值用Sobel算子进行边缘检测。

BW=edge(I,’sobel’,thresh)根据所指定的敏感度阈值thresh，用Sobel算子进行边缘检测，它忽略了所有小于阈值的边缘。当thresh为空时，自动选择阈值。

BW=edge (I, ‘sobel’, thresh,direction)根据所指定的敏感度阈值thresh，在所指定的方向direction上，用Sobel算子进行边缘检测。Direction可取的字符串值为horizontal(水平方向)、vertical(垂直方向)或both(两个方向)。

[BW, thresh]=edge (I,’sobel’…)返回阈值。

f=imread('1.jpg');
f=rgb2gray(f);%转化成灰度图
f=im2double(f);%函数im2double 将其值归一化到0～1之间
%使用垂直Sobcl箅子．自动选择阈值
[VSFAT Threshold]=edge(f, 'sobel','vertical'); %边缘探测
figure,imshow(f),title(' 原始图像，');%显示原始图像
figure,imshow(VSFAT),title( '垂直图像边缘检测');
%显示边缘探测图像
%使用水平和垂直Sobel算子，自动选择阈值
SFST=edge(f,'sobel',Threshold);
figure,imshow(SFST),title('水平和垂直图像边缘检测');
%显示边缘探测图像
%使用指定45度角Sobel算子滤波器，指定阂值
s45=[-2 -1 0;-1 0 1;0 1 2];
SFST45=imfilter(f,s45,'replicate');%功能：对任意类型数组或多维图像进行滤波。
SFST45=SFST45>=Threshold;
figure,imshow(SFST45),title('45度角图像边缘检测') ;

3. Prewitt算子

Prewitt于1970年左右提出了Prewitt算子。由上面对Sobel算子的推导，同时也得出了Prewitt算子。Prewitt算子和Sobel算子的方程完全一样，不同的是只是a和b的系数不同。以下是Prewitt的两个算子：

edge函数实现的语法格式如下：

BW=edge(I, ‘prewitt’)

BW=edge (I, ‘prewitt’,thresh)

BW=edge (I, ‘prewitt’, thresh,direction)

[BW, thresh]=edge (I,’prewitt’…)

BW=edge(I, ‘prewitt’)自动选择阈值用Prewitt算子进行边缘检测。

BW=edge(I,’prewitt’,thresh)根据所指定的敏感度阈值thresh，用Prewitt算子进行边缘检测，它忽略了所有小于阈值的边缘。当thresh为空时，自动选择阈值。

BW=edge (I, ‘prewitt’, thresh,direction)根据所指定的敏感度阈值thresh，在所指定的方向direction上，用Prewitt 算子进行边缘检测。Direction可取的字符串值为horizontal(水平方向)、vertical(垂直方向)或both(两个方向)。

[BW, thresh]=edge (I, ‘prewitt’…)返回阈值。

I = imread('1.BMP');

BW1 = edge(I,'prewitt',0.04);             % 0.04为梯度阈值

figure(1);

imshow(I);

figure(2);

imshow(BW1);

4. LOG 算子

它是一个二阶算子，将在边缘处产生一个陡峭的零交叉。拉普拉斯算子是一个线性的、移不变的算子，它的传递函数在频域空间的圆点是零，因此经拉普拉斯滤波过的图像具有零平均灰度。LOG算子先用高斯低通滤波器将图像进行预先平滑，然后用拉普拉斯算子找出图像中的陡峭边缘，最后用零灰度值进行二值化产生闭合的、连通的轮廓，消除了所有内部点。

edge函数实现的语法格式如下：

BW=edge(I, ‘log’)

BW=edge (I, ‘log’,thresh)

BW=edge (I, ‘log’, thresh,sigma)

[BW, thresh]=edge (I, ‘log’…)

BW=edge(I, ‘log’)自动选择阈值用LOG算子进行边缘检测。

BW=edge(I,’log’,thresh)根据所指定的敏感度阈值thresh，用LOG算子进行边缘检测，它忽略了所有小于阈值的边缘。当thresh为空时，自动选择阈值。当指定thresh为0时，输出图像具有闭合的轮廓，因为其中包含了输入图像中的所有零交叉点。

BW=edge (I, ‘log’, thresh,sigma)根据所指定的敏感度阈值thresh和标准偏差sigma，用LOG算子进行边缘检测，默认时sigma等于2，滤波器是n×n维的。其中

n=ceil(sigma×3)×2+1

[BW, thresh]=edge (I, ‘log’…)返回阈值。

a = imread('1.png');
sigma=0.15;    %sigma=0.5（参数可以调不同的试一下）;
[m,n]=size(a);
e=repmat(logical(uint8(0)),m,n);
rr=2:m-1;
cc=2:n-1;
fsize=ceil(sigma*7)*6+1;   %fsize=ceil(sigma*3)*2+1;
op=fspecial('log',fsize,sigma);
op=op-sum(op(:))/prod(size(op));

5、Canny边缘检测

Canny算子检测边缘的方法是寻找图像梯度的局部极大值，梯度是用高斯滤波器的导数计算的。Canny方法使用两个阈值来分别检测强边缘和弱边缘，而且仅当弱边缘与强边缘相连时，弱边缘才会包含在输出中。因此，此方法不容易受噪声的干扰，能够检测到真正的弱边缘。

BW=edge(I, ‘canny’)

BW=edge (I, ‘canny’,thresh)

BW=edge (I, ‘canny’, thresh,sigma)

[BW, thresh]=edge (I, ‘canny’…)

BW=edge(I, ‘canny’)自动选择阈值用Canny算子进行边缘检测。

BW=edge(I,’canny’,thresh)根据所指定的敏感度阈值thresh，用Canny算子进行边缘检测，thresh是一个含两个元素的矢量，第一个元素是低阈值，第二个元素是高阈值；如果只给thresh指定一个值，则此值作为高阈值，而0.4×thresh作为低阈值；当thresh为空时，自动选择低阈值和高阈值。

BW=edge (I, ‘canny’, thresh,sigma)根据所指定的敏感度阈值thresh和标准偏差sigma，用Canny算子进行边缘检测，默认时sigma等于1，滤波器的尺寸sigma自动选择。

[BW, thresh]=edge (I, ‘canny’…)返回含两个元素的阈值矢量。

function can()

 

    I = imread('un.bmp');
   
    gray = rgb2gray(I);
   
    a = im2single(gray);
   
    [m,n] = size(a);
    % 用于输出的边界位图
    e = false(m,n);
   
    % Magic numbers
    GaussianDieOff = .0001; 
    PercentOfPixelsNotEdges = .7; % 用于阀值选择
    ThresholdRatio = .4;          % 低阀值相对高阀值的比值
    sigma = 1; %设置sigma
    thresh = [];
   
    % 设计滤波器 - a gaussian和它的导数
    pw = 1:30; % possible widths
    ssq = sigma^2;
    width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');%find函数很给力...
    if isempty(width)
    width = 1;  % the user entered a really small sigma
    end
   
    t = (-width:width);
    gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);     % 高斯一维滤波
   
    % Find the directional derivative of 2D Gaussian (along X-axis)
    % Since the result is symmetric along X, we can get the derivative along
    % Y-axis simply by transposing the result for X direction.
    [x,y]=meshgrid(-width:width,-width:width);
    dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);%二维高斯方向导数
   
    % Convolve the filters with the image in each direction
    % The canny edge detector first requires convolution with
    % 2D gaussian, and then with the derivitave of a gaussian.
    % Since gaussian filter is separable, for smoothing, we can use
    % two 1D convolutions in order to achieve the effect of convolving
    % with 2D Gaussian.  We convolve along rows and then columns.
   
    %用一阶高斯滤波器平滑图像
    aSmooth=imfilter(a,gau,'conv','replicate');   % run the filter accross rows
    aSmooth=imfilter(aSmooth,gau','conv','replicate'); % and then accross columns
   
    %应用方向导数
    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;   % normalize
    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 < high < 1.';
      error(eid,'%s',msg);
     end
    end
   
    % The next step is to do the non-maximum supression. 
    % We will accrue indices which specify ON pixels in strong edgemap
    % The array e will become the weak edge map.
    idxStrong = []; 
    for dir = 1:4
      idxLocMax = cannyFindLocalMaxima(dir,ax,ay,mag);
      idxWeak = idxLocMax(mag(idxLocMax) > lowThresh);
      e(idxWeak)=1;
      idxStrong = [idxStrong; idxWeak(mag(idxWeak) > highThresh)];
    end
   
    if ~isempty(idxStrong) % result is all zeros if idxStrong is empty
      rstrong = rem(idxStrong-1, m)+1;
      cstrong = floor((idxStrong-1)/m)+1;
      e = bwselect(e, cstrong, rstrong, 8);
      e = bwmorph(e, 'thin', 1);  % Thin double (or triple) pixel wide contours
    end
    
    imshow(e);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Local Function : cannyFindLocalMaxima
%
function idxLocMax = cannyFindLocalMaxima(direction,ix,iy,mag)

 [m,n] = size(mag);
 
 % Find the indices of all points whose gradient (specified by the
 % vector (ix,iy)) is going in the direction we're looking at. 
   idx = find((iy<=0 & ix>-iy)  | (iy>=0 & ix<-iy));
 
   % Exclude the exterior pixels
 if ~isempty(idx)
   v = mod(idx,m);
   extIdx = find(v==1 | v==0 | idx<=m | (idx>(n-1)*m));
   idx(extIdx) = [];
 end
 
 ixv = ix(idx); 
 iyv = iy(idx);  
 gradmag = mag(idx);
 
 % Do the linear interpolations for the interior pixels
  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);

6、Susan边缘检测

Susan准则用一个圆形模板扫描图像，若模板内其他任意像素的灰度值与模板中心元素(核)的灰度值的差小于一定阈值，就认为该点与核具有相同(或相近)的灰度值，满足这样条件的像素组成的区域称为核值相似区(UnivalueSegment AssimilatingNucle2us，USAN)。把图像中的每个像素与具有相近灰度值的局部区域联系是Susan准则的基础。

具体检测时，用圆形模板扫描整个图像，比较模板内每一个像素与中心像素的灰度值，并给定阈值来判别该像素是否属于USAN区域。如下式：

式中，c(r,r₀)为模板内属于USAN区域的像素的判别函数，I(r₀)是模板中心像素(核)的灰度值，I(r)为模板内其他任意像素的灰度值，t是灰度差门限。它影响检测到角点的个数。t减少，获得图像中更多精细的变化，从而给出相对较多信息，则图像中某一点的USAN区域大小可由下式表示：

USAN区域包含了图像局部许多重要的结构信息，它的大小反映了图像局部特征的强度，模板中的白色区域即为USAN区。在平坦区域USAN区最大，如模板e所示；在边缘处USAN 区大小降为一半，如模板b所示；而在角点附近USAN区变得更小，如模板a所示。由此可得到SUSAN 提取边缘和角点算法的基本原理：即在边缘和角点处的USAN 区最小。因此，可根据USAN区的大小和矩阵特性来检测图像边缘与角点等特征的位置和方向信息。当模板完全处于背景或目标中时，USAN区域最大，当模板移向目标边缘时，USAN区域逐渐变小；当模板中心处于角点时，USAN区域很小。得到每个像素对应的USAN区域大小之后，由上式可知，对于3×3大小的模板，USAN区的最大值n_max =8。根据试验分析可知，在实际噪声图像中如果核心点在边缘附近，n的值一般不会大于3n_max/4。因此，定义几何阈值g= 3n_max/4，利用下式产生边缘初始响应：

其中，g为几何门限，得到的边缘初始响应值大小符合SUSAN原理，即USAN区域越小，初始边缘响应就越大。

function image_out =susan(im,threshold)

d = length(size(im));

if d==3

image=double(rgb2gray(im));

elseif d==2

image=double(im);

end

mask = ([ 0 0 1 1 1 0 0 ;0 1 1 1 11 0;1 1 1 1 1 1 1;1 1 1 1 1 1 1;1 1 1 1 1 1 1;0 1 1 1 1 1 0;0 0 1 11 0 0]);

R=zeros(size(image));

nmax = 3*37/4;

[a b]=size(image);

new=zeros(a+7,b+7);

[c d]=size(new);

new(4:c-4,4:d-4)=image;

for i=4:c-4

for j=4:d-4

current_image = new(i-3:i+3,j-3:j+3);

current_masked_image = mask.*current_image;

current_thresholded =susan_threshold(current_masked_image,threshold);

g=sum(current_thresholded(:));

if nmax

7、小波边缘检测

除了以上介绍的边缘检测方法，其实在目前的边缘检测算法的研究过程中，小波技术在边缘检测中也成为一种主要的研究方向。这里简单地介绍一下。

基于内容的方法是近来研究的热点，利用小波进行边缘检测也是一种方法，它还可用于特征点的检测，Mexican Hatwavelet尺度交互方法最初由Manjunath等提出。通过识别发生在同一图像不同尺度版本显著强度变化来决定特征点。该方法应用MexicanHat wavelet的两个不同尺度到同一图像并且计算两个尺度图像的尺度交互图像，尺度差图像的局部极值决定了特征点。MexicanHat wavelet，也叫做Marrwavelet，具有旋转不变性，因为它有圆形对称频域响应。位于x的Mexican Hat wavelet定义为

其中尺度差图像由下列式子得到代表尺度为i的Mexican Hatwavelet在位置x的响应。是一个规范常数，代表不同尺度i和j的尺度差，的局部极值就决定了潜在的特征点集。这些点中强度超过一定阈值的就被认为是特征点。

clear all;
load wbarb; %小波变换边缘提取程序
I = ind2gray(X,map);%检索图转成灰度图
imshow(I);
I1 = imadjust(I,stretchlim(I),[0,1]);%调整图像的像素值，可以改变对比度和颜色
figure;
imshow(I1);
[N,M] = size(I);
h = [0.125,0.375,0.375,0.125];
g = [0.5,-0.5];
delta = [1,0,0];
J = 3;
a(1:N,1:M,1,1:J+1) = 0;
dx(1:N,1:M,1,1:J+1) = 0;
dy(1:N,1:M,1,1:J+1) = 0;
d(1:N,1:M,1,1:J+1) = 0;
a(:,:,1,1) = conv2(h,h,I,'same'); %二维卷积
dx(:,:,1,1) = conv2(delta,g,I,'same');
dy(:,:,1,1) = conv2(g,delta,I,'same');
x = dx(:,:,1,1);
y = dy(:,:,1,1);
d(:,:,1,1) = sqrt(x.^2+y.^2);
I1 = imadjust(d(:,:,1,1),stretchlim(d(:,:,1,1)),[0 1]);figure;imshow(I1);
lh = length(h);
lg = length(g);
for j = 1:J+1
lhj = 2^j*(lh-1)+1;
lgj = 2^j*(lg-1)+1;
hj(1:lhj)=0;
gj(1:lgj)=0;
for n = 1:lh
hj(2^j*(n-1)+1)=h(n);
end
for n = 1:lg
gj(2^j*(n-1)+1)=g(n);
end
a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same');
dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same');
x = dx(:,:,1,j+1);
y = dy(:,:,1,j+1);
dj(:,:,1,j+1) = sqrt(x.^2+y.^2);
I1 = imadjust(dj(:,:,1,j+1),stretchlim(dj(:,:,1,j+1)),[0 1]);figure;imshow(I1);
end

8、用水线阈值法分割图像

m = imread('afmsurf.tif');figure, imshow(afm);

se = strel('disk', 15);

Itop = imtophat(afm, se); % 高帽变换

Ibot = imbothat(afm, se); % 低帽变换

figure, imshow(Itop, []);   % 高帽变换，体现原始图像的灰度峰值

figure, imshow(Ibot, []);   % 低帽变换，体现原始图像的灰度谷值

Ienhance = imsubtract(imadd(Itop, afm), Ibot);% 高帽图像与低帽图像相减，增强图像

figure, imshow(Ienhance);

Iec = imcomplement(Ienhance); % 进一步增强图像

Iemin = imextendedmin(Iec, 20); figure,imshow(Iemin) % 搜索Iec中的谷值

Iimpose = imimposemin(Iec, Iemin);

wat = watershed(Iimpose); % 分水岭分割

rgb = label2rgb(wat); figure, imshow(rgb); % 用不同的颜色表示分割出的不同区域