图像拼接算法的基本原理

转自：http://blog.csdn.net/yuyin86/article/details/6690423

全景视频是一种利用360 度全景图象建立虚拟环境的新方法。全景图象是通过将普通照相机拍照到的边界部分重叠的图象进行拼接而创建的。可以利用图象重叠部分对应像素的相似性, 通过采用一种行之有效的拼接算法, 使得到的图象无缝平滑。

来自研学论坛 Walkfarer和SCQ的帖子：

http://bbs.matwav.com/post/view?bid=6&id=371051&sty=3&age=0&tpg=1&ppg=1#371051

图像拼接是计算机视觉中的重要分支，它是将两幅以上的具有部分重叠的图像进行无缝拼接从而得到较高分辨率或宽视角的图像。目前成形算法原理大致如下：
1、频率域：（不甚清楚）
利用"相位相关法"完成两幅图像的平移估计（只能精确到像素级）。
2、空间域：
（1）基于特征的方法:找出特征点，进行匹配。
a.使用HARRIS角点,然后使用灰度相关加松弛匹配找到对应点,如果两幅图像重叠区域较大,且透视变形较小,可以考虑使用这种方法。

b.通过特征点的梯度方向等信息，确定一组最佳的特征匹配，利用这一组数据给出两幅图像间变换矩阵的估计初值，再利用递归算法找到最终的精确变换关系。在拍摄图片的相机的旋转与缩放不是很大的情况下可以实现较好效果。
（2）和基于光差的方法（方法精确但收敛慢）。应用最广泛的一类拼接算法是柱面与球面图像的拼接，经过球面与柱面变换后，问题就归结为确定每幅图像的平移量。

以上是我阅读了一些文章总结出来的，贻笑大方了。希望各位修正或添加之。

csq***
我来补充吧.
1 频率域： 一般是用fourier的相位相关，可以估计出频移，旋转，缩放。
频移没什么好讨论的，比较简单， 估计缩放和旋转一般是变到极坐标系去做，无数的paper讨论这个问题，前面讨论“怎样判断两幅图像有没有重叠”的帖子，我给了一篇paper，今年ieee tran ip的。 还不错。频率域的有人做sub-pixel. 这个我也说过， university of centrl florida 的那个什么faroon (名字记不清了)， 写了一篇这样的，但好像只能处理频移的sub-pixel. 我还发email问过他有没有处理旋转的sub-pixel. 他没有回。 ft.
2.空域：
1 。楼主遗漏了基于intensity的方法，不过确实用的不多了。
2。 特征的方法：
看来最多的是 point-rgistration. 不过还是有做用边缘啊那些其他特征去拼的人 。
下面我就介绍一下point发.
Point-registration:
1．经典的harris point, 现在有不少修正版 ，因为harris 当时对参数的选择并没有给出很好的建议，所以参数选择比较烦，我以前就试过很多参数.
参考文献：A Mathematical Comparison of Point Detectors
Evaluation of Interest Point Detectors
找到特征点那就要去匹配拉。
一般是先初步估计一下，剔出差太多的匹配对。 可以用intensity的各种方法，但是这就摆出一个问题： 光照。 很麻烦，两幅图的光照差很多的话，有可能根本弄不出来匹配的点，但你从图上明显可以看到很多对都是匹配的。这个等会讨论。
2．然后进一步估计匹配，RANSAC用的最多，也有其他的方法， 像paper: MLESAC- A new robust estimator with application to estimating image geometry.
估计出匹配对，然后就要算那个乱七八糟的矩阵。 又是很多方法可以来做。一般是各种优化算法像 LM之类的。
3．算出矩阵，然后把一个变到另外一个的坐标系，就是融合的问题了。咋个无缝拼接是个问题。 Szeliski的方法用的多，paper我忘了，等会那篇review里面可以查到。

没有考虑的问题：
1． 假设是perfect 的各种关系，像频移，旋转，缩放，仿射，实际上并不那么简单，实际数码照的片子关系很复杂。
2． 光照。刚才说了，大问题。在估计匹配点，矩阵，甚至最后融合都会引入不少麻烦。

SIFT 方法好，可以解决一些问题，

经典的一篇综述
Image alignment and stitching- A tutorial

http://hi.baidu.com/simonyuee/blog/item/24961f6dfaa543fa431694bd.html

---

