Python计算机视觉编程 第三章 图像到图像的映射

    本章讲解图像之间的变换，以及一些计算变换的实用方法。这些变换可用于图像扭曲变形和图像配准。

3.1 单应性变换

    单应性变换是一个人平面内的点映射到另一个平面内的二维投影变换。本质上，单应性变换$H$，按照下面的方程映射二维中的点（齐次坐标意义下）：
$$\left[\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ z^{{}^{\prime }}\end{array}\right]=\left[\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]或X^{{}^{\prime }}=HX$$
创建homography.py文件，并添加下列代码：

from numpy import *

def normalize(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]))))

    进行点和变换的处理时，我们会按照列优先的原则存储这些点，因此，n个二维点集将会存储为齐次坐标意义下的一个3*n数组。
    在投影变换中，有一些特别重要的变换。比如：
仿射变换：
$$\left[\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{a}_1& a_2& t_x\\ a_3& a_4& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]或X^{{}^{\prime }}=\left[\begin{array}{cc}AA& tt\\ 00& 1\end{array}\right]X$$
    保持了$w=1$，不具有投影变换所具有的强大变形能力。仿射变换包含了一个可逆矩阵$AA$和一个平移向量$tt=[t_x,t_y]$。仿射变换可用于图像扭曲。
相似变换：
$$\left[\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}\right]=\left[\begin{array}{ccc}scos(\theta )& -ssin(\theta )& t_x\\ sin(\theta )& scos(\theta )& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]或X^{{}^{\prime }}=\left[\begin{array}{cc}sAA& tt\\ 00& 1\end{array}\right]X$$
    是一个包含尺度变化的二维刚体变换。其中向量$s$指定了变换的尺度，$R$是角度为$\theta$的旋转矩阵，$tt=[t_x,t_y]$在这里也是一个平移向量。相似变换可以用于图像配准。
刚体变换
    如果相似变换中$s=1$，那么该变换能够保持距离不变，此时变换为刚体变换。
    旋转和平移变换，点与点之间的距离不变。

3.1.1直接线性变换算法

    DLT（Direct Linear Transformation，直接线性变换）是给定4个或者更多对应点对矩阵，来计算单应性矩阵$HH$的算法。将单应性矩阵$HH$作用在对应点对上，重新写出该方程，可以得到：
$$\left[\begin{array}{ccccccccc}-x_1& -y_1& -1& 0& 0& 0& x_1{x}_1^{{}^{\prime }}& y_1{x}_1^{{}^{\prime }}& x_1^{{}^{\prime }}\\ 0& 0& 0& -x_1& -y_1& -1& x_1{y}_1^{{}^{\prime }}& y_1{y}_1^{{}^{\prime }}& y_1^{{}^{\prime }}\\ -x_2& -y_2& -1& 0& 0& 0& x_2{x}_2^{{}^{\prime }}& y_2{x}_2^{{}^{\prime }}& x_2^{{}^{\prime }}\\ 0& 0& 0& -x_2& -y_2& -1& x_2{y}_2^{{}^{\prime }}& y_2{y}_2^{{}^{\prime }}& y_2^{{}^{\prime }}\\ & ...& & ...& & ...& & ...& \end{array}\right]\left[\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\\ h_9\end{array}\right]=00或者AAhh=00$$
    其中$AA$是一个具有对应点对二倍数量行数的矩阵。
DLT算法推导过程：
    在图像拼接中，得到了两张图像的特征匹配，两个点集分别记作$X$和$X_{{}^{\prime }}$。用单应性变换来拟合二者的关系，表达式为：
$$c\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=HH\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$
    其中，$(u,v,1{)}^T$是$X$中特征点的坐标，$(x,y,1{)}^T$是$X^{{}^{\prime }}$中特征点的坐标，$HH$是单应性矩阵，代表它们之间的变换关系。$HH$是个3×3的矩阵，有8个自由度，所以待求未知参数有8个。
$$HH=\left[\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}\right]$$
    通过矩阵乘法将上述表达式展开，并且消去$c$，得到：
$$-h_{1}x+-h_{2}y-h_3+(h_{7}x+h_{8}y+h_9)u=0$$
$$-h_{4}x+-h_{5}y-h_6+(h_{7}x+h_{8}y+h_9)v=0$$
    整理为：$A_{i}h=0$
    其中，$A=\left[\begin{array}{ccccccccc}-x& -y& -1& 0& 0& 0& ux& uy& u\\ 0& 0& 0& -x& -y& -1& vx& vy& v\end{array}\right]$
          $h=\left[\begin{array}{ccccccccc}{h}_1& h_2& h_3& h_4& h_5& h_6& h_7& h_8& h_9\end{array}\right]$
    由未知变量的个数可知，求解出H HH至少需要4对匹配点。通常情况下为了得到更稳定的结果，会用到多于4对的特征匹配。所以，这个方程会变成超定的，可以将最小二乘解作为最后的解。
    将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用SVD（Singular Value Decomposition，奇异值分解）算法找到$HH$的最小二乘法。
在homography.py文件中添加下列代码：

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]

    对这些点进行归一化操作，其均值为0，方差为1，因为算法的稳定性取决于坐标的表示情况和部分数值计算的问题，所以归一化操作非常重要。接下来使用对应点来构造矩阵A。最小二乘解即为矩阵SVD分解后所得矩阵V的最后一行。该行经过变形后得到矩阵H。然后对这个矩阵进行处理和归一化，返回输出。

3.1.2仿射变换

    一种二维坐标到二维坐标之间的线性变换（相同平面），它保持了二维图形的“平直性”（直线经过变换之后依然是直线）和“平行性”（二维图形之间的相对位置关系保持不变，平行线依然是平行线，且直线上点的位置顺序不变），但是角度会改变。任意的仿射变换都能表示为乘以一个矩阵(线性变换)，再加上一个向量 (平移) 的形式。单应性变换有8个自由度，仿射有6个自由度， 因此需要三个对应点对来估计矩阵H，通过将最后两个元素设置为0，即h7=h8=0，仿射变换可以用上面的DLT算法估计得到。
    这里我们使用不同的方法来计算单应性矩阵$H$。
在homograph.py文件中添加下列代码：

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著的Multiple View Geometry in computer,scond eidtion所示，
    #创建矩阵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] 

3.2图像扭曲

    对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。扭曲操作可以使用SciPy工具包中的ndimage包来简单完成。
运行代码：

from scipy import ndimage
from pylab import *
from PIL import Image

im = array(Image.open('JMU1.jpg').convert('L'))
H = array([[1.4,0.05,-100],[0.05,1.5,-100],[0,0,1]])
im2 = ndimage.affine_transform(im,H[:2,:2],(H[0,2],H[1,2]))

figure()
gray()

subplot(121)
imshow(im)
subplot(122)
imshow(im2)
show()

结果：

分析：
    命令ndimage.affine_transform(im,A,b,size)，上面程序使用线性变换H[:2,:2]和两个平移向量H[0,2],H[1,2]来对图像块应用仿射变换。参数size可以用来指定输出图像的大小，这里默认输出图像设置为和原始图像同样大小。可以看到，输出图像结果中丢失的像素用零来填充。

3.2.1图像中的图像

    仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对其。
在homograph.py中添加下列代码：

from scipy import ndimage
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 = 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

测试代码：

from pylab import *
from PIL import Image
import homograph

#仿射扭曲im1到im2的例子
im1 = array(Image.open('girlfriend.jpg').convert('L'))
im2 = array(Image.open('pavilion.jpg').convert('L'))

gray()
subplot(131)
imshow(im1)
axis('equal')
axis('off')
subplot(132)
imshow(im2)
axis('equal')
axis('off')

# 选定一些目标点
tp = array([[854, 1118, 1162, 845],[833, 838, 1212, 1207], [1, 1, 1, 1]])

im3 = homograph.image_in_image(im1, im2, tp)
subplot(133)
imshow(im3)
axis('equal')
axis('off')
show()

结果：

