计算机视觉（三）

这次博客的主题是全景拼接，这个是建立在前面学过的特征提取、特征描述子、特征匹配等。

虽然图片是反应三维世界，但实质是二维的，因此图像拼接的可以简化理解为2D图像的变换、叠加过程。

图像之间经过特征匹配后的变换、叠加过程包括了位移、旋转、尺度、仿射、透视。

单应性变换

以上对图像的操作都涉及单应性变换，单应性变换是将一个平面内的点映射到另一个平面内的二维投影变换，对于平面中的点，用齐次坐标可以有效表示，即单应性变换矩阵H：

单应性矩阵仅依赖尺度定义，所以，单应性矩阵具有 8 个独立的自由度。我们通常使用 w=1 来归
一化点，这样，点就具有唯一的图像坐标 x 和 y。以上，我们就可以仅通过一个矩阵来进行图像的变换。

归一化和转换齐次坐标方程：

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

直接线性变换

一个完全射影变换具有 8 个自由度。每个对应点对可以写出两个方程，分别对应于 x 和 y 坐标：比如平移存在两个未知数，仅需要一个点对。因此，计算单应性矩阵 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)
    
    # --to points--
    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))    
    
    # decondition
    H = dot(linalg.inv(C2),dot(H,C1))
    
    # normalize and return
    return H / H[2,2]