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)); k7=(a1(i+1,j-1)-u1(i,j))*(a2(m+1,n-1)-u2(m,n)); k8=(a1(i+1,j)-u1(i,j))*(a2(m+1,n)-u2(m,n)); k9=(a1(i+1,j+1)-u1(i,j))*(a2(m+1,n+1)-u2(m,n)); num=k1+k2+k3+k4+k5+k6+k7+k8+k9;

den=sqrt(A(i,j)*B(m,n)); ncc=num/den; if ncc>max max=ncc;p=m;q=n; end

end

res2(p,q)=1;

end

%做不出来可以先搜索一个点 res=res2; return

file3------------------------------------------------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %特征点的匹配，主要应用 harris 角点的检测,match单向匹配 函数 %适合于有白边的图像，因为，在加窗滤波的时候，没有限定范围的哈，尽量保证角点不在边上 clc,clear all; a1=imread('qiao1.bmp'); a2=imread('qiao2.bmp'); %double的作用很大的啊 subplot(1,2,1); imshow(a1);subplot(1,2,2);imshow(a2); title('原来的图像');

[result1,cnt1,r1,c1]=harris(a1);%角点检测，得到原始的焦点位置图result [result2,cnt2,r2,c2]=harris(a2); figure;subplot(1,2,1);imshow(a1);hold on;plot(c1,r1,'g.'); subplot(1,2,2);imshow(a2);hold on;plot(c2,r2,'g.'); title('图1，2的角点图');

res2=match(a1,cnt1,r1,c1,a2,cnt2,r2,c2);%从result1中开始搜索在result2中可能达到的 [r22,c22]=find(res2==1);