分析：
    标记物的坐标tp是用其次坐标意义下的坐标表示的。改变tp = array([[900, 1300, 1300, 900],[400, 400, 800, 800], [1, 1, 1, 1]]) 中的坐标位置，可以将图像放置在不同的位置。

3.2.2图像配准

    图像配准是对图像进行变换，是变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准，也可以是非严格配准。其目的在于比较或融合针对同一对象在不同条件下获取的图像，例如图像会来自不同的采集设备，取自不同的时间，不同的拍摄视角等等，有时也需要用到针对不同对象的图像配准问题。具体地说，对于一组图像数据集中的两幅图像，通过寻找一种空间变换把一幅图像映射到另一幅图像，使得两图中对应于空间同一位置的点一一对应起来，从而达到信息融合的目的。
方法分类：
    1、基于灰度和模板的： 这类方法直接采用相关运算等方式计算相关值来寻求最佳匹配位置，模板匹配（Blocking Matching）是根据已知模板图像到另一幅图像中寻找与模板图像相似的子图像。基于灰度的匹配算法也称作相关匹配算法，用空间二维滑动模板进行匹配，不同匹配算法主要体现在相关准则的选择方面。
    常用的算法： 平均绝对差算法（MAD）、绝对误差和算法（SAD）、误差平方和算法（SSD）、平均误差平方和算法（MSD）、归一化积相关算法（NCC）、序贯相似性检测算法（SSDA）、hadamard变换算法（SATD）、局部灰度值编码算法、PIU。
    2、基于特征的匹配方法： 首先提取图像的特征，再生成特征描述子，最后根据描述子的相似程度对两幅图像的特征之间进行匹配。图像的特征主要可以分为点、线（边缘）、区域（面）等特征，也可以分为局部特征和全局特征。区域（面）特征提取比较麻烦，耗时，因此主要用点特征和边缘特征。
    点特征匹配包括： Harris、Moravec、KLT、Harr-like 、HOG 、LBP、SIFT、SURF、BRIEF、SUSAN、FAST 、CENSUS、FREAK、BRISK、ORB、光流法、A-KAZE等。
    边缘特征包括： LoG算子、Robert算子、Sobel算子、Prewitt算子、Canny算子等。
    3、基于域变换的方法： 采用相位相关（傅里叶-梅林变换）、沃尔什变换、小波等方法。

应用： 目标检测、模型重建、运动估计、特征匹配、肿瘤检测、病变定位、血管造影、地质勘探、航空侦察

配准算法的一般步骤：
    1、特征提取：特征提取是指分别提取两幅图像中共有的图像特征，这种特征是出现在两幅图像中对比列、旋转、平移等变换保持一致性的特征，如线交叉点、物体边缘角点、虚圆闭区域的中心等可提取的特征。特征包括：点、线和面三类。点特征是最常用的一种图像特征，包括物体边缘点、角点、线交叉点等；根据各特征点的兴趣值将特征点分成几个等级。对不同的目的，特征点的提取应各有不同。
    2、特征匹配：特征匹配分两步。a. 对特征作描述。现有的主要特征描述子：SIFT特征描述子，SURF特征描述子，对比度直方图（CCH）,DAISY特征描述子，矩方法。b. 利用相似度准则进行特征匹配。常用的相似性测度准则有如欧式距离，马氏距离，Hausdorff距离等。
    3、估计变换模型：空间变换模型是所有配准技术中需要考虑的一个重要因素，各种配准技术都要建立自己的变换模型，变换空间的选取与图像的变形特性有关。常用的空间变换模型有：刚体变换、仿射变换、投影变换、非线性变换。
    4、图像重采样及变换：在得到两幅图像的变换参数后，要将输入图像做相应参数的变换，使之与参考图像处于同一坐标系下，则矫正后的输入图像与参考图像可用作后续的图像融合、目标变化检测处理或图像镶嵌；涉及输入图像变换后所得点坐标不一定为整像素数，则应进行插值处理。常用的插值算法有最近领域法，双线性插值法和立方卷积插值法。

