图像基础

图像基础

图像基础（image fundamentals）解释了在任何模型看到图像之前，数字图像是如何被表示、形成和预处理的。本文件涵盖像素、色彩空间（RGB、HSV、YCbCr、LAB）、针孔相机模型、卷积、边缘检测（Sobel、Canny）、直方图以及特征描述子（SIFT、ORB），这是低层视觉的工具箱。

• 数字图像（digital image）是一个二维数字网格。网格中的每个单元格叫做一个像素（pixel，picture element），其值代表亮度或颜色。灰度图是一张二维矩阵，每个像素存一个亮度值，对 8 位图像而言通常从 0（黑）到 255（白）。

• 彩色图像把它扩展到三个通道。在 RGB 色彩空间中，每个像素存储三个值：红、绿、蓝的强度。

• 彩色图像是一个形状为 (高度, 宽度, 3) 的三维张量（矩阵）。把这三个通道以不同强度混合，就能产生可见光谱中的全部颜色。

• **位深（bit depth）**决定每个通道能表示多少个不同的强度级别。

• 8 位图像每个通道有 2^8 = 256 个级别，共可表示 256^3 \approx 1670 万种颜色。16 位图像每个通道有 65,536 个级别，用于医疗影像和 HDR 摄影这类对细微强度差别非常敏感的场景。

• RGB 对显示器很方便，但其他色彩空间更适合不同的任务。

• HSV（Hue Saturation Value，色相、饱和度、明度）把颜色信息和亮度分开。色相是纯色（在色环上 0-360 度），饱和度是颜色的鲜艳程度（0 = 灰色，1 = 纯色），明度就是亮度。HSV 非常适合做基于颜色的分割，因为你可以不管光照条件，只对色相做阈值。检测"红色物体"在 HSV 里比在 RGB 里容易得多。

• YCbCr 把亮度（Y，感知到的明度）和色度（Cb、Cr，色差信号）分开。这是 JPEG 压缩和视频编码使用的色彩空间。人眼对亮度比对颜色更敏感，所以色度可以用更低的分辨率来存（色度子采样），几乎察觉不到损失。

• LAB（CIELAB）被设计成这样：两种颜色之间的数值距离对应于人感知上的差异。在 LAB 空间里等距的步长，在人眼看来也是等距的。L 通道是亮度，A 从绿到红，B 从蓝到黄。当你需要"感知均匀"的颜色比较时，就用 LAB。

• 成像（image formation）描述三维场景如何变成二维图像。最简单的模型是针孔相机（pinhole camera）：来自场景的光线穿过一个小孔，投射到它后面的感光面上。世界坐标里的一个点 (X, Y, Z) 投影到像素坐标 (u, v)：

\begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = \frac{1}{Z} \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} X \\ Y \\ Z \end{bmatrix}

• 这个 3x3 矩阵就是内参矩阵（intrinsic matrix） K。它编码了相机的内部属性：焦距 f_x, f_y（镜头汇聚光线的能力有多强）和主点（principal point）(c_x, c_y)（光轴与感光面相交的位置，通常接近图像中心）。对于给定的相机和镜头组合，这些值是固定的。

• **外参（extrinsic parameters）**描述相机在世界中的位置：一个旋转矩阵 R（3x3，来自第 2 章）和一个平移向量 t（3x1）。它们一起把世界坐标变换为相机坐标。完整的投影是：

• 其中 \mathbf{P} = [X, Y, Z, 1]^T 是齐次坐标下的三维点，\mathbf{p} = [u, v, 1]^T 是投影后的像素。[R \mid t] 矩阵是 3x4，把旋转和平移并排放在一起。这些全都是第 2 章里的线性代数。

• 真实的镜头会引入畸变（distortion）。

  • **径向畸变（radial distortion）**把直线弯成曲线（桶形畸变让图像向外鼓；枕形畸变把它向内挤）。
    **切向畸变（tangential distortion）**则是因为镜头没有完美地平行于感光面。

• 相机标定（camera calibration）从已知图案（比如棋盘格）的图像中同时估计内参和畸变系数，然后对图像做校正（去畸变）。

• **空间滤波（spatial filtering）**是经典图像处理的基础。滤波器（filter，也叫核 kernel）是一个小矩阵（通常是 3x3 或 5x5），它在图像上滑动。在每个位置，把滤波器的值和重叠的图像小块逐元素相乘再求和，就得到一个输出像素。这就是二维卷积（2D convolution），和驱动 CNN（见文件 02）的是同一个运算，只不过这里的滤波器权重是手工设计的，而不是学出来的。

