SUSAN角点检测算法，及其Matlab实现

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

结果如图所示：

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