§1.2.8

图像旋转、缩放、剪切的矩阵形式？逆变换 + 双线性插值的反向映射做法？

手写练习

—用 numpy + 双线性插值实现 warp_affine，与 cv2.warpAffine 对比

核心概念

仿射变换（Affine Transformation）是一种二维坐标间的线性变换，外加一个平移。它保持了图像的“平直性”（直线变换后仍是直线）和“平行性”（平行线变换后仍是平行线）。旋转（Rotation）、缩放（Scaling）、剪切（Shearing）都是仿射变换的特例。任何仿射变换都可以用一个 2x3 的矩阵来表示，对图像进行操作。

为了避免变换后图像出现孔洞或像素重叠，实际应用中采用“反向映射”（Backward Mapping）或“逆变换”（Inverse Warping）的策略。即遍历目标图像的每一个像素，通过逆变换矩阵计算它在源图像中的（非整数）坐标，再通过插值算法（如双线性插值）得到该点的像素值。

原理与推导

1. 仿射变换矩阵

一个二维点 $(x, y)$ 经过仿射变换后得到新点 $(x', y')$ ，可以用矩阵运算表示：

\begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix} t_x \\ t_y \end{bmatrix}

其中，左边的 2x2 矩阵负责线性变换（旋转、缩放、剪切），右边的向量负责平移。

为了将线性和平移统一到一个矩阵运算中，我们引入齐次坐标 (Homogeneous Coordinates)。二维点 $(x, y)$ 表示为 $(x, y, 1)$ 。变换矩阵扩展为 3x3：

\begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix} a & b & t_x \\ c & d & t_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}

由于最后一行总是 [0, 0, 1]，因此在实践中通常用一个 2x3 的矩阵 $M = \begin{bmatrix} a & b & t_x \\ c & d & t_y \end{bmatrix}$ 来表示。

各种基本变换的矩阵形式（以原点为中心）：

缩放 (Scaling)：沿 x 轴缩放 $s_x$ 倍，沿 y 轴缩放 $s_y$ 倍。

M_{scale} = \begin{bmatrix} s_x & 0 & 0 \\ 0 & s_y & 0 \\ 0 & 0 & 1 \end{bmatrix}

旋转 (Rotation)：逆时针旋转 $\theta$ 角度。

M_{rotate} = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}

剪切 (Shearing)：沿 x 轴剪切 $sh_x$ ，沿 y 轴剪切 $sh_y$ 。

M_{shear} = \begin{bmatrix} 1 & sh_x & 0 \\ sh_y & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}

平移 (Translation)：沿 x 轴平移 $t_x$ ，沿 y 轴平移 $t_y$ 。

M_{translate} = \begin{bmatrix} 1 & 0 & t_x \\ 0 & 1 & t_y \\ 0 & 0 & 1 \end{bmatrix}

组合变换：复杂的仿射变换可以通过矩阵连乘实现，例如先缩放，再旋转，最后平移，其变换矩阵为 $M = M_{translate} \cdot M_{rotate} \cdot M_{scale}$ 。注意矩阵乘法不满足交换律，变换顺序至关重要。

2. 反向映射与双线性插值

动机：

正向映射：遍历源图像每个像素 $(u, v)$ ，计算其在目标图像中的位置 $(x', y') = M \cdot (u, v, 1)^T$ 。这会导致两个问题：1) $(x', y')$ 可能是浮点数，取整会造成精度损失；2) 放大时，目标图像上会出现“孔洞”（没有源像素映射过来）；缩小时，多个源像素可能映射到同一个目标像素，造成信息丢失或“混叠”。
反向映射：遍历目标图像的每个像素 $(x, y)$ $(x, y)$ ，通过逆变换矩阵 $M^{-1}$ 找到它在源图像中对应的位置 $(u', v') = M^{-1} \cdot (x, y, 1)^T$ $(u^{'}, v^{'}) = M^{- 1} \cdot (x, y, 1)^{T}$ 。
- 优点：保证目标图像的每个像素都被赋值，不会产生孔洞。
- 挑战：计算出的源坐标 $(u', v')$ 通常是浮点数，而图像像素值定义在整数格点上。

双线性插值 (Bilinear Interpolation)：为了解决浮点坐标的像素值问题，我们采用双线性插值。假设我们要在源图像中找到浮点坐标 $(u, v)$ 处的像素值。

