图像基础

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

• 数字图像是一个二维数字网格。网格中的每个单元格是一个像素（图像元素），其值表示强度或颜色。灰度图像是一个单一的二维矩阵，其中每个像素包含一个亮度值，对于 8 位图像，通常范围从 0（黑色）到 255（白色）。

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

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

• 位深度决定每个通道可以表示的离散强度级别数量。

• 8 位图像每个通道有 \(2^8 = 256\) 个级别，总共 \(256^3 \approx 1670\) 万种可能的颜色。16 位图像每个通道有 65,536 个级别，用于医学成像和高动态范围摄影等对精细强度差异敏感的场景。

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

• HSV（色调、饱和度、明度）将颜色信息与亮度分离。色调是纯色（在色环上 0-360 度），饱和度是颜色的鲜艳程度（0 = 灰色，1 = 纯色），明度是亮度。HSV 适合基于颜色的分割，因为你可以仅根据色调设定阈值，而无需考虑光照条件。在 HSV 中检测"红色物体"比在 RGB 中容易得多。

• YCbCr 将亮度（Y，感知亮度）与色度（Cb、Cr，颜色差异信号）分离。这是 JPEG 压缩和视频编解码器中使用的色彩空间。人眼对亮度比对颜色更敏感，因此色度可以以较低分辨率存储（色度子采样）而几乎不产生感知损失。

• LAB（CIELAB）的设计目标是使两种颜色之间的数值距离对应于感知差异。在 LAB 空间中相等的步长对人眼观察者来说看起来也是相等的。L 通道是明度，A 从绿色到红色，B 从蓝色到黄色。当需要感知均匀的颜色比较时，使用 LAB。

• 图像形成描述三维场景如何变成二维图像。最简单的模型是针孔相机：来自场景的光线通过一个小孔投射到其后的传感器平面上。世界坐标系中的点 \((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 矩阵是内参矩阵 \(K\)。它编码了相机的内部属性：焦距 \(f_x, f_y\)（透镜会聚光线的强度）和主点 \((c_x, c_y)\)（光轴与传感器的交点，通常靠近图像中心）。对于给定的相机和镜头组合，这些参数是固定的。

• 外参描述相机在世界中的位置：一个旋转矩阵 \(R\)（3x3，来自第 02 章）和一个平移向量 \(t\)（3x1）。它们共同将世界坐标转换为相机坐标。完整的投影是：

\[\mathbf{p} = K [R \mid t] \mathbf{P}\]

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

• 真实镜头会引入畸变。

  • 径向畸变使直线弯曲成曲线（桶形畸变使图像向外凸出；枕形畸变使其向内收缩）。 切向畸变源于镜头未与传感器完全平行。

• 相机标定通过拍摄已知图案（如棋盘格）的图像来估计内参和畸变系数，然后校正（去畸变）图像。

• 空间滤波是经典图像处理的基础。一个滤波器（或卷积核）是一个小矩阵（通常为 3x3 或 5x5），它在图像上滑动。在每个位置，滤波器的值与重叠的图像块逐元素相乘并求和，产生一个输出像素。这就是二维卷积，与驱动 CNN（文件 02）的运算相同，但这里的滤波器权重是手工设计而非学习得到的。

\[(\text{图像} * K)[i,j] = \sum_{m} \sum_{n} \text{图像}[i+m, j+n] \cdot K[m, n]\]

• 这是第 06 章中一维卷积的二维扩展。滤波器决定了该运算检测的内容：不同的滤波器检测不同的特征。

• 模糊通过对相邻像素取平均来平滑图像。盒式滤波器对所有相邻像素赋予相同的权重。

• 高斯滤波器通过二维高斯函数（第 05 章）对相邻像素加权，给相邻像素更大的权重，给远处的像素更小的权重。高斯模糊是最常见的平滑操作，由 \(\sigma\) 参数化：\(\sigma\) 越大，平滑程度越高。

• 中值滤波用邻域的中值代替每个像素，而非加权平均。它在去除椒盐噪声（随机的黑白像素）方面特别有效，同时保留边缘，因为中值对异常值具有鲁棒性（如第 04 章所讨论的）。

• 边缘检测识别像素强度急剧变化的边界。边缘承载了图像中的大部分结构信息；仅凭边缘就可以识别物体。

• 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)\) 共同描述每个像素处的边缘强度和方向。这是第 03 章中梯度在图像域的对应概念。

• Canny 边缘检测器是边缘检测的黄金标准。它包含四个步骤：

  1. 使用高斯滤波器平滑图像以减少噪声
  2. 计算梯度幅值和方向（使用 Sobel）
  3. 非极大值抑制：仅保留沿梯度方向为局部最大值的像素，细化边缘
  4. 滞后阈值处理：使用两个阈值（高阈值和低阈值）。高于高阈值的像素是确定边缘。介于两个阈值之间的像素仅当连接到确定边缘时才被视为边缘。低于低阈值的像素被舍弃。

• Canny 中的双阈值使其比单阈值更鲁棒：强边缘始终被保留，弱边缘仅当属于连续边缘结构时才被保留。

• 频域分析揭示了在空间域难以看到的模式。二维傅里叶变换（扩展自第 03 章的一维版本）将图像分解为不同频率和方向的正弦模式之和：

\[F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \cdot e^{-j2\pi(ux/M + vy/N)}\]

• 低频对应平滑、缓慢变化的区域（天空、墙壁）。高频对应锐利变化（边缘、纹理、噪声）。幅度谱显示每个频率上存在多少能量，相位谱编码了空间排列信息。

• 低通滤波去除高频，从而平滑图像（相当于空间域的高斯模糊）。高通滤波去除低频，从而强调边缘和细节。带通滤波只保留一定范围的频率，用于纹理分析。

• 在实践中，对于大尺寸滤波器，频域滤波可能比空间卷积更快，因为空间域中的卷积等价于频域中的逐元素乘法（卷积定理）。这直接联系到第 03 章中的傅里叶变换性质。

• 直方图总结像素强度的分布。直方图统计每个强度值有多少像素（对于 8 位图像为 0-255）。这是第 04 章中的频率分布应用于像素值。

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

• 直方图均衡化将直方图拉伸以覆盖整个强度范围，从而改善对比度。其思路是找到一个映射，使像素强度的累积分布函数（CDF）近似为线性。这是第 04 章中 CDF 概念的直接应用。

• Otsu 方法自动找到将图像分割为前景和背景的最佳阈值。它尝试每个可能的阈值，并选择使类内方差最小（或等价地，使类间方差最大）的阈值。这是第 04 章中方差概念应用于像素强度群体的体现。

• 特征提取识别图像中可用于匹配、识别和三维重建的独特点或区域。好的特征应具有可重复性（在不同视角下能被再次找到）、独特性（可与其他特征区分）和计算高效性。

• 角点检测寻找图像强度在多个方向上显著变化的点。平滑区域在任何方向上的变化都很小。边缘在一个方向上有变化。角点在至少两个方向上都有变化，使其在局部是唯一的，因此是可靠的标志点。

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

\[ 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\) 的特征值（来自第 02 章）告诉你特征的类型：

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

• Harris 不显式计算特征值，而是使用角点响应函数：\(R = \det(M) - k \cdot (\text{tr}(M))^2\)，其中 \(\det(M) = \lambda_1 \lambda_2\) 且 \(\text{tr}(M) = \lambda_1 + \lambda_2\)（均来自第 02 章）。\(R\) 为正且较大时表示角点。常数 \(k\) 通常为 0.04-0.06。

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

• 斑点检测寻找与周围环境不同的区域。与角点（属于点特征）不同，斑点具有特征尺寸。

• SIFT（尺度不变特征变换，Lowe，2004）在多个尺度上检测斑点，并构建对旋转、尺度具有不变性，对光照变化具有部分不变性的描述子。它的工作原理是：

  1. 使用逐渐增大 \(\sigma\) 的高斯模糊构建尺度空间（见下文）
  2. 在尺度间的 Gaussian 差分（DoG）中寻找极值点
  3. 精炼关键点位置，去除低对比度点和边缘响应
  4. 基于局部梯度方向分配主方向
  5. 从关键点周围 16x16 块中的梯度直方图构建 128 维描述子

• SURF（加速稳健特征）使用盒式滤波器和积分图像近似 SIFT 以实现更快的计算。ORB（定向 FAST 和旋转 BRIEF）是一个快速、开源的替代方案，它将 FAST 角点检测器与 BRIEF 二进制描述子结合，并增加了旋转不变性。

• HOG（方向梯度直方图）描述子将图像划分为小单元格，计算每个单元格内梯度方向的直方图，并在单元格块间进行归一化。HOG 捕捉边缘方向的分布，这对物体形状具有高度信息量。在深度学习之前，HOG + SVM（第 06 章）是行人检测和物体识别的主流方法。

• 图像金字塔以多种分辨率表示图像。

  • 高斯金字塔通过重复模糊和下采样（分辨率减半）构建。每一层都是原始图像的粗略版本。
  • 拉普拉斯金字塔存储连续高斯层之间的差异，捕捉每一步下采样丢失的细节。拉普拉斯金字塔是可逆的：你可以从中重建原始图像。

• 尺度空间形式化了物体存在于不同尺度这一概念。一棵树是一个大斑点；树上的一片叶子是一个小斑点。要同时检测两者，你需要跨尺度搜索。图像的尺度空间是通过将图像与逐渐增大 \(\sigma\) 的高斯函数卷积得到的图像族：

\[L(x, y, \sigma) = G(x, y, \sigma) * I(x, y)\]

• 其中 \(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
   
   # Create a synthetic test image with distinct colours
   H, W = 128, 256
   img = np.zeros((H, W, 3), dtype=np.uint8)
   img[:, :64] = [255, 50, 50]     # red
   img[:, 64:128] = [50, 255, 50]  # green
   img[:, 128:192] = [50, 50, 255] # blue
   img[:, 192:] = [255, 255, 50]   # yellow
   
   # Add a brightness gradient
   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
   
   # Manual RGB to HSV conversion
   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
   
       # Hue
       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):
       """2D convolution (valid mode) from scratch."""
       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
   
   # Create a test image: white rectangle on dark background
   img = jnp.zeros((64, 64))
   img = img.at[15:50, 20:45].set(1.0)
   # Add some noise
   key = jax.random.PRNGKey(42)
   img = img + jax.random.normal(key, img.shape) * 0.05
   
   # Sobel filters
   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)
   
   # Gaussian blur kernel (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()
   
   # Apply filters
   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
   
   # Create a low-contrast image (values clustered in a narrow range)
   key = __import__('jax').random.PRNGKey(42)
   img = __import__('jax').random.uniform(key, (128, 128)) * 0.3 + 0.3  # values in [0.3, 0.6]
   
   def histogram_equalise(img, n_bins=256):
       """Histogram equalisation for a grayscale image."""
       # Quantise to bins
       bins = jnp.linspace(0, 1, n_bins + 1)
       hist = jnp.histogram(img, bins=bins)[0]
   
       # Compute CDF
       cdf = jnp.cumsum(hist)
       cdf_normalised = (cdf - cdf.min()) / (cdf.max() - cdf.min())
   
       # Map each pixel through the 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 corner detection from scratch."""
       # Compute gradients with 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)
   
       # Pad image for valid convolution to preserve size
       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))
   
       # Structure tensor components
       Ixx = Ix * Ix
       Iyy = Iy * Iy
       Ixy = Ix * Iy
   
       # Gaussian smoothing of structure tensor (approximate with window sum)
       w = 3  # window half-size
       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)
   
       # Threshold
       corners = R > threshold * R.max()
       return R, corners
   
   # Test image: checkerboard pattern (lots of 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()