计算机视觉-图像到图像的映射

---

单应性变换

单应性变换是将一个平面内的点映射到另一个平面内的二维投影映射。单应性变换具有很强的实用性，比如图像配准、图像纠正和纹理扭曲，以及创建全景图像。本质上，单应性变换H，在齐次坐标意义下，按照下面的方程映射二维中的点：
$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}\right]=\left[\begin{array}{ccc}{h}_{00}& h_{01}& h_{02}\\ h_{10}& h_{11}& h_{12}\\ h_{20}& h_{21}& h_{22}\end{array}\right]\left[\begin{array}{c}x\\ y\\ w\end{array}\right]$$
或
$$x^{\prime }=Hx$$
在这里，点的齐次坐标是依赖于其尺度定义的，所以x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点，而对于单映矩阵H，将其乘上任意的一个数其表示的意义都相同，我们通常把$H_{22}$的位置变换为1。因此，不论单应性矩阵怎么变化，始终具有8个自由度。
上述矩阵计算结果可以得出三个矩阵，将w’对应的等式带入x’和y’可以得到：

直接线性变换算法

DLT（直接线性变换）是给定4个或更多对应点对矩阵，来计算单应性矩阵H的算法，其是由上述x’和y’推导出来的：

或者Ah=0。将这些对应点对方程的系数堆叠到一个矩阵中，我们可以使用SVD（奇异值分解）算法找到H的最小二乘解。

---

仿射变换

$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array}\right]=\left[\begin{array}{ccc}{a}_1& a_2& t_x\\ a_3& a_4& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$
相对于原始原始变换，保持了w=1，不具有投影变换的变形能力，其包括一个可逆矩阵A和一个平移变量$t=[t_x,t_y]$，可以运用于图像扭曲。
由于仿射变换具有6个自由度，因此我们需要三个对应点对来估计矩阵H，仿射变换可以用上面的DLT算法估计得出。

图像扭曲

扭曲操作可以使用Scipy工具包中的ndimage包来实现。

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

使用如上所示的一个线性变换A和一个平移变量b来对图像块应用仿射变换，size可以用来指定输出图像的大小，实例演示：

	# 图像扭曲1
	im = array(Image.open('data/q1.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(im2)
    show()

我们也可以精确的把一张图像嵌入到另一幅图像内：

	# 图像扭曲2
    im1 = array(Image.open('data/q3.jpg').convert('L'))
    im2 = array(Image.open('data/q7.jpg').convert('L'))
    # tp = array([[2110,2433,2433,2105],[1038,1030,1258,1262],[1,1,1,1]])  # 1(2110,1038)、 2(2433,1030)、3(2433,1258)、4(2105,1262)
    # tp = array([ [1038, 1030, 1258, 1262],[2110, 2433, 2433, 2105], [1, 1, 1, 1]])  # 图像是倒着拍		  的，所以会呈现x和y相反的情况
    tp = array([[1030, 1258, 1262, 1038], [2433, 2433, 2105, 2110], [1, 1, 1, 1]])
    im3=image_in_image(im1,im2,tp)
    figure()
    gray()
    imshow(im3)
    axis('equal')
    axis('off')
    show()

通过和原图对比我们可以看到图像的中央部分已经映射了另一幅图像，这其中有关于图像横纵坐标系的选取问题，我是用横屏拍摄的照片，那么，相机拍摄出的图像在设备上就会被识别成旋转过的图像，所以坐的选取应当取决于拍摄时旋转后的坐标系对应的坐标，而不是设备上显示的坐标。

关于仿射变换还有一个很有用的技巧。对于三个点，仿射变换可以将一幅图像进行扭曲，使这三对对应点对可以完美地匹配上，这是因为仿射变换具有6个自由度，三个对应点对可以给出6个约束条件。所以，如果想要仿射变换的结果更加精确和称心，可以把图像分为两个三角形，然后对其分别仿射变换。实现：

# 图像扭曲3
    im1 = array(Image.open('data/q3.jpg').convert('L'))
    im2 = array(Image.open('data/q7.jpg').convert('L'))
    m,n=im1.shape[:2]
    fp=array([[0,m,n,0],[0,0,n,n],[1,1,1,1]])
    tp = array([[1030, 1258, 1262, 1038], [2433, 2433, 2105, 2110], [1, 1, 1, 1]])
    tp2=tp[:,:3]
    fp2=fp[:,:3]

    H=Haffine_from_points(tp2,fp2)
    im1_t=ndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])
    alpha=alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
    im3=(1-alpha)*im2+alpha*im1_t

    tp2=tp[:,[0,2,3]]
    fp2=fp[:,[0,2,3]]
    H=Haffine_from_points(tp2,fp2)
    im1_t=ndimage.affine_transform(im1,H[:2,:2],(H[0,2],H[1,2]),im2.shape[:2])

    alpha=alpha_for_triangle(tp2,im2.shape[0],im2.shape[1])
    im4=(1-alpha)*im3+alpha*im1_t

    figure()
    gray()
    subplot(121)
    imshow(im3)
    axis('equal')
    axis('off')
    subplot(122)
    imshow(im4)
    axis('equal')
    axis('off')
    show()

分段仿射扭曲

给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。为了三角化这些点，我们经常使用狄洛克三角剖分方法。样例：

# 狄洛克三角剖分
    x, y = array(random.standard_normal((2, 100)))
    tri=Delaunay(np.c_[x,y]).simplices
    figure()
    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]]
        plot(x[t_ext], y[t_ext], 'r')
    plot(x, y, '*')
    axis('off')
    show()

