图像拼接

---

前言

经过前两次实验，我们了解了一种局部图像描述方法SIFT，图像到图像的多种变换方式。接下来，要通过结合两者，来实现图像拼接

---

RANSAC

RANSAC是“RANdom SAmple Consensus”（随即一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC的基本思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。

拼接

1.简单实现图像拼接

import cv2
import numpy as np
from matplotlib import pyplot as plt
import time
import cv_read

# 计时开始
start_time = time.time()

#看看效果
MIN = 10
img1 = cv2.imread('pic/1.jpg')  # query
img2 = cv2.imread('pic/2.jpg')  # train
# SIFT特征点描述
sift = cv2.xfeatures2d.SIFT_create()
# 可以改为SURF
# surf=cv2.xfeatures2d.SURF_create(10000, nOctaves=4, extended=False, upright=True)
kp1, descrip_1 = sift.detectAndCompute(img1, None)
kp2, descrip_2 = sift.detectAndCompute(img2, None)
# FLANN算法 快速fast·近似approximate·近邻nearest neighbors·算法库library
# K-D树算法 k-d TreeAlgorithm
FLANN_INDEX_KDTREE = 0
# 参数algorithm用来指定匹配所使用的算法 参数tree用来指定算法中树的数量
indexParams = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
# 参数SearchParams用来指定递归遍历的次数。值越高结果越准确，耗时也越长。
searchParams = dict(checks=50)
flann = cv2.FlannBasedMatcher(indexParams, searchParams)
match = flann.knnMatch(descrip_1, descrip_2, k=2)  # 最近邻算法 k nearest neighbor
# 建立一个列表保存优秀的特征索引
good = []
for i, (m, n) in enumerate(match):
    if (m.distance < 0.75 * n.distance):  # 欧氏距离比率 匹配准则
        good.append(m)

if len(good) > MIN:
# 图二作为目标图像
    src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)  # 行转列,queryIdx查询点的索引升序
    ano_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)  # trainIdx是被查询点的索引
    H, mask = cv2.findHomography(src_pts, ano_pts, cv2.RANSAC, 5.0) # 得到单应性矩阵H
    warpImg = cv2.warpPerspective(img2, np.linalg.inv(H),
                                  (img1.shape[1] + img2.shape[1], img1.shape[0])) # 将单应性变换应用到图像img2上
    direct = warpImg.copy()
    direct[0:img1.shape[0], 0:img1.shape[1]] = img1

    simple = time.time()

    cv2.namedWindow('Result', cv2.WINDOW_NORMAL)
    cv2.imshow('Result',warpImg)

    rows, cols = img1.shape[:2]
    for col in range(0, cols):
        if img1[:, col].any() and warpImg[:, col].any():  # 开始重叠的最左端
            left = col
            break
    for col in range(cols - 1, 0, -1):
        if img1[:, col].any() and warpImg[:, col].any():  # 重叠的最右一列
            right = col
            break

    res = np.zeros([rows, cols, 3], np.uint8)
    for row in range(0, rows):
        for col in range(0, cols):
            if not img1[row, col].any():  # 如果没有原图，用旋转的填充
                res[row, col] = warpImg[row, col]
            elif not warpImg[row, col].any():  # 如果没有Warp图，用原图替代
                res[row, col] = img1[row, col]
            else:
                srcImgLen = float(abs(col - left))
                testImgLen = float(abs(col - right))
                alpha = srcImgLen / (srcImgLen + testImgLen)
                res[row, col] = np.clip(img1[row, col] * (1 - alpha) + warpImg[row, col] * alpha, 0, 255)

    warpImg[0:img1.shape[0], 0:img1.shape[1]] = res
    final = time.time()
    img3 = cv2.cvtColor(direct, cv2.COLOR_BGR2RGB)
    plt.imshow(img3, ), plt.show()
    img4 = cv2.cvtColor(warpImg, cv2.COLOR_BGR2RGB)
    plt.imshow(img4, ), plt.show()
    print('simple stich cost %f' % (simple - start_time))
    print('\ntotal cost %f' % (final - start_time))
    cv2.imwrite(r'pic\simplepanorma.png', direct)
    cv2.imwrite(r'pic\bestpanorma.png', warpImg)

else:
    print('not enough matches!')

实验用的照片：


运行结果如下：

2.全景图像拼接

图像拼接主函数：

from numpy import *
from matplotlib.pyplot import *
from PIL import Image
import homography
import warp
from PCV.localdescriptors import sift


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

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

# 将匹配转换成齐次坐标点的函数
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


# 估计单应性矩阵
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