找到包围 $(u, v)$ 的四个最近的整数坐标像素点：
- $Q_{11} = (u_1, v_1) = (\lfloor u \rfloor, \lfloor v \rfloor)$
- $Q_{21} = (u_2, v_1) = (\lfloor u \rfloor + 1, \lfloor v \rfloor)$
- $Q_{12} = (u_1, v_2) = (\lfloor u \rfloor, \lfloor v \rfloor + 1)$
- $Q_{22} = (u_2, v_2) = (\lfloor u \rfloor + 1, \lfloor v \rfloor + 1)$ 令它们对应的像素值分别为 $I(Q_{11}), I(Q_{21}), I(Q_{12}), I(Q_{22})$ 。
计算 $(u, v)$ 与左上角点 $(\lfloor u \rfloor, \lfloor v \rfloor)$ 的距离比例：
- $dx = u - \lfloor u \rfloor$
- $dy = v - \lfloor v \rfloor$
几何直观：先在 x 方向进行两次线性插值，再在 y 方向进行一次线性插值。
- 在 $v_1$ 行上插值： $I(R_1) = (1 - dx) \cdot I(Q_{11}) + dx \cdot I(Q_{21})$
- 在 $v_2$ 行上插值： $I(R_2) = (1 - dx) \cdot I(Q_{12}) + dx \cdot I(Q_{22})$
- 在 $R_1$ 和 $R_2$ 之间沿 y 方向插值得到最终结果： $I(u, v) = (1 - dy) \cdot I(R_1) + dy \cdot I(R_2)$
合并公式：
$I(u, v) = (1-dx)(1-dy)I(Q_{11}) + dx(1-dy)I(Q_{21}) + (1-dx)dyI(Q_{12}) + dx \cdot dy \cdot I(Q_{22})$
这可以看作是对四个邻近像素值的加权平均，权重与面积成反比。

3. 算法复杂度

对于一个大小为 $W_{out} \times H_{out}$ 的目标图像：

时间复杂度: $O(W_{out} \cdot H_{out})$ 。因为需要遍历目标图像的每个像素，而每个像素的反向映射和双线性插值操作都是常数时间 $O(1)$ 。
空间复杂度: $O(W_{out} \cdot H_{out})$ 。需要额外的空间来存储生成的目标图像。

代码实现

下面的代码使用 NumPy 从零开始实现一个仿射变换函数 warp_affine_numpy，它内部使用反向映射和双线性插值。然后将其结果与 OpenCV 高度优化的 cv2.warpAffine 函数进行对比。

