【计算机视觉】Harris角点检测

目录

一、角点是什么

二、基于SSD（平方差和）的角点检测

三、基于结构张量的检测

3.1 Shi-Tomasi角点检测

3.2 Harris角点检测

相关代码实现

（1）基于SSD的角点检测算法

(2) Shi-Tomasi角点检测

(3) Harris角点检测

（4）opencv实现

---

一、角点是什么

首先，我们要明白计算机是如何“看待”角点的。想象一下，用一个很小的方框（这个框框我们后面通常称之为“窗口”）在图像上移动。

• 平坦区域 (Flat Region)：如果窗口在一个颜色均匀的区域（比如蓝天、白墙），无论朝哪个方向移动这个窗口，窗口里的内容基本没什么变化。

• 边缘 (Edge)：如果窗口在一条线上（比如建筑的轮廓），沿着线的方向移动，窗口内容变化不大；但如果垂直于线的方向移动，内容就会有剧烈变化。

• 角点 (Corner)：如果窗口在一个角上，那么无论朝哪个方向移动这个窗口，窗口里的内容都会发生显著的变化。

二、基于SSD（平方差和）的角点检测

思路：它先在图像上取一个像素点，以这个点为中心画一个“窗口”，然后，它将这个窗口向周围的8个方向（上、下、左、右、左上、右上、左下、右下）各平移一小步（比如一个像素）。每平移一次，就计算一下“平移后的窗口”和“原始窗口”之间对应像素的差异。怎么算差异呢？就是计算平方差之和 (Sum of Squared Differences, SSD)。

这个值越大，说明两个窗口的差异越大。

在8个方向的SSD值中，我们只关心最小的那个。为什么是最小的？因为如果连最小的SSD值都很大，就说明这个点往任何方向移动，变化都很大，那它很可能就是个角点！

优点：原理超级简单，容易理解和实现。

缺点：计算量巨大！ 因为对于图像中的每一个像素，都要进行8次窗口平移和SSD计算。如果图像很大，会非常耗时。

三、基于结构张量的检测

（1）基础：图像梯度

梯度简单说就是图像中像素灰度值的变化率。

• Ix​：像素在 x 方向（水平）上的梯度。如果一个地方颜色从左到右变化很快，Ix​的绝对值就很大。

• Iy​：像素在 y 方向（垂直）上的梯度。如果一个地方颜色从上到下变化很快，Iy​的绝对值就很大。

（2）结构张量M

现在，我们把一个小窗口内的所有像素的梯度信息汇总起来，构建一个2x2的矩阵 M，这个 M 就是结构张量。

• w(x,y) 是一个加权函数，通常是高斯函数（Gaussian window）。使用高斯加权是为了让窗口中心的梯度贡献更大，边缘的贡献更小，从而提高检测的鲁棒性。

这个矩阵 M 非常神奇，它概括了窗口内所有梯度变化的方向和强度。通过分析这个矩阵，我们就能知道这个区域是平坦、是边缘还是角点。

怎么分析呢？答案是：看 M 的特征值 (Eigenvalues)，我们记作 λ1​ 和 λ2​。

（3）用特征值解读图像

特征值在线性代数中代表了矩阵在特定方向上变换的“伸缩”程度。在这里，λ1​ 和 λ2​ 的大小直接反映了窗口内梯度在两个主要方向上的强度。

• 平坦区域：窗口内梯度很小，所以 λ1​ 和 λ2​ 都很小。

• 边缘：只有一个方向梯度大（垂直于边缘的方向），另一个方向梯度小（沿着边缘的方向）。所以一个特征值大，另一个特征值小（比如 λ1​≫λ2​）。

• 角点：两个方向的梯度都很大。所以两个特征值都很大（λ1​ 和 λ2​ 都大）。

3.1 Shi-Tomasi角点检测

Shi-Tomasi 的作者认为，一个好的角点，应该是两个特征值都很大。最直接的判断方法就是：看两个特征值中较小的那一个。如果连较小的那个都很大，那另一个肯定也很大。

所以，它的角点响应函数 (Response Function) R 非常简单：

                                             

只要 R 大于某个阈值，就认为这个点是角点。

3.2 Harris角点检测

Harris 和 Stephens 提出了一个巧妙的响应函数，它不需要直接计算特征值，因为在当时计算特征值是很耗费计算资源的。

他们利用了矩阵的两个重要性质：行列式 (Determinant) 和 迹 (Trace)。

                                        

Harris 的角点响应函数 R 定义为：

                       

其中 k 是一个经验常数，通常取 0.04 到 0.06 之间。

好那么这个公式是怎么工作的？

• 角点：λ1​ 和 λ2​ 都很大，所以 det(M)=λ1​λ2​ 会是一个很大的正数。trace(M) 也很大，但因为 k 很小，所以 det(M) 占主导，最终的 R 是一个很大的正数。

• 边缘：一个特征值远大于另一个（λ1​≫λ2​≈0）。此时 det(M)=λ1​λ2​≈0。而 trace(M)=λ1​+λ2​≈λ1​。所以 R≈0−k(λ1​)2，会是一个很大的负数。