配准算法的主要方法
1、 基于点：
    a. SIFT算法
    SIFT特征匹配算法包括两个阶段：SIFT特征的生成与SIFT特征向量的匹配。
    SIFT特征向量的生成算法：
    （1）尺度空间极值检测，以初步确定关键点位置和所在尺度。
    （2）拟和三维二次函数精确确定位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点。
    （3）利用关键点领域像素的梯度方向分布特性为每个关键点指定参数方向，使算子具备旋转不变性。
    （4）生成SIFT特征向量。
    SIFT特征向量的匹配：
    对图像1中的某个关键点，找到其与图像2中欧式距离最近的前两个关键点的距离NN和SCN，如果NN/SCN小于某个比例阈值，则接受这一对匹配点。
    b.ASIFT算法
    通过原始图像来模拟得到场景在各个视角下的图像，再对这些得到的图像提起SIFT特征点，然后进行匹配。其放射性要好于SIFT，具有全仿射不变性。
    c.SUFR算法
    特征点提取：
    （1）计算原图像的积分图像
    （2）用不同尺寸的框状滤波器来计算不同阶以及不同层上的每个点图像点的行列式。一般计算4阶4层
    （3）在3维(x,y,S）尺度空间中，在每个3x3x3的局部区域里，进行非最大值抑制。只有比邻近的26个点的响应值都大的点才被选为兴趣点。
    特征匹配：
    对图像1中的某个关键点，找到其与图像2中欧式距离最近的前两个关键点的距离NN和SCN，如果NN/SCN小于某个比例阈值，则接受这一对匹配点。
2、基于边缘：
    基本思想：用边缘检测算子提取出边缘
    边缘匹配：
    （1）根据边缘的相似度
    （2）提取边缘上的控制点，如曲率比较大的点等，然后用这些点来进行匹配。
    （3）将边缘拟合成直线，然后匹配直线。
    估计变换参数：用边缘上的控制点或直线的端点，中点等等。

3.3创建全景图

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

3.3.1RANSAC

    RANSAC (随机一致性采样)是一种迭代算法，该算法从一组包含“外点(outlier)”的观测数据中估计数学模型的参数。“外点”指观测数据中的无效数据，通常为噪声或错误数据，比如图像匹配中的误匹配点和曲线拟合中的离群点。与“外点”相对应的是“内点(inlier)”，即用来估计模型参数的有效数据。因此，RANSAC也是一种“外点”检测算法。此外，RANSAC算法是一种非确定算法，它只能在一定概率下产生可信的结果，当迭代次数增加时，准确的概率也会增加。RANSAC算法是用来找到正确模型来拟合带有噪声数据的迭代方法。

算法基本思想和流程：
    RANSAC是通过反复选择数据集去估计出模型，一直迭代到估计出认为比较好的模型。具体的实现步骤可以分为以下几步：
    1、随机采样K个点，K是求解模型参数的最少点个数(对于直线拟合来说就是两个点（两点确定一条直线），对于拟合圆就是3个点)
    2、使用采样出的k个点估计模型参数，得到数据模型；
    3、将所有数据带入这个模型，计算出“内点”的数目；(累加在一定误差范围内的适合当前迭代得出的模型的数据)；
    4、比较当前模型和之前推出的最好的模型的“内点“的数量，记录最大“内点”数的模型参数和“内点”数；
    5、重复1-4步，直到迭代结束或者当前模型已经足够好了(“内点数目大于一定数量”)。

3.3.2稳健的单应性矩阵估计

    在任何模型中都可以使用RANSAC模块，在使用RANSAC模块时，只需要在相应Python类中实现fit()和get_error()方法，正确使用ransac.py即可。这里使用可能的对应点集来自动找到用于全景图像的单应性矩阵。
    总而言之，就是平面的单应性被定义为一个平面到另外一个平面的投影映射。通过ransac算法来求解单应性矩阵，调用homegraphy.py中相应函数：fit() 是计算选取的4个对应的单应性矩阵。get_error() 是对所有的对应计算单应性矩阵，然后对每个变换后的点，返回相应的误差。

3.3.2拼接图像

    估计出图像间的单应性矩阵（使用RANSAC算法），现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面。一种方法是创建一个很大的图像，比如图像中全部填充0，使其和中心图像平行，然后将所有的图像扭曲到上面，由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充0，以便为扭曲的图像腾出空间。

全景图像拼接步骤：
    1、特征点的匹配：
    两张图像要能拼接在一起成为一张图像，就需要这两张图像中存在有重合的部分。通过这些重合的部分使用sift特征点匹配的算法，来寻找到重合部分的特征点。需要注意的是，虽然sift算法比Harris角点的效果更好，但是也会出现错误点，并非完美的匹配方法。
合成全景图的第一步是提取并且匹配所有素材图片的局部特征点。
    什么是特征点？普遍来讲，一张图片所包含的特征点通常就是周围含有较大信息量的点，而仅通过这些富有特征的局部，基本就可以推测出整张图片。比如说物体的棱角、夜景闪耀的星星，或是图片里的图案和花纹。
    2、图像的匹配：
    在找到特征点对之后，因为上文提到sift并非所有匹配点都是正确的，这里我们用到了RANSAC这个方法。这个方法的作用是找到一个合理的模型来描述正确的数据并且尽量忽视噪点的影响。
    接下来就是， 找到所有匹配(也就是重叠)的图片部分，连接所有图片之后就可以形成一个基本的全景图了。因为每张图片有可能和其他每张图片有重叠部分，所以匹配全部图片需要差不多匹配图片个数的平方次。不过实际上每两张图片之间只需要那么几个相对精准匹配的点就可以估算出这两张图像里的几何关系。
    3、全景图矫直：
    矫正拍摄图片时相机的相对3D旋转，主要原因是拍摄图片时相机很可能并不在同一水平线上，并且存在不同程度的倾斜，略过这一步可能导致全景图变成波浪形状。
    4、图像均衡补偿：
    全局平衡所有图片的光照和色调。
    5、图像频段融合：
    步骤4之后仍然会存在图像之间衔接边缘、晕影效果（图像的外围部分的亮度或饱和度比中心区域低）、视差效果（因为相机透镜移动导致）。

    单应性矩阵是要实现图像拼接重要方法。 它表示了两张图象之间的对应特征点的变换关系。有了这个关系，我们就可以实现图像正确的拼接到另一张图像上。
实验代码：

import matplotlib.pyplot as plt
from numpy import *
from PIL import Image
import os
from PCV.geometry import homography, warp
from PCV.localdescriptors import sift


# 批量修改图像尺寸
path = "D:\\CV\\test\\"
for maindir, subdir,file_name_list in os.walk(path):
    print(file_name_list)
    for file_name in file_name_list:
        image=os.path.join(maindir,file_name) #获取每张图片的路径
        file=Image.open(image)
        out=file.resize((1000,900),Image.Resampling.LANCZOS)  #以高质量修改图片尺寸为（1000.900）
        out.save(image)                            #以同名保存到原路径

# 设置数据文件夹的路径
featname = [path + str(i + 1) + '.sift' for i in range(5)]
imname = [path + str(i + 1) + '.jpg' for i in range(5)]

# 提取特征并匹配使用sift算法
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])

# 可视化匹配
for i in range(4):
    im1 = array(Image.open(imname[i]))
    im2 = array(Image.open(imname[i + 1]))
    plt.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]  # 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 的单应性矩阵

# 扭曲图像
delta = 2000  # 用于填充和平移

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)

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

结果：





原图像：

分析：
    可以看到拼接结果还是不错的，主要还是角度偏差比较小，特征点匹配较多。我又尝试过较大角度的拼接，实验结果不是很理想。
    在实验过程中要注意以下两个问题：
    1、图像要能够拼接起来需要在对准同一个物体进行拍摄，从右到左的顺序，使手机保持在同一水平线上拍摄，能够检测出相应的点。
    2、拼接图像大小会影响的程序的速度，所以在程序中增加了一段代码用于修改图像尺寸。