如何用Numpy正确实现Harris Corner Detector的矩阵求和计算

你对问题的判断是正确的，Harris算子中的M矩阵需要对指定窗口内的Ixx、Ixy、Iyy做求和计算，原有代码仅使用了当前点的单像素值，不符合算法逻辑。

方案1：修改原有循环逻辑（易理解） #

仅需要调整双重循环内的窗口取值与求和逻辑即可，由于你已经提前对图像做了边缘填充，直接取对应窗口区域不会出现下标越界问题：

for y in range(output.shape[0]):
    for x in range(output.shape[1]):
        # 提取以(y, x)为中心的window_size大小窗口
        win_Ixx = Ixx[y:y+window_size, x:x+window_size]
        win_Ixy = Ixy[y:y+window_size, x:x+window_size]
        win_Iyy = Iyy[y:y+window_size, x:x+window_size]
        # 对窗口内值求和构造M矩阵
        m = np.array([
            [np.sum(win_Ixx), np.sum(win_Ixy)],
            [np.sum(win_Ixy), np.sum(win_Iyy)]
        ])
        output[y,x] = np.linalg.det(m) - k*(np.trace(m)**2)

方案2：向量化实现（高性能，推荐） #

嵌套循环在处理大尺寸图像时效率极低，可以用scipy的uniform_filter实现快速窗口求和，完全避免Python层循环，速度提升数个量级：

from scipy.ndimage import uniform_filter

padded_size = window_size//2
padded_image = np.pad(image, padded_size, 'edge')
Ix, Iy = partial_derivative(padded_image)
Ixx = Ix**2
Ixy = Iy*Ix
Iyy = Iy**2

# 窗口求和：均值滤波结果乘以窗口像素总数
sxx = uniform_filter(Ixx, size=window_size) * window_size * window_size
sxy = uniform_filter(Ixy, size=window_size) * window_size * window_size
syy = uniform_filter(Iyy, size=window_size) * window_size * window_size

# 裁剪掉padding部分，匹配原图尺寸
sxx = sxx[padded_size:-padded_size, padded_size:-padded_size]
sxy = sxy[padded_size:-padded_size, padded_size:-padded_size]
syy = syy[padded_size:-padded_size, padded_size:-padded_size]

# 向量化计算所有像素的角点响应值
det = sxx * syy - sxy ** 2
trace = sxx + syy
output = det - k * (trace ** 2)

return output

注意事项 #

• 如果不能引入scipy依赖，也可以用numpy的np.lib.stride_tricks.as_strided构造滑动窗口后求和，实现复杂度稍高，需要注意维度与内存对齐问题
• 窗口大小通常取奇数（3、5、7等），保证中心像素位置准确

内容的提问来源于stack exchange，提问作者johnywalker