引言
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个像素的圆形窗口,如图:
阈值的分析
通过上面对a、b、c、d、e等几个圆形模板的USAN值的分析,当模板的中心位于角点处时,USAN的值最小。
下面简单叙述下利用SUSAN算子检测角点的步骤:
- 利用圆形模板遍历图像,计算每点处的USAN值。
- 设置一阈值g,一般取值为1/2(Max(n), 也即取值为USAN最大值的一半,进行阈值化,得到角点响应。
- 使用非极大值抑制来寻找角点。
- 计算USAN区域的重心,然后计算重心和模板中心的距离,如果距离较小则不是正确的角点;
- 判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素,如果属于那么这个模板中心的点就是角点。
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代码:
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;
else
nimg->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;
}
测试输出结果
补充事项
- 情况一:在点B的处的USAN区域如右边的白色区域所示,呈带状分布在中间,并且是单像素的带状。此时计算出的角点相应量会比较大,会误认为是角点,而实际上是边缘点。解决方法:计算出USAN区域的重心,如果重心与圆心的距离大于一定的阈值才认为是角点。
- 情况二:多个噪声点或是复杂的结构。例如单点噪声,其计算出来的角点量很大,会误认为是角点。解决方法:计算出USAN区域的重心,从中心指向USAN区域重心的直线上的所有像素都必须是USAN的一部分,否则认为不是角点。