Harris Corner Detector 原理及编程实现

原理：

灰度变化率有函数如下：

其中的W(x,y)函数如下

有函数f(x+u,y+v)，则其泰勒展开为：

对图片pitch来说，套用刚刚的展开，我们想要计算patch和相邻patch的距离就是：

设

下面，让我们回到最开始的公式，并对中括号里的内容做相应展开：

我们把W(x,y)和矩阵结合起来，用一个新的符号M来表示：

最开始公式就简化成了以下的形式：

以上就是Harris的公式推导。

Harris intuition：

Harris采用了一种新的角点判定方法。矩阵M的两个特征向量l1和l2与矩阵M的主曲率成正比。Harris利用l1, l2来表征变化最快和最慢的两个方向.若两个都很大就是角点,一个大一个小就是边缘,两个都小就是在变化缓慢的图像区域.

　　　　　　　　　　　　　　　　　　　　　　　　　　　　

图1- 4用矩阵Ｍ的特征向量分类图像像素点

但是解特征向量需要比较多的计算量,且两个特征值的和等于矩阵M的迹，两个特征值的积等于矩阵M的行列式。所以用(1-4）式来判定角点质量。(k常取0.04-0.06)

R在两个特征值平面上的梯度表示如下：

Harris算法实现

Step 1:对每一像素点计算相关矩阵M。

Step 2:计算每像素点的Harris 角点响应。　

Step 3.在w*w范围内寻找极大值点，若Harris 角点响应大于阀值，则视为角点。

Harris算子对灰度的平移是不变的,因为只有差分，对旋转也有不变性，但是对尺度很敏感,在一个尺度下是角点, 在在另一个尺度下可能就不是了.

编程实现：

fx = [-1 0 1;-1 0 1;-1 0 1];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法) 
%fx = [-2 -1 0 1 2];                 % x方向梯度算子(用于Harris角点提取算法) 
Ix = filter2(fx,ori_im);              % x方向滤波 
fy = [-1 -1 -1;0 0 0;1 1 1];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法) 
% fy = [-2;-1;0;1;2];                 % y方向梯度算子(用于Harris角点提取算法) 
Iy = filter2(fy,ori_im);              % y方向滤波 

3)考虑到图像一般情况下的噪声影响，采用高斯滤波去除噪声点。

Ix2 = Ix.^2; 
Iy2 = Iy.^2; 
Ixy = Ix.*Iy; 
clear Ix; 
clear Iy; 
h= fspecial('gaussian',[7 7],2);      % 产生7*7的高斯窗函数，sigma=2 
 
Ix2 = filter2(h,Ix2); 
Iy2 = filter2(h,Iy2); 
Ixy = filter2(h,Ixy); 

4)计算角点的准则函数R（即用一个值来判断该点来衡量这个点是否是角点）,并标记角点（R(i,j)>0.01*Rmax,且R(i,j)为3x3邻域局部最大值）。

M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
R(i,j) = det(M)-k*(trace(M))^2; % 计算R

【可以通过改变准则函数的计算来改进harris算法，上式中的k一般取0.04~0.06】

[height,width] = size(ori_im); 
result = zeros(height,width);         % 纪录角点位置，角点处result的值为1 
R = zeros(height,width); 

　　

Rmax = 0;                              % 图像中最大的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)-0.06*(trace(M))^2;          
        if R(i,j) > Rmax 
            Rmax = R(i,j); 
        end 
    end
end
 
cnt = 0; %角点个数
for i = 2:height-1 
    for j = 2:width-1 
        % 进行非极大抑制，窗口大小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
 
[posc, posr] = find(result == 1); 
disp(cnt);                 % 显示角点个数 
imshow(ori_im); 
hold on; 
plot(posr,posc,'r+'); 

以上信息来源：

详细原理及公式推导http://www.cse.psu.edu/~rcollins/CSE486/lecture06.pdf

算法精讲及总结http://blog.163.com/zhaowei0425@126/blog/static/47586030201132611115984/

编程实现http://www.cnblogs.com/blue-lg/archive/2011/12/17/2291139.html