• 这是第 6 章一维卷积的二维推广。滤波器决定了运算能检测到什么：不同的滤波器检测不同的特征。

• **模糊（blurring）**通过把相邻像素取平均来平滑图像。**盒式滤波器（box filter）**给所有邻居相同的权重。

• **高斯滤波器（Gaussian filter）**按二维高斯（第 5 章）给邻居加权，离得近的像素权重大、离得远的权重小。高斯模糊是最常见的平滑操作，由参数 \sigma 控制：\sigma 越大，平滑得越厉害。

• **中值滤波（median filtering）**用邻域的中值替换每个像素，而不是加权平均。它在去除椒盐噪声（随机的黑白像素）的同时还能保留边缘，因为中值对离群值很鲁棒（正如第 4 章所讨论的）。

• **边缘检测（edge detection）**找出像素强度剧烈变化的边界。边缘承载了图像中大部分的结构信息；光凭边缘你就能认出物体。

• Sobel 算子用两个 3x3 滤波器来估计水平和垂直方向的梯度：

G_x = \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix}, \quad G_y = \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix}

• 用 G_x 对图像卷积得到水平梯度（对垂直边缘响应强），G_y 给出垂直梯度（对水平边缘响应强）。

• 梯度幅值 \sqrt{G_x^2 + G_y^2} 和方向 \arctan(G_y / G_x) 一起描述了每个像素处的边缘强度和朝向。这是第 3 章梯度的图像域对应物。

• Canny 边缘检测器是边缘检测的金标准。它分四步：

  1. 用高斯滤波器平滑图像以降噪
  2. 计算梯度幅值和方向（用 Sobel）
  3. 非极大值抑制（non-maximum suppression）：只保留沿梯度方向上的局部最大像素，把边缘变细
  4. 滞后阈值化（hysteresis thresholding）：用两个阈值（高和低）。高于高阈值的像素肯定是边缘；介于两个阈值之间的像素只有在和确定边缘连通时才算边缘；低于低阈值的像素被丢弃。

• Canny 里的两个阈值让它比单一阈值更鲁棒：强边缘总是被保留，而弱边缘只有在属于连续边缘结构时才被保留。

• **频域（frequency domain）**分析能揭示空间域里很难看出的模式。二维傅里叶变换（2D Fourier transform）（第 3 章一维版本的推广）把图像分解成不同频率和朝向的二维正弦图案之和：

• 低频对应平滑、缓慢变化的区域（天空、墙面）。高频对应剧烈的过渡（边缘、纹理、噪声）。**幅值谱（magnitude spectrum）**显示每个频率上有多少能量，**相位谱（phase spectrum）**则编码空间排列。

• **低通滤波（low-pass filtering）**去掉高频，使图像变平滑（空间域上等价于高斯模糊）。**高通滤波（high-pass filtering）**去掉低频，强调边缘和细节。**带通滤波（band-pass filtering）**只保留一段频率，对纹理分析很有用。

• 在实践中，当滤波器很大时，在频域里做滤波可能比空间卷积更快，因为空间域的卷积等价于频域的逐元素相乘（卷积定理，convolution theorem）。这直接呼应了第 3 章的傅里叶变换性质。

• **直方图（histogram）**汇总像素强度的分布。直方图统计有多少像素取每个强度值（8 位图像为 0-255）。它就是把第 4 章的频率分布用在了像素值上。

• 偏暗的图像直方图集中在左边（低值）。偏亮的图像集中在右边。低对比度图像的直方图很窄。高对比度图像的直方图宽而分散。

• **直方图均衡化（histogram equalisation）**把直方图拉伸到覆盖整个强度范围，从而提升对比度。思路是找一个映射，让像素强度的累积分布函数（CDF）近似线性。这是第 4 章 CDF 概念的直接应用。

• Otsu 方法自动找出把图像分成前景和背景的最佳阈值。它尝试所有可能的阈值，挑出使类内方差最小的那个（等价地，使类间方差最大）。这还是第 4 章的方差概念，只是用在了像素强度总体上。

