关于Harris角点检测的实现

harris角点检测算法非常有名，网上也有一堆源代码，这里自己贴下自己找到的一些资料，以备将来查看。

1.matlab中有内置的corner函数可以实现harris焦点的提取
2.网上也有很多博文中贴出了源代码，这里我们都做了测试，至于网上源代码中 的这段，研究了好久都没有研究出来，为啥这个高斯函数的一阶微分会是这个结果呢，不知道呀2333333333333
如果有大神知道，请劳烦告知一下，谢谢了

fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)

3.自己参考大神们的源码自己也写了一个附在其中，其中也遇到一些问题，比如imfilter， conv2， filter2之间的关系，简单来讲conv2是卷积操作，filter2是相关操作，imfilter理论上都能实现以上两者的效果，在我们放出的代码中的test部分就是对matlab论坛中一个帖子的测试代码了。

% -----------------【harris corner detect】--------------
% author  :   zhyh2010
% date    :   20150511
% target  :   harris corner detect
% --------------------------【end】----------------------


function main
    % clc
    % close all
    % clear all

    %harris_matlab_inside;
%     test
    myharris_1;
%     harris_refer1;
%     harris_refer2;
end

% -----------------【matlab 内置corner函数进行角点检测】--------------
% method ：  corner  只支持单通道
% --------------------------【end】----------------------
function harris_matlab_inside
    src = imread('1.jpg');
    figure, imshow(src)
    title('original pic')

    src_corner = corner(rgb2gray(src), 'Harris', 1000);
    figure, imshow(src);
    hold on
    plot(src_corner(:, 1), src_corner(:, 2), 'b*')
    hold off
end

% -----------------【matlab conv2 filter2 imfilter 区别】--------------
% 检验博文  http://www.ilovematlab.cn/thread-293710-1-1.html
% --------------------------【end】----------------------
function test
    A=[4 3 1 2;0 1 1 3;5 2 0 0]
    B=[1 2 3;0 -1 2;1 1 0]        % kernel
    C = conv2(A, B, 'full')
    D = filter2(B, A, 'full')
    E_1 = imfilter(A, B, 0, 'full', 'conv')
    E_2 = imfilter(A, B, 0, 'full', 'corr')
end


% -----------------【myharris_1】--------------
% 根据原理自行实现 harris 角点检测
% --------------------------【end】----------------------
function myharris_1
     src = imread('1.jpg');
     img = rgb2gray(src);

     h_x = [-2 -1 0 1 2];
%      Ix = imfilter(img, h_x, 0, 'full', 'corr');
% %      Ix = imfilter(img, h_x, 0, 'same', 'corr');
%      Ix = Ix(2:end-1, 2:end -1);
    Ix = filter2(h_x, img);

     h_y = [-2 -1 0 1 2]';
% %      Iy = imfilter(img, h_y, 0, 'same', 'corr');
%      
%      Iy = imfilter(img, h_y, 0, 'full', 'corr');
%      Iy = Iy(2:end-1, 2:end -1);
    Iy = filter2(h_y, img);

     Ix2 = Ix.^2;
     Iy2 = Iy.^2;
     Ixy = Ix.*Iy;

     kernel_size = 3;
     kernel_sigma = 2;
     h = fspecial('gaussian', kernel_size, kernel_sigma);

     Ix2 = filter2(h, Ix2);
     Iy2 = filter2(h, Iy2);
     Ixy = filter2(h, Ixy);

     pic_width = size(img, 2);
     pic_height = size(img, 1);

     alpha = 0.05;
     Judge_standard = zeros(size(img));
     result = zeros(size(img));
     for i = 1 : pic_height
         for j = 1 : pic_width
            M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
            Judge_standard(i, j) = det(M) - alpha * trace(M)^2;
         end
     end


     MaxJS = max(Judge_standard(:));
     R = Judge_standard;
     threshold = 0.1 * MaxJS;
     for i = 2 : pic_height - 1
         for j = 2 : pic_width - 1
             if (Judge_standard(i, j) > threshold) ...
                   && (bigger_around(Judge_standard, i, j))
                 result(i, j) = 1;                
             end
         end
     end

    [posx, posy] = find(result == 1);
    fprintf('---------------find %d corners---------------\n', length(posx));
    figure, imshow(src)
    hold on
    plot(posy, posx, 'b*')
    hold off

