SUSAN算子——边缘检测

引言

SUSAN算子很好听的一个名字，其实SUSAN算子除了名字好听外，她还很实用，而且也好用，SUSAN的全名是：Smallest Univalue Segment Assimilating Nucleus)。它是一种很有特色高效的边缘和角点检测算子，它不仅可以检测图像目标的边界点，而且能够较Robust地检测目标的角点。并且具有结构保留的降噪功能。

基本SUSAN原理

先借助如图所示来解释检测的原理，其中图片是白色背景，有一个颜色比较暗淡的矩形(dark area)。用一个园形模板在图像上移动，若模板内的像素灰度与模板中心的像素（被称为核Nucleus）灰度值小于一定的阈值，则认为该点与核Nucleus具有相同的灰度，满足该条件的像素组成的区域就称为USAN（Univalue Segment Assimilating Nucleus）。在图片上有5个圆形区域。圆形区域表示的是掩码区域。把圆形区域内的每一个位置的像素值与圆心处的像素值相比较，那么圆中的的像素可以分为两类，一类是像素值与圆心处的像素值相近的，另一类是像素值与圆心的处的像素值相差比较大的。

如果将模板中各个像素的灰度都与模板中心的核像素的灰度进行比较，那么就会发现总有一部分模板区域和灰度与核像素的灰度相同或相似，这部分区域可以称为USAN(Univalue Segment Assimilating Nuclues).USAN区域包含很多与图像结构有关的信息。利用这种区域的尺寸、重心、二阶矩的分析，可以得到图像中的角点，边缘等信息。从上图所示，当核像素处在图像中的灰度一致区域时，USAN的面积会达到最大。第e个模板就是属于这种情况。

SUSAN（Smallest Univalue Segment Assimilating Nuclues）进行角点检测时，遵循了常规的思路：使用一个窗口在图像上逐点滑动，在每一个位置上计算出一个角点量，再进行局部极大值抑制，得到最终的角点。我们这里使用的窗口是圆形窗口，最小的窗口是3*3的，本文中使用的是37个像素的圆形窗口，如图：

阈值的分析

在角点检测中，有两种类型的阈值，一种用来约束角点的数量，另一种用来约束角点的质量。当然，一个阈值不能完全做到只影响质量或数量，只是会有一个侧重点。那么，SUSAN中的两个阈值t和g在角点检测中起到一个什么样的作用呢？很明显，阈值g是角点质量的。尽管也会影响数量，但是相对来说更侧重于影响质量（角点的形状）。例如，g值减小，那么SUSAN会更加侧重于检测到更加“尖锐”的角点。所以，可以更加自己的实际需求来确定阈值g。而阈值t，是角点的数量。当t减小时，会检测到更多的角点。所以，阈值t可以在不影响角点质量的情况下，控制检测到的角点的数量。在大多数情况下，设t为25比较合适，如果图像的对比度比较低，可以修改t值以适应变化。
通过上面对a、b、c、d、e等几个圆形模板的USAN值的分析，当模板的中心位于角点处时，USAN的值最小。
下面简单叙述下利用SUSAN算子检测角点的步骤：

