注明:这是一个转载的程序,是我见过的最好的Harris角点的讲解程序了,拿出来分享,原作者找不到了,希望见谅
Harris角点检测算法编程步骤及示例演示
也不说那么多废话了,要介绍啥背景意义之类的,角点检测,顾名思义,就是检测角点,最简单的就是两条线的交点了,还有比如下国际象棋的棋盘格子的交点之类的,反正就是检测这些点。
简单将Harris角点检测算法的思想说下,就是拿一个小窗在图像中移动,通过考察这个小窗口内图像灰度的平均变换值来确定角点。
(1)如果窗口内区域图像的灰度值恒定,那么所有不同方向的偏移几乎不发生变化;
(2)如果窗口跨越一条边,那么沿着这条边的偏移几乎不发生变化, 但是与边垂直的偏移会发生很大的变化;
(3)如果窗口包含一个孤立的点或者角点,那么所有不同方向的偏移会发生很大的变化。
下面给出具体数学推导:
设图像窗口平移量为(u,v),产生的灰度变化为E(u,v),
有E(u,v)=sum[w(x,y)[I(x+u,y+v)-I(x,y)]^2],其中w(x,y)为窗口函数,I(x+u,y+v)为平移后的灰度值,I(x,y)为平移前的灰度值。
有泰勒公式展开可得:
I(x+u,y+v)=I(x,y)+Ix*u+Iy*v+O(u^2,v^2);
Ix,Iy分别为偏微分,在图像中为图像的方向导数.
因此E(u,v)=sum[w(x,y) [Ix*u+Iy*v+O(u^2,v^2)]^2],
可以近似得到E(u,v)=sum[w(x,y) [Ix*u+Iy*v]^2],即
E(u,v)=[u,v][Ix^2,Ix*Iy;Ix*Iy,Iy^2][u,v]T
令M=[Ix^2,Ix*Iy;Ix*Iy,Iy^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.利用水平,竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M中的四个元素的值。
M=[Ix^2,Ix*Iy;Ix*Iy,Iy^2]
2.对M的四个元素进行高斯平滑滤波,为的是消除一些不必要的孤立点和凸起,得到新的矩阵M。
3.接下来利用M计算对应每个像素的角点响应函数R,即:
R=det(M)-k(trace(M))^2;
也可以使用改进的R:
R=[Ix^2*Iy^2-(Ix*Iy)^2]/(Ix^2+Iy^2);里面没有随意给定的参数k,取值应当比第一个令人满意。
4.在矩阵R中,同时满足R(i,j)大于一定阈值threshold和R(i,j)是某领域内的局部极大值,则被认为是角点。
下面给出程序代码:还不错,自己一步步查资料写的,还算详细。
function main
clc;
clear;
close all;
%The coner points extractions using Harrismethods;
frame=imread('*.jpg');
figure(1);
imshow(frame);
%调用harris角点检测子函数;
%输入参数解释:frame为输入图像,7为高斯滤波窗口大小,2为均方差sigma的值,
%0,04为推荐的k值,winsize为极大抑制窗口的大小且给参数时为奇数;
%输出参数解释:posX为检测到角点X坐标,posY为检测到角点Y坐标,
%cnt为检测到角点的个数,Out_Image为输出图像;
[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,sigma,k,winsize)
%ImageData: gracyscale image of input
%GaussWindow: The sizes of Gauss window
%sigma:The variance
%default value
%winsize为极大抑制窗口大小
Out_Image=frame;
ImageData=frame;
ImageData= double(ImageData(:,:,2)); %ImageData数据矩阵式一个三通道的,我们角点标记只需选择一个可以;
%ImageData=im2bw(ImageData,0.5); %或者将3通道的的图像转换成2值化的图像,完成提取;
%算法解释:
%1:利用水平,竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M矩阵中四个元素的值;
%M=[Ix*Ix,Ix*Iy;Ix*Iy,Iy*Iy]
orig_image=ImageData;
fx=[-2,-1,0,1,2];
Ix=filter2(fx,orig_image);
fy=[-2;-1;0;1;2];
Iy=filter2(fy,orig_image);
Ix2=Ix.*Ix;
Iy2=Iy.*Iy;
Ixy=Ix.*Iy;
%2:对M的四个元素进行高斯平滑滤波,得到新的矩阵M;
%滤波平滑,消除突出点,得到新的矩阵M;
h=fspecial('gaussian',[GaussWindow,GaussWindow],sigma); %建立滤波算子
Ix2=filter2(h,Ix2); %filter2是用h滤波器放在Ix2移动进行模板滤波
Iy2=filter2(h,Iy2); %消除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
M=[Ix2(i,j),Ixy(i,j);Ixy(i,j),Iy2(i,j)]; %%自相关矩阵
R(i,j)=det(M)-0.04*(trace(M))^2; %% 计算R值,det()求一个方阵的行列式(Determinant);trace()求方阵的迹,即该方阵对角线上元素之和;
if R(i,j)>Rmax
Rmax=R(i,j);
end
end
end
%winsize为非极大抑制窗口
winr=(winsize-1)/2; %the radius of the neighborhood
istart=winr+1;
jstart=winr+1;
iend=height-winr;
jend=width-winr;
cnt=0;
for i=istart:iend
forj=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
以下这幅图就是角点检测的一个实际例子:
当然这幅图像是用手机随便照的,如果做标定肯定是不能满足的。下面这幅是使用标准的黑白棋盘格时的角点检测值: