图像拼接matlab

clcclearall%读入图片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');functionpoints=myHarris(pic)%功能：寻找Harris角点%输入：RGB图像或gray图%输出：角点所在的行、纵的N×2矩阵iflength(size(pic))==3pic=rgb2gray(pic);endpic=double(pic);hx=[-101];Ix=filter2(hx,pic);hy=[-1;0;1];Iy=filter2(hy,pic);Ix2=Ix.*Ix;Iy2=Iy.*Iy;Ixy=Ix.*Iy;h=fspecial('gaussian',[77],2);Ix2=filter2(h,Ix2);Iy2=filter2(h,Iy2);Ixy=filter2(h,Ixy);[heigth,width]=size(pic);alpha=0.06;R=zeros(heigth,width);fori=1:heigthforj=1:widthM=[Ix2(i,j)Ixy(i,j);Ixy(i,j)Iy2(i,j)];R(i,j)=det(M)-alpha*(trace(M)^2);endendRmax=max(max(R));pMap=zeros(heigth,width);fori=2:heigth-1forj=2:width-1ifR(i,j)0.01*Rmaxtm=R(i-1:i+1,j-1:j+1);tm(2,2)=0;ifR(i,j)tmpMap(i,j)=1;endendendend[row,col]=find(pMap==1);points=[row,col];functiondrawHarrisCorner(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');holdonplot(X1,Y1,'b*');plot(X2+dif,Y2,'b*');set(gcf,'Color','w');functiondes=myHarrisCornerDescription(pic,points)%功能：Harris角点特征描述%输入：%pic：原图像%points：角点位置%输出：%des：8×N的角点特征描述矩阵iflength(size(pic))==3pic=rgb2gray(pic);endlen=length(points);des=zeros(8,len);fork=1:lenp=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));endfunctionmatchs=myMatch(des1,des2)%功能：特征点双向匹配%输入：%des1、des2：特征点描述信息构成的矩阵%输出：%matchs：匹配的特征点对应关系len1=length(des1);len2=length(des2);match1=zeros(len1,2);cor1=zeros(1,len2);fori=1:len1d1=des1(:,i);forj=1:len2d2=des2(:,j);cor1(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));end[~,indx]=max(cor1);match1(i,:)=[i,indx];endmatch2=zeros(len2,2);cor2=zeros(1,len1);fori=1:len2d2=des2(:,i);forj=1:len1d1=des1(:,j);cor2(j)=(d1'*d2)/sqrt((d1'*d1)*(d2'*d2));end[~,indx]=max(cor2);match2(i,:)=[indx,i];endmatchs=[];fori=1:length(match1)forj=1:length(match2)ifmatch1(i,:)==match2(j,:)matchs=[matchs;match1(i,:)];endendendfunctiondrawLinedCorner(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');holdonfork=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');endset(gcf,'Color','w');function[newLoc1,newLoc2]=pointsSelect(loc1,loc2)%功能：筛选匹配特征点对，获取高精度的控制点%输入：%loc1、loc2：粗匹配特征点位置%输出：%newLoc1、newLoc2：精匹配控制点位置slope=(loc2(:,1)-loc1(:,1))./(loc2(:,2)-loc1(:,2));fork=1:3slope=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,:);endnewLoc1=loc1;newLoc2=loc2;functionim=picMatched(pic1,newLoc1,pic2,newLoc2)%功能：获取拼接之后的图片%输入：%pic1、pic2：待拼接图片%newLoc1、newLoc2：变换控制点矩阵%输出：%im：拼接后的图片iflength(size(pic1))==2pic1=cat(3,pic1,pic1,pic1);endiflength(size(pic2))==2pic2=cat(3,pic2,pic2,pic2);endSZ=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);fori=1:size(pic2,1)forj=1:size(pic2,2)tmp=T*[j;i;1];nx=round(tmp(1))+cntrX;ny=round(tmp(2))+cntrY;ifnx=1&&nx=SZ&&ny=1&&ny=SZim(ny,nx,:)=pic2(i,j,:);endendendim=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));fork=1:3ta=im(:,:,k);tb=tpic1(:,:,k);ta(re==0)=tb(re==0);im(:,:,k)=ta;endcleartatbretpic1im=getPicture(im,SZ);im=uint8(im);iflength(size(pic1))==2im=rgb2gray(im);endfunctionim=getPicture(pic,SZ)%功能：获取图像有用区域%输入：%pic：拼接图像%SZ：预定图像尺寸%输出：%im：有用区域图像iflength(size(pic))==2pic=cat(3,pic,pic,pic);endk=1;whilekSZifany(any(pic(k,:,:)))breakendk=k+1;endceil=k;%上边界k=SZ;whilek0ifany(any(pic(k,:,:)))breakendk=k-1;endbottom=k;%下边界k=1;whilekSZifany(any(pic(:,k,:)))breakendk=k+1;endleft=k;%左边界k=SZ;whilek0ifany(any(pic(:,k,:)))breakendk=k-1;endright=k;%右边界%%获取图像im=pic(ceil:bottom,left:right,:);