局部特征-SURF特征描述子

SURF 特征描述子

小象实战讲义 · 基于深度学习的计算机视觉

在上一节我们学习了经典的 SIFT 特征，它以其优异的尺度、旋转不变性而闻名，但其计算复杂度也相对较高。本节我们将介绍 SIFT 的一个著名改进版本——SURF（Speeded-Up Robust Features）。SURF 在保持 SIFT 大部分鲁棒性的前提下，通过一系列巧妙的简化，显著提升了计算速度，是传统特征描述子发展历程中的一个重要里程碑。学完本节，你将理解 SURF 的核心加速思想，并掌握其构建特征描述子的完整流程。

💡 核心导读

• SURF 的定位：作为 SIFT 的加速改进版，旨在保持鲁棒性的同时大幅提升计算效率。
• 核心加速技术：引入 Haar-like 特征和积分图，将复杂的卷积运算简化为快速的加减运算。
• Hessian 矩阵近似：用盒式滤波器（Box Filter）近似高斯二阶微分，是 SURF 速度提升的关键。
• 构建流程：包括构建 Hessian 矩阵尺度空间、特征点定位、确定主方向以及生成 64 维描述子。
• 性能对比：理解 SURF 相比 SIFT 在速度、亮度与模糊鲁棒性上的优势，以及在尺度和旋转不变性上的不足。

预备知识：Haar-like 特征与积分图

在深入 SURF 之前，我们需要理解其加速运算的两大基石：Haar-like 特征和积分图。

Haar-like 特征

Haar-like 特征是一种用于描述图像局部灰度变化的简单特征。它通过定义一系列包含白色和黑色矩形的特征模板来工作。特征模板的类型包括边缘特征、线性特征、中心特征和对角线特征等。

特征值定义为白色矩形区域内所有像素值的和，减去黑色矩形区域内所有像素值的和。这个差值反映了该区域内的灰度变化情况。

import numpy as np

# 示例：手动计算一个简单的 Haar-like 特征（边缘特征）
# 假设我们有一个 4x4 的图像块
image_patch = np.array([
    [10, 12, 15, 18],
    [11, 13, 16, 19],
    [ 9, 14, 20, 22],
    [ 8, 15, 21, 24]
], dtype=np.float32)

# 定义一个 2x4 的 Haar-like 模板：左半边为白(+1)，右半边为黑(-1)
# 模板覆盖整个图像块的上两行
white_sum = np.sum(image_patch[0:2, 0:2])
black_sum = np.sum(image_patch[0:2, 2:4])
haar_value = white_sum - black_sum

print(f"白色区域(左两列)像素和: {white_sum}")
print(f"黑色区域(右两列)像素和: {black_sum}")
print(f"Haar-like 特征值: {haar_value}")

白色区域(左两列)像素和: 46.0
黑色区域(右两列)像素和: 68.0
Haar-like 特征值: -22.0

积分图（Integral Image）

直接计算 Haar-like 特征需要反复对图像区域求和，效率低下。积分图技术应运而生，它能够实现快速矩形区域求和。

积分图 I(x, y) 定义为从图像左上角 (0,0) 到点 (x,y) 所构成的矩形区域内所有像素值的和： $$I(x, y) = \sum_{x' \le x, y' \le y} i(x', y')$$ 其中 i(x', y') 是原始图像的像素值。

一旦计算出整幅图像的积分图，任意矩形区域 (x1, y1, x2, y2)（其中 (x1,y1) 为左上角，(x2,y2) 为右下角）的像素和可以通过四次加减运算得到： $$\text{sum} = I(x2, y2) - I(x1-1, y2) - I(x2, y1-1) + I(x1-1, y1-1)$$ 注意处理边界情况（当 x1-1 或 y1-1 小于 0 时，对应项为 0）。

import numpy as np

# 计算一个小图像的积分图
def compute_integral_image(img):
    h, w = img.shape
    integral = np.zeros((h+1, w+1), dtype=np.float32) # 多一行一列方便边界处理
    for i in range(1, h+1):
        row_sum = 0.0
        for j in range(1, w+1):
            row_sum += img[i-1, j-1]
            integral[i, j] = integral[i-1, j] + row_sum
    return integral

# 使用积分图快速计算矩形区域和
def sum_region(integral, x1, y1, x2, y2):
    # x1,y1: 矩形左上角 (0-based)
    # x2,y2: 矩形右下角 (0-based)
    # integral 是 (h+1, w+1) 的，索引要加1
    return (integral[y2+1, x2+1] - integral[y1, x2+1] -
            integral[y2+1, x1] + integral[y1, x1])

