图像配准经典算法SIFT

clear all; I1=imread('F:\乱七八糟\car-avi-001.bmp'); %read the first frame I2=imread('F:\乱七八糟\car-avi-002.bmp'); %read the second frame dm=7; %给图像扩边，每个边都扩dm大小。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % |-----|-------------------------------|------| % | B | C | D | % |-----|-------------------------------|------| % | | | | % | | | | % | E | A | F | ----> 图像 % | | | | % | | | | % | | | | % |-----|-------------------------------|------| % | G | H | I | % |-----|-------------------------------|------| %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I1=double(I1); %I1 代表第一幅图像 I2=double(I2); %I2 代表第一幅图像，作为当前帧 [rownum colnum] = size(I1); %给I1扩边，并赋值 % 学生自己添加。扩边后的图像矩阵用II表示。 % ................................ A=I1; B=zeros(dm,dm); for m=1:dm for n=1:dm B(m,n)=128; end end D=B;G=B;I=B; C=zeros(dm,colnum); for n=1:colnum C(:,n)=I1(1,n); end E=zeros(rownum,dm); for m=1:rownum E(m,:)=I1(m,1); end F=zeros(rownum,dm); for m=1:rownum F(m,:)=I1(m,colnum); end H=zeros(dm,colnum); for n=1:colnum H(:,n)=I1(rownum,n); end %imgrownumber = rownum + 2*dm; %imgcolnumber = colnum + 2*dm; II=[B,C,D; %拼接 E,A,F; G,H,I]; %下面进行全搜索算法 blocksize=16; rowblocks =rownum/blocksize; colblocks =colnum/blocksize; A=99999999999999999999; %为了找到最小的均方误差，A用于设定一个很大的初值 Eij=0; xrecord=zeros(32,32); %xrecord，yrecord这两个数组的下标表示存放的匹配块的块号，而其值表示即运动矢量 yrecord=zeros(32,32); %diff=zeros(256,256); %这幅图像的大小为256*256，diff用于存放像素差值 diff=zeros(rownum,colnum); tic count=0; for x=0:(rowblocks-1) %x表示行中第几个子块 row=x*blocksize; for y=0:(colblocks-1) %y表示列中第几个子块 col=y*blocksize; for p=-dm:dm for q=-dm:dm %（p，q）表示（x，y）对应子块在前一帧中所搜索的偏离距离，即运动矢量 Eij=0; %因为图像II是对图像I1扩边的结果，所以在计算Eij时，需要由dm的偏移。如果不扩边，Eij的计算应该为： %Eij=sum(sum((I2(row+1:row+blocksize,col+1:col+blocksize)-II(row+p+1:row+dm+p+blocksize,col+q+1:col+dm+q+blocksize)).^2))/(blocksize^2); Eij=sum(sum((I2(row+1:row+blocksize,col+1:col+blocksize)-II(row+dm+p+1:row+dm+p+blocksize,col+dm+q+1:col+dm+q+blocksize)).^2))/(blocksize^2); count=count+1; if Eij<A A=Eij; DFD(x*16+y+1)=A; xrecord(x+1,y+1)=p; yrecord(x+1,y+1)=q; end end end A=999999999999999999; for mx=1:blocksize for ny=1:blocksize diff(row+mx,col+ny)=I2(row+mx,col+ny)-II(row+mx+dm+xrecord(x+1,y+1),col+ny+dm+yrecord(x+1,y+1)); end end end end toc figure,imshow(I1,[]); title('the first frame'); figure,imshow(I2,[]); title('the second frame'); IIII=I2-I1; figure,imshow(IIII,[]); title('帧间差值'); figure,imshow(diff,[]); title('利用全搜索算法块匹配后的帧间差'); for x=0:(rowblocks-1) %x表示行中第几个子块 row=x*blocksize; for y=0:(colblocks-1) %y表示列中第几个子块 col=y*blocksize; %由于没有量化与反量化，因此，恢复出的III跟I2的内容完全一样 III(row+1:row+blocksize,col+1:col+blocksize)=II(row+dm+xrecord(x+1,y+1)+1:row+dm+xrecord(x+1,y+1)+blocksize,col+dm+yrecord(x+1,y+1)+1:col+dm+yrecord(x+1,y+1)+blocksize)+diff(row+1:row+blocksize,col+1:col+blocksize); end end figure,imshow(III,[]); title('恢复后的第二帧图像'); ERR=diff; figure,imshow(ERR,[]); title('利用全搜索算法块匹配后的帧间差DFD'); % 学生自己编程统计DFD分布，并画出曲线 % ............................rownum colnum figure,plot(DFD); title('全搜索法帧间差DFD分布'); %figure,plot(DFD(1,:));