最佳缝合线算法有助于消除鬼影，以得到较好的视觉效果，这里转载该博主的文章，方便自己查阅学习，若有不妥，请联系本人。

原文：http://blog.csdn.net/wd1603926823/article/details/49536691

理论根据《图像拼接的改进算法》，据说这个算法可以消除重叠部分运动物体的重影和鬼影问题，所以就编下试试看，反正之前编的那种很老的取平均值法融合、渐入渐出（基于距离）融合、以及改进的三角函数权重融合都只是适合静态图像融合 有重影和鬼影问题 不适合有运动物体的图像融合，所以还是要最佳缝合线算法：看论文上就四步 很清晰也很好懂 结果自己写的时候才发现看起来很容易的也许编起来没那么容易 之前 想得太简单了。

缝合线算法：

function D=bestlinefusion(A,B)
%《图像拼接的改进算法》最佳缝合线算法 图像融合
%先根据之前得到的H矩阵计算重叠区域Rect
[H,W,k]=size(A);
rdata1=-118;
L=W+1+rdata1;
R=W;
n=R-L+1;
%计算得到重叠区域的差值图像 其实我不懂计算差值图像干嘛 只要计算重叠区域的大小就好了 为什么还要计算差值图 后面又没用到
Rect=zeros(H,n);
for i=1:H
for j=L:R
Rect(i,j-L+1)=A(i,j)-B(i,j-L+1);
end
end
Rect=uint8(Rect);%这句要不要呢？
%最终融合图的大小
rdata1=-118;
rdata2=3;
Y=2W+rdata1+1;
D=zeros(H,Y);
%放路径的矩阵
path=zeros(H,n);
%放强度值 每条路径的强度值strength=color^2+geometry
color=zeros(1,n);
geometry=zeros(1,n);
strength1=zeros(1,n);
strength2=zeros(1,n);
%计算第一行即初始化的强度值
for j=L:R
y=j-L+1;
color(y)=A(1,j)-B(1,y);
if(y==1)
Bxdao=B(1,y+1)+2B(2,y+1);
Bydao=B(2,y)+2B(2,y+1);
Aydao=2A(2,j-1)+A(2,j)+2A(2,j+1);
Axdao=A(1,j+1)+2A(2,j+1)-A(1,j-1)-2A(2,j-1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength1(y)=color(y)^2+geometry(y);
path(1,y)=y;
continue;
end
if(jR)
Axdao=A(1,j-1)-2A(2,j-1);
Aydao=2A(2,j-1)+A(2,j);
Bxdao=B(1,y+1)+2B(2,y+1)-B(1,y-1)-2B(2,y-1);
Bydao=2B(2,y-1)+B(2,y)+2B(2,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength1(y)=color(y)^2+geometry(y);
path(1,y)=y;
continue;
end
Axdao=A(1,j+1)+2A(2,j+1)-A(1,j-1)-2A(2,j-1);
Bxdao=B(1,y+1)+2B(2,y+1)-B(1,y-1)-2B(2,y-1);
Aydao=2A(2,j-1)+A(2,j)+2A(2,j+1);
Bydao=2B(2,y-1)+B(2,y)+2B(2,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength1(y)=color(y)^2+geometry(y);
path(1,y)=y;
end
color=zeros(1,n);
geometry=zeros(1,n);
small=0;
%开始扩展 向下一行 从第二行到倒数第二行 最后一行单独拿出来 像第一行一样 因为它的结构差值geometry不好算
for i=2:H-1
%先把下一行的强度值全部计算出来 到时候需要比较哪三个就拿出哪三个
for j=L:R
x=i;
y=j-L+1;
color(y)=A(i,j)-B(x,y);
if(y1)
Axdao=2A(i-1,j+1)+A(i,j+1)+2A(i+1,j+1)-2A(i-1,j-1)-A(i,j-1)-2A(i+1,j-1);
Bxdao=2B(x-1,y+1)+B(x,y+1)+2B(x+1,y+1);
Aydao=-2A(i-1,j-1)-A(i-1,j)-2A(i-1,j+1)+2A(i+1,j-1)+A(i+1,j)+2A(i+1,j+1);
Bydao=-B(x-1,y)-2B(x-1,y+1)+B(x+1,y)+2B(x+1,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength2(y)=color(y)^2+geometry(y);
continue;
end
if(j==R)
Axdao=-2A(i-1,j-1)-A(i,j-1)-2A(i+1,j-1);
Bxdao=2B(x-1,y+1)+B(x,y+1)+2B(x+1,y+1)-2B(x-1,y-1)-B(x,y-1)-2B(x+1,y-1);
Aydao=-2A(i-1,j-1)-A(i-1,j)+2A(i+1,j-1)+A(i+1,j);
Bydao=-2B(x-1,y-1)-B(x-1,y)-2B(x-1,y+1)+2B(x+1,y-1)+B(x+1,y)+2B(x+1,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength2(y)=color(y)^2+geometry(y);
continue;
end
Axdao=2A(i-1,j+1)+A(i,j+1)+2A(i+1,j+1)-2A(i-1,j-1)-A(i,j-1)-2A(i+1,j-1);
Bxdao=2B(x-1,y+1)+B(x,y+1)+2B(x+1,y+1)-2B(x-1,y-1)-B(x,y-1)-2B(x+1,y-1);
Aydao=-2A(i-1,j-1)-A(i-1,j)-2A(i-1,j+1)+2A(i+1,j-1)+A(i+1,j)+2A(i+1,j+1);
Bydao=-2B(x-1,y-1)-B(x-1,y)-2B(x-1,y+1)+2B(x+1,y-1)+B(x+1,y)+2B(x+1,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength2(y)=color(y)^2+geometry(y);
end
for j=1:n
if(path(i-1,j)==1)
if(strength2(1)<strength2(2))
strength1(j)=strength1(j)+strength2(1);
path(i,j)=1;
else
strength1(j)=strength1(j)+strength2(2);
path(i,j)=2;
end
else
if(path(i-1,j)n)
if(strength2(n-1)<strength2(n))
strength1(j)=strength1(j)+strength2(n-1);
path(i,j)=n-1;
else
strength1(j)=strength1(j)+strength2(n);
path(i,j)=n;
end
else
if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)))
if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)+1))
small=strength2(path(i-1,j)-1);
path(i,j)=path(i-1,j)-1;
else
small=strength2(path(i-1,j)+1);
path(i,j)=path(i-1,j)+1;
end
else
if(strength2(path(i-1,j))<strength2(path(i-1,j)+1))
small=strength2(path(i-1,j));
path(i,j)=path(i-1,j);
else
small=strength2(path(i-1,j)+1);
path(i,j)=path(i-1,j)+1;
end
end
strength1(j)=strength1(j)+small;
end
end
small=0;
end
strength2=zeros(1,n);
color=zeros(1,n);
geometry=zeros(1,n);
end
%单独计算最后一行
i=H;
for j=L:R
x=i;
y=j-L+1;
color(y)=A(i,j)-B(x,y);
if(y1)
Axdao=2A(i-1,j+1)+A(i,j+1)-2A(i-1,j-1)-A(i,j-1);
Aydao=-2A(i-1,j-1)-A(i-1,j)-2A(i-1,j+1);
Bxdao=2B(x-1,y+1)+B(x,y+1);
Bydao=-B(x-1,y)-2B(x-1,y+1);
continue;
end
if(j==R)
Bxdao=2B(x-1,y+1)+B(x,y+1)-2B(x-1,y-1)-B(x,y-1);
Bydao=-2B(x-1,y-1)-B(x-1,y)-2B(x-1,y+1);
Axdao=-2A(i-1,j-1)-A(i,j-1);
Aydao=-2A(i-1,j-1)-A(i-1,j);
continue;
end
Axdao=2A(i-1,j+1)+A(i,j+1)-2A(i-1,j-1)-A(i,j-1);
Bxdao=2B(x-1,y+1)+B(x,y+1)-2B(x-1,y-1)-B(x,y-1);
Aydao=-2A(i-1,j-1)-A(i-1,j)-2A(i-1,j+1);
Bydao=-2B(x-1,y-1)-B(x-1,y)-2B(x-1,y+1);
geometry(y)=(Axdao-Bxdao)(Aydao-Bydao);
strength2(y)=color(y)^2+geometry(y);
end
for j=1:n
if(path(i-1,j)==1)
if(strength2(1)<strength2(2))
strength1(j)=strength1(j)+strength2(1);
path(i,j)=1;
else
strength1(j)=strength1(j)+strength2(2);
path(i,j)=2;
end
else
if(path(i-1,j)==n)
if(strength2(n-1)<strength2(n))
strength1(j)=strength1(j)+strength2(n-1);
path(i,j)=n-1;
else
strength1(j)=strength1(j)+strength2(n);
path(i,j)=n;
end
else
if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)))
if(strength2(path(i-1,j)-1)<strength2(path(i-1,j)+1))
small=strength2(path(i-1,j)-1);
path(i,j)=path(i-1,j)-1;
else
small=strength2(path(i-1,j)+1);
path(i,j)=path(i-1,j)+1;
end
else
if(strength2(path(i-1,j))<strength2(path(i-1,j)+1))
small=strength2(path(i-1,j));
path(i,j)=path(i-1,j);
else
small=strength2(path(i-1,j)+1);
path(i,j)=path(i-1,j)+1;
end
end
strength1(j)=strength1(j)+small;
end
end
small=0;
end
%比较strength1里放的每条路径的强度值的总和 谁最小 就选path中对应的那一列的路径
[minzhi,minth]=min(strength1);
mypath=path(:,minth);
%mypath放的就是最佳缝合线选出的路径 这条路径坐标是参考图A 右边是目标图B
for i=1:H
for j=1:mypath(i)+L-1
D(i,j,1)=A(i,j,1);
D(i,j,2)=A(i,j,2);
D(i,j,3)=A(i,j,3);
end
for j=mypath(i)+L-1:Y-1
x=i;
y=j-L+1;
D(i,j,1)=B(x,y,1);
D(i,j,2)=B(x,y,2);
D(i,j,3)=B(x,y,3);
end
end

