图像特征之SUSAN角点

由于局部梯度的方法对噪声影响比较敏感而且计算量大，英国学者Smith和Brady提出一种基于形态学的角点特征检测方法。这种方法是一种基于灰度的特征点获取方法，适用于图像中的边缘检测，角点检测且计算速度快，适用于实时图像处理。

SUSAN角点检测的原理是：用一个固定半径的圆形模板在图像上滑动，该模板中心像素点称为核。若模板内其他点的灰度值与核的灰度值之差小于某一个阈值，则认为该点与核具有相似的灰度，所以满足这样条件的像素组成的区域称为核值相似区(Univalue Segment Assimilating Nucleus, USAN)。通过这种方式生成的USAN区域包含了该像素点的图像结构信息。USAN区域越大，模板内相似像素点较多，一般为图像平滑区域，USAN区域大小中等，一般为图像边缘区域，USAN区域很小，一般为角点。通过计算USAN区域的大小，就可以知道该点是否为角点。

自己按照SUSAN算子的理论编写了一段代码，是在opencv2.3.1+vs2008上实现的，若有错误，还望指正：

#include "stdafx.h"     #include <opencv2/opencv.hpp>  #include "highgui.h"  #include <math.h>  typedef  unsigned long uint32;typedef  unsigned int  uint16;typedef  unsigned char uint8;#define THRESHOLD  20#define RADIUS  3IplImage *src_gray1, *src_gray2, *src_gray3;    IplImage* src_img, *dst_img;CvMat  *NumR0,*Rr0;int MaxNumR0=0;void AllocateImage(IplImage* I)   //给图像分配大小    {        CvSize sz   = cvGetSize(I);        dst_img   = cvCreateImage( sz, IPL_DEPTH_8U, 1);cvSetZero(dst_img);    src_gray1 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    //原图的三个通道        src_gray2 = cvCreateImage( sz, IPL_DEPTH_8U, 1);        src_gray3 = cvCreateImage( sz, IPL_DEPTH_8U, 1);  NumR0     = cvCreateMat(sz.height,sz.width,CV_16SC1);Rr0       = cvCreateMat(sz.height,sz.width,CV_16SC1);cvSetZero(NumR0);cvSetZero(Rr0);    }  void DeallocateImage(){cvReleaseImage(&src_img); cvReleaseImage(&dst_img);     cvReleaseImage(&src_gray1);      cvReleaseImage(&src_gray2);      cvReleaseImage(&src_gray3);  cvReleaseMat(&NumR0);} void SUSAN_check(IplImage* I, IplImage* dst, int r)                    //SUSAN角点检测函数{int i,j,x,y;int num,MaxNum=0,g,CValue;int R0Value,temp;for( i=r; i<(I->height-r); i++ )               //计算n(r0),模板内图像USAN的像元数量{for( j=r; j<(I->width-r); j++ ){R0Value=cvGetReal2D(I,i,j);num=0;for( y=i-r; (y-i<=r)&&(y-i>=-r);y++ ){for( x=j-r;(x-j<=r)&&(x-j>=-r);x++ ){if(((x-j)*(x-j)+(y-i)*(y-i))<=r*r&&(x!=j||y!=i))                //中心点圆形邻域内{temp=cvGetReal2D(I,y,x);if( abs(R0Value-temp) <= THRESHOLD  )num++;}}}cvSetReal2D(NumR0,i,j,num);if(num>MaxNum)  MaxNum = num;}}g=3*MaxNum/4;for( i=r; i<(I->height-r); i++ )               //计算R(r0),USAN特征图像{for( j=r; j<(I->width-r); j++ ){temp = cvGetReal2D(NumR0,i,j);if( temp<g )cvSetReal2D(Rr0,i,j,(g-temp));}}for( i=r; i<(I->height-r); i++ )             //非极大抑制{for( j=r; j<(I->width-r); j++ ){CValue = cvGetReal2D(Rr0,i,j);       //USAN特征图像中心点值num=0;if(CValue!=0){for( y=i-r; (y-i<=r)&&(y-i>=-r);y++ )         //查看周围半径r圆形邻域是否为极大值{for( x=j-r;(x-j<=r)&&(x-j>=-r);x++ ){if(((x-j)*(x-j)+(y-i)*(y-i))<=r*r&&(x!=j||y!=i))                //中心点圆形邻域内{temp = cvGetReal2D(Rr0,y,x);if(j==27&&i==129){MaxNum=MaxNum;}if(CValue>temp)num++;}}}if(num==MaxNum){cvCircle(src_img,cvPoint(j,i),3,cvScalar(0,0,255),1,8,0);//cvSet2D(src_img,i,j,cvScalar(0,0,255));*(dst->imageData+i*dst->widthStep+j) =255;}}}}}int _tmain(int argc, _TCHAR* argv[]){//src_img  = cvLoadImage("角点.bmp");//src_img  = cvLoadImage("5.2.09.pgm");src_img  = cvLoadImage("7.1.02.pgm");AllocateImage(src_img);cvSplit( src_img, src_gray1, src_gray2, src_gray3, 0); cvNamedWindow("my picture",CV_WINDOW_AUTOSIZE); cvNamedWindow("my dst",CV_WINDOW_AUTOSIZE);  cvShowImage("my picture",src_img); SUSAN_check(src_gray1, dst_img, RADIUS) ;cvShowImage("my dst",src_img); cvWaitKey(0);      DeallocateImage();cvDestroyWindow("my picture");  cvDestroyWindow("my dst"); return 0;}