# 示例图像
small_img = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
], dtype=np.float32)

integral_img = compute_integral_image(small_img)
print("积分图 (带零边界):")
print(integral_img)

# 计算矩形区域 (1,1) 到 (2,2) 的和 (即像素5,6,8,9)
rect_sum = sum_region(integral_img, 1, 1, 2, 2)
print(f"\n矩形区域(1,1)-(2,2)的像素和 (使用积分图): {rect_sum}")
print(f"直接求和验证: {np.sum(small_img[1:3, 1:3])}")

积分图 (带零边界):
[[ 0.  0.  0.  0.]
 [ 0.  1.  3.  6.]
 [ 0.  5. 12. 21.]
 [ 0. 12. 27. 45.]]

矩形区域(1,1)-(2,2)的像素和 (使用积分图): 28.0
直接求和验证: 28.0

有了积分图，无论 Haar-like 模板有多大、多复杂，计算其特征值都变得异常快速，只需几次查表与加减运算。这正是 SURF 能够加速的基础。

SURF 的核心思想与构建流程

SURF 由 Herbert Bay 等人在 2006 年提出，是对 SIFT 的改进，目标是在保持鲁棒性的同时将速度提升数倍。

1. 构建尺度空间与特征点检测

SIFT 使用高斯差分（DoG）来近似高斯拉普拉斯（LoG），并在图像金字塔的不同尺度上寻找极值点。SURF 则采用了不同的策略。

使用 Hessian 矩阵行列式 SURF 的核心是使用 Hessian 矩阵 的行列式来检测斑点状的特征点。对于图像 I(x,y)，在尺度 σ 下的 Hessian 矩阵 H(x, y, σ) 定义为： $$H(x, y, σ) = \begin{bmatrix} L_{xx}(x, y, σ) & L_{xy}(x, y, σ) \\ L_{xy}(x, y, σ) & L_{yy}(x, y, σ) \end{bmatrix}$$ 其中 L_{xx}, L_{xy}, L_{yy} 是图像经过高斯滤波后（尺度为 σ）的二阶偏导数。

Hessian 矩阵的行列式 det(H) 可以用来判断该点是否为极值点（斑点）：

• det(H) > 0：该点是极值点（如果 L_{xx} < 0 则为极大值点，反之为极小值点）。
• det(H) < 0：该点是鞍点。
• det(H) = 0：无法判断。

因此，SURF 通过计算并寻找 det(H) 的极大值点来定位特征点。

盒式滤波器近似 直接计算高斯二阶导数卷积仍然很慢。SURF 的关键创新在于用盒式滤波器（Box Filter） 来近似这些二阶高斯微分模板。

• L_{xx} 被近似为三段竖条模板（左右为正、中间为负，检测水平方向的变化）。
• L_{yy} 被近似为三段横条模板（上下为正、中间为负，检测垂直方向的变化）。
• L_{xy} 被近似为四个象限正负交替的棋盘状模板，对斜向结构敏感。

这些模板的值只有 -2, -1, 0, 1, 2 等简单的整数。更重要的是，对这些模板的卷积操作，可以完全转化为对积分图的加减运算，速度极快。

由于近似会引入误差，实际计算的行列式公式中引入了一个权重系数 0.9： $$det(H_{approx}) = D_{xx}D_{yy} - (0.9 D_{xy})^2$$ 其中 D_{xx}, D_{xy}, D_{yy} 是用盒式滤波器近似计算得到的响应值。

尺度空间的构建 与 SIFT 不断下采样图像不同，SURF 采用了一种更高效的方式构建尺度空间：

1. 图像尺寸保持不变。
2. 逐渐增大盒式滤波器的尺寸。滤波器尺寸与尺度 σ 成正比。
3. 将不同尺寸滤波器得到的响应图堆叠起来，形成尺度空间。
4. 在尺度空间（三维空间：x, y, 滤波器尺寸）中寻找 det(H_approx) 的极值点，作为候选特征点。

2. 特征点定位与筛选

初步得到的极值点还需要进一步精炼：

1. 三维插值：在尺度空间和图像空间中进行插值，以获得亚像素精度的特征点位置和尺度。
2. 去除低对比度点：剔除 |det(H_approx)| 值过小的点，它们对噪声敏感。
3. 边缘响应抑制：SURF 不需要 SIFT 式的迹/行列式比值测试。det(H) 本身就是斑点响应度量，边缘处两个特征值异号会导致 det(H)<0，因此设定 det(H) 阈值即可同时剔除低对比度点与边缘响应点。

3. 确定特征点主方向

为了保证旋转不变性，需要为每个特征点分配一个主方向。SURF 的方法同样利用了 Haar 小波（Haar-like 特征的一种形式）。

1. 以特征点为中心，计算半径为 6s 的圆形邻域（s 是特征点所在的尺度值）。
2. 在这个圆形区域内，用边长为 4s 的 Haar 小波模板（分别计算 x 方向和 y 方向的响应）对每个像素进行计算。Haar 小波模板的响应计算同样得益于积分图，非常快。
3. 用一个滑动扇形窗口（例如 60° 扇形）遍历整个圆形区域。对于每个扇形方向，累加扇形内所有点的 Haar 小波响应向量 (dx, dy)。
4. 将各个扇形方向的响应向量相加，得到该方向的总响应向量。模长最大的那个扇形方向，即为该特征点的主方向。

import numpy as np
import math

# 模拟主方向确定过程（概念演示，非完整实现）
# 假设我们有一个特征点周围的梯度方向（由Haar小波响应模拟）
np.random.seed(42)
num_samples = 100
# 模拟梯度方向：大部分集中在 30 度（pi/6 弧度）附近
mean_angle = math.pi / 6
angles = np.random.normal(loc=mean_angle, scale=0.3, size=num_samples)
angles = np.mod(angles, 2 * math.pi) # 把角度归一化到 [0, 2pi)，而不是直接截断
# 模拟梯度幅值
magnitudes = np.random.uniform(0.5, 2.0, size=num_samples)

# 将梯度分解为x,y分量
dx = magnitudes * np.cos(angles)
dy = magnitudes * np.sin(angles)

# 滑动扇形统计（简化版：分成12个30度扇形）
num_bins = 12
bin_width = 2 * math.pi / num_bins
bin_responses = np.zeros((num_bins, 2)) # 存储每个扇形的(dx_sum, dy_sum)

for dxi, dyi in zip(dx, dy):
    # 计算该梯度向量的角度
    angle = math.atan2(dyi, dxi)
    if angle < 0:
        angle += 2 * math.pi
    # 确定所属扇形
    bin_idx = int(angle // bin_width)
    bin_responses[bin_idx, 0] += dxi
    bin_responses[bin_idx, 1] += dyi

# 计算每个扇形响应向量的模长
bin_magnitudes = np.linalg.norm(bin_responses, axis=1)
# 找到模长最大的扇形
dominant_bin = np.argmax(bin_magnitudes)
dominant_angle = (dominant_bin + 0.5) * bin_width # 取扇形中心角度

print(f"模拟的梯度样本数: {num_samples}")
print(f"每个扇形的响应模长: {bin_magnitudes}")
print(f"主方向所在的扇形索引: {dominant_bin}")
print(f"主方向角度（弧度）: {dominant_angle:.3f}")
print(f"主方向角度（度）: {math.degrees(dominant_angle):.1f}")

模拟的梯度样本数: 100
每个扇形的响应模长: [59.37200873 56.46805236  1.44845875  0.          0.          0.
  0.          0.          0.          0.          0.          4.56506744]
主方向所在的扇形索引: 0
主方向角度（弧度）: 0.262
主方向角度（度）: 15.0

（注：样本角度的均值是30度，恰好落在第0、1个扇形的边界上；本演示用固定30°分箱，量化后最大响应落在第0个扇形，其中心角为15°，故出现15°偏差。真实 SURF 采用连续滑动的60°扇形窗取最长响应向量，并不存在这种固定分箱量化；SURF 旋转不变性弱于 SIFT 的主要原因是 Haar 小波近似与方向估计精度不及 SIFT 的8方向梯度直方图法。）

4. 生成 SURF 描述子

确定主方向后，接下来要生成一个对该特征点周围区域进行描述的向量，即 SURF 描述子。

1. 构造描述区域：以特征点为中心，沿着上一步确定的主方向，构造一个边长为 20s 的正方形区域。
2. 划分子区域：将该正方形区域均匀划分为 4x4 = 16 个子区域。
3. 子区域统计：在每个子区域内（通常包含 5x5 个采样点），进行如下操作：
   • 使用尺寸为 2s 的 Haar 小波模板，计算每个采样点在 x 方向（相对于主方向）和 y 方向的响应值 dx 和 dy。
   • 为了增加对光照变化的鲁棒性，同时统计响应值的绝对值 |dx| 和 |dy|。
   • 对子区域内所有采样点的这4个值（dx, dy, |dx|, |dy|）分别求和。
4. 形成描述向量：每个子区域贡献4个值，16个子区域共得到 16 * 4 = 64 维的特征向量。这就是 SURF-64 描述子。
   • 为了获得更好的区分度，通常还会对这个64维向量进行归一化（例如 L2 归一化），以减弱光照变化的影响。

SURF 描述子生成示意图（文本示意）：
特征点
  |
  |-- 主方向
  |
  v
+------------------+ 边长为 20s
| 子区域1 | 子区域2 | ... |
|----------|----------| ... | 划分为 4x4 网格
| 子区域5 | 子区域6 | ... |
|----------|----------| ... |
|   ...   |   ...   | ... |
+------------------+

每个子区域 (5x5个采样点)：
  对每个点计算 Haar 响应 (dx, dy)
  统计：sum(dx), sum(dy), sum(|dx|), sum(|dy|)
   => 每个子区域输出4个值。

16个子区域 => 64维特征向量。

SURF 与 SIFT 的对比及特点

性能对比

• 速度：SURF 的核心优势。通过积分图和盒式滤波器，其计算速度通常比 SIFT 快约 3 倍，在某些实现中甚至更快。
• 尺度不变性：SURF 通过改变滤波器尺寸来构建尺度空间，其尺度不变性略逊于 SIFT 的多尺度金字塔下采样方法。在图像缩放比例非常大时，SURF 的稳定性可能下降。
• 旋转不变性：SURF 使用 60° 扇形滑动窗口确定主方向，其精度低于 SIFT 的梯度方向直方图法。因此，在图像发生大角度旋转时，SURF 的匹配性能下降比 SIFT 更明显。
• 鲁棒性：
  • 亮度变化：SURF 描述子对亮度变化具有很好的鲁棒性，部分原因是其使用了响应值的绝对值统计以及向量归一化。
  • 模糊：在一些测试中，SURF 对图像模糊的鲁棒性甚至优于 SIFT。
  • 视角变化：两者性能相近。

历史定位与现状

SURF 是传统手工设计特征时代的杰出代表，它在 SIFT 的基础上做出了极佳的工程优化，在需要实时性或计算资源有限的场景中（如2010年代的移动设备、嵌入式视觉）曾得到广泛应用。

然而，随着深度学习的崛起，尤其是基于卷积神经网络（CNN）的特征提取方法（如用预训练 CNN 的中间层输出作为通用特征）在性能上全面超越了 SURF、SIFT 等手工特征，并且得益于 GPU 加速，其速度也不再是瓶颈。因此，在今天的主流计算机视觉研究和产品开发中，SURF 已不再是首选特征。

尽管如此，理解 SURF 的设计思想——如何通过巧妙的近似和数据结构（积分图）来极大化算法效率——对于算法工程师的思维训练仍然具有很高的价值。在一些对传统库依赖较重或极度轻量化的边缘场景中，偶尔还能见到它的身影。

📝 动手练一练

1. 概念辨析：SURF 在构建尺度空间时，是保持图像大小不变而改变滤波器尺寸，这与 SIFT 的方法有何本质不同？这种设计带来的最大好处是什么？
2. 计算思考：假设一个特征点的尺度 s = 2，用于生成描述子的正方形区域边长是 20s。如果将该区域划分为 4x4 个子区域，每个子区域在 x 方向上应该有多少个像素（假设边长为像素数）？如果每个子区域内有 5x5 个采样点，这些采样点是密集覆盖还是间隔采样？

本章小结

本节我们深入学习了 SIFT 的加速改进版——SURF 特征描述子。SURF 通过引入 Haar-like 特征和积分图实现快速区域求和，并用盒式滤波器近似 Hessian 矩阵计算，将复杂的卷积运算转化为加减运算，从而实现了数倍于 SIFT 的计算速度。我们详细剖析了其构建流程：基于 Hessian 矩阵行列式检测特征点、利用 Haar 小波响应确定主方向、最后在方向对齐的区域内统计小波响应生成 64 维描述子。尽管在尺度和旋转不变性上略逊于 SIFT，但其在速度、亮度及模糊鲁棒性上的优势，使其成为传统特征提取时代的一个重要里程碑。

📋 行动清单

• 理解积分图的原理：亲手编写代码计算一个小图像的积分图，并用它快速计算任意矩形区域的和，体会其“一次计算，多次快速查询”的精髓。
• 对比 SIFT 与 SURF：在纸上或思维导图中列出 SIFT 和 SURF 在尺度空间构建、特征点检测算子、主方向确定和描述子生成四个关键步骤上的异同点。
• 思考算法设计：回顾 SURF 如何用“近似”（盒式滤波器）和“预计算”（积分图）来换取速度。尝试思考在你自己遇到的工程问题中，是否有可以借鉴这种“以精度换效率”或“以空间换时间”思想的地方。

—— 小象教研组