• **特征提取（feature extraction）**识别图像中可用于匹配、识别和三维重建的独特点或区域。好的特征应当可重复（在另一个视角下还能再找到）、有区分度（能和其他特征区分开），而且计算高效。

• **角点检测（corner detection）**找的是在多个方向上图像强度都显著变化的点。平滑区域在任何方向上变化都很小。边缘在一个方向上有变化。角点在至少两个方向上有变化，因此它在局部是独一无二的，是可靠的路标。

• Harris 角点检测器分析每个像素处的结构张量（structure tensor，也叫二阶矩矩阵）：

M = \sum_{(x,y) \in W} w(x,y) \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix}

• 其中 I_x 和 I_y 是图像梯度（用 Sobel 计算），W 是局部窗口，w 是高斯加权函数。M 的特征值（来自第 2 章）告诉你特征的类型：

  • 两个特征值都小：平坦区域（没有特征）
  • 一个大一个小：边缘
  • 两个都大：角点

• Harris 不直接算特征值，而是用一个角点响应函数：R = \det(M) - k \cdot (\text{trace}(M))^2，其中 \det(M) = \lambda_1 \lambda_2、\text{trace}(M) = \lambda_1 + \lambda_2（都来自第 2 章）。R 是个大正数就表明是角点。常数 k 通常取 0.04-0.06。

• Shi-Tomasi 检测器把它简化为 R = \min(\lambda_1, \lambda_2)，直接检查较小的特征值是否足够大。它在实践中稍微稳定一些。

• **斑点检测（blob detection）**找的是和周围环境不同的区域。和角点（点特征）不同，斑点有特征性的大小。

• SIFT（Scale-Invariant Feature Transform，尺度不变特征变换，Lowe，2004）在多个尺度上检测斑点，并构造一个对旋转、尺度不变，对光照变化也部分不变的描述子。它的流程是：

  1. 用逐级增大的 \sigma 做高斯模糊，构建尺度空间（scale space）（见下文）
  2. 在不同尺度的差分高斯（DoG）中找极值
  3. 精炼关键点位置，去除低对比度点和边缘响应
  4. 根据局部梯度方向赋予一个主方向
  5. 在关键点周围的 16x16 小块里用梯度直方图构建一个 128 维描述子

• SURF（Speeded-Up Robust Features）用盒式滤波器和积分图像来近似 SIFT，计算更快。ORB（Oriented FAST and Rotated BRIEF）是一个快速、开源的替代方案，它把 FAST 角点检测器和 BRIEF 二值描述子结合起来，并加上了旋转不变性。

• HOG（Histogram of Orientated Gradients，方向梯度直方图）描述子把图像分成小单元格，在每个单元格内计算梯度方向的直方图，并在单元格块之间做归一化。HOG 捕捉的是边缘朝向的分布，这对物体形状的信息量很大。在深度学习出现之前，HOG + SVM（第 6 章）是行人检测和物体识别的主流方法。

• **图像金字塔（image pyramids）**用多个分辨率来表示一张图像。

  • **高斯金字塔（Gaussian pyramid）**通过反复模糊和下采样（把分辨率减半）来构建。每一层都是原图的更粗糙版本。
  • **拉普拉斯金字塔（Laplacian pyramid）**存储相邻高斯层之间的差值，捕捉每次下采样时丢失的细节。拉普拉斯金字塔是可逆的：你可以从它重建原图。

• **尺度空间（scale space）**把"物体存在于不同尺度"这一想法形式化。一棵树是一个大斑点；树上的一片叶子是一个小斑点。要同时检测到两者，你需要在多个尺度上搜索。图像的尺度空间是用逐级增大的 \sigma 做高斯卷积得到的一族图像：

• 其中 G 是标准差为 \sigma 的二维高斯。在多个尺度上都能持续存在的特征，更可能是真实的结构而不是噪声。尺度空间是 SIFT 的理论基础，也是整个现代计算机视觉里多尺度处理（包括目标检测中的特征金字塔网络，文件 03）的基础。

编程练习（使用 CoLab 或 notebook）

1. 加载一张图像，把它转换到不同的色彩空间（RGB、HSV、LAB），并可视化各个通道。观察颜色信息在不同空间里是如何分布的。

import jax.numpy as jnp
import matplotlib.pyplot as plt
from PIL import Image
import numpy as np

