序言

??本章讲解了图像之间的变换，以及一些计算变换的应用。这些应用可以用于图像扭曲变换和图像配准。

运行环境：

1. python3.8
2. PyCharm2020.3

一、单应性变换

??单应性变换 是将一个平面内的点映射到另一个平面内的二维投影变换。在这里，平面是指图像或者三维中的平面表示。单应性变换具有很强的实用性，比如图像配准，图像纠正和纹理扭曲，以及创建全景图像，我们将频繁的使用单应性变换。本质上，单应性变换H，按照下面的方程映射二维中的点（齐次坐标意义下）：$$?\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ w^{{}^{\prime }}\end{array}?=?\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& h_9\end{array}??\begin{array}{c}x\\ y\\ z\end{array}?$$或者$$x^{{}^{\prime }}=Hx$$对于图像平面内（甚至是三维中的点，后面我们会介绍到）的点，齐次坐标是个非常有用的表示方式。点的齐次坐标是依赖于其尺度定义的，所以，x=[x,y,w]=[ax,ay,aw]=[x/w,y/w,1]都表示同一个二维点。因此，单应性矩阵H也仅依赖尺度定义，所以，单应性矩阵具有8个独立的自由度。我们通常使用w=1来归一化点，这样，点具有唯一的图像坐标x和y。这个额外的坐标是的我们可以简单地使用一个矩阵来表示变换。

??我们创建 homography.py 文件。用下面的函数可以实现对点进行归一化和转换齐次坐标的功能：

from numpy import *

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数组。这种格式使得矩阵乘法和点的变换操作更加容易。对于其他的例子，比如对于聚类和分类的特征，我们将使用典型的行数组来存储数据。

??在这些投影变换中，有一些特别重要的变换。比如，仿射变换：$$?\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}?=?\begin{array}{ccc}{a}_1& a_2& t_x\\ a_3& a_4& t_y\\ 0& 0& 1\end{array}??\begin{array}{c}x\\ y\\ 1\end{array}?$$或$$x^{{}^{\prime }}=\left[\begin{array}{cc}A& t\\ 0& 1\end{array}\right]x$$保持了w=1，不具有投影变换所具有的强大变形能力，反射变换包括一个可逆矩阵A和一个平移向量t=[t_x,t_y]。仿射变换可以用于很多应用，比如 图像扭曲。

相似变换：$$?\begin{array}{c}{x}^{{}^{\prime }}\\ y^{{}^{\prime }}\\ 1\end{array}?=?\begin{array}{ccc}s\cos (\theta )& ?s\sin (\theta )& t_x\\ s\sin (\theta )& s\cos (\theta )& t_y\\ 0& 0& 1\end{array}??\begin{array}{c}x\\ y\\ 1\end{array}?$$或$$x^{{}^{\prime }}=\left[\begin{array}{cc}sR& t\\ 0& 1\end{array}\right]x$$是一个包含尺度变化的二维刚体变换。上式中的向量s指定了变换的尺度，R是角度为θ的旋转矩阵，t=[t_x,t_y]在这里也是一个平移向量。如果s=1，那么该变换能够保持距离不变。此时，变换称为刚体变换。相似变换可以用于很多应用，比如图像配准。

??下面让我们来一起探讨如何设计用于估计单应性矩阵的算法，然后看一下使用仿射变换进行图像扭曲，使用相似变换进行图像匹配，以及使用完全投影变换进行创建全景图像的一些例子：

1.1 直接线性变换算法

??单应性矩阵可以有两幅图像（或者平面）中对应点对计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要4个对应点对。

??DLT（Direct Linear Transformation，直接线性变换）是给定4个点或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点上，重新写出该方程，我们可以得到下面的方程：$$?\begin{array}{ccccccccc}?x_1& ?y_1& ?1& 0& 0& 0& x_1{x}_1^{{}^{\prime }}& y_1{x}_1^{{}^{\prime }}& x_1^{{}^{\prime }}\\ 0& 0& 0& ?x_1& ?y_1& ?1& x_1{y}_1{}^{\prime }& y_1{y}_1{}^{\prime }& y_1^{{}^{\prime }}\\ ?x_2& ?y_2& ?1& 0& 0& 0& x_2{x}_2^{{}^{\prime }}& y_2{x}_2^{{}^{\prime }}& x_2^{{}^{\prime }}\\ 0& 0& 0& ?x_2& ?y_2& ?1& x_2{y}_2^{{}^{\prime }}& y_2{y}_2^{{}^{\prime }}& y_2^{{}^{\prime }}\\ & ...& & ...& & ...& & ...& \end{array}??\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\\ h_9\end{array}?=0$$或者Ah=0，其中A是一个具有对应点对二倍数量行数的矩阵。将这些对应点对方程的系数堆叠到一个矩阵红，我们可以使用SVD算法找到H的最小二乘解。下面是算法代码，我们也将其加入到 homography.py 文件里：

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

1.2 仿射变换

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

??下面的函数使用对应点来计算放射变换矩阵，我们继续将其添加到 homography.py 文件中。

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]

??同样的，类似于DLT算法，这些点需要经过预处理和去处理化操作。

