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