• 平坦区域：两个特征值都接近0，所以 det(M) 和 trace(M) 都接近0，最终的 R 也接近0。

因此，只需要计算 R 值，如果它是一个远大于0的数，那它就是一个角点。这就避免了开方和求解特征值，在计算上更高效。

后处理：非极大值抑制 (NMS)

无论是哪种方法，在计算出响应值 R 之后，一个角点周围的像素点很可能也有比较高的响应值。这会导致检测到的角点“扎堆”。

非极大值抑制 (Non-Maximum Suppression, NMS) 就是解决这个问题的关键步骤。它的逻辑是：

1. 在一个小邻域内（比如 7x7 的窗口）。

2. 找到这个邻域里响应值最大的那个点。

3. 只保留这个最大的点作为角点，把邻域里其他所有的点都“抑制”掉（即认为它们不是角点）。

通过这个步骤，就可以把一堆密集的响应点精简成一个最能代表该位置的角点，得到干净、准确的结果。

相关代码实现

（1）基于SSD的角点检测算法

import cv2import numpy as npimport osimport timedef non_max_suppression(response, window_size):h, w = response.shape#创建一个和输入响应图response大小一样，内容全为0的矩阵suppressed#suppressed用于存放抑制后的结果suppressed = np.zeros_like(response, dtype=np.float32)pad = window_size // 2for r in range(pad, h - pad):for c in range(pad, w - pad):window = response[r - pad:r + pad + 1, c - pad:c + pad + 1]if response[r, c] >= window.max():suppressed[r, c] = response[r, c]return suppresseddef corner_detector_ssd(image, window_size=5, threshold=0.01, nms_size=15):gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)gray = np.float32(gray)h, w = gray.shape#创建一个与灰度图等大的矩阵，用来存储每个像素的“角点响应值”。#响应值越高，代表该点越可能是角点response = np.zeros_like(gray)pad = window_size // 2for r in range(pad + 1, h - pad - 1):print(f"  正在处理行: {r}/{h-pad-1}", end='\r')for c in range(pad + 1, w - pad - 1):min_ssd = float('inf')shifts = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]original_window = gray[r - pad:r + pad + 1, c - pad:c + pad + 1]for dr, dc in shifts:shifted_window = gray[r + dr - pad : r + dr + pad + 1, c + dc - pad : c + dc + pad + 1]ssd = np.sum((original_window - shifted_window) ** 2)if ssd < min_ssd:min_ssd = ssdresponse[r, c] = min_ssdresponse_thresh = response > (response.max() * threshold)final_response = non_max_suppression(response * response_thresh, window_size=nms_size)corners = np.argwhere(final_response > 0)return cornersif __name__ == '__main__':IMAGE_PATH = 'test.png'OUTPUT_DIR = 'result'OUTPUT_FILENAME = 'method_1.png'if not os.path.exists(OUTPUT_DIR):os.makedirs(OUTPUT_DIR)image = cv2.imread(IMAGE_PATH)if image is None:print(f"错误: 无法加载图片 '{IMAGE_PATH}'")else:start_time = time.time()corners = corner_detector_ssd(image.copy())for r_coord, c_coord in corners:cv2.circle(image, (c_coord, r_coord), 3, (0, 0, 255), -1)output_path = os.path.join(OUTPUT_DIR, OUTPUT_FILENAME)cv2.imwrite(output_path, image)end_time = time.time()print(f"方法1:检测到 {len(corners)} 个角点。耗时: {end_time - start_time:.2f} 秒。")print(f"结果已保存至 {output_path}")

(2) Shi-Tomasi角点检测

import cv2import numpy as npimport scipy.signalimport osimport timedef non_max_suppression(response, window_size):h, w = response.shapesuppressed = np.zeros_like(response, dtype=np.float32)pad = window_size // 2for r in range(pad, h - pad):for c in range(pad, w - pad):window = response[r - pad:r + pad + 1, c - pad:c + pad + 1]if response[r, c] >= window.max():suppressed[r, c] = response[r, c]return suppresseddef corner_detector_shi_tomasi(image, window_size=5, threshold=0.02):#图像预处理gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)gray = np.float32(gray) / 255.0h, w = gray.shape#计算图像梯度sobel_kernel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])sobel_kernel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])Ix = scipy.signal.convolve2d(gray, sobel_kernel_x, mode='same', boundary='symm')Iy = scipy.signal.convolve2d(gray, sobel_kernel_y, mode='same', boundary='symm')#构建结构张量矩阵Ixx = Ix ** 2Iyy = Iy ** 2Ixy = Ix * Iy#高斯模糊 (GaussianBlur) 对 Ixx, Iyy, Ixy 进行平滑处理Sxx = cv2.GaussianBlur(Ixx, (window_size, window_size), 1.5)Syy = cv2.GaussianBlur(Iyy, (window_size, window_size), 1.5)Sxy = cv2.GaussianBlur(Ixy, (window_size, window_size), 1.5)response = np.zeros_like(gray, dtype=np.float32)for r in range(h):for c in range(w):M = np.array([[Sxx[r, c], Sxy[r, c]], [Sxy[r, c], Syy[r, c]]])eigenvalues, _ = np.linalg.eig(M)response[r, c] = np.min(eigenvalues)response_thresh = response > (response.max() * threshold)final_response = non_max_suppression(response * response_thresh, window_size=15)corners = np.argwhere(final_response > 0)return cornersif __name__ == '__main__':IMAGE_PATH = 'test.png'OUTPUT_DIR = 'result'OUTPUT_FILENAME = 'method_2.png'if not os.path.exists(OUTPUT_DIR):os.makedirs(OUTPUT_DIR)image = cv2.imread(IMAGE_PATH)if image is None:print(f"错误: 无法加载图片 '{IMAGE_PATH}'")else:start_time = time.time()corners = corner_detector_shi_tomasi(image.copy())for r_coord, c_coord in corners:cv2.circle(image, (c_coord, r_coord), 3, (0, 0, 255), -1)output_path = os.path.join(OUTPUT_DIR, OUTPUT_FILENAME)cv2.imwrite(output_path, image)end_time = time.time()print(f"方法2：检测到 {len(corners)} 个角点。耗时: {end_time - start_time:.2f} 秒。")print(f"结果已保存至 {output_path}")