利用圆形模板遍历图像，计算每点处的USAN值。
设置一阈值g，一般取值为1/2(Max(n)，也即取值为USAN最大值的一半，进行阈值化，得到角点响应。
使用非极大值抑制来寻找角点。

通过上面的方式得到的角点，存在很大伪角点。为了去除伪角点，SUSAN算子可以由以下方法实现：

计算USAN区域的重心，然后计算重心和模板中心的距离，如果距离较小则不是正确的角点；
判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素，如果属于那么这个模板中心的点就是角点。

综上所述，我们基本知道SUSAN角点检测算法。

SUSAN角点检测代码

1.#include <stdio.h>
2.#include <cv.h>
3.#include <highgui.h>
4.
5.#define max_corners 300
6.int main( int argc, char** argv )
7.{
8.    int cornerCount=max_corners;
9.    CvPoint2D32f corners[max_corners];
10.    double qualityLevel = 0.05;
11.    double minDistance = 5;
12.    IplImage *srcImage = 0, *grayImage = 0, *corners1 = 0, *corners2 = 0;
13.    int i;
14.    CvScalar color = CV_RGB(255,0,0);
15.    cvNamedWindow( "image", 2);
16.    srcImage = cvLoadImage("test.jpg", 1);
17.    grayImage = cvCreateImage(cvGetSize(srcImage), IPL_DEPTH_8U, 1);
18.
19.    cvCvtColor(srcImage, grayImage, CV_BGR2GRAY);
20.    corners1= cvCreateImage(cvGetSize(srcImage), IPL_DEPTH_32F, 1);
21.    corners2= cvCreateImage(cvGetSize(srcImage),IPL_DEPTH_32F, 1);
22.    cvGoodFeaturesToTrack (grayImage, corners1, corners2, corners,&cornerCount,
23.                               qualityLevel, minDistance, 0);
24.    printf("num corners found: %d\n", cornerCount);
25.    if(cornerCount>0)
26.    {
27.        for (i = 0; i < cornerCount; ++i)
28.        {
29.                   cvCircle(srcImage,cvPoint((int)(corners[i].x),(int)(corners[i].y)),
30.                                 3, color, 1, CV_AA, 0);
31.        }
32.    }
33.    cvShowImage("image", srcImage);
34.    cvWaitKey(0);
35.    cvReleaseImage(&srcImage);
36.    cvReleaseImage(&grayImage);
37.    cvReleaseImage(&corners1);
38.    cvReleaseImage(&corners2);
39.
40.    return 0;
41.}

#include <stdio.h>
#include <cv.h>
#include <highgui.h>#define max_corners 300
int main( int argc, char** argv )
{int cornerCount=max_corners;CvPoint2D32f corners[max_corners];double qualityLevel = 0.05;double minDistance = 5;IplImage *srcImage = 0, *grayImage = 0, *corners1 = 0, *corners2 = 0;int i;CvScalar color = CV_RGB(255,0,0);cvNamedWindow( "image", 2); srcImage = cvLoadImage("test.jpg", 1);grayImage = cvCreateImage(cvGetSize(srcImage), IPL_DEPTH_8U, 1);cvCvtColor(srcImage, grayImage, CV_BGR2GRAY);corners1= cvCreateImage(cvGetSize(srcImage), IPL_DEPTH_32F, 1);corners2= cvCreateImage(cvGetSize(srcImage),IPL_DEPTH_32F, 1);cvGoodFeaturesToTrack (grayImage, corners1, corners2, corners,&cornerCount,qualityLevel, minDistance, 0);printf("num corners found: %d\n", cornerCount);if(cornerCount>0){for (i = 0; i < cornerCount; ++i){cvCircle(srcImage,cvPoint((int)(corners[i].x),(int)(corners[i].y)),3, color, 1, CV_AA, 0);}}cvShowImage("image", srcImage);cvWaitKey(0);cvReleaseImage(&srcImage);cvReleaseImage(&grayImage);cvReleaseImage(&corners1);cvReleaseImage(&corners2);return 0;
}

测试结果输出

源图像

SUSAN算子处理输出结果图像

在些，我们给出SUSAN角点测测另一种应用。

检测动态匹配轮廓的SUSAN代码：

1.#include<math.h>
2.#include<stdlib.h>
3.#include<stdio.h>
4.#include "cv.h"
5.#include "highgui.h"
6.
7.int main( int argc, char** argv )
8.{
9.    int height ,width ,step ,channels ;
10.    int i,j,k,same ,max,min,sum;
11.    float thresh;
12.    uchar *data0,*data1 ;
13.    //char *filename="result.bmp";
14.    IplImage* Img, *nimg;
15.
16.    Img = cvLoadImage( "test.jpg",0);
17.    cvNamedWindow( "Image", 2);
18.    cvShowImage( "Image", Img );
19.
20.    nimg = cvCreateImage(cvGetSize(Img),8,1);
21.
22.    height    = Img->height;
23.    width     = Img->width;
24.    step = Img->widthStep/sizeof(uchar);
25.    channels = Img->nChannels;
26.    data0   = (uchar*)Img->imageData;
27.    data1 =    (uchar*)nimg->imageData;
28.
29.    printf("Processing a %d X %d image with %d channels\n",width,height,channels);
30.    int OffSetX[37] =
31.    { -1, 0, 1,
32.    -2,-1, 0, 1, 2,
33.    -3,-2,-1, 0, 1, 2, 3,
34.    -3,-2,-1, 0, 1, 2, 3,
35.    -3,-2,-1, 0, 1, 2, 3,
36.    -2,-1, 0, 1, 2,
37.    -1, 0, 1 };
38.    int OffSetY[37] =
39.    {
40.        -3,-3,-3,
41.        -2,-2,-2,-2,-2,
42.        -1,-1,-1,-1,-1,-1,-1,
43.        0, 0, 0, 0, 0, 0, 0,
44.        1, 1, 1, 1, 1, 1, 1,
45.        2, 2, 2, 2, 2,
46.        3, 3, 3
47.    };
48.
49.    max = min = data0[0];
50.
51.    for(i=3;i<height-3;i++)
52.        for(j=3;j<width-3;j++)
53.        {
54.            same =0;
55.            sum = 0;
56.            for(k=0;k<37;k++)
57.            {
58.                sum+=data0[(i+OffSetY[k])*step+(j+OffSetX[k])];
59.                thresh = (float)sum/37;
60.                float data_fabs;
61.                data_fabs= (float)(data0[(i+OffSetY[k])*step+(j+OffSetX[k])]-data0[i*step+j]);
62.
63.                if(fabs( data_fabs)<=thresh)
64.                    same++;
65.            }
66.
67.            if(same<18)
68.                nimg->imageData[i*step+j] = 255;
69.            else
70.                nimg->imageData[i*step+j]  = 0;
71.
72.            printf("same = %d\n", same);
73.
74.        }
75.
76.        cvNamedWindow( "Image", 2);
77.        cvShowImage( "Image", nimg );
78.        cvWaitKey(0);
79.        cvDestroyWindow( "Image" );
80.        cvReleaseImage( &Img );
81.        cvReleaseImage( &nimg );
82.        return 0;
83.}

#include<math.h>
#include<stdlib.h>
#include<stdio.h>
#include "cv.h"
#include "highgui.h"int main( int argc, char** argv )
{int height ,width ,step ,channels ;int i,j,k,same ,max,min,sum;float thresh;uchar *data0,*data1 ;//char *filename="result.bmp";IplImage* Img, *nimg; Img = cvLoadImage( "test.jpg",0);cvNamedWindow( "Image", 2); cvShowImage( "Image", Img ); nimg = cvCreateImage(cvGetSize(Img),8,1);height    = Img->height;width     = Img->width;step = Img->widthStep/sizeof(uchar);channels = Img->nChannels;data0   = (uchar*)Img->imageData;data1 =    (uchar*)nimg->imageData;printf("Processing a %d X %d image with %d channels\n",width,height,channels);int OffSetX[37] = { -1, 0, 1,-2,-1, 0, 1, 2,-3,-2,-1, 0, 1, 2, 3,-3,-2,-1, 0, 1, 2, 3,-3,-2,-1, 0, 1, 2, 3,-2,-1, 0, 1, 2,-1, 0, 1 };int OffSetY[37] = { -3,-3,-3,-2,-2,-2,-2,-2,-1,-1,-1,-1,-1,-1,-1,0, 0, 0, 0, 0, 0, 0,1, 1, 1, 1, 1, 1, 1,2, 2, 2, 2, 2,3, 3, 3 };max = min = data0[0];for(i=3;i<height-3;i++)for(j=3;j<width-3;j++){same =0;sum = 0;for(k=0;k<37;k++){sum+=data0[(i+OffSetY[k])*step+(j+OffSetX[k])];thresh = (float)sum/37;float data_fabs;data_fabs= (float)(data0[(i+OffSetY[k])*step+(j+OffSetX[k])]-data0[i*step+j]);if(fabs( data_fabs)<=thresh)same++;}if(same<18)nimg->imageData[i*step+j] = 255;elsenimg->imageData[i*step+j]  = 0;printf("same = %d\n", same);}cvNamedWindow( "Image", 2);cvShowImage( "Image", nimg ); cvWaitKey(0); cvDestroyWindow( "Image" );cvReleaseImage( &Img ); cvReleaseImage( &nimg );return 0;
}

测试输出结果

源图像

输出结果

补充事项

SUSAN算子是一个原理简单、易于了解的算子。由于其指数基于对周边象素的灰度比较，完全不涉及梯度的运算，因此其抗噪声能力很强，运算量也比较小；同时，SUSAN算子还是一个各向同性的算子；最后，通过控制参数t和g，可以根据具体情况很容易地对不同对比度、不同形状的图像通过设置恰当的t和g进行控制。比如图像的对比度较大，则可选取较大的t值，而图像的对比度较小，则可选取较小的t值。当SUSAN进行角点检测时，在某些情况会产生误报：

情况一：在点B的处的USAN区域如右边的白色区域所示，呈带状分布在中间，并且是单像素的带状。此时计算出的角点相应量会比较大，会误认为是角点，而实际上是边缘点。解决方法：计算出USAN区域的重心，如果重心与圆心的距离大于一定的阈值才认为是角点。
情况二：多个噪声点或是复杂的结构。例如单点噪声，其计算出来的角点量很大，会误认为是角点。解决方法：计算出USAN区域的重心，从中心指向USAN区域重心的直线上的所有像素都必须是USAN的一部分，否则认为不是角点。

总之，SUSAN算子是一个非常难得的算子，不仅具有很好的边缘检测性能；而且对角点检测也具有很好的效果。

参考文献

[1]SUSAN Low Level Image Processing.