D=uint8(D);

还是以那两幅图为例：

结果：

差值图像Rect：

可是最终得到的融合后的图像D怎么是下面这个鬼样子：

眼啊！这哪里最佳了？？！！

貌似理论部分没编错啊 哪里有问题呢 还是运行最佳缝合线本来就应该是这样？为了验证下 在另一篇也讲最佳缝合线算法的论文里 《全景图像拼接关键技术研究》下载了两张图片，看找出来最佳缝合线是否和这个作者的近似或者一样

原图：

通过放进OPENCV求匹配对 有32对匹配然后将这些匹配对输出给MATLAB计算H变换矩阵 我用了30对坐标 得到H矩阵 再进行H矩阵变换后：

然后用下渐入渐出融合 结果：

用下改进的三角函数权重的融合结果:

果然改进的肉眼看不出来效果 但它们两个都有重影 比较严重的重影！用下刚刚编的最佳缝合线算法看看 这是差值图像:

然后这是最佳缝合线得到的D：

可以看到用最佳缝合线后没有重影了 那棵树那里以及那棵大树旁边两棵树没有重影了 可以和渐入渐出以及改进的那个对比下 这个已经没有重影了 原来最佳缝合线还是有用的喔

我想把这幅图中找到的最佳缝合线画出来 ，其实我放大后看得到那条最佳缝，我想想怎么画出来：