(3) Harris角点检测

import cv2
import numpy as np
import os
import timedef non_max_suppression(response, window_size=5):h, w = response.shapesuppressed = np.zeros_like(response, dtype=np.float32)pad = window_size // 2for r in range(pad, h - pad):for c in range(pad, w - pad):window = response[r - pad:r + pad + 1, c - pad:c + pad + 1]if response[r, c] == window.max():suppressed[r, c] = response[r, c]return suppresseddef corner_detector_harris(image, window_size=5, k=0.04, threshold=0.01, nms_window=7):gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)gray = np.float32(gray) / 255.0Ix = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3)Iy = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3)Ixx = Ix * IxIyy = Iy * IyIxy = Ix * IySxx = cv2.GaussianBlur(Ixx, (window_size, window_size), 1.5)Syy = cv2.GaussianBlur(Iyy, (window_size, window_size), 1.5)Sxy = cv2.GaussianBlur(Ixy, (window_size, window_size), 1.5)det_M = (Sxx * Syy) - (Sxy ** 2)trace_M = Sxx + Syyresponse = det_M - k * (trace_M ** 2)response_thresh = response > (response.max() * threshold)final_response = non_max_suppression(response * response_thresh, window_size=nms_window)corners = np.argwhere(final_response > 0)corners = [(int(c), int(r)) for r, c in corners]  return cornersif __name__ == '__main__':IMAGE_PATH = 'test.png'OUTPUT_DIR = 'result'OUTPUT_FILENAME = 'method_3.png'if not os.path.exists(OUTPUT_DIR):os.makedirs(OUTPUT_DIR)image = cv2.imread(IMAGE_PATH)if image is None:print(f"错误: 无法加载图片 '{IMAGE_PATH}'")else:start_time = time.time()corners = corner_detector_harris(image.copy(), window_size=5, k=0.04, threshold=0.05, nms_window=7)for x, y in corners:cv2.circle(image, (x, y), 3, (0, 0, 255), -1)output_path = os.path.join(OUTPUT_DIR, OUTPUT_FILENAME)cv2.imwrite(output_path, image)end_time = time.time()print(f"方法3：检测到 {len(corners)} 个角点。耗时: {end_time - start_time:.2f} 秒。")print(f"结果已保存至 {output_path}")

（4）opencv实现

import cv2
import numpy as np
import os
import timedef corner_detector_opencv_nms(image, block_size=7, ksize=5, k=0.04, threshold=0.05, nms_size=15):gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)gray = np.float32(gray)dst = cv2.cornerHarris(gray, block_size, ksize, k)dst = cv2.dilate(dst, None)mask = dst > threshold * dst.max()corners = []h, w = dst.shapefor y in range(h):for x in range(w):if mask[y, x]:y0 = max(0, y - nms_size // 2)y1 = min(h, y + nms_size // 2 + 1)x0 = max(0, x - nms_size // 2)x1 = min(w, x + nms_size // 2 + 1)if dst[y, x] == np.max(dst[y0:y1, x0:x1]):corners.append((x, y))return cornersif __name__ == '__main__':IMAGE_PATH = 'test.png'OUTPUT_DIR = 'result'OUTPUT_FILENAME = 'method_4.png'if not os.path.exists(OUTPUT_DIR):os.makedirs(OUTPUT_DIR)image = cv2.imread(IMAGE_PATH)if image is None:print(f"错误: 无法加载图片 '{IMAGE_PATH}'")else:start_time = time.time()corners = corner_detector_opencv_nms(image.copy())for x, y in corners:cv2.circle(image, (x, y), 3, (0, 0, 255), -1)output_path = os.path.join(OUTPUT_DIR, OUTPUT_FILENAME)cv2.imwrite(output_path, image)end_time = time.time()print(f"方法4：检测到 {len(corners)} 个角点。耗时: {end_time - start_time:.2f} 秒。")print(f"结果已保存至 {output_path}")