# 构造一张颜色分明的合成测试图
H, W = 128, 256
img = np.zeros((H, W, 3), dtype=np.uint8)
img[:, :64] = [255, 50, 50]     # 红色
img[:, 64:128] = [50, 255, 50]  # 绿色
img[:, 128:192] = [50, 50, 255] # 蓝色
img[:, 192:] = [255, 255, 50]   # 黄色

# 加一个亮度渐变
for y in range(H):
    scale = 0.3 + 0.7 * y / H
    img[y] = (img[y] * scale).astype(np.uint8)

img_jnp = jnp.array(img, dtype=jnp.float32) / 255.0

# 手动实现 RGB 到 HSV 的转换
def rgb_to_hsv(rgb):
    r, g, b = rgb[..., 0], rgb[..., 1], rgb[..., 2]
    maxc = jnp.max(rgb, axis=-1)
    minc = jnp.min(rgb, axis=-1)
    diff = maxc - minc + 1e-7

    # 色相
    h = jnp.where(maxc == minc, 0.0,
        jnp.where(maxc == r, 60 * ((g - b) / diff % 6),
        jnp.where(maxc == g, 60 * ((b - r) / diff + 2),
                              60 * ((r - g) / diff + 4))))
    s = jnp.where(maxc < 1e-7, 0.0, diff / maxc)
    v = maxc
    return jnp.stack([h / 360, s, v], axis=-1)

hsv = rgb_to_hsv(img_jnp)

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
for i, (ch, name) in enumerate(zip([img_jnp[...,0], img_jnp[...,1], img_jnp[...,2]],
                                     ['Red', 'Green', 'Blue'])):
    axes[0, i].imshow(ch, cmap='gray', vmin=0, vmax=1)
    axes[0, i].set_title(f'RGB: {name}'); axes[0, i].axis('off')

for i, (ch, name) in enumerate(zip([hsv[...,0], hsv[...,1], hsv[...,2]],
                                     ['Hue', 'Saturation', 'Value'])):
    axes[1, i].imshow(ch, cmap='gray', vmin=0, vmax=1)
    axes[1, i].set_title(f'HSV: {name}'); axes[1, i].axis('off')

plt.suptitle('RGB vs HSV Channels')
plt.tight_layout(); plt.show()
​

2. 从零开始用二维卷积实现 Sobel 边缘检测和高斯模糊。把它们应用到一张图像上并比较结果。

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

def conv2d(image, kernel):
    """从零实现的二维卷积（valid 模式）。"""
    H, W = image.shape
    kH, kW = kernel.shape
    out_h, out_w = H - kH + 1, W - kW + 1
    output = jnp.zeros((out_h, out_w))
    for i in range(out_h):
        for j in range(out_w):
            patch = image[i:i+kH, j:j+kW]
            output = output.at[i, j].set(jnp.sum(patch * kernel))
    return output

# 测试图：暗背景上的白色矩形
img = jnp.zeros((64, 64))
img = img.at[15:50, 20:45].set(1.0)
# 加点噪声
key = jax.random.PRNGKey(42)
img = img + jax.random.normal(key, img.shape) * 0.05

# Sobel 滤波器
sobel_x = jnp.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=jnp.float32)
sobel_y = jnp.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=jnp.float32)

# 高斯模糊核（5x5，sigma=1）
ax = jnp.arange(-2, 3, dtype=jnp.float32)
xx, yy = jnp.meshgrid(ax, ax)
gaussian = jnp.exp(-(xx**2 + yy**2) / (2 * 1.0**2))
gaussian = gaussian / gaussian.sum()

# 应用滤波器
gx = conv2d(img, sobel_x)
gy = conv2d(img, sobel_y)
edges = jnp.sqrt(gx**2 + gy**2)
blurred = conv2d(img, gaussian)

fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for ax, data, title in zip(axes,
    [img, edges, blurred, gx],
    ['Original', 'Edge Magnitude', 'Gaussian Blur', 'Horizontal Gradient']):
    ax.imshow(data, cmap='gray')
    ax.set_title(title); ax.axis('off')
plt.tight_layout(); plt.show()
​

3. 从零实现直方图均衡化，并应用到一张低对比度的灰度图上。比较前后的直方图。

import jax.numpy as jnp
import matplotlib.pyplot as plt

# 构造一张低对比度图像（值集中在一个窄区间）
key = __import__('jax').random.PRNGKey(42)
img = __import__('jax').random.uniform(key, (128, 128)) * 0.3 + 0.3  # 值落在 [0.3, 0.6]