# 扭曲图像
# img delta
# delta = 100  # 用于填充和平移 for padding and translation
# libraryimg delta
delta = 1000

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

单应性矩阵：

from numpy import *


# 1.单应性变换
def normallize(points):
    """在齐次坐标意义下，对点集进行归一化，是最后一行为1"""
    for row in points:
        row /= points[-1]
    return points


def make_homog(points):
    """将点集（dim×n的数组）转换为齐次坐标表示"""
    return vstack((points, ones((1, points.shape[1]))))


# 1.1直接线性变换算法
def H_from_points(fp, tp):
    """使用线性DLT方法，计算单应性矩阵H，使fp映射到tp。点自动进行归一化"""

    if fp.shape != tp.shape:
        raise RuntimeError("number of points do not match")

    # 对点进行归一化（对数值计算很重要）
    # --- 映射起始点 ---
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1)) + 1e-9
    C1 = diag([1 / maxstd, 1 / maxstd, 1])
    C1[0][2] = -m[0] / maxstd
    C1[1][2] = -m[1] / maxstd
    fp = dot(C1, fp)

    # --- 映射对应点 ---
    m = mean(tp[:2], axis=1)
    maxstd = max(std(tp[:2], axis=1)) + 1e-9
    C2 = diag([1 / maxstd, 1 / maxstd, 1])
    C2[0][2] = -m[0] / maxstd
    C2[1][2] = -m[1] / maxstd
    tp = dot(C2, tp)

    # 创建用于线性方法的矩阵，对于每个对应对，在矩阵中会出现两行数值
    nbr_correspondences = fp.shape[1]
    A = zeros((2 * nbr_correspondences, 9))
    for i in range(nbr_correspondences):
        A[2 * i] = [-fp[0][i], -fp[1][i], -1, 0, 0, 0,
                    tp[0][i] * fp[0][i], tp[0][i] * fp[1][i], tp[0][i]]
        A[2 * i + 1] = [0, 0, 0, -fp[0][i], -fp[1][i], -1,
                        tp[1][i] * fp[0][i], tp[1][i] * fp[1][i], tp[1][i]]

    U, S, V = linalg.svd(A)
    H = V[8].reshape((3, 3))

    # 反归一化
    H = dot(linalg.inv(C2), dot(H, C1))

    # 归一化，然后返回
    return H / H[2, 2]


# 1.2 仿射变换，计算H矩阵
def Haffine_from_points(fp, tp):
    """计算H仿射变换，使得tp是fp经过仿射变换H得到的"""

    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')

    # 对点进行归一化（对数值计算很重要）
    # --- 映射起始点 ---
    m = mean(fp[:2], axis=1)
    maxstd = max(std(fp[:2], axis=1)) + 1e-9
    C1 = diag([1 / maxstd, 1 / maxstd, 1])
    C1[0][2] = -m[0] / maxstd
    C1[1][2] = -m[1] / maxstd
    fp_cond = dot(C1, fp)

    # --- 映射对应点 ---
    m = mean(tp[:2], axis=1)
    C2 = C1.copy()  # 两个点集，必须都进行相同的缩放
    C2[0][2] = -m[0] / maxstd
    C2[1][2] = -m[1] / maxstd
    tp_cond = dot(C2, tp)

    # 因为归一化后点的均值为0，所以平移量为0
    A = concatenate((fp_cond[:2], tp_cond[:2]), axis=0)
    U, S, V = linalg.svd(A.T)

    # 如Hartley和Zisserman著的Multiplr View Geometry In Computer,Scond Edition所示，
    # 创建矩阵B和C
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]

    tmp2 = concatenate((dot(C, linalg.pinv(B)), zeros((2, 1))), axis=1)
    H = vstack((tmp2, [0, 0, 1]))

    # 反归一化
    H = dot(linalg.inv(C2), dot(H, C1))

    return H / H[2, 2]


# RANSAC算法求解单应性矩阵
class RanSacModel(object):
    """
    用于测试单应性矩阵的类，其中单应性矩阵是由网站www.spicy.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))


def H_from_ransac(fp, tp, model, maxiter=1000, match_theshpld=10):
    """
    使用RANSAC稳健性估计点对应间的单应性矩阵H（ransac.py为从
    www.scipy.org/Cookbook/RANSAC下载的版本）
    """

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

    import ransac

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

    # 计算H，并返回
    H, ransac_data = ransac.ransac(data.T, model, 4, maxiter, match_theshpld, 10,
                                   return_all=True)

    return H, ransac_data['inliers']

RANSAC:

import numpy
import scipy # use numpy if scipy unavailable
import scipy.linalg # use numpy if scipy unavailable

## Copyright (c) 2004-2007, Andrew D. Straw. All rights reserved.

## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:

##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.

##     * Redistributions in binary form must reproduce the above
##       copyright notice, this list of conditions and the following
##       disclaimer in the documentation and/or other materials provided
##       with the distribution.

##     * Neither the name of the Andrew D. Straw nor the names of its
##       contributors may be used to endorse or promote products derived
##       from this software without specific prior written permission.

## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

def ransac(data,model,n,k,t,d,debug=False,return_all=False):
    """fit model parameters to data using the RANSAC algorithm
    
