程序设计与计算思维¶

Computer Programming and Computational Thinking

第 5 讲：图像变换 III¶

2025—2026学年度春季学期

清华大学 韩文弢

内容回顾¶

  • Python 语言
    • 元组
    • 函数（参数默认值、关键字参数、Lambda 表达式、函数对象）
  • 变换
    • 线性变换与非线性变换
    • 线性变换与矩阵

准备工作¶

In [1]:
from urllib.request import urlretrieve
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets

# 设置作图时的中文字体，避免乱码
plt.rcParams['font.sans-serif'] = ['Noto Sans CJK SC']

img_path, _ = urlretrieve('https://pacman.cs.tsinghua.edu.cn/~hanwentao/cpct/05/longcorgi.png')
img = plt.imread(img_path)[:, :, :3]
plt.imshow(img)
Out[1]:
<matplotlib.image.AxesImage at 0x10aa99be0>
No description has been provided for this image

最终效果¶

In [2]:
black = np.array([0, 0, 0])
white = np.array([1, 1, 1])
green = np.array([0, 1, 0])
red = np.array([1, 0, 0])

from math import sin, cos, hypot, atan2

# 恒等
identity = lambda x, y: (x, y)

# 伸缩
scalex = lambda alpha: lambda x, y: (alpha*x, y)
scaley = lambda alpha: lambda x, y: (x, alpha*y)
scale = lambda alpha: lambda x, y: (alpha*x, alpha*y)

# 翻转
swap = lambda x, y: (y, x)
flipy = lambda x, y: (x, -y)

# 旋转
rotate = lambda theta: lambda x, y: (cos(theta)*x + sin(theta)*y, -sin(theta)*x + cos(theta)*y)

# 剪切
shear = lambda alpha: lambda x, y: (x + alpha*y, y)

# 一般线性变换
lin = lambda a, b, c, d: lambda x, y: (a*x + b*y, c*x + d*y)

# 用矩阵表示
def lin_mat(A):
    a = A[0, 0]
    b = A[0, 1]
    c = A[1, 0]
    d = A[1, 1]
    return lambda x, y: lin(a, b, c, d)

# 平移，是仿射变换，但不是线性变换
translate = lambda alpha, beta: lambda x, y: (x + alpha, y + beta)

# 非线性剪切
nonlin_shear = lambda alpha: lambda x, y: (x, y + alpha*x**2)

# 扭曲
warp = lambda alpha: lambda x, y: rotate(alpha*hypot(x, y))(x, y)

# 极坐标转直角坐标
xy = lambda rho, theta: (rho*cos(theta), rho*sin(theta))

# 直角坐标转极坐标
rho_theta = lambda x, y: (hypot(x, y), atan2(y, x))


def transform_image(img, transform_func, pixels=400): 
    """
    对图像应用变换
    
    参数:
    - img: 输入图像 (NumPy 数组)
    - transform_func: 变换函数 (x, y) -> (new_x, new_y)
    - pixels: 输出图像的像素数
    """
    rows, cols = img.shape[:2]
    m = max(cols, rows)
    
    # 创建输出图像
    output = np.zeros((pixels, pixels, img.shape[2]), dtype=img.dtype)
    
    for i in range(pixels):
        for j in range(pixels):
            # 将像素坐标转换为 [-1, 1] 范围
            x, y = scale(2/pixels)(*flipy(*swap(*translate(-pixels/2, -pixels/2)(i, j))))
            
            # 应用变换的逆（因为我们是从输出像素找输入像素）
            new_x, new_y = transform_func(x, y)
            
            # 将坐标转换回像素索引
            src_i, src_j = translate(rows/2, cols/2)(*swap(*flipy(*scale(m/2)(new_x, new_y))))
            src_i = int(src_i)
            src_j = int(src_j)
            
            if 0 <= src_i < rows and 0 <= src_j < cols:
                output[i, j] = img[src_i, src_j]
            else:
                output[i, j] = black  # 黑色背景
    
    return output


def interactive_transform(a=1.0, b=0.0, c=0.0, d=1.0):
    """
    交互式变换演示
    """
    transformed = transform_image(img, lin(a, b, c, d), pixels=400)
    
    plt.figure(figsize=(8, 8))
    plt.imshow(transformed)
    plt.title(f'变换矩阵: a={a:.1f}, b={b:.1f}, c={c:.1f}, d={d:.1f}')
    plt.axis('off')
    plt.show()