def histogram_equalise(img, n_bins=256):
    """灰度图的直方图均衡化。"""
    # 量化到若干个 bin
    bins = jnp.linspace(0, 1, n_bins + 1)
    hist = jnp.histogram(img, bins=bins)[0]

    # 计算 CDF
    cdf = jnp.cumsum(hist)
    cdf_normalised = (cdf - cdf.min()) / (cdf.max() - cdf.min())

    # 通过 CDF 映射每个像素
    indices = jnp.clip((img * n_bins).astype(jnp.int32), 0, n_bins - 1)
    equalised = cdf_normalised[indices]
    return equalised

eq_img = histogram_equalise(img)

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes[0, 0].imshow(img, cmap='gray', vmin=0, vmax=1)
axes[0, 0].set_title('Original (Low Contrast)'); axes[0, 0].axis('off')
axes[0, 1].imshow(eq_img, cmap='gray', vmin=0, vmax=1)
axes[0, 1].set_title('After Histogram Equalisation'); axes[0, 1].axis('off')

axes[1, 0].hist(img.ravel(), bins=64, color='#3498db', alpha=0.8)
axes[1, 0].set_title('Histogram Before'); axes[1, 0].set_xlim(0, 1)
axes[1, 1].hist(eq_img.ravel(), bins=64, color='#e74c3c', alpha=0.8)
axes[1, 1].set_title('Histogram After'); axes[1, 1].set_xlim(0, 1)

plt.tight_layout(); plt.show()
​

4. 从零实现 Harris 角点检测器。在一张简单图像中检测角点并可视化。

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

def harris_corners(img, k=0.05, threshold=0.01):
    """从零实现的 Harris 角点检测。"""
    # 用 Sobel 计算梯度
    sobel_x = jnp.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=jnp.float32)
    sobel_y = jnp.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=jnp.float32)

    # 对图像做 padding，保证 valid 卷积后尺寸不变
    img_pad = jnp.pad(img, 1, mode='edge')
    H, W = img.shape

    Ix = jnp.zeros_like(img)
    Iy = jnp.zeros_like(img)
    for i in range(H):
        for j in range(W):
            patch = img_pad[i:i+3, j:j+3]
            Ix = Ix.at[i, j].set(jnp.sum(patch * sobel_x))
            Iy = Iy.at[i, j].set(jnp.sum(patch * sobel_y))

    # 结构张量的各分量
    Ixx = Ix * Ix
    Iyy = Iy * Iy
    Ixy = Ix * Iy

    # 对结构张量做高斯平滑（这里用窗口求和近似）
    w = 3  # 窗口半径
    R = jnp.zeros_like(img)
    pad_xx = jnp.pad(Ixx, w, mode='constant')
    pad_yy = jnp.pad(Iyy, w, mode='constant')
    pad_xy = jnp.pad(Ixy, w, mode='constant')

    for i in range(H):
        for j in range(W):
            sxx = jnp.sum(pad_xx[i:i+2*w+1, j:j+2*w+1])
            syy = jnp.sum(pad_yy[i:i+2*w+1, j:j+2*w+1])
            sxy = jnp.sum(pad_xy[i:i+2*w+1, j:j+2*w+1])
            det = sxx * syy - sxy * sxy
            trace = sxx + syy
            R = R.at[i, j].set(det - k * trace * trace)

    # 阈值
    corners = R > threshold * R.max()
    return R, corners

# 测试图：棋盘格图案（有大量角点）
block = 16
n = 4
checker = jnp.zeros((block * n, block * n))
for i in range(n):
    for j in range(n):
        if (i + j) % 2 == 0:
            checker = checker.at[i*block:(i+1)*block, j*block:(j+1)*block].set(1.0)

R, corners = harris_corners(checker)
cy, cx = jnp.where(corners)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
axes[0].imshow(checker, cmap='gray')
axes[0].set_title('Checkerboard'); axes[0].axis('off')
axes[1].imshow(R, cmap='hot')
axes[1].set_title('Harris Response'); axes[1].axis('off')
axes[2].imshow(checker, cmap='gray')
axes[2].scatter(cx, cy, c='#e74c3c', s=15, marker='x')
axes[2].set_title(f'Detected Corners ({len(cx)})'); axes[2].axis('off')
plt.tight_layout(); plt.show()
​