计算单应性矩阵 python_计算机视觉编程系列-图像到图像的映射

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

3.1 单应性变换

单应性变换是一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表面。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像。本质上，单应性变换

，按照下面的方程映射二位中的点：

对于图像平面内的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以x = [x, y, w]= [αx, αy, αw]= [x/w, y/w, 1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以单应性矩阵具有8个独立的自由度。我们通常使用w = 1来归一化点，这样点具有唯一的图像坐标x和y。这个额外的坐标使得我们可以简单的使用一个矩阵来表示变换。下面的函数可以用来实现对点归一化和转换齐次坐标的功能。

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

def make_homog(points):
    return vastack((points,ones((1,points,shape[1]))))

进行点和变换处理时，我们会按照列优先的原则存储这些点。因此，n个二维点集将会存储为齐次坐标意义下的一个3 *n数组。这种格式使得矩阵乘法和点的变换操作更加容易。对于其他的例子，比如对于聚类和分类特征，我们将使用典型的行数组来存储数据。在这些投影的变换中，有一些特别重要的变换。比如，仿射变换：

保持了w = 1，不具有投影变换所具有的强大变形能力。仿射变换包含一个可逆矩阵A和一个平移向量t = [tx , ty].仿射变换可以用于很多应用，比如图像扭曲。相似变换：

是一个包含尺度变化的二维刚体变换。上式中的向量s 指定了变换的尺度，R是角度为

的旋转矩阵t = [tx , ty]在这里也是一个平移向量。如果s=1，那么该变换能够保持距离不变。此时，变换为刚体变换。相似变换可以用于很多应用比如图像配准。下面让我们一起探讨如何设计用于估计单应性矩阵的算法，然后看一下使用仿射变换进行图像扭曲，使用相似变换进行图像匹配，以及使用完全投影变换进行创建全景图像的例子。

3.1.1 直接线性变换算法

单应性矩阵可以由两幅图像中对应点计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点束，每个对应点可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要4个对应点对。

DLT（Direct Linear Transformation）是给定4个或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点对上，重新写出该方程，我们可以得到下面的方程：

或者Ah = 0，其中A是一个具有对应点对二倍数行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用SVD奇异值分解算法找到H的最小二乘解。下面是该算法的代码：

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[0][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 仿射变换

由于仿射变换具有6个自由度，因此我们需要三个对应点对来估计矩阵H。通过将最后两个元素设置为0，即h7 = h8 = 0,仿射变换可以用上面的DLT算法估计得出。

这里我们将使用不同方法来计算单应性矩阵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,fp)
    
    #因为归一化后点的均值为0&#