Harris、SUSAN角点检测Matlab源代码

　　角点检测主要运用于运动检测、图像匹配、视频跟踪、三维建模、目标识别。常用的几种基于灰度变化的算法有K.R，Harris，KLT，SUSAN，MORAVEC。

　　　　　　

角点检测=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪灰度图⎧⎩⎨梯度模版混合二值图轮廓图

　 由于以下缺点，Moverac渐渐被Harris取代　　
　

Moravec角点检测缺陷⎧⎩⎨强度值计算各向异性窗口函数混有噪声相似判别太简单min(SSD)

1.Harris角点检测　　

　　Harris角点检测的思路较为清晰，其原理主要是找取双边变化最大的点，主要在于理解其响应值（评价灰度变化情况）函数的构造。当然虽说不难，但也是找了一大堆资料才看明白的。原本是打算自己写一遍，但发现有几篇博客里面介绍的已经不能再详细了，以下是传送门：http://www.cnblogs.com/ronny/p/4009425.html
　　以下介绍一下Matlab实现的Harris角点检测，具体步骤都标注好了，大家可以对照着博客中的原理进行理解、
　　

clear; % 读取图像 grayImage = imread('C:\Users\Administrator\Desktop\3.bmp');     % 转化为灰度图像%grayImage=rgb2gray(srcImage);% 求取图像宽高[ImageHeight,ImageWidth]=size(grayImage);% 显示原始灰度图%imshow(ori_im);% 方法1：x方向梯度算子(用于Harris角点提取算法) fx = [-2 -1 0 1 2];% 方法2：x方向高斯卷积核(高斯尺度空间的多尺度优化)%fx = [5 0 -5;8 0 -8;5 0 -5]; % x方向滤波微分Ix = filter2(fx,grayImage);            % 显示x方向滤波图%figure;imshow(Ix);% 方法1：y方向梯度算子(用于Harris角点提取算法) fy = [-2;-1;0;1;2];      % 方法2：y方向高斯卷积核(高斯尺度空间的多尺度优化) %fy = [5 8 5;0 0 0;-5 -8 -5];% y方向滤波微分Iy = filter2(fy,grayImage);% 显示y方向滤波图figure;imshow(Iy);%（相关参数说明见harris理论，文中前面有链接）Ix2 = Ix.^2; Iy2 = Iy.^2; Ixy = Ix.*Iy; % 产生7*7的高斯窗函数，sigma=2,用于窗口的高斯平滑 w= fspecial('gaussian',[7 7],2);      Ix2 = filter2(w,Ix2); Iy2 = filter2(w,Iy2); Ixy = filter2(w,Ixy); % 纪录角点位置，角点处值为1 corner = zeros(ImageHeight,ImageWidth);% 图像中最大的响应值 Rmax = 0;                               % 计算各点的响应值R = zeros(ImageHeight,ImageWidth); for i = 1:ImageHeight     for j = 1:ImageWidth         M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];                     R(i,j) = det(M)-0.06*(trace(M))^2;                              if R(i,j) > Rmax             Rmax = R(i,j);         end;     end; end; % 角点个数cnt = 0; % 进行非极大抑制，窗口大小3*3 for i = 2:ImageHeight-1     for j = 2:ImageWidth-1         if ( R(i,j) > 0.01*Rmax && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1))             corner(i,j) = 1;             cnt = cnt+1;         end;     end; end; [upix, vpix] = find(corner == 1); %角点个数cnt      %绘制角点figure;imshow(grayImage)hold on;plot(vpix,upix,'r.');      

得到结果如下图

2.SUSAN角点检测

SUSAN算法容易把边缘也提取出来，但对噪声的容忍度比较强。
原理参考http://blog.csdn.net/tostq/article/details/49305615

clear;clc;% 读取图像Image=imread('C:\Users\Administrator\Desktop\1.bmp');% 转化为灰度图像%Image=rgb2gray(Image);% 显示图像%imshow(Image);% 获取图像高宽（行烈）[ImageHeight,ImageWidth]=size(Image);% 这一步没太大必要%Image=double(Image);% 判断灰度相近的阈值threshold=45;  % 当前像素和窗体内像素差别在t以下的个数，即相似的个数usan=[];% 计算以像素为中心的窗体内包含的% 包含37个像素的圆窗口，面积为12*pi=37，因此是以sqrt（12）为半径的原% 没有在外围扩展图像，最终图像会缩小for i=4:ImageHeight-3            for j=4:ImageWidth-3         %从原图中截取7*7的区域再在其中挑选圆窗        tmp=Image(i-3:i+3,j-3:j+3);          %c表示灰度值相近的程度，越大越相近        c=0;        for p=1:7           for q=1:7               %在7*7的区域中选取圆窗包含的像素                if (p-4)^2+(q-4)^2<=12                     %usan(k)=usan(k)+exp(-(((img(i,j)-tmp(p,q))/t)^6));                    %判断灰度是否相近，t是自己设置的                    if abs(Image(i,j)-tmp(p,q))<threshold                          c=c+1;                    end                end           end        end        usan=[usan c];   endend%相当于进一步调整阈值，在threshold的基础上进一步减少角点个数g=2*max(usan)/3;for i=1:length(usan)   if usan(i)<g        usan(i)=g-usan(i);   else       usan(i)=0;   endend% 由于usan是一维的，所以要重新变换为二维，对应图像位置imgn=reshape(usan,[ImageWidth-6,ImageHeight-6])';%figure;%imshow(imgn)%非极大抑制[m n]=size(imgn);re=zeros(m,n);for i=2:m-1   for j=2:n-1         if imgn(i,j)>max([max(imgn(i-1,j-1:j+1)) imgn(i,j-1) imgn(i,j+1) max(imgn(i+1,j-1:j+1))]);            re(i,j)=1;        else            re(i,j)=0;        end   endendfigure;imshow(Image)hold on;[x,y]=find(re==1);plot(y,x,'*')

“`c

效果图如下