python cv 全景拼接

python cv 全景拼接

一、全景拼接介绍

在同一位置（即图像的照相机位置相同）拍摄的两幅或者多幅图像是单应性相关的。我们经常使用该约束将很多图像缝补，来拼成一个大的图像来创建全景图像。

1.1 流程：

1.根据给定图像/集，实现特征匹配。
2.通过匹配特征计算图像之间的变换结构。
3.利用图像变换结构，实现图像映射。
4.针对叠加后的图像，采用APAP之类的算法，对齐特征点。(图像配准)
5.通过图割方法，自动选取拼接缝。
6.根据multi-band blending策略实现融合。

1.2 sift匹配与RANSAC

SIFT算法的描述子稳健性很强，比Harris角点要来得精确，但是它的匹配正确率也并不是百分百的，会受到一些噪声点的干扰，有时就会因为不同地方有类似的图案导致匹配错误。那么RANSAC算法便是用来找到正确模型来拟合带有噪声数据的迭代方法。RANSAC通过反复取样，也就是从整个观测数据中随机抽一些数据估算模型参数之后看和所有数据误差有多大，然后取误差最小视为最好以及分离内群与离群数据。基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。
举个简单的例子，直线的拟合便是RANSAC的一个标准化体现：

在所有的点中随机选择两个点。
根据这两个点作一条直线。
设定一个阈值，计算在这条线上的点的数量，记为inliners。
根据最大的inliners那条线进行后续计算。
同理，RANSAC算法可以应用到其它模块中，例如用于图像变换的单应性矩阵的计算。
在拼接的过程中，通过将响速和单应矩阵H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，我们可以决定应该将该图像填补到左边还是右边。当该图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边”情况中，需要在单应矩阵中加入平移。

1.3 图像配准

图像配准是对图像进行变换，使变换后的图像能够在很好的拼接在上一张图片的坐标系。
为了能够进行图像对比和更精细的图像分析，图像配准是一步非常重要的操作。

1.4图割方法

最大流最小割算法原理，
1.最小割问题
一个有向图，并有一个源顶点（source vertex）和目标顶点（target vertex）.边的权值为正,又称之为容量（capacity）。如下图

一个st-cut(简称割cut)会把有向图的顶点分成两个不相交的集合，其中s在一个集合中，t在另外一个集合中。
这个割的容量（capacity of the cut）就是A到B所有边的容量和。注意这里不包含B到A的。最小割问题就是要找到割容量最小的情况。
2.最大流问题
跟mincut问题类似，maxflow要处理的情况也是一个有向图，并有一个原顶点（source vertex）和目标（target vertex），边的权值为正,又称之为容量（capacity）。
（1）初始化，所有边的flow都初始化为0。
（2）沿着增广路径增加flow。增广路径是一条从s到t的无向路径，但也有些条件，可以经过没有满容量的前向路径（s到t）或者是不为空的反向路径(t->s)。

1.5 图像融合

图像拼接之后可以发现，在拼接的交界处有明显的衔接痕迹，存在边缘效应，因为光照色泽的原因使得图片交界处的过渡很糟糕，所以需要特定的处理解决这种不自然。那么这时候可以采用blending方法。multi-band blending是目前图像融和方面比较好的方法。
原理：
1.建立两幅图像的拉普拉斯金字塔
2.求高斯金字塔（掩模金字塔-为了拼接左右两幅图像）因为其具有尺度不变性

3.进行拼接blendLapPyrs() ; 在每一层上将左右laplacian图像直接拼起来得结果金字塔resultLapPyr
4.重建图像: 从最高层结果图将左右laplacian图像拼成的resultLapPyr金字塔中每一层，从上到下插值放大并和下一层相加，即得blend图像结果（reconstructImgFromLapPyramid）

二、代码

from pylab import *
from numpy import *
from PIL import Image

# If you have PCV installed, these imports should work
from PCV.geometry import homography, warp
from PCV.localdescriptors import sift

"""
This is the panorama example from section 3.3.
"""

# set paths to data folder
featname = ['ch3/ima3/'+str(i+1)+'.sift' for i in range(5)] 
imname = ['ch3/ima3/'+str(i+1)+'.jpg' for i in range(5)]

# extract features and match
l = {}
d = {}
for i in range(5): 
    sift.process_image(imname[i],featname[i])
    l[i],d[i] = sift.read_features_from_file(featname[i])

matches = {}
for i in range(4):
    matches[i] = sift.match(d[i+1],d[i])

# visualize the matches (Figure 3-11 in the book)
for i in range(4):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i+1]))
    figure()
    sift.plot_matches(im2,im1,l[i+1],l[i],matches[i],show_below=True)


# function to convert the matches to hom. points
def convert_points(j):
    ndx = matches[j].nonzero()[0]
    fp = homography.make_homog(l[j+1][ndx,:2].T) 
    ndx2 = [int(matches[j][i]) for i in ndx]
    tp = homography.make_homog(l[j][ndx2,:2].T) 
    
    # switch x and y - TODO this should move elsewhere
    fp = vstack([fp[1],fp[0],fp[2]])
    tp = vstack([tp[1],tp[0],tp[2]])
    return fp,tp


# estimate the homographies
model = homography.RansacModel() 

fp,tp = convert_points(1)
H_12 = homography.H_from_ransac(fp,tp,model)[0] #im 1 to 2 

fp,tp = convert_points(0)
H_01 = homography.H_from_ransac(fp,tp,model)[0] #im 0 to 1 

tp,fp = convert_points(2) #NB: reverse order
H_32 = homography.H_from_ransac(fp,tp,model)[0] #im 3 to 2 

tp,fp = convert_points(3) #NB: reverse order
H_43 = homography.H_from_ransac(fp,tp,model)[0] #im 4 to 3    


# warp the images
delta = 2000 # for padding and translation

im1 = array(Image.open(imname[1]), "uint8")
im2 = array(Image.open(imname[2]), "uint8")
im_12 = warp.panorama(H_12,im1,im2,delta,delta)

im1 = array(Image.open(imname[0]), "f")
im_02 = warp.panorama(dot(H_12,H_01),im1,im_12,delta,delta)

im1 = array(Image.open(imname[3]), "f")
im_32 = warp.panorama(H_32,im1,im_02,delta,delta)

im1 = array(Image.open(imname[4]), "f")
im_42 = warp.panorama(dot(H_32,H_43),im1,im_32,delta,2*delta)


figure()
imshow(array(im_42, "uint8"))
axis('off')
show()

运行效果