Python计算机视觉编程第三章

最新推荐文章于 2024-09-15 08:52:10 发布

asdaasddsa

最新推荐文章于 2024-09-15 08:52:10 发布

阅读量704

点赞数 28

文章标签：计算机视觉开发语言

本文链接：https://blog.csdn.net/zxsdss/article/details/140879451

版权

一、单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。平面指的是图像或者三维中的平面表面。单应性矩阵H依赖于尺度定义。其映射点的方程为：

$\begin{bmatrix}x'\\y'\\w'\end{bmatrix}=\begin{bmatrix}h_1&h_2&h_3\\h_4&h_5&h_6\\h_7&h_8&h_9\end{bmatrix}\begin{bmatrix}x\\y\\w\end{bmatrix}$ 或者 $\mathbf{x}^{\prime}=H\mathbf{x}$

其中实现对点进行归一化和转换齐次坐标的代码为：

import numpy as np
def normalize(points):
    for row in points:
        row /= points[-1]
    return points

def make_homog(points):
    return np.vstack([points, np.ones(points.shape[0])])

仿射变换的方程为：

$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}a_1&a_2&t_x\\a_3&a_4&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}$ 或者 $x^{'}=\begin{bmatrix}A&t\\0&1\end{bmatrix}x$

相似变换的方程为：

$\begin{bmatrix}x'\\y'\\1\end{bmatrix}=\begin{bmatrix}s\cos(\theta)&-s\sin(\theta)&t_x\\s\sin(\theta)&s\cos(\theta)&t_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}$ 或者 $x'=\begin{bmatrix}sR&t\\0&1\end{bmatrix}x$

其为一个包含尺度变化的二维刚体变换。

1.直接线性变换算法

单应性矩阵可以由两幅图像(或平面)中对应点对计算出来。计算单应性矩阵日需要4个对应点对。DLT（直接线性变换）是给定4个或者更多对应点对应矩阵，来计算单应性矩阵H的算法。将单应性矩阵 H 作用在对应点对上，可以得到Ah=0的方程或：

$\begin{bmatrix}-x_1&-y_1&-1&0&0&0&x_1x_1^{\prime}&y_1x_1^{\prime}&x_1^{\prime}\\0&0&0&-x_1&-y_1&-1&x_1y_1^{\prime}&y_1y_1^{\prime}&y_1^{\prime}\\-x_2&-y_2&-1&0&0&0&x_2x_2^{\prime}&y_2x_2^{\prime}&x_2^{\prime}\\0&0&0&-x_2&-y_2&-1&x_2y_2^{\prime}&y_2y_2^{\prime}&y_2^{\prime}\\&\vdots&&\vdots&\vdots&\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\\h_4\\h_5\\h_6\\h_7\\h_8\\h_9\end{bmatrix}=\mathbf{0}$

可以使用SVD（奇异值分解）算法找到H的最小二乘解。其实现的代码为：

def H_from_points(fp, tp):
    if fp.shape != tp.shape:
        raise RuntimeError("number of points don't match")
    m = np.mean(fp[:2],axis=1)
    maxstd = max(np.std(fp[:2],axis=1))+1e-9
    C1 = np.diag([1/maxstd,1/maxstd,1])
    C1[0][2] = -m[0]/maxstd
    C1[1][2] = -m[1]/maxstd
    fp = np.dot(C1,fp)

    m = np.mean(tp[:2],axis=1)
    maxstd = max(np.std(tp[:2],axis=1))+1e-9
    C2 = np.diag([1/maxstd,1/maxstd,1])
    C2[0][2] = -m[0]/maxstd
    C2[1][2] = -m[1]/maxstd
    tp = np.dot(C2,tp)

    nbr_correspondences = fp.shape[1]
    A = np.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 = np.linalg.svd(A)
    H = V[8].reshape((3,3))
    H = np.dot(np.linalg.inv(C2),np.dot(H,C1))
    return H / H[2,2]

2.仿射变换

由于仿射变换具有6个自由度，因此需要三个对应点对来估计矩阵。通过将最后两个元素设置为0，即 $h_7$ = $h_8$ =0，仿射变换可以用上面的DLT算法估计得出。其实现代码为：

