创建全景图

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

RANSAC

RANSAC 是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是 用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间 的单应性矩阵，RANSAC 基本的思想是，数据中包含正确的点和噪声点，合理的模 型应该能够在描述正确数据点的同时摒弃噪声点。 RANSAC 的标准例子：用一条直线拟合带有噪声数据的点集。简单的最小二乘在该 例子中可能会失效，但是 RANSAC 能够挑选出正确的点，然后获取能够正确拟合 的直线。该算法只选择了和直线模型一致的数据点，成功 地找到了正确的解。 

我们在任何模型中都可以使用 RANSAC 模块。在使用 RANSAC 模块时，我们只需 要在相应 Python 类中实现 fit() 和 get_error() 方法，剩下就是正确地使用 ransac.py。 我们这里使用可能的对应点集来自动找到用于全景图像的单应性矩阵。图 3-11 所示 为使用 SIFT 特征自动找到匹配对应。这可以通过运行下面的命令来实现：

import sift

featname = ['Univ'+str(i+1)+'.sift' for i in range(5)]

imname = ['Univ'+str(i+1)+'.jpg' for i in range(5)]

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])

显然，并不是所有图像中的对应点对都是正确的。实际上，SIFT 是具有很强稳健性 的描述子，能够比其他描述子，例如图像块相关的 Harris 角点，产生更少的错误的 匹配。但是该方法仍然远非完美。

我们使用 RANSAC 算法来求解单应性矩阵，首先需要将下面模型类添加到 homography.py 文件中：

class RansacModel(object):

""" 用于测试单应性矩阵的类，其中单应性矩阵是由网站 http://www.scipy.org/Cookbook/RANSAC 上 的 ransac.py 计算出来的 """

def __init__(self,debug=False):

self.debug = debug

def fit(self, data)：

""" 计算选取的 4 个对应的单应性矩阵 """ # 将其转置，来调用 H_from_points() 计算单应性矩阵 data = data.T # 映射的起始点

fp = data[:3,:4]

# 映射的目标点

tp = data[3:,:4]

# 计算单应性矩阵，然后返回

return H_from_points(fp,tp)

def get_error( self, data, H):

""" 对所有的对应计算单应性矩阵，然后对每个变换后的点，返回相应的误差 """

data = data.T

# 映射的起始点

fp = data[:3]

# 映射的目标点

tp = data[3:]

# 变换

fp fp_transformed = dot(H,fp)

# 归一化齐次坐标

for i in range(3): fp_transformed[i] /= fp_transformed[2]

# 返回每个点的误差

return sqrt( sum((tp-fp_transformed)**2,axis=0) )

可以看到，这个类包含 fit() 方法。

该方法仅仅接受由 ransac.py 选择的４个对应点 对（data 中的前４个点对），然后拟合一个单应性矩阵。记住，４个点对是计算单 应性矩阵所需的最少数目。由于 get_error() 方法对每个对应点对使用该单应性矩 阵，然后返回相应的平方距离之和，因此 RANSAC 算法能够判定哪些点对是正确 的，哪些是错误的。在实际中，我们需要在距离上使用一个阈值来决定哪些单应性 矩阵是合理的。为了方便使用，将下面的函数添加到 homography.py 文件中：

def H_from_ransac(fp,tp,model,maxiter=1000,match_theshold=10):

""" 使用 RANSAC 稳健性估计点对应间的单应性矩阵 H

# 输入：齐次坐标表示的点 fp，tp（3×n 的数组）"""

import ransac

# 对应点组 data = vstack((fp,tp))

# 计算 H，并返回 H,

ransac_data = ransac.ransac(data.T,model,4,maxiter,match_theshold,10, return_all=True) return H,

ransac_data['inliers']

该函数同样允许提供阈值和最小期望的点对数目。最重要的参数是最大迭代次数： 程序退出太早可能得到一个坏解；迭代次数太多会占用太多时间。函数的返回结果 为单应性矩阵和对应该单应性矩阵的正确点对。 类似于下面的操作，你可以将 RANSAC 算法应用于对应点对上： # 将匹配转换成齐次坐标点的函数

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)

return fp,tp

# 估计单应性矩阵