python

 1import cv2
 2import numpy as np
 3import matplotlib.pyplot as plt
 4
 5def warp_affine_numpy(src_img, M, dsize):
 6    """
 7    使用 Numpy 和双线性插值实现仿射变换，模拟 cv2.warpAffine。
 8    
 9    Args:
10        src_img (np.ndarray): 源图像.
11        M (np.ndarray): 2x3 的仿射变换矩阵.
12        dsize (tuple): 输出图像的尺寸 (width, height).
13    
14    Returns:
15        np.ndarray: 变换后的图像.
16    """
17    H_src, W_src = src_img.shape[:2]
18    W_out, H_out = dsize
19    
20    # 创建目标图像
21    dst_img = np.zeros((H_out, W_out, src_img.shape[2]), dtype=src_img.dtype)
22    
23    # 为什么这样做: 为了进行反向映射，我们需要变换矩阵 M 的逆。
24    # M 是 2x3，无法直接求逆。我们将其扩展为 3x3 的齐次矩阵再求逆。
25    M_homo = np.vstack([M, [0, 0, 1]])
26    M_inv = np.linalg.inv(M_homo)
27    
28    # 为什么这样做: 创建目标图像的网格坐标。这是反向映射的第一步，
29    # 我们要为目标图像的每一个像素 (x_dst, y_dst) 找到它在源图像中的对应位置。
30    x_dst, y_dst = np.meshgrid(np.arange(W_out), np.arange(H_out))
31    
32    # 为什么这样做: 将网格坐标转换为齐次坐标形式 (x, y, 1) 的列向量。
33    # 这样我们就可以通过一次矩阵乘法，为所有目标像素计算出其在源图像的坐标。
34    # 形状: (3, H_out * W_out)
35    dst_coords = np.vstack((x_dst.ravel(), y_dst.ravel(), np.ones(H_out * W_out)))
36    
37    # 为什么这样做: 应用逆变换，计算出所有目标像素在源图像中的对应坐标。
38    # src_coords = M_inv @ dst_coords
39    # 形状: (3, H_out * W_out)
40    src_coords_homo = M_inv @ dst_coords
41    
42    # 为什么这样做: 从齐次坐标转换回2D坐标，并处理透视除法（虽然仿射变换的w总是1）。
43    # 然后将坐标 reshape 回图像的 HxW 形状。
44    u_src = (src_coords_homo[0, :] / src_coords_homo[2, :]).reshape(H_out, W_out)
45    v_src = (src_coords_homo[1, :] / src_coords_homo[2, :]).reshape(H_out, W_out)
46
47    # --- 双线性插值 ---
48    # 为什么这样做: 获取插值所需的四个最近邻像素的整数坐标。
49    u1 = np.floor(u_src).astype(int)
50    v1 = np.floor(v_src).astype(int)
51    u2 = u1 + 1
52    v2 = v1 + 1
53    
54    # 为什么这样做: 计算插值权重。dx, dy 是浮点坐标与左上角整数坐标的距离比例。
55    dx = u_src - u1
56    dy = v_src - v1
57    # 需要扩展维度以匹配多通道图像的计算
58    dx = np.expand_dims(dx, axis=-1)
59    dy = np.expand_dims(dy, axis=-1)
60
61    # 为什么这样做: 边界检查。将超出源图像范围的坐标裁剪到边界内，
62    # 防止索引错误。这是一种 'border replicate' 的处理方式。
63    u1 = np.clip(u1, 0, W_src - 1)
64    u2 = np.clip(u2, 0, W_src - 1)
65    v1 = np.clip(v1, 0, H_src - 1)
66    v2 = np.clip(v2, 0, H_src - 1)
67    
68    # 为什么这样做: 获取四个邻近像素的像素值。
69    # NumPy的高级索引允许我们用坐标数组一次性取出所有需要的像素值。
70    I11 = src_img[v1, u1]
71    I21 = src_img[v1, u2]
72    I12 = src_img[v2, u1]
73    I22 = src_img[v2, u2]
74    
75    # 为什么这样做: 应用双线性插值公式。这是核心计算步骤。
76    # (1-dx)(1-dy)I11 + dx(1-dy)I21 + (1-dx)dyI12 + dx*dy*I22
77    w1 = (1 - dx) * (1 - dy)
78    w2 = dx * (1 - dy)
79    w3 = (1 - dx) * dy
80    w4 = dx * dy
81    
82    interpolated_values = w1 * I11 + w2 * I21 + w3 * I12 + w4 * I22
83    
84    # 为什么这样做: 创建一个掩码，标记那些映射到源图像外部的像素。
85    # 原始的 u_src, v_src 坐标在 0 到 W_src-1/H_src-1 范围之外。
86    mask = (u_src < 0) | (u_src >= W_src) | (v_src < 0) | (v_src >= H_src)
87    
88    # 为什么这样做: 将插值结果填充到目标图像中，并将边界外的像素设置为0（黑色）。
89    # 这模拟了 cv2.warpAffine 的 BORDER_CONSTANT 行为。
90    dst_img = interpolated_values.astype(src_img.dtype)
91    dst_img[mask] = 0
92    
93    return dst_img
94
95if __name__ == '__main__':
96    # 1. 创建一个测试图像
97    img = cv2.imread('lena.jpg') # 请准备一张名为 lena.jpg 的图片
98    if img is None:
99        print("警告: lena.jpg 未找到。将创建一个棋盘格图像用于测试。")
100        img = np.zeros((256, 256, 3), dtype=np.uint8)
101        img[::16, :, :] = 255
102        img[:, ::16, :] = 255
103        img[8:248:16, 8:248:16, :] = 128
104        
105    H, W = img.shape[:2]
106    dsize = (W, H)
107    
108    # 2. 定义一个仿射变换矩阵 (旋转30度，缩小0.8倍)
109    center = (W / 2, H / 2)
110    angle = 30
111    scale = 0.8
112    M = cv2.getRotationMatrix2D(center, angle, scale)
113    
114    # 3. 使用我们自己的 Numpy 实现
115    img_warped_numpy = warp_affine_numpy(img, M, dsize)
116    
117    # 4. 使用 OpenCV 的实现
118    img_warped_cv2 = cv2.warpAffine(img, M, dsize)
119    
120    # 5. 对比结果
121    plt.figure(figsize=(15, 5))
122    plt.subplot(1, 3, 1)
123    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
124    plt.title("原始图像")
125    plt.axis('off')
126    
127    plt.subplot(1, 3, 2)
128    plt.imshow(cv2.cvtColor(img_warped_numpy, cv2.COLOR_BGR2RGB))
129    plt.title("Numpy 实现 (warp_affine_numpy)")
130    plt.axis('off')
131    
132    plt.subplot(1, 3, 3)
133    plt.imshow(cv2.cvtColor(img_warped_cv2, cv2.COLOR_BGR2RGB))
134    plt.title("OpenCV 实现 (cv2.warpAffine)")
135    plt.axis('off')
136    
137    plt.tight_layout()
138    plt.show()
139
140    # 计算差异
141    diff = np.mean(np.abs(img_warped_numpy.astype(float) - img_warped_cv2.astype(float)))
142    print(f"Numpy 实现和 OpenCV 实现的平均像素绝对差值: {diff:.4f}")
143    print("注：微小的差异是由于浮点数精度和边界处理的细微不同导致的，属于正常现象。")

