SUSAN边缘检测算法，及其Matlab和OpenCV实现

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.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;
}

结果如图所示：