# 创建交互式控件
interact = widgets.interactive(
    interactive_transform,
    a=widgets.FloatSlider(value=1.0, min=-1.5, max=1.5, step=0.1, description='a:'),
    b=widgets.FloatSlider(value=0.0, min=-1.5, max=1.5, step=0.1, description='b:'),
    c=widgets.FloatSlider(value=0.0, min=-1.5, max=1.5, step=0.1, description='c:'),
    d=widgets.FloatSlider(value=1.0, min=-1.5, max=1.5, step=0.1, description='d:'),
)

display(interact)

interactive(children=(FloatSlider(value=1.0, description='a:', max=1.5, min=-1.5), FloatSlider(value=0.0, desc…

对象变换与坐标变换¶

  • 对象变换是指通过改变对象的坐标来实现变换，如点 $(x, y) \to (2x, 2y)$ 是进行了 2 倍的放大
  • 坐标变换是指通过改变坐标轴来实现变换，例如 $x$ 和 $y$ 坐标轴都缩小一半，点的坐标也会发生上述变换

这两种变换的效果可能相同（点的坐标发生同样的变化），但是含义不同。

为了做想要的变换时有一致的效果，我们将图像的坐标系变换到 $[-1, 1]^2$ 的区域里再进行变换。

image.png

完整的过程¶

完整的变换过程如下所示，注意我们是从变换后的图像中取每个像素，然后计算对应变换前的图像中哪个像素：

image.png

逆变换¶

如果 $f$ 是从 2 维向量到 2 维向量的函数，我们定义 $f$ 的逆，记为 $f^{-1}$，具有“撤销” $f$ 效果的性质，即

$$f(f^{-1}(v)) = v$$

和

$$f^{-1}(f(v)) = v$$

这个等式可能对所有 $v$ 成立，也可能对某些 $v$ 成立。

示例：放大和缩小¶

In [3]:
# 验证缩放变换的逆
v = np.random.rand(2)

# scale(2) 的逆是 scale(0.5)
T = scale(0.5)
T_inv = scale(2)

# 应用 T 然后应用 T_inv
result = T_inv(*T(*v))

print(f'原始向量: {v}')
print(f'变换后: {result}')
print(f'相等？ {np.allclose(v, result)}')

原始向量: [0.82284975 0.77955379]
变换后: (np.float64(0.822849746691098), np.float64(0.7795537925491557))
相等？ True
In [4]:
# 验证旋转变换的逆
v = np.random.rand(2)

# rotate(θ) 的逆是 rotate(-θ)
theta = np.radians(30)
T = rotate(theta)
T_inv = rotate(-theta)

result = T_inv(*T(*v))

print(f'原始向量: {v}')
print(f'变换后: {result}')
print(f'相等？ {np.allclose(v, result)}')

原始向量: [0.11806597 0.27107506]
变换后: (np.float64(0.1180659687326121), np.float64(0.2710750580133623))
相等？ True

求解逆变换¶

逆变换到底是在做什么？

再次考虑缩放。假设将输入向量 $v$ 缩放 2 得到输出向量 $u$：

$$u = 2 v$$

现在假设想进行反向操作。如果给定 $u$，如何找到 $v$？

在这个特定情况下，有 $v = \frac{1}{2} u$。

考虑更加一般的情况，如果有一个线性变换，可以写成

$$u = A \, v$$

其中 $A$ 是矩阵。

如果给定 $u$ 并想反向找到 $v$，需要求解线性方程组。

通常可以（什么时候不可以？）解这个方程组找到一个新矩阵 $B$ 使得

$$v = B \, u$$

即 $B$ 撤销了 $A$ 的效果。那么，

$$v = (B \, A) v$$

因此 $B A$ 必须是单位矩阵，称 $B$ 为 $A$ 的矩阵逆，并写成

$$B = A^{-1}$$

对于 $2 \times 2$ 的矩阵，可以写出矩阵逆的公式。对于更大规模的矩阵，通常需要计算机运行算法来找到逆。

In [5]:
# 验证矩阵逆
A = np.random.randn(2, 2)
v = np.random.rand(2)

# 使用矩阵逆
A_inv = np.linalg.inv(A)

# 应用 A 然后应用 A_inv
result = A_inv @ (A @ v)

print(f'原始向量: {v}')
print(f'变换后: {result}')
print(f'相等？ {np.allclose(v, result)}')

原始向量: [0.83213606 0.74802034]
变换后: [0.83213606 0.74802034]
相等？ True
In [6]:
# 验证 (A*B)^(-1) = B^(-1) * A^(-1)
A = np.random.randn(2, 2)
B = np.random.randn(2, 2)

left = np.linalg.inv(A @ B)
right = np.linalg.inv(B) @ np.linalg.inv(A)

print(f'(A*B)^(-1):\n{left}')
print(f'B^(-1) * A^(-1):\n{right}')
print(f'相等？ {np.allclose(left, right)}')

(A*B)^(-1):
[[-0.81035038 -2.08689412]
 [ 0.52604683  0.31146011]]
B^(-1) * A^(-1):
[[-0.81035038 -2.08689412]
 [ 0.52604683  0.31146011]]
相等？ True

2x2 矩阵逆的公式¶

如果 $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$，则

$$A^{-1} = \frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix}$$

In [7]:
# 手动计算 2x2 矩阵逆
def inv_2x2_manual(A):
    a, b = A[0, 0], A[0, 1]
    c, d = A[1, 0], A[1, 1]
    det = a * d - b * c
    if abs(det) < 1e-10:
        raise ValueError("矩阵不可逆")
    return np.array([[d, -b], [-c, a]]) / det

A = np.random.randn(2, 2)
print(f'NumPy 逆:\n{np.linalg.inv(A)}')
print(f'手动计算逆:\n{inv_2x2_manual(A)}')
print(f'相等？ {np.allclose(np.linalg.inv(A), inv_2x2_manual(A))}')

NumPy 逆:
[[ 0.72936789 -1.0798709 ]
 [-0.11181915  2.1705302 ]]
手动计算逆:
[[ 0.72936789 -1.0798709 ]
 [-0.11181915  2.1705302 ]]
相等？ True

非线性变换的逆变换¶

如果有一个非线性变换 $T$，能求逆吗？换句话说，如果 $u = T(v)$，能否解这个方程找到用 $v$ 关于 $u$ 的函数表示？

通常这是一个困难的问题，有时可以解析地做到，但通常不能。

尽管如此，有一些数值方法有时可以解这些方程，例如牛顿法。在 scipy 中可以使用 fsolve 函数来进行数值求解。

In [8]:
from scipy.optimize import fsolve

# 使用数值方法求非线性变换的逆
def inverse_numerical(f, y, x0=(0.0, 0.0)):
    """
    数值求解非线性变换的逆
    f: 变换函数 (x, y) -> (u, v)
    y: 目标输出 (u, v)
    x0: 初始猜测
    """
    def equations(x):
        u, v = f(x[0], x[1])
        return [u - y[0], v - y[1]]
    
    solution = fsolve(equations, x0)
    return solution[0], solution[1]

# 测试：非线性剪切的逆
alpha = 0.5
original_point = (0.5, 0.5)
transformed = nonlin_shear(alpha)(*original_point)
# 创建逆函数
inverse_nonlin_shear = lambda alpha: lambda x, y: inverse_numerical(nonlin_shear(alpha), (x, y))
recovered = inverse_nonlin_shear(alpha)(*transformed)
# 等价于 recovered = inverse_numerical(nonlin_shear(alpha), transformed)

print(f'原始点: {original_point}')
print(f'变换后: {transformed}')
print(f'恢复的点: {recovered}')

原始点: (0.5, 0.5)
变换后: (0.5, 0.625)
恢复的点: (np.float64(0.5), np.float64(0.5))

图像变换的完整流程¶

下面实现完整的图像变换（反向生成）：

In [9]:
def transform_image(img, transform_func, pixels=400): 
    """
    对图像应用变换
    
    参数:
    - img: 输入图像 (NumPy 数组)
    - transform_func: 变换函数 (x, y) -> (new_x, new_y)
    - pixels: 输出图像的像素数
    """
    rows, cols = img.shape[:2]
    m = max(cols, rows)
    
    # 创建输出图像
    output = np.zeros((pixels, pixels, img.shape[2]), dtype=img.dtype)
    
    for i in range(pixels):
        for j in range(pixels):
            # 将像素坐标转换为 [-1, 1] 范围
            x, y = scale(2/pixels)(*flipy(*swap(*translate(-pixels/2, -pixels/2)(i, j))))
            
            # 应用变换的逆（因为我们是从输出像素找输入像素）
            new_x, new_y = transform_func(x, y)
            
            # 将坐标转换回像素索引
            src_i, src_j = translate(rows/2, cols/2)(*swap(*flipy(*scale(m/2)(new_x, new_y))))
            src_i = int(src_i)
            src_j = int(src_j)
            
            if 0 <= src_i < rows and 0 <= src_j < cols:
                output[i, j] = img[src_i, src_j]
            else:
                output[i, j] = black  # 黑色背景
    
    return output
In [10]:
# 示例：应用旋转变换
theta = np.radians(30)
rotated_img = transform_image(img, rotate(-theta))

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(img)
axes[0].set_title('原始图像')
axes[0].axis('off')
axes[1].imshow(rotated_img)
axes[1].set_title(f'旋转 {np.degrees(theta):.1f}°')
axes[1].axis('off')
Out[10]:
(np.float64(-0.5), np.float64(399.5), np.float64(399.5), np.float64(-0.5))
No description has been provided for this image
In [11]:
# 示例：应用剪切变换
alpha = 0.5
sheared_img = transform_image(img, shear(-alpha))

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(img)
axes[0].set_title('原始图像')
axes[0].axis('off')
axes[1].imshow(sheared_img)
axes[1].set_title(f'剪切 (α={alpha})')
axes[1].axis('off')
Out[11]:
(np.float64(-0.5), np.float64(399.5), np.float64(399.5), np.float64(-0.5))
No description has been provided for this image
In [12]:
# 示例：应用扭曲变换
alpha = 1.0
warped_img = transform_image(img, warp(-alpha))

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(img)
axes[0].set_title('原始图像')
axes[0].axis('off')
axes[1].imshow(warped_img)
axes[1].set_title(f'扭曲变换 (α={alpha})')
axes[1].axis('off')
Out[12]:
(np.float64(-0.5), np.float64(399.5), np.float64(399.5), np.float64(-0.5))
No description has been provided for this image
In [13]:
# 示例：应用非线性剪切
alpha = 0.3
nonlin_sheared_img = transform_image(img, nonlin_shear(-alpha))

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(img)
axes[0].set_title('原始图像')
axes[0].axis('off')
axes[1].imshow(nonlin_sheared_img)
axes[1].set_title(f'非线性剪切 (α={alpha})')
axes[1].axis('off')
Out[13]:
(np.float64(-0.5), np.float64(399.5), np.float64(399.5), np.float64(-0.5))
No description has been provided for this image

正向处理的问题¶

前面我们用的方法是反向变换：对于输出图像的每个像素，计算它在原图中对应的位置。这样可以确保输出图像的每个像素都有值。

如果用正向变换：对于原图的每个像素，计算它变换后的位置。这种方法可能会产生问题：

  • 间隙：输出图像的某些像素可能没有被覆盖
  • 碰撞：多个原图像素可能映射到同一个输出像素

为了让正向变换和逆向变换可比，我们使用：

  • 正向变换：使用变换 $T$，将原图坐标映射到输出坐标
  • 逆向变换：使用 $T$ 的逆变换 $T^{-1}$，从输出坐标追溯到原图坐标
In [14]:
# 定义变换的逆变换

def inverse_rotate(theta):
    """旋转的逆：旋转 -theta"""
    return rotate(-theta)

def inverse_scale(alpha):
    """缩放的逆：缩放 1/alpha"""
    return scale(1.0 / alpha)

def inverse_shear(alpha):
    """剪切的逆：剪切 -alpha"""
    return shear(-alpha)

def inverse_translate(alpha, beta):
    """平移的逆：平移 (-alpha, -beta)"""
    return translate(-alpha, -beta)

def inverse_nonlin_shear(alpha):
    """非线性剪切的逆：需要数值求解"""
    def inv_f(x, y):
        # 原变换: (x, y) -> (x, y + alpha * x^2)
        # 逆变换: 已知 (X, Y)，求 (x, y)
        # X = x, Y = y + alpha * x^2
        # 所以 x = X, y = Y - alpha * X^2
        return np.array([x, y - alpha * x**2])
    return inv_f
In [15]:
# 正向变换函数
def transform_image_forward(img, transform_func, pixels=400):
    """
    使用正向变换变换图像
    transform_func: 变换函数 T，将原图坐标映射到输出坐标
    注意：正向变换可能产生间隙和碰撞问题
    """
    rows, cols = img.shape[:2]
    m = max(rows, cols)
    
    # 创建输出图像，初始化为黑色
    result = np.zeros((pixels, pixels, img.shape[2]), dtype=img.dtype)
    
    # 创建计数矩阵，用于处理碰撞
    count = np.zeros((pixels, pixels))
    
    # 遍历原图的每个像素
    for i_src in range(rows):
        for j_src in range(cols):
            # 将原图像素索引转换为 xy 坐标
            x_src, y_src = scale(2/m)(*flipy(*swap(*translate(-rows/2, -cols/2)(i_src, j_src))))
            
            # 应用正向变换 T
            x_dst, y_dst = transform_func(x_src, y_src)
            
            # 将变换后的 xy 坐标转换为输出图像的像素索引
            i_dst, j_dst = translate(pixels/2, pixels/2)(*swap(*flipy(*scale(pixels/2)(x_dst, y_dst))))
            i_dst, j_dst = int(i_dst), int(j_dst)
            
            # 检查边界
            if 0 <= i_dst < pixels and 0 <= j_dst < pixels:
                result[i_dst, j_dst] = img[i_src, j_src]
                count[i_dst, j_dst] += 1
    
    return result, count
In [16]:
# 比较正向变换和逆向变换
def compare_forward_backward(img, T, T_inv, T_name, pixels=400):
    """
    比较正向变换和逆向变换的效果
    T: 正向变换函数
    T_inv: T 的逆变换函数（用于逆向变换方法）
    """
    # 逆向变换：使用 T 的逆
    img_backward = transform_image(img, T_inv, pixels)
    
    # 正向变换：使用 T
    img_forward, count = transform_image_forward(img, T, pixels)
    
    # 绘制结果
    fig, axes = plt.subplots(1, 4, figsize=(16, 4))
    
    axes[0].imshow(img)
    axes[0].set_title('原始图像')
    axes[0].axis('off')
    
    axes[1].imshow(img_backward)
    axes[1].set_title(f'逆向变换 ($T^{{-1}}$)')
    axes[1].axis('off')
    
    axes[2].imshow(img_forward)
    axes[2].set_title(f'正向变换 ($T$)')
    axes[2].axis('off')
    
    # 显示覆盖次数（碰撞情况）
    axes[3].imshow(count, cmap='hot')
    axes[3].set_title('像素覆盖次数（正向）')
    axes[3].axis('off')
    
    plt.suptitle(f'变换: {T_name}', fontsize=14, y=1.02)
    plt.tight_layout()
    plt.show()
    
    # 统计信息
    print(f"正向变换统计:")
    print(f"  未被覆盖的像素: {np.sum(count == 0)} / {pixels*pixels} ({100*np.sum(count == 0)/(pixels*pixels):.1f}%)")
    print(f"  被覆盖1次的像素: {np.sum(count == 1)}")
    print(f"  被覆盖多次的像素: {np.sum(count > 1)}")
In [17]:
# 测试几种变换的正向和逆向效果

# 1. 旋转
theta = np.radians(30)  # 30度
compare_forward_backward(img, rotate(theta), inverse_rotate(theta), 'rotate(30°)')
No description has been provided for this image 
正向变换统计:
  未被覆盖的像素: 56917 / 160000 (35.6%)
  被覆盖1次的像素: 178
  被覆盖多次的像素: 102905
In [18]:
# 2. 缩放
alpha_scale = 4
compare_forward_backward(img, scale(alpha_scale), inverse_scale(alpha_scale), 'scale(4)')
No description has been provided for this image 
正向变换统计:
  未被覆盖的像素: 102400 / 160000 (64.0%)
  被覆盖1次的像素: 57600
  被覆盖多次的像素: 0
In [19]:
# 3. 组合变换：旋转 + 缩放

# 正向变换 T = scale(4) ∘ rotate(30°)
def compose_forward(x, y):
    r = rotate(np.pi/6)
    s = scale(4)
    return s(*r(x, y))

# 逆变换 T^{-1} = rotate(-30°) ∘ scale(0.25)
def compose_inverse(x, y):
    s_inv = scale(0.25)
    r_inv = rotate(-np.pi/6)
    return r_inv(*s_inv(x, y))

compare_forward_backward(img, compose_forward, compose_inverse, 'scale(4)·rotate(30°)')
No description has been provided for this image 
正向变换统计:
  未被覆盖的像素: 102135 / 160000 (63.8%)
  被覆盖1次的像素: 57840
  被覆盖多次的像素: 25

正向变换和反向变换的比较¶

特性 反向变换 正向变换
实现方式 从输出像素追溯到原图 从原图像素映射到输出
间隙问题 无间隙 存在间隙（黑色像素）
碰撞问题 无碰撞 多个像素可能映射到同一位置
图像质量 高，边缘平滑 较低，需要额外处理
计算复杂度 需要求逆变换 直接应用变换
适用场景 图像处理的主流方法 特定场景（如粒子系统）

结论：对于图像变换，反向变换往往是更好的选择，因为它保证了输出图像的每个像素都有确定的值，避免了间隙和碰撞问题。

本讲小结¶

  • 坐标变换：从输出图像反向映射到输入图像
  • 逆变换：撤销变换效果的变换
    • 线性变换的逆：矩阵逆
    • 非线性变换的逆：数值方法