L=W+1+rdata1;
imshow(D)
hold on;
for i=1:137
D(i,L+mypath(i),1)=255;
D(i,L+mypath(i),2)=255;
D(i,L+mypath(i),3)=255;
end
hold off;

这条白线就是用上面的最佳缝合线程序找到的最佳缝合线 感觉好像一道白色的闪电 上面看不太出来 用黑色画下：

这样就看得清了这就是找到的最佳缝合线 的确是避开了那些重影物体 我是这样来写的 其实就是和上面的H矩阵不一样 这个计算出来平移距离是76 竖直方向其实有5个像素的平移量 但我懒得改上面的了 直接只改了rdata1=-76 所以和那个论文中找的最佳缝合线有点点差别 不过没多大关系 因为我指望着像《图像拼接的改进算法》里说的那样找到最佳缝合线后还要进行多分辨率拼接的 就像之前我编渐入渐出时初步得到H矩阵变换后的结果也是有点没对齐的，但是融合时用渐入渐出后它就完完全全对齐了 这个也是一样的 我猜想后续的多分辨率拼接也可以做到这样的效果 虽然现在得到的最佳缝合线有点没对齐的样子 因为我没用竖直方向的平移量rdata2=5，其实加上一样的 但太麻烦了因为我再看一次程序 然后重新计算 我下次再去改！现在这样就行了 竖直方向没对齐的我就依靠后续的多分辨率算法好了 但当竖直方向的平移量rdata2很大时以及有旋转关系时 也要重新计算 不能只算一个rdata1 绝对不能这样 会影响后续多分辨率的效果！