SUSAN角点检测算法，及其Matlab实现_susan算子检测角点的步骤

1、SUSAN角点检测算法步骤
（1）在图像上放置一个37个像素的圆形模板，模板在图像上滑动，依次比较模板内各个像素点的灰度与模板核的灰度，判断是否属于USAN区域。判别函数如下：
其中，${\vec{r}}_0$.表示二维图像中核心点的位置，$\vec{r}$表示模板中其他点的位置，$I({\vec{r}}_0)$表示图像在${\vec{r}}_0$处的像素值，$I(\vec{r})$表示图像在$\vec{r}$处的像素值，$t$是表示亮度插值的一个门限值，决定了USAN区域各点之间最大的亮度差值，$c$是一个用来比较输出的函数。也可以用平滑的县来代替这种直接的分割方式，这样可以获得更稳定而敏感的结果，虽然计算复杂，但是可以通过查表来获得较快的速度。公式如下：
$$c(\vec{r},{\vec{r}}_0)=e^{-{\left(\frac{I(\vec{r})-I({\vec{r}}_0)}{t}\right)}^6}$$
（2）统计圆形模板中和核心点有相似亮度值的像素值个数$n(r_0)$$$n({\vec{r}}_0)=\sum _{\vec{r}\in D({\vec{r}}_0)}c(\vec{r},{\vec{r}}_0)$$其中，$D({\vec{r}}_0)$是以${\vec{r}}_0$为中心的圆形模板区域，$n$是USAN区域中像素的个数。
（3）使用如下角点响应函数，若某个像素点的USAN值小于某一特定阈值，则该点被认为是初始角点。将$n$同一个固定的阈值$g$比较（一般设置为最大与中心相似点数$n_{\max }$的0.5倍左右，这是与边缘检测不同的地方），初始的边缘响应可以用下面的等式计算：

角点误报检查：
约束1：计算出USAN区域的重心，如果重心与圆心的距离大于一定的阈值才认为是角点。
约束2：计算出USAN区域的重心，从中心指向USAN区域重心直线上所有像素都必须是USAN的一部分，否则认为不是角点。
（4）对初始角点进行非极值抑制来求得最后的角点。
2、Matlab实现

%% SUSAN角点检测算法
clear 
close all
clc

img=imread('board.jpg');
[m,n,c]=size(img);
if c>1
   img_gray=rgb2gray(img);
else
   img_gray=img;
end
img_gray=double(img_gray);
t=25;   % 阈值
g=37/2+1;
R=zeros(m,n);
%% 1、圆形模板在图像上滑动，依次比较模板内各个像素点的灰度与模板核的灰度，判断是否属于USAN区域
for i=4:m-3
   for j=4:n-3
      tmp=img_gray(i-3:i+3,j-3:j+3);
      usan=0;  % USAN计数
      centre_x=0;  % USAN区域的重心
      centre_y=0;
      %% 2、统计圆形模板中和核心点有相似亮度值的像素值个数
      for p=1:7
         for q=1:7
            if (p-4)^2+(q-4)^2<12
               if abs(img_gray(i,j)-tmp(p,q))<t
                  centre_x=centre_x+p;  
                  centre_y=centre_y+q;
                  usan=usan+1; 
               end
            end
         end
      end
      centre_x=centre_x/usan;
      centre_y=centre_y/usan;
      %% 模板中心与USAN区域重心约束
      if sqrt((centre_x-4)^2+(centre_y-4)^2)<0.2  
         usan=37; 
         %% 模板中心与USAN区域重心连线上的像素约束
         k=(centre_y-4)/(centre_x-4);
         b=(centre_x*4-4*centre_y)/(centre_x-4);
         if centre_x<4
            centre_x=ceil(centre_x);    % 向上取整
         else
            centre_x=floor(centre_x);   % 向下取整 
         end
         if centre_y<4
            centre_y=ceil(centre_y);    % 向上取整
         else
            centre_y=floor(centre_y);   % 向下取整 
         end        
         for p=min(centre_x,4):max(centre_x,4) 
            for q=min(centre_y,4):max(centre_y,4) 
               if (abs(k*p-q)<1) && (abs(img_gray(i,j)-tmp(p,q))>t) 
                  usan=37; 
                  break;
               end
            end
         end
      end
      %% 3、计算角点响应，确定初始角点
      if usan<g
         R(i,j)=g-usan; 
      end
   end
end
%% 4、非极大值抑制
s=5;    % 邻域大小
Re=zeros(m,n);
Rmax=0.02*max(R,[],'all');
for i=4+(s+1)/2:m-3-(s-1)/2
   for j=4+(s+1)/2:n-3-(s-1)/2
      if (R(i,j)==max(R(i-(s-1)/2:i+(s-1)/2,j-(s-1)/2:j+(s-1)/2),[],'all')) && (R(i,j)>0) 
         Re(i,j)=1;
      end
   end
end
[x,y]=find(Re==1);
figure(1)
imshow(img);
hold on
plot(y,x,'r+','MarkerSize',20);
hold off
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

结果如图所示：

局部角点检测结果（符合算法预期）如图所示：