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的迹,两个特征值的积等于矩阵M的行列式。所以用(1-4)式来判定角点质量。(k常取0.04-0.06)

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




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

阅读更多

没有更多推荐了,返回首页