图像到图像的映射

单应性变换

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

对于图像平面内的点， 齐次坐标 是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，
x=[x,y,w]=[αx,αy,αw]=[x/w,y/w,1]
都表示同一个二维点。因此，单应性矩阵 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]))))

直接线性变换算法

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

或者 Ah=0，其中 A 是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用 SVD（Singular Value Decomposition，
奇异值分解）算法找到 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[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]

仿射变换

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

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 Edition 所示，
 # 创建矩阵 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]

图像扭曲

对图像块应用仿射变换，称为图像扭曲（或者仿射扭曲）。扭曲操作可以使用
SciPy 工具包中的 ndimage 包来简单完成。命令：

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

使用如上所示的一个线性变换 A 和一个平移向量 b 来对图像块应用仿射变换。选项
参数 size 可以用来指定输出图像的大小。默认输出图像设置为和原始图像同样大
小。例如：

from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image

from scipy import ndimage
im = array(Image.open('cat1.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()
imshow(im)
figure()
imshow(im2)
show()