This implementation written from pseudocode found at
http://en.wikipedia.org/w/index.php?title=RANSAC&oldid=116358182

{{{
Given:
    data - a set of observed data points
    model - a model that can be fitted to data points
    n - the minimum number of data values required to fit the model
    k - the maximum number of iterations allowed in the algorithm
    t - a threshold value for determining when a data point fits a model
    d - the number of close data values required to assert that a model fits well to data
Return:
    bestfit - model parameters which best fit the data (or nil if no good model is found)
iterations = 0
bestfit = nil
besterr = something really large
while iterations < k {
    maybeinliers = n randomly selected values from data
    maybemodel = model parameters fitted to maybeinliers
    alsoinliers = empty set
    for every point in data not in maybeinliers {
        if point fits maybemodel with an error smaller than t
             add point to alsoinliers
    }
    if the number of elements in alsoinliers is > d {
        % this implies that we may have found a good model
        % now test how good it is
        bettermodel = model parameters fitted to all points in maybeinliers and alsoinliers
        thiserr = a measure of how well model fits these points
        if thiserr < besterr {
            bestfit = bettermodel
            besterr = thiserr
        }
    }
    increment iterations
}
return bestfit
}}}
"""
    iterations = 0
    bestfit = None
    besterr = numpy.inf
    best_inlier_idxs = None
    while iterations < k:
        maybe_idxs, test_idxs = random_partition(n,data.shape[0])
        maybeinliers = data[maybe_idxs,:]
        test_points = data[test_idxs]
        maybemodel = model.fit(maybeinliers)
        test_err = model.get_error( test_points, maybemodel)
        also_idxs = test_idxs[test_err < t] # select indices of rows with accepted points
        alsoinliers = data[also_idxs,:]
        if debug:
            print('test_err.min()',test_err.min())
            print('test_err.max()',test_err.max())
            print('numpy.mean(test_err)',numpy.mean(test_err))
            print('iteration %d:len(alsoinliers) = %d'%(
                iterations,len(alsoinliers)))
        if len(alsoinliers) > d:
            betterdata = numpy.concatenate( (maybeinliers, alsoinliers) )
            bettermodel = model.fit(betterdata)
            better_errs = model.get_error( betterdata, bettermodel)
            thiserr = numpy.mean( better_errs )
            if thiserr < besterr:
                bestfit = bettermodel
                besterr = thiserr
                best_inlier_idxs = numpy.concatenate( (maybe_idxs, also_idxs) )
        iterations+=1
    if bestfit is None:
        raise ValueError("did not meet fit acceptance criteria")
    if return_all:
        return bestfit, {'inliers':best_inlier_idxs}
    else:
        return bestfit

def random_partition(n,n_data):
    """return n random rows of data (and also the other len(data)-n rows)"""
    all_idxs = numpy.arange( n_data )
    numpy.random.shuffle(all_idxs)
    idxs1 = all_idxs[:n]
    idxs2 = all_idxs[n:]
    return idxs1, idxs2

class LinearLeastSquaresModel:
    """linear system solved using linear least squares

    This class serves as an example that fulfills the model interface
    needed by the ransac() function.
    
    """
    def __init__(self,input_columns,output_columns,debug=False):
        self.input_columns = input_columns
        self.output_columns = output_columns
        self.debug = debug
    def fit(self, data):
        A = numpy.vstack([data[:,i] for i in self.input_columns]).T
        B = numpy.vstack([data[:,i] for i in self.output_columns]).T
        x,resids,rank,s = scipy.linalg.lstsq(A,B)
        return x
    def get_error( self, data, model):
        A = numpy.vstack([data[:,i] for i in self.input_columns]).T
        B = numpy.vstack([data[:,i] for i in self.output_columns]).T
        B_fit = scipy.dot(A,model)
        err_per_point = numpy.sum((B-B_fit)**2,axis=1) # sum squared error per row
        return err_per_point
        
def test():
    # generate perfect input data

    n_samples = 500
    n_inputs = 1
    n_outputs = 1
    A_exact = 20*numpy.random.random((n_samples,n_inputs) )
    perfect_fit = 60*numpy.random.normal(size=(n_inputs,n_outputs) ) # the model
    B_exact = scipy.dot(A_exact,perfect_fit)
    assert B_exact.shape == (n_samples,n_outputs)

    # add a little gaussian noise (linear least squares alone should handle this well)
    A_noisy = A_exact + numpy.random.normal(size=A_exact.shape )
    B_noisy = B_exact + numpy.random.normal(size=B_exact.shape )

    if 1:
        # add some outliers
        n_outliers = 100
        all_idxs = numpy.arange( A_noisy.shape[0] )
        numpy.random.shuffle(all_idxs)
        outlier_idxs = all_idxs[:n_outliers]
        non_outlier_idxs = all_idxs[n_outliers:]
        A_noisy[outlier_idxs] =  20*numpy.random.random((n_outliers,n_inputs) )
        B_noisy[outlier_idxs] = 50*numpy.random.normal(size=(n_outliers,n_outputs) )

    # setup model

    all_data = numpy.hstack( (A_noisy,B_noisy) )
    input_columns = range(n_inputs) # the first columns of the array
    output_columns = [n_inputs+i for i in range(n_outputs)] # the last columns of the array
    debug = False
    model = LinearLeastSquaresModel(input_columns,output_columns,debug=debug)

    linear_fit,resids,rank,s = scipy.linalg.lstsq(all_data[:,input_columns],
                                                  all_data[:,output_columns])

    # run RANSAC algorithm
    ransac_fit, ransac_data = ransac(all_data,model,
                                     50, 1000, 7e3, 300, # misc. parameters
                                     debug=debug,return_all=True)
    if 1:
        import pylab

        sort_idxs = numpy.argsort(A_exact[:,0])
        A_col0_sorted = A_exact[sort_idxs] # maintain as rank-2 array

        if 1:
            pylab.plot( A_noisy[:,0], B_noisy[:,0], 'k.', label='data' )
            pylab.plot( A_noisy[ransac_data['inliers'],0], B_noisy[ransac_data['inliers'],0], 'bx', label='RANSAC data' )
        else:
            pylab.plot( A_noisy[non_outlier_idxs,0], B_noisy[non_outlier_idxs,0], 'k.', label='noisy data' )
            pylab.plot( A_noisy[outlier_idxs,0], B_noisy[outlier_idxs,0], 'r.', label='outlier data' )
        pylab.plot( A_col0_sorted[:,0],
                    numpy.dot(A_col0_sorted,ransac_fit)[:,0],
                    label='RANSAC fit' )
        pylab.plot( A_col0_sorted[:,0],
                    numpy.dot(A_col0_sorted,perfect_fit)[:,0],
                    label='exact system' )
        pylab.plot( A_col0_sorted[:,0],
                    numpy.dot(A_col0_sorted,linear_fit)[:,0],
                    label='linear fit' )
        pylab.legend()
        pylab.show()

if __name__=='__main__':
    test()
    

图像扭曲:

from numpy import *
import numpy as np
from scipy import ndimage
from scipy.spatial import Delaunay
from matplotlib.pyplot import *
from PIL import Image
import matplotlib.pyplot as plt
import homography

# 将图像放到图像中
def image_in_image(im1, im2, tp):
    """使用仿射变换将im1放置在im2上，使im1图像的角和tp尽可能的靠近
        tp是齐次表示的，并且是按照从左上角逆时针计算的"""

    # 扭曲的点
    m, n = im1.shape[:2]
    fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])

    # 计算仿射变换，并且将其应用于图像im1中
    # 计算H矩阵
    H = homography.Haffine_from_points(tp, fp)
    im1_t = ndimage.affine_transform(im1, H[:2, :2],
                                     (H[0, 2], H[1, 2]), im2.shape[:2])
    alpha = (im1_t > 0)

    return (1 - alpha) * im2 + alpha * im1_t


# 三角形的alpha图像
def alpha_for_triangle(points, m, n):
    """对于带有由points定义角点的三角形，创建大小为（m,n）的alpha图
    （在归一化的齐次坐标意义下）"""

    alpha = zeros((m, n))

    points = np.array(points, dtype=np.int)

    for i in range(min(points[0]), max(points[0])):
        for j in range(min(points[1]), max(points[1])):
            x = linalg.solve(points, [i, j, 1])
            if min(x) > 0:
                alpha[i, j] = 1
    return alpha


# 三角剖分的函数
def triangulate_points(x, y):
    """二维点的 Delaunay 三角剖分"""
    tri = Delaunay(np.c_[x, y]).simplices
    return tri


def pw_affine(fromim, toim, fp, tp, tri):
    """ 从一幅图像中扭曲矩形图像块
    fromim= 将要扭曲的图像
    toim= 目标图像
    fp= 齐次坐标表示下，扭曲前的点
    tp= 齐次坐标表示下，扭曲后的点
    tri= 三角剖分 """

    im = toim.copy()

    # 检查图像是灰度图像还是彩色图象
    is_color = len(fromim.shape) == 3

    # 创建扭曲后的图像（如果需要对彩色图像的每个颜色通道进行迭代操作，那么有必要这样做）
    im_t = zeros(im.shape, 'u8')

    for t in tri:
        # 计算仿射变换
        H = homography.Haffine_from_points(tp[:, t], fp[:, t])

        if is_color:
            for col in range(fromim.shape[2]):
                im_t[:, :, col] = ndimage.affine_transform(
                    fromim[:, :, col], H[:2, :2], (H[0, 2], H[1, 2]), im.shape[:2])
                # im1_t = ndimage.affine_transform(im1, H[:2, :2],
                #                                  (H[0, 2], H[1, 2]), im2.shape[:2])
        else:
            im_t = ndimage.affine_transform(
                fromim, H[:2, :2], (H[0, 2], H[1, 2]), im.shape[:2])
        # 三角形的 alpha
        alpha = alpha_for_triangle(tp[:, t], im.shape[0], im.shape[1])

        # 将三角形加入到图像中
        im[alpha > 0] = im_t[alpha > 0]
    return im


# 绘制三角形
def plot_mesh(x, y, tri):
    """ 绘制三角形 """
    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]]  # 将第一个点加入到最后
        plot(x[t_ext], y[t_ext], 'r')

# 图像拼接
def panorama(H, fromim, toim, padding=2400, delta=2400):
    """
    使用单应性矩阵H（使用RANSAC健壮性估计得出），协调两幅图像，创建水平全景图形。结果
    为一副和toim具有相同高度的图像。padding指定填充像素的数目，delta指定额外的平移量
    """

    # 检查图像是灰度图像，还是彩色图像
    is_color = len(fromim.shape) == 3

    # 用于geometric_transform() 的单应性变换
    def transf(p):
        p2 = dot(H, [p[0], p[1], 1])
        return (p2[0] / p2[2], p2[1] / p2[2])

    if H[1, 2] < 0:  # fromim在右边
        print("warp - right")
        # 变换fromim
        if is_color:
            # 在目标图像的右边填充0
            toim_t = hstack((toim, zeros((toim.shape[0], padding, 3))))
            fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的右边填充0
            toim_t = hstack((toim, zeros((toim.shape[0], padding))))
            fromim_t = ndimage.geometric_transform(fromim, transf,
                                                   (toim.shape[0], toim.shape[1] + padding))
    else:
        print("warp - left")
        # 为了补偿填充效果，在左边加入平移量
        H_delta = array([[1, 0, 0], [0, 1, -delta], [0, 0, 1]])
        H = dot(H, H_delta)
        # fromim 变换
        if is_color:
            # 在图像左边填充0
            toim_t = hstack((zeros((toim.shape[0], padding, 3)), toim))
            fromim_t = zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
            for col in range(3):
                fromim_t[:, :, col] = ndimage.geometric_transform(fromim[:, :, col],
                                                                  transf, (toim.shape[0], toim.shape[1] + padding))
        else:
            # 在目标图像的左边填充0
            toim_t = hstack((zeros((toim.shape[0], padding)), toim))
            fromim_t = ndimage.geometric_transform(fromim, transf,
                                                   (toim.shape[0], toim.shape[1] + padding))

    # 协调后返回（将fromim放置在toim上）
    if is_color:
        # 所有非黑色像素
        alpha = ((fromim_t[:, :, 0] * fromim_t[:, :, 1] * fromim_t[:, :, 2]) > 0)
        for col in range(3):
            toim_t[:, :, col] = fromim_t[:, :, col] * alpha + toim_t[:, :, col] * (1 - alpha)
    else:
        alpha = (fromim_t > 0)
        toim_t = fromim_t * alpha + toim_t * (1 - alpha)
    return toim_t

实验用图：





结果：





坏的结果：

原因不明，为什么会把最左的图片拼接到最右？

时间紧尚未有答案，之后还会持续更新