def Haffine_from_points(fp, tp):
    if fp.shape != tp.shape:
        raise RuntimeError('number of points do not match')
    m = np.mean(fp[:2], axis=1)
    maxstd = max(np.std(fp[:2], axis=1)) + 1e-9
    C1 = np.diag([1 / maxstd, 1 / maxstd, 1])
    C1[0][2] = -m[0] / maxstd
    C1[1][2] = -m[1] / maxstd
    fp_cond = np.dot(C1, fp)
 
    m = np.mean(tp[:2], axis=1)
    C2 = C1.copy()
    C2[0][2] = -m[0] / maxstd
    C2[1][2] = -m[1] / maxstd
    tp_cond = np.dot(C2, tp)

    A = np.concatenate((fp_cond[:2],tp_cond[:2]), axis=0)
    U,S,V = np.linalg.svd(A.T)
    tmp = V[:2].T
    B = tmp[:2]
    C = tmp[2:4]
 
    tmp2 = np.concatenate((np.dot(C,np.linalg.pinv(B)), np.zeros((2,1))), axis=1)
    H = np.vstack((tmp2,[0,0,1]))
    H = np.dot(np.linalg.inv(C2),np.dot(H,C1))
    return H / H[2,2]

类似于 DLT算法，这些点需要经过预处理和去处理化操作。

二、图像扭曲

对图像块应用仿射变换，我们将其称为图像扭曲。扭曲操作可以通过SciPy工具包的ndimage包来实现。使用其中的如下的命令：

transformed_im = ndimage.affine_transform(im,A,b,size)

from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
from scipy import ndimage

img = np.array(Image.open(r'D:\test\pil.png').convert('L'))
H = np.array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
img2 = ndimage.affine_transform(img,H[:2,:2],(H[0,2],H[1,2]))

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(121), plt.imshow(img, plt.cm.gray), plt.title('原始图像'), plt.axis('off')
plt.subplot(122), plt.imshow(img2, plt.cm.gray), plt.title('扭曲图像'), plt.axis('off')
plt.show()

实验结果为：

输出图像结果中丢失的像素用零来填充。

1.图像中的图像

仿射扭曲的一个简单例子是将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。

def image_in_image(img,img2,tp):
    m,n = img.shape[:2]
    fp = np.array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
    H = Haffine_from_points(fp, tp)
    img_t = ndimage.affine_transform(img,H[:2,:2],(H[0,2],H[1,2]),img2.shape[:2])
    alpha = (img_t > 0)
    return (1 - alpha) * img2 + alpha * img_t

该代码的功能是将扭曲的图像和第二幅图像融合，创建出 alpha 图像。总体实现代码为：

from PIL import Image
from matplotlib import pyplot as plt
import numpy as np
from scipy import ndimage

img2 = np.array(Image.open(r'D:\test\pil.png').convert('L'))
img = np.array(Image.open(r'D:\test\kanni.png').convert('L'))
tp =  np.array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
img3 = image_in_image(img,img2,tp)


plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(131), plt.imshow(img, plt.cm.gray), plt.title('原始图像1'), plt.axis('off')
plt.subplot(132), plt.imshow(img2, plt.cm.gray), plt.title('原始图像2'), plt.axis('off')
plt.subplot(133), plt.imshow(img3, plt.cm.gray), plt.title('生成图像'), plt.axis('off')
plt.show()

其结果为：

其中使用了两个前面所提到的函数来实现这一功能。其次，通过修改代码中tp的值可以修改标记图片的位置，坐标tp是用齐次坐标意义下的坐标表示的。

对于三个点，仿射变换可以将一幅图像进行扭曲，使这三对对应点对可以完美地匹配上。因为仿射变换具有6个自由度，三个对应点对可以给出6个约束条件。所以，如果打算使用仿射变换将图像放置到另一张图像上，可以将图像分成两个三角形，然后对它们分别进行扭曲图像操作。其实现代码为：

def alpha_for_triangle(points,m,n):
    alpha = np.zeros((m,n))
    for i in range(min(points[0]),max(points[0])): 
        for j in range(min(points[1]),max(points[1])):
            x = np.linalg.solve(points,[i,j,1])
            if min(x) > 0: # 所有系数都大于零
                alpha[i,j] = 1
    return alpha