工程实践

数据增强: 仿射变换是计算机视觉，特别是深度学习中，最常用的数据增强手段。通过对训练图像进行随机的旋转、缩放、剪切、平移，可以极大地扩充数据集，提高模型的泛化能力。torchvision.transforms 和 albumentations 等库提供了高效的实现。
插值方法选择:
- cv2.INTER_NEAREST (最近邻插值): 最快，但效果最差，图像会出现锯齿。适用于对掩码（mask）或标签图像进行变换，因为不会产生新的像素值。
- cv2.INTER_LINEAR (双线性插值): 默认选项，速度和效果的良好平衡。适用于大多数场景。
- cv2.INTER_CUBIC (双三次插值): 考虑16个邻域像素，效果更平滑，细节保留更好，但计算量更大。
- cv2.INTER_LANCZOS4 (Lanczos插值): 考虑64个邻域像素，效果更锐利，但可能在边缘产生振铃效应，计算最慢。
性能优化: 在实际项目中，应优先使用高度优化的库如 OpenCV, Pillow-SIMD。在深度学习训练流程中，应使用 GPU 加速的变换，如 PyTorch 的 torch.nn.functional.grid_sample，其原理与本节描述的反向映射+插值完全一致，但利用了 GPU 的并行计算能力，速度极快。
变换中心: cv2.getRotationMatrix2D 允许指定旋转中心。要围绕点 $(c_x, c_y)$ 旋转，实际执行的变换是：先将图像平移 $(-c_x, -c_y)$ ，再绕原点旋转，最后平移回 $(c_x, c_y)$ 。
边界处理: cv2.warpAffine 的 borderMode 参数可以控制映射到源图像外部的像素如何填充，常用选项有 BORDER_CONSTANT (填充常数，默认为0), BORDER_REPLICATE (复制边界像素), BORDER_REFLECT (反射)。

常见误区与边界情况

误区：正向映射 vs 反向映射: 初学者常会想到用正向映射，但这会导致前述的“孔洞”和“混叠”问题。面试中必须清晰地解释为什么工业界和学术界标准实现都采用反向映射。
误区：坐标系混淆: 图像处理中 (row, col) 坐标系（NumPy 风格）和几何中的 (x, y) 坐标系（OpenCV 风格）常常混用。x 对应 col，y 对应 row。在实现时必须保持坐标系定义的一致性，否则结果会错乱。
误区：矩阵乘法顺序: 变换的组合顺序会影响最终结果。例如，“先旋转再平移”和“先平移再旋转”得到的变换矩阵和结果图像都不同。
边界情况：变换矩阵不可逆: 如果仿射变换矩阵 M 是奇异的（例如，缩放因子为0，将图像压缩到一条线或一个点），其3x3齐次矩阵将不可逆。此时 np.linalg.inv 会抛出异常。
边界情况：数值精度: 自己实现的 NumPy 版本和 OpenCV 版本可能存在微小的像素值差异（通常小于1）。这是由于浮点数计算的精度误差、插值算法的内部实现细节差异导致的，属于正常现象。
面试追问:
- "如何对 bounding box 进行仿射变换？": 对 bounding box 的四个角点应用仿射变换，然后取变换后四个新点的 min/max 坐标，形成一个新的、与坐标轴平行的外包矩形。注意，旋转后的物体其 AABB（Axis-Aligned Bounding Box）会变大。
- "如果变换不是仿射变换，比如透视变换呢？": 原理相同，仍然使用反向映射。只是变换矩阵变为 3x3 的单应性矩阵（Homography），逆变换公式更复杂（涉及透视除法）。PyTorch 的 grid_sample 和 OpenCV 的 warpPerspective 就是处理这种情况。
- "如何在 GPU 上实现这个过程？": 解释 torch.nn.functional.grid_sample 的工作原理：它接收一个输入特征图和一个“采样网格”（sampling grid）。这个网格的每个 (x, y) 位置都存储了源特征图上的一个采样坐标。grid_sample 内部会完成反向查找和双线性（或最近邻）插值，整个过程在 GPU 上并行完成，效率极高。