harris corner原理

原理：

灰度变化率有函数如下：

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

对公式中中括号中的部分Taylor展开并取一阶式，得到

设

则矩阵形式

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

　　　　　　　　　　　　　　　　　　　　　　　　　　　　

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

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

Harris算法总结

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

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

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

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

编程实现：

[objc]  view plain copy

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

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

[objc]  view plain copy

1. Ix2 = Ix.^2;   
2. Iy2 = Iy.^2;   
3. Ixy = Ix.*Iy;   
4. clear Ix;   
5. clear Iy;   
6. h= fspecial('gaussian',[7 7],2);      % 产生7*7的高斯窗函数，sigma=2   
7.    
8. Ix2 = filter2(h,Ix2);   
9. Iy2 = filter2(h,Iy2);   
10. 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】

[objc]  view plain copy

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

　　

[objc]  view plain copy

1. Rmax = 0;                              % 图像中最大的R值   
2. for i = 1:height   
3.     for j = 1:width   
4.         M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];         
5.         R(i,j) = det(M)-0.06*(trace(M))^2;            
6.         if R(i,j) > Rmax   
7.             Rmax = R(i,j);   
8.         end   
9.     end  
10. end  
11.    
12. cnt = 0; %角点个数  
13. for i = 2:height-1   
14.     for j = 2:width-1   
15.         % 进行非极大抑制，窗口大小3*3   
16.         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)   
17.             result(i,j) = 1;   
18.             cnt = cnt+1;   
19.         end  
20.     end  
21. end  
22.    
23. [posc, posr] = find(result == 1);   
24. disp(cnt);                 % 显示角点个数   
25. imshow(ori_im);   
26. hold on;   
27. 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