自定义函数都在下面 clc

clear all % 读入图片

pic1=imread('lena1.jpg'); pic2=imread('lena2.jpg'); % Harris角点检测

points1=myHarris(pic1); points2=myHarris(pic2); % 画出Harris角点 figure(1)

drawHarrisCorner(pic1,points1,pic2,points2); % 角点特征描述

des1=myHarrisCornerDescription(pic1,points1); des2=myHarrisCornerDescription(pic2,points2); % 角点粗匹配

matchs=myMatch(des1,des2); % 获取各自出匹配角点位置

matchedPoints1=points1(matchs(:,1),:); matchedPoints2=points2(matchs(:,2),:); % 粗匹配角点连线 figure(2)

drawLinedCorner(pic1,matchedPoints1,pic2,matchedPoints2); % 角点精匹配

[newLoc1,newLoc2]=pointsSelect(matchedPoints1,matchedPoints2); % 精匹配角点连线 figure(3)

drawLinedCorner(pic1,newLoc1,pic2,newLoc2); % 图像拼接

im=picMatched(pic1,newLoc1,pic2,newLoc2); % 显示拼接图像 figure(4) imshow(im);

set(gcf,'Color','w');

function points=myHarris(pic) % 功能：寻找Harris角点 % 输入：RGB图像或gray图

% 输出：角点所在的行、纵的N×2矩阵 if length(size(pic))==3 pic=rgb2gray(pic); end

pic=double(pic); hx=[-1 0 1];

Ix=filter2(hx,pic);

hy=[-1;0;1]; Iy=filter2(hy,pic); Ix2=Ix.*Ix; Iy2=Iy.*Iy; Ixy=Ix.*Iy;

h=fspecial('gaussian',[7 7],2); Ix2=filter2(h,Ix2); Iy2=filter2(h,Iy2); Ixy=filter2(h,Ixy);

[heigth,width]=size(pic); alpha=0.06;

R=zeros(heigth,width); for i=1:heigth for j=1:width

M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; R(i,j)=det(M)-alpha*(trace(M)^2); end end

Rmax=max(max(R));

pMap=zeros(heigth,width); for i=2:heigth-1 for j=2:width-1

if R(i,j)>0.01*Rmax tm=R(i-1:i+1,j-1:j+1); tm(2,2)=0; if R(i,j)>tm pMap(i,j)=1; end end end end

[row,col]=find(pMap==1); points=[row,col];

function drawHarrisCorner(pic1,points1,pic2,points2) % 功能：画出Harris角点的连接 % 输入：

% pic1、pic2：待拼接的图像

% points1、points2：Harris角点位置 X1=points1(:,2); Y1=points1(:,1); X2=points2(:,2); Y2=points2(:,1); dif=size(pic1,2);

imshowpair(pic1,pic2,'montage');

hold on

plot(X1,Y1,'b*'); plot(X2+dif,Y2,'b*'); set(gcf,'Color','w');

function des=myHarrisCornerDescription(pic,points) % 功能：Harris角点特征描述 % 输入：

% pic：原图像

% points：角点位置 % 输出：

% des：8×N的角点特征描述矩阵 if length(size(pic))==3 pic=rgb2gray(pic); end

len=length(points); des=zeros(8,len); for k=1:len

p=points(k,:);

pc=pic(p(1),p(2));

des(1,k)=pic(p(1)-1,p(2)-1)-pc; des(2,k)=pic(p(1),p(2)-1)-pc; des(3,k)=pic(p(1)+1,p(2)-1)-pc; des(4,k)=pic(p(1)+1,p(2))-pc; des(5,k)=pic(p(1)+1,p(2)+1)-pc; des(6,k)=pic(p(1),p(2)+1)-pc; des(7,k)=pic(p(1)-1,p(2)+1)-pc; des(8,k)=pic(p(1)-1,p(2))-pc; des(:,k)=des(:,k)/sum(des(:,k)); end

function matchs=myMatch(des1,des2) % 功能：特征点双向匹配 % 输入：

% des1、des2：特征点描述信息构成的矩阵 % 输出：

% matchs：匹配的特征点对应关系 len1=length(des1); len2=length(des2); match1=zeros(len1,2); cor1=zeros(1,len2); for i=1:len1 d1=des1(:,i); for j=1:len2 d2=des2(:,j);

cor1(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));