二、图像扭曲

??对图像块应用仿射变换，我们将其称为图像扭曲（或者 仿射扭曲）。该操作不仅经常在计算机图形学中，而且经常出现在计算机视觉算法总。扭曲的操作可以使用SciPy工具包中的ndimage包来简单完成。命令：

transform_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

# 解决中文乱码
rcParams['font.sans-serif'] = 'SimHei'
rcParams['axes.unicode_minus'] = False

im = array(Image.open('image/image1.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)
title("原始图像")
axis('off')
subplot(122)
imshow(im2)
title("扭曲后的图像")
axis('off')
show()

运行结果：

??从运行结果可以看出，原始图像（左）和扭曲后的图像（右）的差别，扭曲后的图像中丢失的像素用零来填充。

2.1 图像中的图像

??仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，是的他们能够和指定的区域或者标记物对齐。

??将函数 image_in_image() 添加到 wary.py 文件中。该函数的输入参数为两幅图像和一个坐标。该坐标为将第一幅图像放置到第二幅图像中的角点坐标：

import homography
from scipy import ndimage
from numpy import *

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 = homography.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

??将扭曲的图像和第二幅图像融合，就创建 alpha 图像。该图像定义了每个像素从各个图像中获取的像素值成分多少。这里基于以下事实：扭曲的图像是在扭曲区域边界之外以 0 来填充的图像，来创建一个二值的 alpha 图像。严格意义上，需要在第一幅图象中的潜在 0 像素上加上一个小的数值，或者合理的处理这些 0 像素。

注意：我们使用的图像坐标是齐次坐标意义下的。

使用仿射变换放置图像

??我们试着使用该函数将一幅图像插入另一幅图像。

import warp
from numpy import *
from matplotlib.pyplot import *
from PIL import Image

# 将im1仿射扭曲到im2的指定位置
im1 = array(Image.open('image/1.jpg').convert('L'))
im2 = array(Image.open('image/image3.jpg').convert('L'))

figure()
gray()
subplot(121)
imshow(im1)
axis('equal')
axis('off')
subplot(122)
imshow(im2)
axis('equal')
axis('off')

# 选定一些目标点
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

# 调用的warp.py的image_in_image函数，从而实现仿射变换
im3 = warp.image_in_image(im1, im2, tp)

figure()
imshow(im3)
axis('equal')
axis('off')
show()

运行结果：

??从运行结果上看，图片将原本右下角的一个公告去给替换了。因为 tp 中的齐次坐标是对应公告牌的坐标，所以将其放置到该位置。那么，

使用ginput函数获取 tp 中的齐次坐标值

??其中，tp 中选取的点的坐标值可以通过手工确定，也可以通过 PyLab 类库中的 ginput() 函数获得。可以使用下面的代码实现点齐次坐标的获取：

from PIL import Image
from pylab import *
import numpy as np

im = array(Image.open('image/image3.jpg'))
imshow(im)

print('Please click 4 point')
points = ginput(4)

# 转整数
int_points = np.int_(points)

# 行列互换操作
new_points = list(map(list, zip(*int_points)))

print('you clicked x:', new_points[1])
print('you clicked y:', new_points[0])
show()

运行结果：

选取点的顺序	运行结果

??然后就可以用 x 的坐标替换掉 tp 中的第一个参数列表，用 y 替换掉 tp 中的第而个参数列表：

# tp = array([[264, 538, 540, 264], [40, 36, 605, 605], [1, 1, 1, 1]])
# 替换成
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

注意：运行上述代码选取点时，要按照上图的给的顺序点击，不然可能会出现图片位置错误、图片旋转等一些问题！因为 image_in_image() 函数中的 tp 是按照从左上角逆时针计算的。

??函数 Haffine_from_points() 会返回给定对应点对的最优仿射变换。在上面的例子中，对应点对为图像和公告牌的角点。如果透视效应比较弱，那么这种方法会返回很好的结果。

??对于三个点，仿射变换可以将一幅图像进行扭曲，使这三对对应点对可以完美地匹配上。这是因为，仿射变换具有 6 个自由度，三个对应点对可以给出 6 个约束条件（对于这三个对应点对，x 和 y 坐标必须都要匹配）。所以，如果你真的打算使用仿射变换将图像放置到公告牌上，可以将图像分成两个三角形，然后对它们分别进行扭曲图像操作。代码如下：

from warp import *
from numpy import *
from matplotlib.pyplot import *
from PIL import Image

im1 = array(Image.open('image/1.jpg').convert('L'))
im2 = array(Image.open('image/image3.jpg').convert('L'))

# 选取 im1 角上的一些点
m, n = im1.shape[:2]
fp = array([[0, m, m, 0], [0, 0, n, n], [1, 1, 1, 1]])

# 选定一些目标点
tp = array([[862, 1019, 1019, 866], [905, 911, 1107, 1096], [1, 1, 1, 1]])

# 第一个三角形
tp2 = tp[:, :3]
fp2 = fp[:, :3]

# 计算 H
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 = 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
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 = 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()

运行结果：

??这里我们简单地为每个三角形创建了 alpha 图像，然后将所有的图像合并起来。该三角形的 alpha 图像可以简单地通过检查像素的坐标是否能够写成三角形顶点坐标的凸组合来计算得出。如果坐标可以表示成这种形式，那么该像素就位于三角形的内部。上面的例子使用了下面的函数 alpha_for_triangle()，将其添加到 warp.py 文件中。

def alpha_for_triangle(points, m, n):
    """对于带有由 points 定义角点的三角形，创建大小为 (m, n) 的alpha 图
    （在归一化的齐次坐标意义下）"""

    alpha = zeros((m, n))
    for i in range(min(points[0]), max(points[0])):
        for j in range(min(points[1]), max(points[1])):
            x = linalg.solve(points, [i, j, 1])
            if min(x) > 0:  # 所有系数都大于零
                alpha[i, j] = 1
    return alpha

??这样显卡可以极其快速地操作上面的代码。Python 语言的处理速度比你的显卡（或者 C/C++ 实现）慢很多，但是对于我们来说已经够用了。

2.2 分段仿射扭曲

??正如上面的例子所示，三角形图像块的仿射扭曲可以完成角点的精确匹配。让我们看一下对应点对集合之间最常用的扭曲方式：分段仿射扭曲。给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这些都是最基本的操作。下面我们来演示一下如何使用 Matplotlib 和 SciPy 来完成该操作。

为了三角化这些点，我们经常使用狄洛克三角剖分方法。在 Matplotlib（但是不在 PyLab 库中）中有狄洛克三角剖分，我们可以用下面的方式使用它：

from numpy import *
from matplotlib.pyplot import *
from scipy.spatial import Delaunay

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

figure()
plot(x, y, '*')
axis('off')
show()

运行结果：

??运行结果显示了一些实例点和三角剖分的结果。狄洛克三角剖分选择一些三角形，使三角剖分中所有三角形的最小角度最大（三角剖分中的边实际上是泰森图的对偶图）。函数 delaunay() 有 4 个输出，其中我们仅需要三角形列表信息（第三个输出）。在 warp.py 文件中创建用于三角剖分的函数：

from scipy.spatial import Delaunay

def triangulate_points(x, y):
    """二维点的 Delaunay 三角剖分"""

    # centers, edges, tri, neighbors = Delaunay(x, y)
    tri = Delaunay(np.c_[x, y]).simplices
    return tri

??函数输出的是一个数组，该数组的每一行包含对应数组 x 和 y 中每个三角形三个点的切片。

??现在将该算法应用于一个例子，在该例子中，在 5×6 的网格上使用 30 个控制点，将一幅图像扭曲到另一幅图像中的非平坦区域。图 3-5b 所示的是将一幅图像扭曲到“turning torso”的表面。目标点是使用 ginput() 函数手工选取出来的，将结果保存在 turningtorso_points.txt 文件中；

??然后我们使用下面简短的脚本将这些操作统一起来：

 # -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.geometry import warp

# 打开图像，并将其扭曲
fromim = array(Image.open('image/sunset_tree.jpg'))
x, y = meshgrid(range(5), range(6))
x = (fromim.shape[1]/4) * x.flatten()
y = (fromim.shape[0]/5) * y.flatten()

# 三角剖分
tri = warp.triangulate_points(x, y)

# 打开图像和目标点
im = array(Image.open('image/turningtorso1.jpg'))
tp = loadtxt('turningtorso1_points.txt') # destination points

figure()
subplot(1, 4, 1)
axis('off')
imshow(im)

# 将点转换成齐次坐标
fp = array(vstack((y, x, ones((1, len(x))))), 'int')
tp = array(vstack((tp[:, 1], tp[:, 0], ones((1, len(tp))))), 'int')

# 扭曲三角形
im = warp.pw_affine(fromim, im, fp, tp, tri)

# 绘制图像
subplot(1, 4, 2)
axis('off')
imshow(fromim)
warp.plot_mesh(fp[1], fp[0], tri)
subplot(1, 4, 3)
axis('off')
imshow(im)
subplot(1, 4, 4)
axis('off')
imshow(im)
warp.plot_mesh(tp[1], tp[0], tri)
show()

??再通过下面的辅助函数（将其添加到 warp.py 文件中）来绘制出图像中的这些三角形：

def plot_mesh(x,y,tri):
    """ 绘制三角形"""
    for t in tri:
        t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
        plot(x[t_ext],y[t_ext],'r')

运行结果：

??使用狄洛克三角剖分标记点进行分段仿射扭曲：第一张为带有标记物的目标图像；第二张为带有三角剖分的图像；第三张为扭曲后的图像；第四张为带有三角剖分的扭曲图像。

如果运行教材上的代码报错：ModuleNotFoundError: No module named ‘matplotlib.delaunay’
将 import matplotlib.delaunay as md 改成 from scipy.spatial import Delaunay
并将 triangulate_points() 函数中的 centers, edges, tri, neighbors = Delaunay(x, y) 改成 tri = Delaunay(np.c_[x,y]).simplices

使用狄洛克三角剖分标记点进行分段仿射扭曲使用的 turningtorso1_points.txt 文件可以从该网站上下载: Programming Computer Vision with Python