图像变换：图像扭曲与图像融合

---

单应性变换

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

如图，两张图片中红色的点叫做对应点对。
单应性变换的目标是通过给定的几个点（通常是4对点）来得到单应性矩阵。
用一个3x3的矩阵来表示单应性，可以写成

考虑第一组对应点 （X1，Y1） 在第一张图像和 （X2，Y2）第二张图像中。然后，Homography H 以下列方式映射它们

坐标变换：

要得到两张图片的H矩阵，就必须至少知道4个相同对应位置的点，通常，这些点对应是通过匹配图像之间的SIFT或Harris等特征自动找到的。

算法如下

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。然后对这个矩阵进行处理和归一化，返回输出。

---

仿射变换

仿射变换也称仿射投影，是指几何中，对一个向量空间进行线性变换并接上一个平移，变换为另一个向量空间。所以，仿射变换其实也就是在讲如何来进行两个向量空间的变换

假设有一个向量空间k：

和另一个向量空间j：

如果我们想要将向量空间由k变为j，可以通过下面的公式来实现：

将上式进行拆分可得

我们再将上式转换为矩阵的乘法

上式表明通过一个两行三列的矩阵M就可以实现两个向量空间之间的转换，即仿射变换具有6个自由度，因此需要三个对应点对来估计矩阵M

算法如下：

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著的Multiplr 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包来简单完成。

源码如下：

from numpy import *
from matplotlib.pyplot import *
from scipy import ndimage
from PIL import Image
im = array(Image.open('raccoon.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]))
gray()
subplot(121)
imshow(im)
axis('off')
subplot(122)
imshow(im2)
axis('off')
show()

实验结果：

如图左边为原始图像，右边为扭曲后生成的图像 ，根据代码，设置了H矩阵值为[[1.4, 0.05, -100], [0.05, 1.5, -100], [0, 0, 1]]，则s cosθ=1.4，-s sinθ=0.05，tx=-100，s sinθ=0.05，scosθ=1.5，ty=-100。从上图可以看到扭曲区域边界之外的地方用0像素来填充，来创建一个二值的alpha图像。

仿射变换之后还是平行四边形

图像融合

通过仿射扭曲原理可以实现图像的融合，即将图像或者图像的一部分放在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。

源码如下：

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
im1 = array(Image.open('raccoon.jpg').convert('L'))
im2 = array(Image.open('dog.jpg').convert('L'))

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

# 选定一些目标点
tp = array([[200, 300, 500, 200], [100, 200, 400, 500], [1, 1, 1, 1]])

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

运行结果

如图，运行结果表明已将左一图像按照指定的位置融入左二的图像中，融合结果如右一所示

通过修改tp = array([[200, 300, 500, 200], [100, 200, 400, 500], [1, 1, 1, 1]])这一行的坐标点可以改变左一映射在左二的位置，其中前四个数字代表四个角点的纵坐标，中间四个数字代表四个角点的横坐标，四个角点的顺序为从左上角开始按照逆时针方向排序，最后四个数字代表四个角点的α通道。