图像拼接matlab.doc

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));
end
[~,indx]=max(cor1);
match1(i,:)=[i,indx];
end
match2=zeros(len2,2);
cor2=zeros(1,len1);
for i=1:len2
d2=des2(:i);
for j=1:len1
d1=des1(:,j);
cor2(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));
end
[~,indx]=max(cor2);
match2(i,:)=[indx,i];
end
matchs=[];
for i=1:length(match1)
for j=1:length(match2)
if match1(i,:)==match2(j,:)
matchs=[matchs;match1(i,:)];
end
end
end
function drawLinedCorner(pic1,loc1,pic2,loc2)
% 功能：画出匹配角点的连接
% 输入：
% pic1、pic2：待拼接的图像
% loc1、loc2：匹配角点位置
X1=loc1(:,2);
Y1=loc1(:,1);
X2=loc2(:,2);
Y2=loc2(:,1);
dif=size(pic1,2);
imshowpair(pic1,pic2,'montage');
hold on
for k=1:length(X1)
plot(X1(k),Y1(k),'b*');
plot(X2(k)+dif,Y2(k),'b*');
line([X1(k),X2(k)+dif],[Y1(k),Y2(k)],'Color','r');
end
set(gcf,'Color','w');
function [newLoc1,newLoc2]=pointsSelect(loc1,loc2)
% 功能：筛选匹配特征点对，获取高精度的控制点
% 输入：
% loc1、loc2：粗匹配特征点位置
% 输出：
% newLoc1、newLoc2：精匹配控制点位置
slope=(loc2(:,1)-loc1(:,1))./(loc2(:,2)-loc1(:,2));
for k=1:3
slope=slope-mean(slope);
len=length(slope);
t=sort(abs(slope));
thresh=t(round(0.5*len));
ind=abs(slope)<=thresh;
slope=slope(ind);
loc1=loc1(ind,:);
loc2=loc2(ind,:);
end
newLoc1=loc1;
newLoc2=loc2;
function im=picMatched(pic1,newLoc1,pic2,newLoc2)
% 功能：获取拼接之后的图片
% 输入：
% pic1、pic2：待拼接图片
% newLoc1、newLoc2：变换控制点矩阵
% 输出：
% im：拼接后的图片
if length(size(pic1))==2
pic1=cat(3,pic1,pic1,pic1);
end
if length(size(pic2))==2
pic2=cat(3,pic2,pic2,pic2);
end
SZ=2000;
X1=newLoc1(:,2);
Y1=newLoc1(:,1);
X2=newLoc2(:,2);
Y2=newLoc2(:,1);
sel=randperm(length(newLoc1),3);
x=X2(sel)';
y=Y2(sel)';
X=X1(sel)';
Y=Y1(sel)';
U=[x;y;ones(1,3)];
V=[X;Y;ones(1,3)];
T=V/U;
cntrX=SZ/2;
cntrY=SZ/2;
im=zeros(SZ,SZ,3);
for i=1:size(pic2,1)
for j=1:size(pic2,2)
tmp=T*[j;i;1];
nx=round(tmp(1))+cntrX;
ny=round(tmp(2))+cntrY;
if nx>=1 && nx<SZ && ny>=1 && ny<=SZ
im(ny,nx,:)=pic2(i,j,:);
end
end
end
im=imresize(im,1,'bicubic');
tpic1=zeros(SZ,SZ,3);
tpic1(1+cntrY:size(pic1,1)+cntrY,1+cntrX:size(pic1,2)+cntrX,:)=pic1;
re=rgb2gray(uint8(im))-rgb2gray(uint8(tpic1));
for k=1:3
ta=im(:,:,k);
tb=tpic1(:,:,k);
ta(re==0)tb(re==0);
im(:,:,k)=ta;
end
clear ta tb re tpic1
im=getPicture(im,SZ);
im=uint8(im);
if length(size(pic1))==2
im=rgb2gray(im);
end
function im=getPicture(pic,SZ)
% 功能：获取图像有用区域
% 输入：
% pic：拼接图像
% SZ：预定图像尺寸
% 输出：
% im：有用区域图像
if length(size(pic))==2
pic=cat(3,pic,picpic);
end
k=1;
while k0
if any(any(pic(k,:,:)))
break
end
k=k-1;
end
bottom=k; %下边界
k=1;
while k0
if any(any(pic(:,k,:)))
break
end
k=k-1;
end
right=k; %右边界
%%获取图像
im=pic(ceil:bottom,left:right,:);