这是一个高中同学Q我让我帮他解决的一个问题，大意就是对飞机拍摄得到的矫正后图像进行拼接。opencv里有相应的实现(cv2.Stitcher_create)，但是我试了一下似乎只能左右拼接，而且我没用过opencv，所以想用matlab做，然后python调用就可以了。
以下为思路：
输入N个图像（为加快速度进行下采样），对每个图像进行特征点检测（getFeatures），获得loc和des，loc为图像每个特征点的坐标。并把每个图像的loc和des放入元组(cell)中。
同样设立一个n*n的元组，用来储存一个图像和另一个图像的匹配映射matched。
循环两次计算元组后(对称矩阵，存储一半就可以)，每次找到元组存储的映射点数最多的个数(matched(i)=j ；若j为0 则表明图像A的第i个点在B无对应点，否则为B的第j个点 A的第i个特征点的坐标为loc(i,1),loc(i,2))，假设match{i,j}中映射点最多，那么就先把第i个图像和第j个图像进行拼接(stitch.m),得到新图像k。i和j再之后不再进行访问，所以 在存储各个图像的元组I、元组loc和元组des中，分别把第n-1个图像的I、loc、des覆盖到i,把第n个覆盖到j，并把第n-1个赋值为新图像的数据。然后n递减。由于我已经晕了，所以二维元组match直接全部重新计算了。
%分别求1-n的loc 和des
%分别求每两个图的match
%对match中的匹配点个数进行排序
%拼接匹配点最多的两个图像 得到新图像 代替原来的图像
%重复 直到只剩一张
%每次迭代 n-1
%每次提出两个放在最后 把合成的放在第n-1
clear
clc
close all
n=4;%一共多少个图片
I1=imread('1.jpg');
I2=imread('2.jpg');
I3=imread('3.jpg');
I4=imread('4.jpg');
I=cell(1,n);%用元祖存储一系列图像
p=0.5;%压缩比例
I{1}=double(rgb2gray(imresize(I1,p)));
I{2}=double(rgb2gray(imresize(I2,p)));
I{3}=double(rgb2gray(imresize(I3,p)));
I{4}=double(rgb2gray(imresize(I4,p)));
des=cell(1,n);
loc=cell(1,n);
for i=1:n
[des{i},loc{i}] = getFeatures(I{i});
end
matched=cell(n,n);%n*n元组
match_num=zeros(n,n);
while n>1
for i =1:n
for j=1:n
if i<j
matched{i,j}=match(des{i},des{j});
match_num(i,j)=sum(sum(matched{i,j}~=0));%统计每个匹配的个数
end
end
end
max_num=0;
max_i=0;
max_j=0;
for i=1:n
for j=1:n
if match_num(i,j)>max_num
max_i=i;
max_j=j;
max_num=match_num(i,j);
end
end
end
%获得匹配点最多的两个图片的下标 max_i和max_j
I_new=stitch(I{max_j},I{max_i},loc{max_j},loc{max_i},matched{max_i,max_j});
figure,imshow(uint8(I_new))
I{max_i}=I{n-1};%n-1和i交换
loc{max_i}=loc{n-1};
des{max_i}=des{n-1};
I{max_j}=I{n};
loc{max_j}=loc{n};
des{max_j}=des{n};
I{n-1}=I_new;
[des{n-1},loc{n-1}]=getFeatures(I{n-1});
n=n-1;
end
(代码中对j,i进行stitch而不是i,j， 是因为拼接的新图像的下标数值大且面积大，最宜用小图像拼接到稍微大的图像)
两个图像拼接的函数stitch.m:
首先获得两个图像的长宽：h1、h2、d1、d2
考虑最大可能，把图像I1分别向四周拓宽max(d1,d2)个像素，并存为result矩阵中。
调用function map=getmap(matched) 对输入的matched进行调整，获得一个x*2维度的map矩阵，x为matched中对应点的个数，第一列为图像I1的点序号，第二列为图像I2的点序号。
之后分别求map中第一列的点和 对应的第二列的点的x坐标的偏移pos_x和y坐标偏移pos_y的均值，也就是I2相对于I1的偏移。
result按如下方式计算：
其中,D为权值，应在不同区域获得不同的取值，本例使用越靠近I1的D取值越大，越靠近I2的D取值越小的相邻元素等比矩阵mask实现。(代码中mask(i,j)越靠近I2权值越大。mask的取值取决于重叠部分的大小以及位移的方向)
其中，mask也即重叠部分的长宽可按如下方法获得：
左上角坐标可以通过遍历I2，在result中找到第一个不为0的点，就是重叠部分的左上角下标I2中下标与result的对应关系为，result(dev_x+i,dev_y+j)=I2(i,j) 其中dev_x/y为相对于result的偏移(之前的pos_x/y为I2相对于I1)：
dev_x=max(h2,d2)+pos_x; %相对于扩展后的I1(result)的
dev_y=max(h2,d2)+pos_y;
右下角直接计算：
mask_i_max=h1-pos_x+1; mask_j_max=d1-pos_y+1;
(此时result只有I1图像)
遍历图二，判定是否该点在重叠范围内，并且该重叠范围内的点的灰度值不为0。若不为0，则由之前的mask计算result，否则，result（i，j）=I2(i,j)。
最后对result裁剪掉全0的行或者列。
(代码中有注释掉的旋转缩放部分，但是好像出错了)
% I1=imread('1.jpg');
% I2=imread('2.jpg');
% I1=double(rgb2gray(I1));
% I2=double(rgb2gray(I2));
% [des1,loc1] = getFeatures(I1);
% [des2,loc2] = getFeatures(I2);
% matched = match(des1,des2);
%drawMatched(matched,I1,I2,loc1,loc2);
function result=stitch(I1,I2,loc1,loc2,matched)
map=getmap(matched);
[h1,d1]=size(I1);
[h2,d2]=size(I2);
result=zeros(h1+2*max(h2,d2),d1+2*max(h2,d2));
result(max(h2,d2)+1:max(h2,d2)+h1,max(h2,d2)+1:max(h2,d2)+d1)=I1;
pos_x=0;
pos_y=0;%平移
for i=1:length(map)
pos_x=pos_x+loc1(map(i,1),1)-loc2(map(i,2),1);
pos_y=pos_y+loc1(map(i,1),2)-loc2(map(i,2),2);
end
pos_x=round(pos_x/length(map));
pos_y=round(pos_y/length(map));
dev_x=max(h2,d2)+pos_x; %相对于扩展后的I1(result)的
dev_y=max(h2,d2)+pos_y;
%result(dev_x+1:dev_x+h2,dev_y+1:dev_y+d2)=I2;
for i =1:h2
for j=1:d2
bre=0;
if result(dev_x+i,dev_y+j)~=0 %在遮罩范围内
mask_i_min=i;%遮罩在I2内的左上角下标
mask_j_min=j;
bre=1;
break
end
end
if bre==1
break
end
end
mask_i_max=h1-pos_x+1;
mask_j_max=d1-pos_y+1;
mask=getmask(pos_x,pos_y,mask_i_max-mask_i_min,mask_j_max-mask_j_min);
for i =1:h2
for j=1:d2
if i>=mask_i_min && i<=mask_i_max && j>=mask_j_min && j<=mask_j_max && result(dev_x+i,dev_y+j)%在遮罩范围内
%result(dev_x+i,dev_y+j)=255;
result(dev_x+i,dev_y+j)=mask(i-mask_i_min+1,j-mask_j_min+1)*I2(i,j)+(1-mask(i-mask_i_min+1,j-mask_j_min+1))*result(dev_x+i,dev_y+j);
else
result(dev_x+i,dev_y+j)=I2(i,j);
end
end
end
%调整
y_min=1;
[x_max,y_max]=size(result);
for x_min=1:x_max
if sum(result(x_min,:)~=0)
break
end
end
for y_min=1:y_max
if sum(result(:,y_min)~=0)
break
end
end
for x_max=x_max:-1:x_min
if sum(result(x_max,:)~=0)
break
end
end
for y_max=y_max:-1:y_min
if sum(result(:,y_max)~=0)
break
end
end
result=result(x_min:x_max,y_min:y_max);
end
%旋转
% theta=0;
% p=0;
% for i=1:length(map)
% for j=1:length(map)
% if i~=j
% theta=theta+(loc1(map(i,1),1)-loc1(map(j,1),1))/(0.01+loc1(map(i,1),2)-loc1(map(j,1),2)) - ...
% (loc2(map(i,2),1)-loc2(map(j,2),1))/(0.01+loc2(map(i,2),2)-loc2(map(j,2),2));
% len1=sqrt((loc1(map(i,1),1)-loc1(map(j,1),1))^2+(loc1(map(i,1),2)-loc1(map(j,1),2))^2);
% len2=sqrt((loc2(map(i,2),1)-loc2(map(j,2),1))^2+(loc2(map(i,2),2)-loc2(map(j,2),2))^2);
% p=p+(len1/len2);
% end
% end
% end
% theta=theta/(length(map)*length(map)-length(map));
% p=p/(length(map)*length(map)-length(map));
% dir=-theta;
% rotate=atan(-theta)*180/pi;
% B=zeros(2*h2,2*d2);
% B(h2+1:2*h2,d2+1:2*d2)=I2;
% B=imrotate(B,rotate);
% B(1:h1,1:d1)=I1;
% imshow(uint8(B));
% figure,
% S=imresize(B,p,'nearest');
% imshow(uint8(S));
% loc1=uint8(loc1);
% loc2=uint8(loc2);
% R=[cos(dir),-sin(dir);sin(dir),cos(dir)];
% for i=1:length(map)
% I1(loc1(map(i,1),1),loc1(map(i,1),2))=0;
% I2(loc2(map(i,2),1),loc2(map(i,2),2))=0;
% end
% subplot(121),imshow((I1))
% subplot(122),imshow((I2))
function map=getmap(matched)%返回 二维矩阵 表示点的对应关系
map=zeros(sum(matched~=0),2);
k=1;
for i =1:length(matched)
if matched(i)~=0
map(k,2)=i;
map(k,1)=matched(i);
k=k+1;
end
end
end
function mask=getmask(pos_x,pos_y,len_x,len_y) %根据偏移值返回一个m*n的mask pox>0 poy>0时候 右下为1 左上为0
if pos_x>1
d_x=0:1/len_x:1;
elseif pos_x<-1
d_x=1:1/(-len_x):0;
else
d_x=ones(1,1);
end
if pos_y>1
d_y=0:1/len_y:1;
elseif pos_y<-1
d_y=1:1/(-len_y):0;
else
d_y=ones(1,1);
end
mask=d_x'*d_y;
end
运行：
在百度地图上随便截了几张图：
运行结果如下：
（压缩比例p=0.5）
注：为提高速度和效率 可以在对图像进行下采样代替图像获得特征点，之后得到的平移大小和缩放比例运算后可以得到原图像的平移大小，在对原图像进行平移融合即可。
源代码与sift代码(一年前计算机视觉课上下的，无法获得来源)已上传至百度网盘：
链接：https://pan.baidu.com/s/1ax12UNZhsrSe2-wPyrnwBQ
提取码：14q0
zip中.m文件:stitch和multstitch为本例代码 run和process为其他项目代码（无用），其他都为sift代码。
2283
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
