harris角点检测与ncc匹配

harris角点检测与ncc匹配

file1:--------------------------------------------------------------------------------------

function [y1,y2,r,c]=harris(X)

% 角点的检测，利用harris 算法

% 输出的是一幅图像

% [result,cnt,r,c]=harris(X)

% clc,clear all;

% filename='qiao1.bmp';

% X= imread('filename.bmp'); % 读取图像

% Info=imfinfo(filename); %这个要习惯用

% % f=rgb2gray(X);

f=X;

%

ori_im=double(f)/255; %unit8转化为64为双精度double64

fx = [-2 -1 0 1 2]; % x方向梯度算子(用于Harris角点提取算法) Ix = filter2(fx,ori_im); % x方向滤波善于使用filter

% fy = [5 8 5;0 0 0;-5 -8 -5]; % 高斯函数一阶微分，y方向(用于改进的Harris 角点提取算法)

fy = [-2;-1;0;1;2]; % y方向梯度算子(用于Harris角点提取算法) Iy = filter2(fy,ori_im); % y方向滤波

Ix2 = Ix.^2;

Iy2 = Iy.^2;

Ixy = Ix.*Iy;

clear Ix;

clear Iy; %消除变量哈

h= fspecial('gaussian',[10 10 ],2); % 产生7*7的高斯窗函数，sigma=2

Ix2 = filter2(h,Ix2);

Iy2 = filter2(h,Iy2);

Ixy = filter2(h,Ixy); %分别进行高斯滤波

height = size(ori_im,1);

width = size(ori_im,2);

result = zeros(height,width); % 纪录角点位置，角点处值为1 ，背景都是黑色的哈

R = zeros(height,width);

Rmax = 0; % 图像中最大的R值以便设置门限

for i = 1:height

for j = 1:width

M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; %2*2的矩阵

R(i,j) = det(M)-0.06*(trace(M))^2; % 计算R ，求得RMAX，看来是整体求得的，角点响应函数

if R(i,j) > Rmax

Rmax = R(i,j);

end;

end;

end;

cnt = 0; %记录点数的

for i = 2:height-1

for j = 2:width-1 % 进行非极大抑制，窗口3*3

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)

result(i,j) = 1;

cnt = cnt+1;

end;

end;

end;

%

% i=1;

% for j=1:height

% for k=1:width

% if result(j,k)==1;

% corners1(i,1)=j;

% corners1(i,2)=k;

% i=i+1;

% end;

% end;

% end;

[posr, posc] = find(result == 1);

% 角点个数

% imshow(ori_im) %和 X的效果是一样的

% hold on;

% plot(posr,posc,'r.');

y1=result;

y2=cnt;

r=posr;c=posc;

return;

file2--------------------------------------------------------------------------------------------------------------------------------

function res=match(a1,cnt1,r1,c1,a2,cnt2,r2,c2)

% res=match(a1,a2)

% 将从a1寻找a2中的最佳匹配点，得到从a2中抽取的res，也就是单向搜索

% [result1,cnt1,r11,c11]=harris(a1);

% [result2,cnt2,r22,c22]=harris(a2);%%%%%%%%%%%%%%%%%%%%%%%%%%%%可以保证想匹配哪些点就匹配哪些

% figure;

% imshow(result1);title('result1角点位置');

% figure;title('result2角点位置');

% imshow(result2);

win=[1/9 1/9 1/9;1/9 1/9 1/9;1/9 1/9 1/9];

u1=filter2(win,a1);

u2=filter2(win,a2); %求均值

%基于点特征的图像配准算法研究山大

a1=double(a1);

a2=double(a2);

A=filter2(win,(a1-u1).^2);%求方差

B=filter2(win,(a2-u2).^2);

[m1,n1]=size(a1);

[m2,n2]=size(a2);

res1=zeros(m1,n1);

res2=zeros(m2,n2); %寻找的匹配的点

for s=1:cnt1

max=0; p=0;q=0;i=r1(s,1);j=c1(s,1); %p.q存放坐标

for v=1:cnt2

m=r2(v,1);n=c2(v,1);

k1=(a1(i-1,j-1)-u1(i,j))*(a2(m-1,n-1)-u2(m,n)); %%%%%%%%%%%用result是找不到什么信息，3*3对于如此稀疏点没有用

k2=(a1(i-1,j)-u1(i,j))*(a2(m-1,n)-u2(m,n));

k3=(a1(i-1,j+1)-u1(i,j))*(a2(m-1,n+1)-u2(m,n)); %用循环要好一些

k4=(a1(i,j-1)-u1(i,j))*(a2(m,n-1)-u2(m,n));

k5=(a1(i,j)-u1(i,j))*(a2(m,n)-u2(m,n));

k6=(a1(i,j+1)-u1(i,j))*(a2(m,n+1)-u2(m,n));