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. Ix2 = filter2(h,Ix2);   
8. Iy2 = filter2(h,Iy2);   
9. 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. cnt = 0; %角點個數  
12. for i = 2:height-1   
13.     for j = 2:width-1   
14.         % 進行非極大抑制，視窗大小3*3   
15.         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)   
16.             result(i,j) = 1;   
17.             cnt = cnt+1;   
18.         end  
19.     end  
20. end  
21. [posc, posr] = find(result == 1);   
22. disp(cnt);                 % 顯示角點個數   
23. imshow(ori_im);   
24. hold on;   
25. plot(posr,posc,'r+');   

以上資訊來源：

詳細原理及公式推導http://www.cse.psu.edu/~rcollins/CSE486/lecture06.pdf

算法精講及總結http://blog.163.com/[email protected]/blog/static/47586030201132611115984/

程式設計實作http://www.cnblogs.com/blue-lg/archive/2011/12/17/2291139.html