3、特征点检测（上）：Harris角点检测原理、数学推导与代码实现

各位同学，欢迎来到特征点检测的第一讲。

图像融合要做出全景大片的效果，第一步是什么？

你得先知道两张图里哪些地方是“同一个位置”。

比如你拍了两张风景照，左边那张的右上角有一座塔，右边那张的左上角也有同一座塔。你得让程序认出这个塔，然后才能把两张图拼起来。这个“认出”的过程，靠的就是特征点检测。

今天咱们聊的，是特征点检测里最经典、最基础的方法——Harris角点检测。我个人觉得，搞懂Harris，后面再学SIFT、SURF这些高级货，会轻松很多。

3.1 什么是角点？为什么它这么重要？

先问个问题：什么样的点适合做特征点？

你想想看，如果我在一张纯白墙上点一个像素点，换个角度拍，这个点还在吗？大概率找不到了，因为周围全是白的，没有任何区分度。

但如果我拍的是窗户的角，或者桌子的棱，换个角度，这个角依然很明显。这种点，就是我们要找的角点。

角点说白了，就是图像里两个不同方向都有明显变化的位置。比如一个“L”形的拐角，你往左挪一点，灰度值变了；往上挪一点，灰度值也变了。这种点在图像里非常稳定，适合用来做匹配。

3.2 Harris角点检测的直觉理解

Harris角点检测的思路其实很朴素。

想象一下，你拿一个小窗口在图像上滑动。这个小窗口可能是3x3或者5x5的像素块。

• 平坦区域：窗口往哪滑，里面的像素值都差不多。灰度变化很小。
• 边缘区域：窗口沿着边缘方向滑，灰度变化不大；但垂直于边缘滑，变化就很大。
• 角点区域：窗口往任何方向滑，灰度变化都很大。

Harris要做的，就是用一个数学公式，把第三种情况给揪出来。

嗯，这里要注意：它不关心你窗口滑了多远，它只关心微小的移动带来的变化。这个“微小”很重要，因为后面我们要用泰勒展开来近似计算。

3.3 数学推导：从直觉到公式

好，咱们来点硬核的。别怕，我会一步步拆开讲。

假设图像灰度函数为 I(x, y)。我们把窗口在 (x, y) 处移动一个微小的量 (u, v)，那么窗口内的灰度变化量 E(u, v) 可以写成：

E(u, v) = Σ w(x, y) * [I(x+u, y+v) - I(x, y)]²

这里的 w(x, y) 是窗口函数，通常用高斯权重，让窗口中心的像素贡献更大。

接下来是关键一步。I(x+u, y+v) 怎么算？用泰勒展开！

I(x+u, y+v) ≈ I(x, y) + u * Ix + v * Iy

其中 Ix 和 Iy 是图像在 x 和 y 方向的梯度。把这个近似代回 E(u, v) 的公式里：

E(u, v) ≈ Σ w(x, y) * [u * Ix + v * Iy]²
        = Σ w(x, y) * (u² * Ix² + 2uv * IxIy + v² * Iy²)

写成矩阵形式会更优雅：

E(u, v) ≈ [u, v] * M * [u, v]ᵀ

其中 M 是一个 2x2 的矩阵，叫做结构张量（Structure Tensor）：

M = Σ w(x, y) * [ Ix²   IxIy ]
                [ IxIy  Iy²  ]

这个 M 矩阵，就是Harris角点检测的核心。它描述了窗口内梯度的分布情况。

3.4 角点响应函数：如何判断是不是角点？

有了 M 矩阵，我们怎么判断一个点是不是角点呢？

Harris 提出了一个巧妙的办法：看 M 矩阵的两个特征值 λ₁ 和 λ₂。

• 两个特征值都很小：平坦区域，不是角点。
• 一个特征值大，一个特征值小：边缘区域，不是角点。
• 两个特征值都很大：角点！

但直接算特征值有点慢。Harris 给出了一个近似公式，叫做角点响应函数 R：

R = det(M) - k * (trace(M))²

其中：

• det(M) = λ₁ * λ₂ = Ix² * Iy² - (IxIy)²
• trace(M) = λ₁ + λ₂ = Ix² + Iy²
• k 是一个经验常数，通常在 0.04 到 0.06 之间。

判断规则很简单：

• R > 0 且较大：角点。
• R < 0：边缘。
• |R| 很小：平坦区域。

3.5 算法流程总结

咱们把整个流程捋一遍，我习惯用下面这张图来理解：

这张图把整个流程串起来了。从计算梯度开始，到最终输出角点坐标，每一步都有它的道理。

3.6 代码实现：从理论到实践

光说不练假把式。咱们用 Python 和 OpenCV 来实现一下。

import cv2
import numpy as np

def harris_corner_detection(image_path, k=0.04, threshold=0.01):
    """
    Harris角点检测实现
    
    参数:
        image_path: 图像路径
        k: Harris角点响应函数中的经验常数
        threshold: 角点筛选阈值（相对于最大R值的比例）
    
    返回:
        标记了角点的图像
    """
    # 1. 读取图像并转为灰度图
    img = cv2.imread(image_path)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    gray = np.float32(gray)
    
    # 2. 计算梯度（Sobel算子）
    Ix = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)
    Iy = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)
    
    # 3. 计算Ix², Iy², IxIy
    Ix2 = Ix * Ix
    Iy2 = Iy * Iy
    Ixy = Ix * Iy
    
    # 4. 高斯加权（构建M矩阵的每个元素）
    sigma = 2.0
    ksize = (5, 5)
    Sxx = cv2.GaussianBlur(Ix2, ksize, sigma)
    Syy = cv2.GaussianBlur(Iy2, ksize, sigma)
    Sxy = cv2.GaussianBlur(Ixy, ksize, sigma)
    
    # 5. 计算角点响应函数R
    det = Sxx * Syy - Sxy * Sxy
    trace = Sxx + Syy
    R = det - k * (trace ** 2)
    
    # 6. 阈值筛选
    R_max = np.max(R)
    if R_max == 0:
        R_max = 1  # 防止除零
    R_normalized = R / R_max
    corners = np.where(R_normalized > threshold)
    
    # 7. 在原图上标记角点
    result_img = img.copy()
    for y, x in zip(corners[0], corners[1]):
        cv2.circle(result_img, (x, y), 3, (0, 0, 255), -1)
    
    return result_img

# 使用示例
# result = harris_corner_detector('building.jpg')
# cv2.imshow('Harris Corners', result)
# cv2.waitKey(0)

3.7 实际项目中的注意事项

最后，分享几个我在项目中踩过的坑：

• 图像噪声问题： Harris 对噪声比较敏感。如果图像质量不好，建议先做一次高斯滤波去噪。
• 尺度问题： Harris 检测到的角点对尺度变化不鲁棒。如果你要拼接的两张图缩放比例不同，Harris 可能找不到对应的点。这时候可以考虑用 SIFT 或 SURF。
• 计算效率： 对于大尺寸图像，Harris 的计算量不小。我一般会先缩小图像，检测到角点后再映射回原图。

好了，Harris角点检测的原理、推导和实现就讲到这里。它虽然简单，但思想非常经典。后面的课程里，我们会基于它继续深入，看看更高级的特征点检测方法是怎么演进的。

---

公众号:蓝海资料掘金营，微信deep3321