Python计算机视觉 第3章-图像到图像的映射
3.1 单应性变换
单应性变换（Homography）是计算机视觉中非常重要的一种几何变换，它用于将一个平面内的点映射到另一个平面内。具体来说，单应性变换可以描述一个图像在摄像机视角变化、平面移动或旋转时，如何从一个视角变换到另一个视角。
这种变换在多个应用场景中非常有用，比如：
- 图像配准：将不同视角或不同时间拍摄的图像对齐，找到它们之间的对应关系。
- 图像校正：修正由于摄像机角度或透视导致的图像扭曲，使图像看起来更平整。
- 纹理扭曲：将一个平面的纹理准确地映射到另一个平面上。
- 全景图像创建：将多个图像拼接成一个大的全景图像。
单应性变换的频繁使用，尤其是在涉及多个视角或需要精确对齐图像的情况下，能够显著提升算法的鲁棒性和精度。在项目中，理解和正确应用单应性变换是处理图像和三维几何信息的关键技能。
单应性变换（Homography）将二维平面上的点映射到另一个平面上的点，在齐次坐标（homogeneous coordinates）下，这种映射可以通过以下方程来表示：
( x ′ y ′ w ′ ) H ⋅ ( x y w ) \begin{pmatrix} x' \\ y' \\ w' \ \end{pmatrix} \mathbf{H} \cdot \begin{pmatrix} x \\ y \\ w \end{pmatrix} x′y′w′ H⋅ xyw
其中，单应性矩阵 H \mathbf{H} H 为：
H = ( h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ) \mathbf{H} = \begin{pmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{pmatrix} H= h11h21h31h12h22h32h13h23h33
经过单应性变换后的目标点的常规二维坐标 ( x ′ , y ′ ) (x', y') (x′,y′) 为：
x ′ = h 11 x + h 12 y + h 13 h 31 x + h 32 y + h 33 x' = \frac{h_{11}x + h_{12}y + h_{13}}{h_{31}x + h_{32}y + h_{33}} x′=h31x+h32y+h33h11x+h12y+h13
y ′ = h 21 x + h 22 y + h 23 h 31 x + h 32 y + h 33 y' = \frac{h_{21}x + h_{22}y + h_{23}}{h_{31}x + h_{32}y + h_{33}} y′=h31x+h32y+h33h21x+h22y+h23
通过这些公式，你可以描述平面间的各种变换，比如旋转、缩放、平移、透视变换等。
3.1.1 直接线性变换算法
单应性矩阵可以由两幅图像（或者平面）中对应点对计算出来。前面已经提到过，一个完全射影变换具有8个自由度。根据对应点约束，每个对应点对可以写出两个方程，分别对应于x和y坐标。因此，计算单应性矩阵H需要４个对应点对。
DLT（Direct Linear Transformation，直接线性变换）是给定４个或者更多对应点对矩阵，来计算单应性矩阵H的算法。将单应性矩阵H作用在对应点对上，重新写出该方程，我们可以得到下面的方程：
[ − x 1 − y 1 − 1 0 0 0 x 1 x 1 ′ y 1 x 1 ′ x 1 ′ 0 0 0 − x 1 − y 1 − 1 x 1 y 1 ′ y 1 y 1 ′ y 1 ′ − x 2 − y 2 − 1 0 0 0 x 2 x 2 ′ y 2 x 2 ′ x 2 ′ 0 0 0 − x 2 − y 2 − 1 x 2 y 2 ′ y 2 y 2 ′ y 2 ′ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ] [ h 1 h 2 h 3 h 4 h 5 h 6 h 7 h 8 h 9 ] = 0 \begin{bmatrix}-x_1&-y_1&-1&0&0&0&x_1x_1^{\prime}&y_1x_1^{\prime}&x_1^{\prime}\\0&0&0&-x_1&-y_1&-1&x_1y_1^{\prime}&y_1y_1^{\prime}&y_1^{\prime}\\-x_2&-y_2&-1&0&0&0&x_2x_2^{\prime}&y_2x_2^{\prime}&x_2^{\prime}\\0&0&0&-x_2&-y_2&-1&x_2y_2^{\prime}&y_2y_2^{\prime}&y_2^{\prime}\\&\vdots&&\vdots&\vdots&\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\\h_4\\h_5\\h_6\\h_7\\h_8\\h_9\end{bmatrix}=\mathbf{0} −x10−x20−y10−y20⋮−10−100−x10−x2⋮0−y10−y2⋮0−10−1⋮x1x1′x1y1′x2x2′x2y2′⋮y1x1′y1y1′y2x2′y2y2′⋮x1′y1′x2′y2′ h1h2h3h4h5h6h7h8h9 =0
或者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]
上面函数的第一步操作是检查点对的两个数组中点的数目是否相同。如果不相同，函数将会抛出异常信息。这对于写出稳健的代码来说非常有用。但是，为了使得代码例子更简单、更容易理解，我们在本书中仅在很少的例子中使用异常处理技巧。
3.1.2 仿射变换
由于仿射变换具有6个自由度，因此我们需要三个对应点对来估计矩阵H。通过将最后两个元素设置为0，即
h
7
=
h
8
=
0
h_7 =h_8=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, Second 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算法，这些点需要经过预处理和去处理化操作。
3.2 图像扭曲
对图像块应用仿射变换，我们将其称为图像扭曲（或者仿射扭曲）。该操作不仅经常应用在计算机图形学中，而且经常出现在计算机视觉算法中。扭曲操作可以使用SciPy工具包中的ndimage包来简单完成。
以下为实验代码：
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from skimage import data, color
# 读取示例图像
image = color.rgb2gray(data.astronaut())
# 定义仿射变换矩阵
# 例如，这里是一个旋转矩阵和一个平移矩阵的组合
affine_matrix = np.array([
[1.2, 0.2, -30], # x轴的缩放和旋转，以及平移
[0.1, 1.2, 20], # y轴的缩放和旋转，以及平移
[0, 0, 1] # 齐次坐标的归一化因子
])
# 对图像应用仿射变换
transformed_image = ndimage.affine_transform(
image,
affine_matrix[:2, :2], # 2x2 仿射矩阵
offset=affine_matrix[:2, 2], # 平移偏移
mode='reflect' # 边界处理模式
)
# 显示原始和变换后的图像
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.title('Original Image')
plt.imshow(image, cmap='gray')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title('Transformed Image')
plt.imshow(transformed_image, cmap='gray')
plt.axis('off')
plt.show()
3.2.1 图像中的图像
仿射扭曲的一个简单例子是，将图像或者图像的一部分放置在另一幅图像中，使得它们能够和指定的区域或者标记物对齐。
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from skimage import io, color
def image_in_image(background, overlay, position):
"""
将图像 overlay 放置到图像 background 中的指定位置。
:param background: 背景图像
:param overlay: 要放置的图像
:param position: 放置 overlay 的坐标 (x, y) 元组
:return: 带有 overlay 的背景图像
"""
# 确保 overlay 图像的尺寸
h, w = overlay.shape[:2]
# 生成仿射变换矩阵，将 overlay 图像的四个角点对齐到背景图像中的指定区域
src_points = np.array([
[0, 0], # overlay 的左上角
[w, 0], # overlay 的右上角
[w, h], # overlay 的右下角
[0, h] # overlay 的左下角
])
dst_points = np.array([
[position[0], position[1]], # 背景图像中放置位置的左上角
[position[0] + w, position[1]], # 背景图像中放置位置的右上角
[position[0] + w, position[1] + h], # 背景图像中放置位置的右下角
[position[0], position[1] + h] # 背景图像中放置位置的左下角
])
# 构建矩阵 A 和向量 b 以求解仿射变换矩阵
A = []
b = []
for i in range(4):
A.append([src_points[i][0], src_points[i][1], 1, 0, 0, 0, -dst_points[i][0] * src_points[i][0],
-dst_points[i][0] * src_points[i][1]])
A.append([0, 0, 0, src_points[i][0], src_points[i][1], 1, -dst_points[i][1] * src_points[i][0],
-dst_points[i][1] * src_points[i][1]])
b.append(dst_points[i][0])
b.append(dst_points[i][1])
A = np.array(A)
b = np.array(b)
# 通过最小二乘法求解仿射变换矩阵
h = np.linalg.lstsq(A, b, rcond=None)[0]
H = np.append(h, [1]).reshape(3, 3)
# 将 overlay 图像进行仿射变换
transformed_overlay = ndimage.affine_transform(
overlay,
H[:2, :2],
offset=H[:2, 2],
output_shape=background.shape,
mode='constant',
cval=0
)
# 合并图像
mask = (transformed_overlay > 0).astype(float)
result = background.copy()
result = result * (1 - mask) + transformed_overlay * mask
return result
# 示例使用
if __name__ == "__main__":
# 读取内置示例图像
background = color.rgb2gray(io.imread('img.png')) # 背景图像
overlay = color.rgb2gray(io.imread('python.png')) # 要放置的图像
position = (100, 100) # 放置位置（x, y）
# 应用函数
result_image = image_in_image(background, overlay, position)
# 显示结果
plt.figure(figsize=(10, 5))
plt.subplot(1, 3, 1)
plt.title('Background Image')
plt.imshow(background, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.title('Overlay Image')
plt.imshow(overlay, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 3)
plt.title('Result Image')
plt.imshow(result_image, cmap='gray')
plt.axis('off')
plt.show()
3.2.2 分段仿射扭曲
对应点对集合之间最常用的扭曲方式：分段仿射扭曲。给定任意图像的标记点，通过将这些点进行三角剖分，然后使用仿射扭曲来扭曲每个三角形，我们可以将图像和另一幅图像的对应标记点扭曲对应。对于任何图形和图像处理库来说，这些都是最基本的操作。
为了三角化这些点，我们经常使用狄洛克三角剖分方法。在Matplotlib（但是不在PyLab 库中）中有狄洛克三角剖分，我们可以用下面的方式使用它：
以下是实验代码：
import matplotlib.pyplot as plt
import numpy as np
from scipy import ndimage
from scipy.spatial import Delaunay
x,y = np.array(np.random.standard_normal((2,100)))
tri = Delaunay(np.c_[x, y]).simplices
plt.figure()
for t in tri:
t_ext = [t[0], t[1], t[2], t[0]] # 将第一个点加入到最后
plt.plot(x[t_ext],y[t_ext],'r')
plt.plot(x,y,'*')
plt.axis('off')
plt.show()
3.2.3 图像配准
图像配准是对图像进行变换，使变换后的图像能够在常见的坐标系中对齐。配准可以是严格配准，也可以是非严格配准。为了能够进行图像对比和更精细的图像分析，图像配准是一步非常重要的操作。
3.3 创建全景图
在同一位置（即图像的照相机位置相同）拍摄的两幅或者多幅图像是单应性相关的（如图3-9所示）。我们经常使用该约束将很多图像缝补起来，拼成一个大的图像来创建全景图像。
3.3.1 RANSAC
RANSAC是“RANdom SAmple Consensus”（随机一致性采样）的缩写。该方法是用来找到正确模型来拟合带有噪声数据的迭代方法。给定一个模型，例如点集之间的单应性矩阵，RANSAC基本的思想是，数据中包含正确的点和噪声点，合理的模型应该能够在描述正确数据点的同时摒弃噪声点。
RANSAC的标准例子：用一条直线拟合带有噪声数据的点集。简单的最小二乘在该例子中可能会失效，但是RANSAC能够挑选出正确的点，然后获取能够正确拟合的直线。
3.3.2 拼接图像
估计出图像间的单应性矩阵（使用RANSAC算法），现在我们需要将所有的图像扭曲到一个公共的图像平面上。通常，这里的公共平面为中心图像平面（否则，需要进行大量变形）。一种方法是创建一个很大的图像，比如图像中全部填充0，使其和中心图像平行，然后将所有的图像扭曲到上面。由于我们所有的图像是由照相机水平旋转拍摄的，因此我们可以使用一个较简单的步骤：将中心图像左边或者右边的区域填充０，以便为扭曲的图像腾出空间。以下为示例代码：
def panorama(H, fromim, toim, padding=2400, delta=2400):
""" 使用单应性矩阵 H（使用 RANSAC 健壮性估计得出），协调两幅图像，创建水平全景图像。结果为一幅和 toim 具有相同高度的图像。padding 指定填充像素的数目，delta 指定额外的平移量。 """
# 检查图像是灰度图像，还是彩色图像
is_color = len(fromim.shape) == 3
# 用于 geometric_transform() 的单应性变换
def transf(p):
p2 = np.dot(H, [p[0], p[1], 1])
return (p2[0] / p2[2], p2[1] / p2[2])
if H[1, 2] < 0: # fromim 在右边
print('warp - right')
# 变换 fromim
if is_color:
# 在目标图像的右边填充 0
toim_t = np.hstack((toim, np.zeros((toim.shape[0], padding, 3))))
fromim_t = np.zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
for col in range(3):
fromim_t[:, :, col] = ndimage.geometric_transform(
fromim[:, :, col], transf, (toim.shape[0], toim.shape[1] + padding))
else:
# 在目标图像的右边填充 0
toim_t = np.hstack((toim, np.zeros((toim.shape[0], padding))))
fromim_t = ndimage.geometric_transform(fromim, transf, (toim.shape[0], toim.shape[1] + padding))
else: # fromim 在左边
print('warp - left')
# 为了补偿填充效果，在左边加入平移量
H_delta = np.array([[1, 0, 0], [0, 1, -delta], [0, 0, 1]])
H = np.dot(H, H_delta)
# 变换 fromim
if is_color:
# 在目标图像的左边填充 0
toim_t = np.hstack((np.zeros((toim.shape[0], padding, 3)), toim))
fromim_t = np.zeros((toim.shape[0], toim.shape[1] + padding, toim.shape[2]))
for col in range(3):
fromim_t[:, :, col] = ndimage.geometric_transform(
fromim[:, :, col], transf, (toim.shape[0], toim.shape[1] + padding))
else:
# 在目标图像的左边填充 0
toim_t = np.hstack((np.zeros((toim.shape[0], padding)), toim))
fromim_t = ndimage.geometric_transform(fromim, transf, (toim.shape[0], toim.shape[1] + padding))
# 协调后返回（将 fromim 放置在 toim 上）
if is_color:
# 所有非黑色像素
alpha = ((fromim_t[:, :, 0] > 0) | (fromim_t[:, :, 1] > 0) | (fromim_t[:, :, 2] > 0))
for col in range(3):
toim_t[:, :, col] = fromim_t[:, :, col] * alpha + toim_t[:, :, col] * ~alpha
else:
alpha = (fromim_t > 0)
toim_t = fromim_t * alpha + toim_t * ~alpha
return toim_t
对于通用的geometric_transform() 函数，我们需要指定能够描述像素到像素间映射的函数。在这个例子中，transf()函数就是该指定的函数。该函数通过将像素和H相乘，然后对齐次坐标进行归一化来实现像素间的映射。通过查看H中的平移量，我们可以决定应该将该图像填补到左边还是右边。当该图像填补到左边时，由于目标图像中点的坐标也变化了，所以在“左边”情况中，需要在单应性矩阵中加入平移。简单起见，我们同样使用0像素的技巧来寻找alpha图。现在在图像中使用该操作，函数如下所示：
# 扭曲图像
delta = 2000 # 用于填充和平移
# 读取图像
im1 = np.array(Image.open(imname[1]))
im2 = np.array(Image.open(imname[2]))
# 图像拼接
im_12 = warp.panorama(H_12, im1, im2, delta, delta)
im1 = np.array(Image.open(imname[0]))
im_02 = warp.panorama(np.dot(H_12, H_01), im1, im_12, delta, delta)
im1 = np.array(Image.open(imname[3]))
im_32 = warp.panorama(H_32, im1, im_02, delta, delta)
im1 = np.array(Image.open(imname[j + 1])) # 确保 imname[j + 1] 是一个有效的索引
im_42 = warp.panorama(np.dot(H_32, H_43), im1, im_32, delta, 2 * delta)