model = homography.RansacModel()

fp,tp = convert_points(1)

H_12 = homography.H_from_ransac(fp,tp,model)[0]

# im1 到 im2 的单应性矩阵

fp,tp = convert_points(0)

H_01 = homography.H_from_ransac(fp,tp,model)[0]

# im0 到 im1 的单应性矩阵

tp,fp = convert_points(2)

# 注意：点是反序的

H_32 = homography.H_from_ransac(fp,tp,model)[0]

# im3 到 im2 的单应性矩阵

tp,fp = convert_points(3)

# 注意：点是反序的

H_43 = homography.H_from_ransac(fp,tp,model)[0]

# im4 到 im3 的单应性矩阵

拼接图像

估计出图像间的单应性矩阵（使用 RANSAC 算法），现在我们需要将所有的图像扭 曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面（否则，需要 进行大量变形）。一种方法是创建一个很大的图像，比如图像中全部填充 0，使其和 中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平 旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域82 ｜ 第 3 章 填充０，以便为扭曲的图像腾出空间。将下面的代码添加到 warp.py 文件中： def panorama(H,fromim,toim,padding=2400,delta=2400): """ 使用单应性矩阵 H（使用 RANSAC 健壮性估计得出），协调两幅图像，创建水平全景图像。

，根据特征点匹配的方式，则利用这些匹配的点来估算单应矩阵（使用上面的RANSAC算法，也就是把其中一张通过个关联性和另一张匹配的方法。通过单应矩阵H，可以将原图像中任意像素点坐标转换为新坐标点，转换后的图像即为适合拼接的结果图像。可以简单分为以下几步：
1.根据给定图像/集，实现特征匹配。
2.通过匹配特征计算图像之间的变换结构。
3.利用图像变换结构，实现图像映射。
4.针对叠加后的图像，采用APAP之类的算法，对齐特征点。(图像配准)
5.通过图割方法，自动选取拼接缝。
6.根据multi-band blending策略实现融合。

3、图像融合

因为相机和光照强度的差异，会造成一幅图像内部，以及图像之间亮度的不均匀，拼接后的图像会出现明暗交替，这样给观察造成极大的不便。 亮度与颜色均衡处理，通常的处理方式是通过相机的光照模型，校正一幅图像内部的光照不均匀性，然后通过相邻两幅图像重叠区域之间的关系，建立相邻两幅图像之间直方图映射表，通过映射表对两幅图像做整体的映射变换，最终达到整体的亮度和颜色的一致性。

实现步骤

读入连续图片并使用SIFT特征查找匹配对应点对
import sift
#sift程序应与运行图片在同一文件夹下
featname = ['C:/Users/Administrator/Documents/python/jmu'+str(i+1)+'.sift' for i in range(5)]
imname = ['C:/Users/Administrator/Documents/python/jmu'+str(i+1)+'.jpg' for i in range(5)]
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])

class RansacModel(object):
""" Class for testing homography fit with ransac.py from
http://www.scipy.org/Cookbook/RANSAC"""
def __init__(self,debug=False):
self.debug = debugdef fit(self, data):
""" 计算选取四个对应的单应性矩阵 """

# 将其转置，来调用H_from_points()计算单应性矩阵
data = data.T
#映射的起始点
fp = data[:3,:4]
# 映射的目标点
tp = data[3:,:4]
#计算单应性矩阵然后返回
return H_from_points(fp,tp)

def get_error( self, data, H):
""" 对所有的对应计算单应性矩阵，然后对每个变换后的点，返回相应的误差 """
data = data.T
#映射的起始点
fp = data[:3]
# 映射的目标点
tp = data[3:]
#变换fp
fp_transformed = dot(H,fp)
#归一化齐次坐标
for i in range(3):
fp_transformed[i] /= fp_transformed[2]
return sqrt( sum((tp-fp_transformed)**2,axis=0) )

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 = ['C:/Users/Administrator/Documents/python/jmu'+str(i+1)+'.sift' for i in range(5)]
imname = ['C:/Users/Administrator/Documents/python/jmu'+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([fp[1],fp[0],fp[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 = 1200 # 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')
savefig("example.png",dpi=300)
show()