# 三角剖分映射
    fromim = array(Image.open('data/q3.jpg'))
    x, y = meshgrid(range(5), range(6))
    x = (fromim.shape[1] / 4) * x.flatten()
    y = (fromim.shape[0] / 5) * y.flatten()

    tri = triangulate_points(x, y)
    im = array(Image.open('data/q7.jpg'))
    tp = loadtxt('turningtorso_points.txt')
    fp = vstack((y, x, ones((1, builtins.len(x)))))
    tp = vstack((tp[:, 1], tp[:, 0], ones((1, builtins.len(tp)))))
    im = pw_affine(fromim, im, fp, tp, tri)

    figure()
    imshow(im)
    plot_mesh(tp[1], tp[0], tri)
    axis('off')
    show()

---

RANSAC

RANSAC是随机一致性采样，该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。RANSAC的基本思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。
下面看一个小例子：

图像拼接

    for i in range(1,6):
        sift.process_image('data/q'+str(i)+'.jpg','q'+str(i)+'.sift')

    featname = ['q' + str(i + 1) + '.sift' for i in range(5)]
    imname = ['data\\q' + str(i + 1) + '.jpg' for i in range(5)]

    # extract features and match
    l = {}
    d = {}
    for i in range(5):
        sift.process_image(imname[i], featname[i])
        l[i], d[i] = sift.read_features_from_file(featname[i])

    matches = {}
    for i in range(4):
        matches[i] = sift.match(d[i + 1], d[i])

    # visualize the matches (Figure 3-11 in the book)
    for i in range(4):
        im1 = array(Image.open(imname[i]))
        im2 = array(Image.open(imname[i + 1]))
        figure()
        sift.plot_matches(im2, im1, l[i + 1], l[i], matches[i], show_below=True)
	    def convert_points(j):
        ndx = matches[j].nonzero()[0]
        fp = make_homog(l[j + 1][ndx, :2].T)
        ndx2 = [int(matches[j][i]) for i in ndx]
        tp = make_homog(l[j][ndx2, :2].T)

        # switch x and y - TODO this should move elsewhere
        fp = vstack([fp[1], fp[0], fp[2]])
        tp = vstack([tp[1], tp[0], tp[2]])
        return fp, tp


    # estimate the homographies
    model = RansacModel()

    fp, tp = convert_points(1)
    H_12 = H_from_ransac(fp, tp, model)[0]  # im 1 to 2

    fp, tp = convert_points(0)
    H_01 = H_from_ransac(fp, tp, model)[0]  # im 0 to 1

    tp, fp = convert_points(2)  # NB: reverse order
    H_32 = H_from_ransac(fp, tp, model)[0]  # im 3 to 2

    tp, fp = convert_points(3)  # NB: reverse order
    H_43 = H_from_ransac(fp, tp, model)[0]  # im 4 to 3

    # warp the images
    delta = 2000  # for padding and translation

    im1 = array(Image.open(imname[1]), "uint8")
    im2 = array(Image.open(imname[2]), "uint8")
    im_12 = panorama(H_12, im1, im2, delta, delta)

    im1 = array(Image.open(imname[0]), "f")
    im_02 = panorama(dot(H_12, H_01), im1, im_12, delta, delta)

    im1 = array(Image.open(imname[3]), "f")
    im_32 = panorama(H_32, im1, im_02, delta, delta)

    im1 = array(Image.open(imname[4]), "f")
    im_42 = panorama(dot(H_32, H_43), im1, im_32, delta, 2 * delta)

    figure()
    imshow(array(im_42, "uint8"))
    axis('off')
    show()

总结

在图像拼接的代码运行时，多张图片会出现视角变换后映射出错的情况，两两拼接比较顺利，这大概率是因为在计算单应性矩阵时对应点有噪声，目前还没有解决这个问题，留待以后。