基于harris角点检测的特征点匹配程序（调试总是有问题）

[uneamie]

function [match_points]=points_matching(I1,I2,c1,r1,c2,r2)
%I1,I2为读入的两幅图像；
%[r1,c1]第一幅图角点的坐标；[r2,c2]第二幅图的
%%%&cross-correlation: CC=(sum(sum(I1(x,y)*I2(x,y))))/(sqrt(sum(sumI1(x,y)))*sqrt(sum(sumI2(x,y))))
I1_s=I1.*I1;
I2_s=I2.*I2;

m1=size(c1,1);
m2=size(c2,1);
[C1,R1]=size(I1);
[C2,R2]=size(I2);
match=zeros(m1,2);
radius=20;
for i=1:m1
if c1(i)>radius&&c1(i)<C1-radius&&r1(i)>radius&&r1(i)<R1-radius
temp1=sum(sum(I1_s(r1(i):r1(i)+radius,c1(i):c1(i)+radius)));
temp1=sqrt(temp1);%sqrt(sum(sumI1(x,y)))
for j=1:m2
max=0;
if c2(j)>radius&&c2(j)<C2-radius&&r2(j)>radius&&r2(j)<R2-radius
temp2=sum(sum(I2_s(r2(j):r2(j)+radius,c2(j):c2(j)+radius)));
（每一次总是在这一步出错，说我超过矩阵大小）
temp2=sqrt(temp2);%sqrt(sum(sumI2(x,y)))
temp3=I1(r1(i):r1(i)+radius,c1(i):c1(i)+radius).*I2(r2(j):r2(j)+radius,c2(j):c2(j)+radius);
temp4=sum(sum(temp3));
CC=temp4/(temp1*temp2);
if CC>max
max=CC;
match(i,1)=i;
match(j,2)=j;
（match这样弄不知道对不对，就是想表示出一个角点它所对应的角点）
end
end
end

end
a=find(match(:,1)>0);
match_points=zeros(size(a),4);
for k=1:size(a)
match_points(k,:)=[r1(a(k)),c1(a(k)),r2(a(k)),c2(a(k))] %in every row, we show x-y of two matched points
end
刚开始写程序，自己调了一个晚上，怎么写都弄不出来……只好来请教了

[uneamie]

自己先顶下~

[x642458]

什么语言的，MATLAB？

[uneamie]

[Quote=引用 2 楼 x642458 的回复:]
什么语言的，MATLAB？

[/Quote]恩，是啊

[x642458]

我自己刚写的，LZ参考下吧
%清屏，所以变量置零
clc,clear all;
%读图
image=imread('6.bmp');
Info=imfinfo('6.bmp'); %读图像信息，并判断是否是灰度图
if Info.BitDepth>8
image=rgb2gray(image);
end
image=double(image)/255;
%邻域平均法滤波
image=filter2(fspecial('average',3),image,'full');
X=[-1 0 1]; %X方向梯度算子
fx = filter2(X,image); % x方向滤波
Y=[-1;0;1]; %Y方向梯度算子
fy = filter2(Y,image); % y方向滤波
fx2=fx.^2;
fy2=fy.^2;
fxy=fx.*fy;
clear fx;
clear fy;

% 产生11*11的高斯窗口函数，sigma=2
g=fspecial('gaussian',[11 11],2);
A = filter2(g,fx2);
B = filter2(g,fy2);
C = filter2(g,fxy);

[height,width]=size(image);
%响应函数
R=zeros(height,width);
% 纪录角点位置，角点处值为1
Result=zeros(height,width);
% 图像中最大的R值，用Rmax*0.05作为阈值T
%系数0.05是实验得出
Rmax=0.0;
for i=1:height
for j=1:width
%对图像上每个点（x,y)计算自相关矩阵M
M=[A(i,j) C(i,j)
C(i,j) B(i,j)];
%对图像上每个点（x,y)计算自相关矩阵M
R(i,j)=det(M)-0.06*(trace(M))^2;
if(R(i,j)>Rmax)
Rmax=R(i,j);
end
end
end
%记录检测的角点数
count=0;
for i1=2:height-1;
for j1=2:width-1
% 在3*3窗口范围内寻找局部最大值点作为检测出的角点
a=[R(i1-1,j1-1) R(i1-1,j1) R(i1-1,j1+1) R(i1,j1-1) R(i1,j1+1) R(i1+1,j1-1) R(i1+1,j1) R(i1+1,j1+1)];
if R(i1,j1)>0.01*Rmax && R(i1,j1)>max(a)
Result(i1,j1)=1;
count=count+1;
end
end
end
% 纪录角点位置
[posr,posc]=find(Result==1);
% 输出角点个数
fprintf('检测到的角点数count=%d\n',count);
imshow(image)
hold on;
%角点处画+
plot(posc,posr,'r+');

[fujilove]

楼上的谢谢了，代码很好，刚运行过，检测的角点的质量很好

[qiaoc123]

楼主的问题解决了吗？