end

% -----------------【myharris_1 接口子函数】--------------
% 本质  3 X 3 局部抑制操作
% isbigger = bigger_around(A, i, j)
%   A  :  存储角点评判标准的值
%   i，j ： 当前坐标 
%   attention ： 2 <= i <= size(A, 1) - 1
%                2 <= j <= size(A, 2) - 1
% --------------------------【end】----------------------
function isbigger = bigger_around(A, i, j)
    assert(i >= 2 && i <= size(A, 1) - 1);
    assert(j >= 2 && j <= size(A, 2) - 1);

    cond = A(i, j) > A(i - 1, j - 1);
    cond = cond && (A(i, j) > A(i - 1, j));
    cond = cond && (A(i, j) > A(i - 1, j + 1));
    cond = cond && (A(i, j) > A(i, j - 1));
%     cond = cond && (A(i, j) > A(i, j));
    cond = cond && (A(i, j) > A(i, j + 1));
    cond = cond && (A(i, j) > A(i + 1, j - 1));
    cond = cond && (A(i, j) > A(i + 1, j));
    cond = cond && (A(i, j) > A(i + 1, j + 1));

    isbigger = cond;
end

% -----------------【别人的harris算法实现1】--------------
% 参考地址： http://blog.sina.com.cn/s/blog_6c73951701017bwi.html
% --------------------------【end】----------------------
function harris_refer1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   Harris角点提取算法                                          %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    clear;
    %filename = 'Lena.jpg';
    filename='1.jpg';
    X = imread(filename);     % 读取图像
    % imshow(X);
    Info = imfinfo(filename); %获取图像相关信息
    if (Info.BitDepth > 8)
        f = rgb2gray(X);
    end
    %《基于特征点的图像配准与拼接技术研究》
    %计算图像亮度f(x,y)在点(x,y)处的梯度-----------------------------------------------
    fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
    ori_im = double(f) / 255;                   %unit8转化为64为双精度double64
    %fx = [-2 -1 0 1 2];                     % x方向梯度算子(用于Harris角点提取算法)
    Ix = filter2(fx, ori_im);                % x方向滤波
    fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
    %fy = [-2; -1; 0; 1; 2];                     % y方向梯度算子(用于Harris角点提取算法)
    Iy = filter2(fy, ori_im);                % y方向滤波
    %构造自相关矩阵---------------------------------------------------------------
    Ix2 = Ix .^ 2;
    Iy2 = Iy .^ 2;
    Ixy = Ix .* Iy;
    clear Ix;
    clear Iy;
    h= fspecial('gaussian', [3 3], 2);        % 产生7*7的高斯窗函数，sigma=2
    Ix2 = filter2(h,Ix2);
    Iy2 = filter2(h,Iy2);
    Ixy = filter2(h,Ixy);
    %提取特征点---------------------------------------------------------------
    height = size(ori_im, 1);
    width = size(ori_im, 2);
    result = zeros(height, width);           % 纪录角点位置，角点处值为1
    R = zeros(height, width);
    Rmax = 0;                              % 图像中最大的R值
    k = 0.06; %k为常系数，经验取值范围为0.04~0.06
    for i = 1 : height
        for j = 1 : width
            M = [Ix2(i, j) Ixy(i, j); Ixy(i, j) Iy2(i, j)];             % auto correlation matrix
            R(i,j) = det(M) - k * (trace(M)) ^ 2;                     % 计算R
            if R(i,j) > Rmax
                Rmax = R(i, j);
            end;
        end;
    end;
    %T = 0.01 * Rmax;%固定阈值，当R(i, j) > T时，则被判定为候选角点
    T = 0.1 * Rmax;%固定阈值，当R(i, j) > T时，则被判定为候选角点

    %在计算完各点的值后，进行局部非极大值抑制-------------------------------------
    cnt = 0;
    for i = 2 : height-1
        for j = 2 : width-1
            % 进行非极大抑制，窗口大小3*3
            if (R(i, j) > T && 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))
                result(i, j) = 1;
                cnt = cnt+1;
            end;
        end;
    end;
    i = 1;
        for j = 1 : height
            for k = 1 : width
                if result(j, k) == 1;
                    corners1(i, 1) = j;
                    corners1(i, 2) = k;
                    i = i + 1;
                end;
            end;
        end;
    [posc, posr] = find(result == 1);
    figure,imshow(ori_im);
    hold on;
    plot(posr, posc, 'r+');

end

% -----------------【别人的harris算法实现2】--------------
% 参考地址： http://blog.csdn.net/makenothing/article/details/12884331
% --------------------------【end】----------------------
function harris_refer2
    %MatLab角点检测程序harris。

    ori_im2=rgb2gray(imread('1.jpg'));      
     %ori_im2=imresize(ori_im2',0.50,'bicubic');  %加上这句图就变成竖着的了    

    fx = [5 0 -5;8 0 -8;5 0 -5];          % % la gaucienne，ver axe x  
    Ix = filter2(fx,ori_im2);              % la convolution vers axe x  
    fy = [5 8 5;0 0 0;-5 -8 -5];          % la gaucienne，ver axe y  
    Iy = filter2(fy,ori_im2);              % la convolution vers axe y  
    Ix2 = Ix.^2;  
    Iy2 = Iy.^2;  
    Ixy = Ix.*Iy;  
    clear Ix;  
    clear Iy;  

    h= fspecial('gaussian',[3 3],2);      % générer une fonction gaussienne，sigma=2  

    Ix2 = filter2(h,Ix2);  
    Iy2 = filter2(h,Iy2);  
    Ixy = filter2(h,Ixy);  

    height = size(ori_im2,1);  
    width = size(ori_im2,2);  
    result = zeros(height,width);         % enregistrer la position du coin  

    R = zeros(height,width);  

    K=0.04;  
    Rmax = 0;                              % chercher la valeur maximale de R  
    for i = 1:height  
        for j = 1:width  
            M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];           
            R(i,j) = det(M)-K*(trace(M))^2;                     % % calcule R  
            if R(i,j) > Rmax  
               Rmax = R(i,j);  
            end;  
        end;  
    end;  

    cnt = 0;  
    for i = 2:height-1  
        for j = 2:width-1  
            % réduire des valuers minimales ，la taille de fenetre 3*3  
            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)  
                result(i,j) = 1;  
                cnt = cnt+1;  
            end;  
        end;  
    end;  

    [posr2, posc2] = find(result == 1);  
    cnt                                      % compter des coins  
    figure  
    imshow(ori_im2);  
    hold on;  
    plot(posc2,posr2,'w*');  
end

% -----------------【别人的harris算法实现3】--------------
% 参考地址： http://blog.csdn.net/makenothing/article/details/12884331
% --------------------------【end】----------------------
function harris_refer3
    %%%Prewitt Operator Corner Detection.m  
%%%时间优化--相邻像素用取差的方法求Harris角点  
 %%   
 %clear;  
 Image = imread('1.jpg');                 % 读取图像  
 Image = im2uint8(rgb2gray(Image));     


dx = [-1 0 1;-1 0 1;-1 0 1];  %dx：横向Prewitt差分模版  
Ix2 = filter2(dx,Image).^2;     
Iy2 = filter2(dx',Image).^2;                                          
Ixy = filter2(dx,Image).*filter2(dx',Image);  

%生成 9*9高斯窗口。窗口越大，探测到的角点越少。  
h= fspecial('gaussian',9,2);  
A = filter2(h,Ix2);       % 用高斯窗口差分Ix2得到A   
B = filter2(h,Iy2);                                   
C = filter2(h,Ixy);                                    
nrow = size(Image,1);                              
ncol = size(Image,2);                               
Corner = zeros(nrow,ncol); %zeros用来产生一个全零矩阵，故矩阵Corner用来保存候选角点位置,初值全零，值为1的点是角点  

 %参数t:点(i,j)八邻域的“相似度”参数，只有中心点与邻域其他八个点的像素值之差在  
 %（-t,+t）之间，才确认它们为相似点，相似点不在候选角点之列  
 t=20;  

 %我并没有全部检测图像每个点，而是除去了边界上boundary个像素，也就是从第8行第8列开始遍历。  
 %因为我们感兴趣的角点并不出现在边界上  
 %个人觉得这一部分是的主要目的是找出可能是角点的点，缩小范围，加快运算速度。  
 %具体思想是如果中心点（i,j）周围8个点中有7、8个点灰度值与之相似，那么该中心点应该处于平坦区域，不可能为角点，  
 %如果中心点（i,j）周围只有1个点或者没有点与之相似，那么该中心点也不可能为角点。  
 boundary=8;  
 for i=boundary:nrow-boundary+1   
    for j=boundary:ncol-boundary+1  
         nlike=0; %相似点个数  
         if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)<Image(i,j)+t   
            nlike=nlike+1;  
         end  
         if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)<Image(i,j)+t    
            nlike=nlike+1;  
         end    
        if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)<Image(i,j)+t    
            nlike=nlike+1;  
         end  
         if nlike>=2 && nlike<=6  
             Corner(i,j)=1;%如果周围有2~6个相似点，那(i,j)就是角点  
         end;  
     end;  
 end;  
CRF = zeros(nrow,ncol);    % CRF用来保存角点响应函数值,初值全零  
 CRFmax = 0;                % 图像中角点响应函数的最大值，作阈值之用   
k=0.05;     
% 计算CRF  
 %工程上常用CRF(i,j) =det(M)/trace(M)计算CRF，那么此时应该将下面第105行的  
 %比例系数k设置大一些，k=0.1对采集的这几幅图像来说是一个比较合理的经验值  
 for i = boundary:nrow-boundary+1  
     for j = boundary:ncol-boundary+1  
     if Corner(i,j)==1  %只关注候选点  
         M = [A(i,j) C(i,j);  
              C(i,j) B(i,j)];        
         CRF(i,j) = det(M)-k*(trace(M))^2;  
         if CRF(i,j) > CRFmax   
            CRFmax = CRF(i,j);      
        end;              
    end  
 end;               
end;    
%CRFmax  
 count = 0;       % 用来记录角点的个数  
 t=0.01;          
% 下面通过一个3*3的窗口来判断当前位置是否为角点  
 for i = boundary:nrow-boundary+1   
for j = boundary:ncol-boundary+1  
         if Corner(i,j)==1  %只关注候选点的八邻域  
             if CRF(i,j) > t*CRFmax && CRF(i,j) >CRF(i-1,j-1) ......%?????为什么要CRF(i,j) > t*CRFmax啊？求大神告知  
                && CRF(i,j) > CRF(i-1,j) && CRF(i,j) > CRF(i-1,j+1) ......  
                && CRF(i,j) > CRF(i,j-1) && CRF(i,j) > CRF(i,j+1) ......  
                && CRF(i,j) > CRF(i+1,j-1) && CRF(i,j) > CRF(i+1,j)......  
                && CRF(i,j) > CRF(i+1,j+1)   
            count=count+1;%这个是角点，count加1  
             else % 如果当前位置（i,j）不是角点，则在Corner(i,j)中删除对该候选角点的记录  
                 Corner(i,j) = 0;       
            end;  
         end;   
end;   
end;   
% disp('角点个数');  
 % disp(count)  
 figure,imshow(Image);      % display Intensity Image  
 hold on;   
% toc(t1)  
 for i=boundary:nrow-boundary+1   
for j=boundary:ncol-boundary+1  
         column_ave=0;  
         row_ave=0;  
         k=0;  
        if Corner(i,j)==1  
             for x=i-3:i+3  %7*7邻域  
                 for y=j-3:j+3  
                     if Corner(x,y)==1  
 % 用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取意义不大  
                         row_ave=row_ave+x;  
                         column_ave=column_ave+y;  
                         k=k+1;  
                     end  
                 end  
             end  
         end  
         if k>0 %周围不止一个角点             
             plot( column_ave/k,row_ave/k ,'g.');  
         end  
 end;   
end;   
end

% http://blog.csdn.net/jinshengtao/article/details/17717731

参考文献
1.imfilter 与conv2 和 filter2 之间的区别

2.harris源码1

3.harris源码2