完整版Harris角点检测算法编程步骤及示例演示

也不说那么多废话了，要介绍啥背景意义之类的，角点检测，顾名思 义，就是检测角点，最简单的就是两条线的交点了，还有比如下国际 象棋的棋盘格子的交点之类的，反正就是检测这些点。

简单将 Harris 角点检测算法的思想说下， 就是拿一个小窗在图像中移 动，通过考察这个小窗口内图像灰度的平均变换值来确定角点。

1) 如果窗口内区域图像的灰度值恒定，那么所有不同方向的偏移

几乎不发生变化；

2) 如果窗口跨越一条边， 那么沿着这条边的偏移几乎不发生变化， 但是与边垂直的偏移会发生很大的变化；

3) 如果窗口包含一个孤立的点或者角点，那么所有不同方向的偏

移会发生很大的变化。

面给出具体数学推导： 设图像窗口平移量为( u,v) ,产生的灰度变化为 E(u,v),

有 E(u,v)二sum[w(x,y)[l(x+u,y+v)-l(x,y)F2棋中 w(x,y)为窗口函数, l(x+u,y+v为平移后的灰度值，l(x,y)为平移前的灰度值。 有泰勒公式展开可得：

l(x+u,y+v)=l(x,y)+lx*u+ly*v+0(uA2,v八2); lx,ly分别为偏微分，在图像中为图像的方向导数

因此 E(u,v)=sum[w(x,y) [lx*u+ly*v+O(u八2,v八2)]八2], 可以近似得到 E(u,v)=sum[w(x,y) [lx*u+ly*v]^2],即

E(u,v)二[u,v][lx八2,lx*ly;lx*ly,ly八2][u,v]

T

令M=[lx八2,lx*ly;lx*ly,ly八2]，因此最后对角点的检测成了对矩阵 M的

特征值的分析了，令 M 其特征值为 x1,x2; 当x1>>x2或者x2>>x1，则检测到的是边缘部分;

当x1,x2都很小，图像窗口在所有移动的方向上移动灰度级都无明显 变化.

当X1,X2都很大时且相当，检测到的是角点。 编程时用x1,x2不方便，因此定义角点响应函数；

R二det(M)-k(trace(M))八2;

其中det(M)为矩阵M的行列式，trace(M)为矩阵M的迹。

面给出更具数学公式实际编程的步骤：

1. 利用水平，竖直差分算子对图像的每个像素进行滤波以求得 lx,ly进而求得M中的四个元素的值。 M二[lx八2,lx*ly;lx*ly,ly八2]

2. 对 M 的四个元素进行高斯平滑滤波，为的是消除一些不必要 的孤立点和凸起，得到新的矩阵 M。

3. 接下来利用M计算对应每个像素的角点响应函数 R,即:

R二det(M)-k(trace(M))八2;

也可以使用改进的 R：

R二[lx八2*ly八2-(lx*lyF2]/(lx八2+ly八2);里面没有随意给定的参数 k，

取值应当比第一个令人满意。

4•在矩阵R中，同时满足R(i,j大于一定阈值threshold和R(i,j)

是某领域内的局部极大值，则被认为是角点。

面给出程序代码：还不错，自己一步步查资料写的，还算详 细。

function main clc; clear; close all;

%The coner points extractions using Harris methods; frame=imread('*.jpg'); figure(1); imshow(frame);

%调用 harris 角点检测子函数 ;

%输入参数解释： frame 为输入图像 ,7 为高斯滤波窗口大小， 2 为均方差 sigma 的值，

%0,04为推荐的 k 值， winsize 为极大抑制窗口的大小且给

参数时为奇数 ;

%输出参数解释：posX为检测到角点X坐标，posY为检测

到角点Y坐标,

%cnt为检测到角点的个数，Out_lmage为输出图像；

[posX,posY,cnt,Out_Image]=conerdetection(frame,7,2,0.04,7) %输出的图像已经是 2 值化的 figure(2); imshow(Out_Image); hold on;

plot(posX,posY,'ro','MarkerSize',15); disp(cnt); end function

[posX,posY,cnt,Out_Image]=conerdetection(frame,GaussWindow,s igma,k,winsize)

%ImageData: gracyscale image of input %GaussWindow: The sizes of Gauss window %sigma:The variance %default value %wi nsize为极大抑制窗口大小 Out_Image=frame; ImageData=frame;

ImageData= double(ImageData(:,:,2));

%ImageData 数据

矩阵式一个三通道的，我们角点标记只需选择一个可以

%ImageData=im2bw(ImageData,0.5);

%或者将 3 通

道的的图像转换成 2 值化的图像，完成提取 ;

%1：利用水平，竖直差分算子对图像的每个像素进行滤波

以求得lx,ly,进而求得M矩阵中四个元素的值；

%M=[Ix*Ix,Ix*Iy;Ix*Iy,Iy*Iy] orig_image=lmageData； fx=[-2,-1,0,1,2]； lx=filter2(fx,orig_image)； fy=[-2；-1；0；1；2]； ly=filter2(fy,orig_image)； lx2=lx.*lx； ly2=ly.*ly； lxy=lx.*ly；

%2:对 M 的四个元素进行高斯平滑滤波，得到新的矩阵 M； %滤波平滑，消除突出点 ,得到新的矩阵 M；

h=fspecial('gaussian',[GaussWindow,GaussWindow],sigma)； %建立滤波算子

lx2=filter2(h,lx2)； %filter2 是用 h 滤波器放在 lx2 移

动

进行模板滤波

ly2=filter2(h,ly2)； %消除y方向上的突兀点

Ixy=filter2(h,Ixy);

%提取前的图像矩阵的预处理 ;

给宽

height=size(orig_image,1);

%返回图像矩阵的行数给

width=size(orig_image,2); result=zeros(height,width);

%返回图像矩阵的列数 % 纪录角点位置 , 角点处

值为 1

R=zeros(height,width);

%创建与图像矩阵大小

相同的零矩阵

Rmax=0;

% 图像中最大的 R

%3:接下来利用 M 计算对应于每个像素的角点响应函数 Cim(即卩 R);

%计算公式为：R二det(M)-k*(trace(M))八2,其中k为一个任

意

数，经验选取 0.04 就可以;

%由于 k 的取值有些太随意，因此，改用此公式，用新的 公式定义 R： R=det(M)/Tr(M);

%即 Cim=R=[Ix*Ix*Iy*Iy-(Ix*Iy)*(Ix*Iy)]/[Ix*Ix+Iy*Iy]; for i=1:height

for j=1:width

%%自相

%% 计算 R

M=[Ix2(i,j),Ixy(i,j);Ixy(i,j),Iy2(i,j)];

关矩阵

R(i,j)二det(M)-0.04*(trace(M))八2;

值,det()求一个方阵的行列式(Determinant);trace()求方阵的迹, 即该方阵对角线上元素之和 ;

if R(i,j)>Rmax

Rmax=R(i,j); end end end

%wi nsize为非极大抑制窗口 neighborhood winr=(winsize-1)/2; istart=winr+1; jstart=winr+1; iend=height-winr; jend=width-winr;

%the radius of the

cnt=0;

for i=istart:iend

for j=jstart:jend

subr=R((i-winr):(i+winr),(j-winr):(j+winr));

取出 winr*winr 这块区域里面的的矩阵 ;

subrmax=max(max(subr));

if(R(i,j)>k*Rmax)&&(R(i,j)==subrmax)

result(i,j)=1; cnt=cnt+1; end end end

[posY,posX]=find(result==1); % %cnt 为检测出来的角点的个数 ; % figure(2);

% imshow(orig_image); % hold on;

% plot(posX,posY,'ro','MarkerSize',15); % disp(cnt); end

以下这幅图就是角点检测的一个实际例子：

当然这幅图像是用手机随便照的，如果做标定肯定是不能满足

的。下面这幅是使用标准的黑白棋盘格时的角点检测值: