计算机视觉：第三章图像到图像映射——单应性变化换（直接线性变换，仿射变换），图像扭曲，图像中的图像

备注：用于测试的图片全部来自福建省厦门集美大学，觉得集美大学长得不错的可以点个赞。

1.单应性变化换

1.1 原理

单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表示。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，我们将频繁的使用单应性变换。本质上，单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：

或者

对于图像平面内（甚至是三维中的点，后面我们会介绍到）的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。我们通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标是的我们可以简单地使用一个矩阵来表示变换。

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

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

如下图所示的平面的两幅图像。红点表示两幅图像中的相同物理点，我们称之为对应点。这里显示了四种不同颜色的四个对应点 - 红色，绿色，黄色和橙色。 一个Homography是一个变换（3×3矩阵），将一个图像中的点映射到另一个图像中的对应点。单应性变换其实就是一个平面到另一个平面的变换关系。

1.2单应性矩阵的理解及直接线性变换算法

在齐次坐标中，假设一点p(xi,yi,1)经过H矩阵的变换变为p‘(xi’,yi’,1)，即 p’ = Hp，通常，对于直接线性变换，一个完全射影变换具有8个自由度，H矩阵有8个自由度，这样至少需要4对特征点对求解。4个特征点对可以建立8个方程。那么对于有n对特征点的情况(超定方程)，解p’ = Hp方程组可以转化为对齐次方程组Ax = 0 的求解。而对 Ax = 0 的求解转化为 min ||Ax||2 的非线性优化问题（超定方程，通过最小二乘拟合得到近似解）。

进一步可得：

这样便可构造系数矩阵：

通过系数矩阵我们可以构造出齐次线性方程组(Ax = 0)：

即：

对于这样的超定方程求解，可以通过最小二乘的方式求解。通过对系数矩阵A求取特征值和特征向量得到。通过以下方式获得最小二乘解：

[V,D] = eig(A’*A) ，其中D是特征值对角矩阵(特征值沿主对角线降序)，V是对应D特征值的特征向量(列向量)组成的特征矩阵,A’表示A的转置。其最小二乘解为V(1)，即系数矩阵A最小特征值对应的特征向量就是超定方程组Ax = 0的最小二乘解，至此，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)) +