img2 = np.array(Image.open(r'D:\test\pil.png').convert('L'))
img = np.array(Image.open(r'D:\test\kanni.png').convert('L'))
tp =  np.array([[264,538,540,264],[40,36,605,605],[1,1,1,1]])
# img3 = image_in_image(img,img2,tp)

m,n = img.shape[:2]
fp = np.array([[0,m,m,0],[0,0,n,n],[1,1,1,1]])
# 第一个三角形
tp2 = tp[:,:3]
fp2 = fp[:,:3]
# 计算 H
H = Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(img,H[:2,:2],
 (H[0,2],H[1,2]),img2.shape[:2]) 
# 三角形的 alpha
alpha = alpha_for_triangle(tp2,img2.shape[0],img2.shape[1])
img3 = (1-alpha)*img2 + alpha*im1_t
# 第二个三角形
tp2 = tp[:,[0,2,3]]
fp2 = fp[:,[0,2,3]]
# 计算 H
H = Haffine_from_points(tp2,fp2)
im1_t = ndimage.affine_transform(img,H[:2,:2],
 (H[0,2],H[1,2]),img2.shape[:2])
# 三角形的 alpha 图像
alpha = alpha_for_triangle(tp2,img2.shape[0],img2.shape[1])
img4 = (1-alpha)*img3 + alpha*im1_t

结果为：

2.分段仿射扭曲

三角形图像块的仿射扭曲可以完成角点的精确匹配。对应点对集合之间最常用的扭曲方式是分段仿射扭曲。可以使用使用 Matplotlib 和 SciPy 来完成这个操作。为了三角化这些点，经常使用狄洛克三角剖分方法。其一个简单的例子为：

import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage
from scipy.spatial import Delaunay

x,y = np.array(np.random.standard_normal((2,100)))
tri = Delaunay(np.c_[x, y]).simplices
plt.figure() 
for t in tri:
    t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
    plt.plot(x[t_ext],y[t_ext],'r')
plt.plot(x,y,'*')
plt.axis('off')
plt.show()

结果为：

其是一个进行三角刨分的结果，狄洛克三角剖分选择一些三角形，使三角剖分中所有三角形的最小角度最大。实际在图像中使用的例子为：

def triangulate_points(x,y):
    tri = Delaunay(np.c_[x, y]).simplices
    return tri

def pw_affine(fromim,toim,fp,tp,tri):
    im = toim.copy()
    is_color = len(fromim.shape) == 3
    im_t = np.zeros(im.shape, 'uint8')
    for t in tri:
        # import pdb;pdb.set_trace()
        H = 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])
        else:
            im_t = ndimage.affine_transform(fromim,H[:2,:2],(H[0,2],H[1,2]),im.shape[:2])
        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]] # 将第一个点加入到最后
        plt.plot(x[t_ext],y[t_ext],'r')

fromim = np.array(Image.open(r'D:\test\pil.png'))
x,y = np.meshgrid(range(5),range(6))
x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()
tri = triangulate_points(x,y)
im = np.array(Image.open(r'D:\test\bianhuan.png'))
tp = np.loadtxt(r'D:\test\bianhuan.txt') # destination points
fp = np.vstack((y,x,np.ones((1,len(x)))))
tp = np.vstack((tp[:,1],tp[:,0],np.ones((1,len(tp)))))
im1 = pw_affine(fromim,im,fp,tp,tri)


plt.rcParams['font.sans-serif'] = ['SimHei']
plt.subplot(141), plt.imshow(im, plt.cm.gray), plt.title('原始图像1'), plt.axis('off')
plt.subplot(142), plt.imshow(fromim, plt.cm.gray), plt.title('原始图像2'), plt.axis('off')
plt.subplot(143), plt.imshow(im1, plt.cm.gray), plt.title('生成图像'), plt.axis('off')
plt.subplot(144), plt.imshow(im1, plt.cm.gray),plot_mesh(tp[1],tp[0],tri), plt.title('生成图像'), plt.axis('off')
plt.show()

结果为：

