1、SUSAN边缘检测计算步骤
（1）在图像上放置一个37个像素的圆形模板，模板在图像上滑动，依次比较模板内各个像素点的灰度与模板核的灰度，判断是否属于USAN区域。判别函数如下：
其中，
r
⃗
0
{{\vec{r}}_{0}}
r
0.表示二维图像中核心点的位置，
r
⃗
\vec{r}
r
表示模板中其他点的位置，
I
(
r
⃗
0
)
I({{\vec{r}}_{0}})
I(r
0)表示图像在
r
⃗
0
{{\vec{r}}_{0}}
r
0处的像素值，
I
(
r
⃗
)
I(\vec{r})
I(r
)表示图像在
r
⃗
\vec{r}
r
处的像素值，
t
t
t是表示亮度插值的一个门限值，决定了USAN区域各点之间最大的亮度差值，
c
c
c是一个用来比较输出的函数。也可以用平滑的县来代替这种直接的分割方式，这样可以获得更稳定而敏感的结果，虽然计算复杂，但是可以通过查表来获得较快的速度。公式如下：
c
(
r
⃗
,
r
⃗
0
)
=
e
−
(
I
(
r
⃗
)
−
I
(
r
⃗
0
)
t
)
6
c(\vec{r},{{\vec{r}}_{0}})={{e}^{-{{\left( \frac{I(\vec{r})-I({{{\vec{r}}}_{0}})}{t} \right)}^{6}}}}
c(r
,r
0)=e−(tI(r
)−I(r
0))6
（2）统计圆形模板中和核心点有相似亮度值的像素值个数
n
(
r
0
)
n({{r}_{0}})
n(r0)
n
(
r
⃗
0
)
=
∑
r
⃗
∈
D
(
r
⃗
0
)
c
(
r
⃗
,
r
⃗
0
)
n({{\vec{r}}_{0}})=\sum\limits_{\vec{r}\in D({{{\vec{r}}}_{0}})}{c(\vec{r},{{{\vec{r}}}_{0}})}
n(r
0)=r
∈D(r
0)∑c(r
,r
0)其中，
D
(
r
⃗
0
)
D({{\vec{r}}_{0}})
D(r
0)是以
r
⃗
0
{{\vec{r}}_{0}}
r
0为中心的圆形模板区域，
n
n
n是USAN区域中像素的个数。
（3）使用如下角点响应函数，若某个像素点的USAN值小于某一特定阈值，则该点被认为是初始角点。将
n
n
n同一个固定的阈值
g
g
g比较（一般设置为最大与中心相似点数
n
max
{{n}_{\max }}
nmax的0.75倍左右），初始的边缘响应可以用下面的等式计算：
（4）对初始边缘进行非极值抑制来求得最后的边缘。（理论上在垂直于边缘的方向上进行非极大值抑制，以细化边缘或拟合亚像素边缘，但这里只进行最简单的非极大值抑制）
2、Matlab实现
%% SUSAN边缘检测算法（不计算边缘方向，只进行最简单的非极大值抑制）
close all
clear
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=45; % 阈值
g=2*37/3;
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);
c=0; % USAN计数
%% 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
c=c+1;
end
end
end
end
%% 3、计算边缘响应，确定初始边缘
if c<g
R(i,j)=g-c;
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)>Rmax)
Re(i,j)=1;
end
end
end
[x,y]=find(Re==1);
figure(1)
imshow(img);
hold on
plot(y,x,'r+','MarkerSize',1);
hold off
结果如图所示：
3、OpenCV实现
// SUSAN边缘检测
#include <iostream>
#include <opencv2\opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat img = imread("board.jpg");
Mat img_gray, tmp, R;
int t = 25; // USAN阈值
if (img.channels() > 1)
cvtColor(img, img_gray, COLOR_BGR2GRAY);
else
img_gray = img.clone();
img_gray.convertTo(img_gray, CV_32FC1);
Size s = img_gray.size();
int m = s.height, n = s.width, usan;
double g = 2 * 37 / 3.0;
float *p1, *p2, *p3;
R = Mat::zeros(s, CV_32FC1);
// 1、圆形模板在图像上滑动，依次比较模板内各个像素点的灰度与模板核的灰度，判断是否属于USAN区域
for (int i = 3 ; i < m - 3; i++)
{
p1 = img_gray.ptr<float>(i);
p2 = R.ptr<float>(i);
for (int j = 3; j < n - 3; j++)
{
tmp = img_gray(Rect(j - 3, i - 3, 7, 7));
usan = 0; // USAN计数
// 2、统计圆形模板中和核心点有相似亮度值的像素值个数
for (int p = 0; p < 7; p++)
{
p3 = tmp.ptr<float>(p);
for (int q = 0; q < 7; q++)
if ((p - 3)*(p - 3) + (q - 3)*(q - 3) < 12)
if (fabsf(p1[j] - p3[q]) < t)
usan = usan + 1;
}
// 3、计算边缘响应，确定初始边缘
if (usan < g)
p2[j] = g - usan;
}
}
// 4、非极大值抑制
int s1 = 5, s2; // 邻域
Mat Re = Mat::zeros(s, CV_32FC1);
double Rmax, rmax;
minMaxIdx(R, NULL, &Rmax, NULL, NULL);
Rmax = 0.02*Rmax;
s2 = (s1 - 1) / 2;
for (int i = 3 + s2; i < m - 3 - s2; i++)
{
p2 = R.ptr<float>(i);
p3 = Re.ptr<float>(i);
for (int j = 3 + s2; j < n - 3 - s2; j++)
{
minMaxIdx(R(Rect(j - s2, i - s2, s1, s1)), NULL, &rmax, NULL, NULL);
if ((p2[j] == rmax) && (p2[j] > Rmax))
{
p3[j] = 255;
drawMarker(img, Point(j, i), Scalar(0, 0, 255), MARKER_CROSS, 10, 1);
}
}
}
namedWindow("img", 2);
imshow("img", img);
waitKey(0);
return 0;
}
结果如图所示：
