Matlab Harris角点匹配

图像课程作业，仅供参考
题目：用Harris算子和归一化互相关实现两幅图的匹配


1、自定义Harris函数，寻找图像的Harris角点坐标
(Harris.m)

% coords:图像角点坐标集合（为了方便后期使用归一化互相关做匹配）
function coords=Harris(I)
    [m,n]=size(I);
    %imshow(I);
    %===================================
    % Step 1:计算相关矩阵M
    %===================================
    tmp=zeros(m+2,n+2);
    tmp(2:m+1,2:n+1)=I; % 补零
    Ix=zeros(m+2,n+2);% x方向梯度
    Iy=zeros(m+2,n+2);% y方向梯度
    Ix(:,2:n)=tmp(:,3:n+1)-tmp(:,1:n-1);
    Iy(2:m,:)=tmp(3:m+1,:)-tmp(1:m-1,:);
    Ix2=Ix(2:m+1,2:n+1).^2;
    Iy2=Iy(2:m+1,2:n+1).^2;
    Ixy=Ix(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);
    h=fspecial('gaussian',[7 7],2);
    Ix2=filter2(h,Ix2);
    Iy2=filter2(h,Iy2);
    Ixy=filter2(h,Ixy);
    %===================================
    % Step 2:计算Harris角点响应
    %===================================
    k=0.05; % k常取0.04-0.06
    Rmax=0;
    R=zeros(m,n);
    for i=1:m
        for j=1:n
            M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
            R(i,j)=det(M)-k*(trace(M))^2;     
            if R(i,j)>Rmax
                Rmax=R(i,j);
            end
        end
    end
    %===================================
    % Step 3:寻找Harris角点
    %===================================
    tmp(2:m+1,2:n+1)=R;
    %res=zeros(m+2,n+2);
    %res(2:m+1,2:n+1)=I;
    coords=[];
    for i=2:m+1
        for j=2:n+1
            % 在3*3邻域内非极大值抑制
            if tmp(i,j)>0.01*Rmax &&...
               tmp(i,j)>tmp(i-1,j-1) && tmp(i,j)>tmp(i-1,j) && tmp(i,j)>tmp(i-1,j+1) &&...
               tmp(i,j)>tmp(i,j-1) && tmp(i,j)>tmp(i,j+1) &&...
               tmp(i,j)>tmp(i+1,j-1) && tmp(i,j)>tmp(i+1,j) && tmp(i,j)>tmp(i+1,j+1)
                    %res(i,j)=255; 
                    coords=[coords,[i-1;j-1]];
            end    
        end
    end
    %coords
    %figure,imshow(mat2gray(res(2:m+1,2:n+1)));
end

2、计算互相关矩阵寻找匹配Harris角点对
(程序入口)

I1=imread('C:\Users\53121\Desktop\boat1.bmp');  
I2=imread('C:\Users\53121\Desktop\boat2.bmp');  
coords1=Harris(I1);  
coords2=Harris(I2);  
%===================================  
% 计算互相关矩阵R  
%===================================  
[~,m]=size(coords1);%获取图像1的角度个数  
[~,n]=size(coords2);%获取图像2的角度个数  
R=zeros(m,n);%互相关矩阵  
[r1,c1]=size(I1);  
[r2,c2]=size(I2);  
r=2;%角点匹配领域半径  
I1_tmp=zeros(r1+2*r,c1+2*r);  
I2_tmp=zeros(r2+2*r,c2+2*r);  
I1_tmp(r+1:r1+r,r+1:c1+r)=I1;  
I2_tmp(r+1:r2+r,r+1:c2+r)=I2;  
for i=1:m  
    x1=coords1(1,i)+r;   
    y1=coords1(2,i)+r;  
    for j=1:n  
        x2=coords2(1,j)+r;  
        y2=coords2(2,j)+r;  
        I1_mean=mean(mean(I1_tmp(x1-r:x1+r,y1-r:y1+r)));  
        I2_mean=mean(mean(I2_tmp(x2-r:x2+r,y2-r:y2+r)));  
        a=I1_tmp(x1-r:x1+r,y1-r:y1+r)-I1_mean;  
        b=I2_tmp(x2-r:x2+r,y2-r:y2+r)-I2_mean;  
        p=sum(sum(a.*b));  
        q=sqrt(sum(sum(a.^2))*sum(sum(b.^2)));  
        R(i,j)=p/q;  
    end  
end  
%===================================  
% 寻找相关性最大的角点对并绘图  
%===================================  
% 构造画布  
canvas=ones(max(r1,r2),c1+c2+100)*255;  
canvas(1:r1,1:c1)=I1;  
canvas(1:r2,c1+101:c1+c2+100)=I2;  
imshow(uint8(canvas));  
hold on  
% Harris角点图  
% for i=1:m  
%     x=coords1(2,i);  
%     y=coords1(1,i);  
%     plot(x,y,'ro');  
% end  
% for i=1:n  
%     x=coords2(2,i);  
%     y=coords2(1,i);  
%     plot(x+c1+100,y,'ro');  
% end  
  
[m,n]=find(R==max(max(R)));% 选取相关系数最大的匹配点  
[number,~]=size(m)  
for i=1:number  
    x1=coords1(2,m(i));  
    x2=coords2(2,n(i));  
    y1=coords1(1,m(i));  
    y2=coords2(1,n(i));  
    line([x1,x2+c1+100],[y1,y2]);%匹配连线  
    plot([x1 x2+c1+100],[y1 y2],'ro');  
end  

效果图