要进行该实验首先需要生成一个对应的目标点的文件，这里使用的是ginput()函数手动选取的30个目标点并保存到txt文本文件中去。代表着使用10个三角形进行变换，对于三角形的选取需要不断地进行优化才能得到书本上的那种效果，这里只是进行一个简单的实验。

3.图像配准

图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。其是进行图像对比和更精细的图像分析重要的一步。其一般的步骤为特征提取、特征匹配、变换模型估计、图像变换、优化和评估。

三、创建全景图

1.RANSAC

RANSAC（随机一致性采样）是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，其思想是合理的模型应该能够在描述正确数据点的同时摒弃噪声点。其中的一个简单的实现的例子为：

import numpy
import scipy # use numpy if scipy unavailable
import scipy.linalg # use numpy if scipy unavailable
def ransac(data,model,n,k,t,d,debug=False,return_all=False):
    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:
    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()

结果为：

其定义了一个基于最小二乘法的方法，适用于二维回归问题。橙色线表示迭代方法找到的最小二乘法参数，该方法成功地忽略了异常值点。

2.稳健的单应性矩阵估计

在任何模型中都可以使用 RANSAC 模块。在使用 RANSAC 模块时，只需要在相应 Python 类中实现 fit() 和 get_error() 方法，剩下就是正确地使用 ransac.py，即上述的示例代码。

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

imname = [r'D:\test\111\0' + str(i + 1) + '.jpg' for i in range(4)]
featname = [r'D:\test\111\0' + str(i + 1) + '.sift' for i in range(4)]
# import pdb;
# pdb.set_trace()
im = [array(Image.open(imname[i]).convert('L')) for i in range(4)]
# 提取特征并匹配使用sift算法
l = {}
d = {}
for i in range(4):
    process_image(imname[i],featname[i])
    l[i], d[i] = read_sift(featname[i])
 
matches = {}
for i in range(3):
    matches[i] = sift.match(d[i+1],d[i])
figure()
gray()
for i in range(3):
    plot_matches(im[i+1], im[i], l[i+1], l[i], matches[i], show_below=True)
    figure()
show()

结果为：

上述的实验为使用SIFT特征自动寻找匹配对应，使用可能的对应点集来自动找到用于全景图像的单应性矩阵。但并不是所有图像中的对应点对都是正确的。SIFT是由很强稳健性的描述子，可以比其他描述子产生更少的错误的匹配。

3.拼接图像

估计出图像间的单应性矩阵后，需要将所有的图像扭曲到一个公共的图像平面上。其中一种方法是创建一个很大的图像，然后将所有图像扭曲到上面。实现代码为：

def panorama(H,fromim,toim,padding=2400,delta=2400):
    is_color = len(fromim.shape) == 3
 
    def transf(p):
        p2 = np.dot(H,[p[0],p[1],1])
        return (p2[0]/p2[2],p2[1]/p2[2])

    if H[1,2]<0:
        print ('warp - right')
        if is_color:
            toim_t = np.hstack((toim,np.zeros((toim.shape[0],padding,3))))
            fromim_t = np.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:
            toim_t = np.hstack((toim,np.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 = np.array([[1,0,0],[0,1,-delta],[0,0,1]])
        H = np.dot(H,H_delta)
        if is_color:
            toim_t = np.hstack((np.zeros((toim.shape[0],padding,3)),toim))
            fromim_t = np.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:
            toim_t = np.hstack((np.zeros((toim.shape[0],padding)),toim))
            fromim_t = ndimage.geometric_transform(fromim,transf,(toim.shape[0],toim.shape[1]+padding))
        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

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 = np.vstack([fp[1], fp[0], fp[2]])
    tp = np.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

delta = 2000
im1 = np.array(Image.open(imname[1]))
im2 = np.array(Image.open(imname[2]))
im_12 = panorama(H_12,im1,im2,delta,delta)
im1 = np.array(Image.open(imname[0]))
im_02 = panorama(np.dot(H_12,H_01),im1,im_12,delta,delta)
im1 = np.array(Image.open(imname[3]))
im_32 = panorama(H_32,im1,im_02,delta,delta)
im1 = np.array(Image.open(imname[4]))
im_42 = panorama(np.dot(H_32,H_43),im1,im_32,delta,2*delta)

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

运行出现如下的错误：