Harris、SUSAN角点检测Matlab源代码

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

　　　　　　

$$\textbf{角点检测}= \begin{cases} \textbf{灰度图} \begin{cases} \textbf{梯度}\\\textbf{模版}\\\textbf{混合} \end{cases}\\\\ \textbf{二值图}\\\\\\\\\textbf{轮廓图} \end{cases}$$

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

$$\textbf{Moravec角点检测缺陷} \begin{cases} \textbf{强度值计算各向异性}\\ \textbf{窗口函数混有噪声}\\ \textbf{相似判别太简单$min(\small{SSD})$} \end{cases}$$

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];
   end
end
%相当于进一步调整阈值，在threshold的基础上进一步减少角点个数
g=2*max(usan)/3;
for i=1:length(usan)
   if usan(i)<g 
       usan(i)=g-usan(i);
   else
       usan(i)=0;
   end
end
% 由于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
   end
end
figure;
imshow(Image)
hold on;
[x,y]=find(re==1);
plot(y,x,'*')

“`c

效果图如下