计算机视觉-实验三（图像映射与全景拼接）

• • 一、全景拼接基本流程
  • 二、图像映射
  • 1.流程
  • 2.类型
  • 3.参数求解
  • 4.仿射变换
  • 5.单应性变化
  • 6.前向映射
  • 7.逆向映射
  • 8.像素插值
  • 三、RANSAC求解单应矩阵
  • RANSAC loop
  • 四、APAP算法
  • 五、图割方法
  • 六、图像融合
  • 七、图像拼接
  • 八、实验过程

1. 针对某个场景拍摄多张/序列图像
2. 计算第二张图像与第一张图像之间的变换关系
3. 将第二张图像叠加到第一张图像的坐标系中
4. 变换后的融合/合成
5. 在多图场景中，重复上述过程

在上述流程当中，我们可以发现图像之间的变换关系是一个关键部分。

1. 针对两张/多张图像提取特征
2. 特征匹配
3. 根据图像变换特点，选取合适的变换结构
4. 根据DLT等方法计算变换结构
5. 采用正向/逆向映射，利用插值方式实现图像映射变换

图像映射类型有平移、旋转、尺度变换、仿射、透视映射、桶装映射等，不同的类型所对应的自由度也不同，如下图所示：

单应性变换即一个平面到另一个平面的映射关系。公式如下所示：

得到最小二乘解，在平移问题中最小二乘解和平均残差得到的结果等价

残差：

平方和代价函数：

矩阵表达式：

对于 f ( x ) f(x) f(x)中每个像素x，根据变换模型计算相应的映射坐标 x ’ = h ( x ) x’=h(x) x’=h(x),并将x像素值赋给 g ( x ’ ) g(x’) g(x’)

对于 g ( x ′ ) g(x^{'}) g(x′)中的每个像素 x ′ x^{'} x′，根据变换模型计算相应的映射坐标 x = h − 1 ( x ′ ) x=h^{-1}(x') x=h−1(x′)

主要插值方法：

• 最近邻插值
• 线性插值
• 双线性插值
• 三线性插值
• …

在直线拟合时，我们需要根据给定的若干二维空间中的点，求出直线 y = a x + b y=ax+b y=ax+b，并使得该直线对空间点的拟合误差最小，但可能存在噪声点对结果造成很大影响，如图所示：

此时我们可以随机选取两个点，并根据选择点构造直线，再给定阈值计算出inliers的数量

APAP算法流程：

1. 通过SIFT得到两幅图像的匹配点对
2. 通过RANSAC剔除外点，得到N对内点
3. 利用DLT和SVD计算全局单应性
4. 将原图划分网格，取网格中心点，计算每个中心点和原图上内点之间的欧氏距离和权重
5. 将权重放到DLT算法的A矩阵中，构建新的W*A矩阵，重新SVD分解，得到当前网格的局部单应性矩阵
6. 遍历每个网格，利用局部单应性矩阵映射到全景画布上，得到APAP变换后的原图
7. 最后进行拼接线的加权融合
   

在寻找拼接缝时，我们希望找到一条不明显的拼接缝，而这就需要通过最大流最小割图割方法进行操作。
最大流：
给定指定有向图，其中有两个特殊的点源S(Sources)和汇T(Sinks)，每条边有指定的容量(Capacity)，求满足条件的从S到T的最大流(MaxFlow)。
最小割：
割是网络中定点的一个划分，它把网络中的所有顶点划分成两个顶点集合S和T，其中源点s∈S,汇点t∈T。记为CUT(S,T)，满足条件的从S到T的最小割（Min cut）。

在进行图像融合时我们通常采用multi-band blending策略进行融合
基本原理：
1.建立两幅图像的拉普拉斯金字塔
2.求高斯金字塔（掩模金字塔-为了拼接左右两幅图像）因为其具有尺度不变性
3.进行拼接blendLapPyrs() ; 在每一层上将左右laplacian图像直接拼起来得结果金字塔resultLapPyr
4.重建图像: 从最高层结果图
将左右laplacian图像拼成的resultLapPyr金字塔中每一层，从上到下插值放大并和下一层相加，即得blend图像结果

流程：

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

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 = ['./panoramicpic3/' + str(i + 1) + '.sift' for i in range(5)]
imname = ['./panoramicpic3/' + 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()

实验结果：