实现 Harris 角点检测器

Posted 2023-03-12

【中文标题】实现 Harris 角点检测器

【英文标题】：Implementing a Harris corner detector

【发布时间】：2011-04-21 04:55:55

【问题描述】：

我正在实施一个 Harris 角检测器以用于教育目的，但我卡在了 harris 响应部分。基本上，我正在做的是：

计算 x 和 y 方向的图像强度梯度 (1) 的模糊输出 根据 (2) 的输出计算 Harris 响应 在 3x3 邻域和阈值输出中抑制 (3) 输出中的非最大值

1 和 2 似乎工作正常；但是，作为 Harris 响应，我得到的值非常小，并且没有任何点达到阈值。输入是标准的户外摄影。

[...]
[Ix, Iy] = intensityGradients(img);
g = fspecial('gaussian');
Ix = imfilter(Ix, g);
Iy = imfilter(Iy, g);
H = harrisResponse(Ix, Iy);
[...]

function K = harrisResponse(Ix, Iy)
    max = 0;
    [sy, sx] = size(Ix);
    K = zeros(sy, sx);
    for i = 1:sx,
        for j = 1:sy,
            H = [Ix(j,i) * Ix(j,i), Ix(j,i) * Iy(j,i)
                Ix(j,i) * Iy(j,i), Iy(j,i) * Iy(j,i)];
            K(j,i) = det(H) / trace(H);
            if K(j,i) > max,
                max = K(j,i);
            end
        end
    end
    max
end

对于示例图片，max 最终为 6.4163e-018，这似乎太低了。

【参考方案1】：

Harris 角点检测中的角点被定义为“区域中的最高值像素”（通常是 3X3 或 5x5），因此您关于没有达到“阈值”的评论对我来说似乎很奇怪。只需收集所有像素值高于其周围5x5 邻域中的所有其他像素的所有像素。

除此之外： 我不是 100% 确定，但我认为你应该有：

K(j,i) = det(H) - lambda*(trace(H)^2) 其中 lambda 是适用于您的情况的正常数（Harris 建议值为 0.04）。

一般来说，过滤输入的唯一明智时机是在此之前：

[Ix, Iy] = intensityGradients(img);

过滤 Ix2、Iy2 和 Ixy 对我来说没有多大意义。

此外，我认为您的示例代码在这里是错误的（函数harrisResponse 有两个或三个输入变量吗？）：

H = harrisResponse(Ix2, Ixy, Iy2);
[...]

function K = harrisResponse(Ix, Iy)

【讨论】：

我已恢复为不再过滤 Ix2 等，因此在 *** 上的副本中留下了一些错误。

问题是我没有将 3x3 正方形中的所有像素相加来找出 Ix2 等；相反，我刚刚使用了相应的像素。在改变 H 以总结所有 9 个像素的所有 Ix2、Ixy 和 Iy2 后，它看起来非常漂亮。

det(H)/trace(H) 是在没有 lambda 的情况下常用的近似值。

我不知道最后一招。好的。看来你自己解决了问题，很好！ （并且古老的技巧仍然有效：向某人解释问题可以帮助您找到解决方案）这是有效的代码吗？

【参考方案2】：

提议的实施效率极低。 让我们在计算梯度之后开始（也可以优化）：

A = Ix.^2;
B = Iy.^2;
C = (Ix.*Iy).^4;
lambda = 0.04;

H = (A.*B - C) - lambda*(A+B).^2;

% if you really need max:
max(H(:))

不需要循环，因为 Matlab 讨厌循环。

【讨论】：

但是为什么要计算 C = (Ix.*Iy).^4 而不是简单的 C = (Ix.*Iy) 呢？

【参考方案3】：

在计算机视觉系统工具箱中有一个函数叫做detectHarrisFeatures。

【参考方案4】：

基本上，Harris 角点检测将有 5 个步骤：

梯度计算 高斯平滑 Harris 度量计算 非最大抑制 阈值

如果你是在MATLAB中实现的，那么很容易理解算法并得到结果。

下面的MATLAB代码可以帮助你解决疑惑：

% Step 1: Compute derivatives of image
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');

% Step 2: Smooth space image derivatives (gaussian filtering)
Ix2 = conv2(Ix .^ 2, g, 'same');
Iy2 = conv2(Iy .^ 2, g, 'same');
Ixy = conv2(Ix .* Iy, g, 'same');

% Step 3: Harris corner measure
harris = (Ix2 .* Iy2 - Ixy .^ 2) ./ (Ix2 + Iy2);

% Step 4: Find local maxima (non maximum suppression)
mx = ordfilt2(harris, size .^ 2, ones(size));

% Step 5: Thresholding
harris = (harris == mx) & (harris > threshold);

【参考方案5】：

我用python实现的解决方案，它对我有用，我希望你能找到你要找的东西

import numpy as np
import matplotlib.pyplot as plt
from PIL.Image import *
from scipy import ndimage

def imap1(im):
    print('testing the picture . . .')
    a = Image.getpixel(im, (0, 0))
    if type(a) == int:
        return im
    else:
        c, l = im.size
        imarr = np.asarray(im)
        neim = np.zeros((l, c))
        for i in range(l):
            for j in range(c):
                t = imarr[i, j]
                ts = sum(t)/len(t)
                neim[i, j] = ts
        return neim

def Harris(im):
    neim = imap1(im)
    imarr = np.asarray(neim, dtype=np.float64)
    ix = ndimage.sobel(imarr, 0)
    iy = ndimage.sobel(imarr, 1)
    ix2 = ix * ix
    iy2 = iy * iy
    ixy = ix * iy
    ix2 = ndimage.gaussian_filter(ix2, sigma=2)
    iy2 = ndimage.gaussian_filter(iy2, sigma=2)
    ixy = ndimage.gaussian_filter(ixy, sigma=2)
    c, l = imarr.shape
    result = np.zeros((c, l))
    r = np.zeros((c, l))
    rmax = 0
    for i in range(c):
        print('loking for corner . . .')
        for j in range(l):
            print('test ',j)
            m = np.array([[ix2[i, j], ixy[i, j]], [ixy[i, j], iy2[i, j]]], dtype=np.float64)
            r[i, j] = np.linalg.det(m) - 0.04 * (np.power(np.trace(m), 2))
            if r[i, j] > rmax:
                rmax = r[i, j]
    for i in range(c - 1):
        print(". .")
        for j in range(l - 1):
            print('loking')
            if r[i, j] > 0.01 * rmax and r[i, j] > r[i-1, j-1] and r[i, j] > r[i-1, j+1]\
                                     and r[i, j] > r[i+1, j-1] and r[i, j] > r[i+1, j+1]:
                result[i, j] = 1

    pc, pr = np.where(result == 1)
    plt.plot(pr, pc, 'r+')
    plt.savefig('harris_test.png')
    plt.imshow(im, 'gray')
    plt.show()
    # plt.imsave('harris_test.png', im, 'gray')

im = open('chess.png')
Harris(im)

【讨论】：

什么是 pc，以及它提供了什么？

@acousticpython pc 和 pr 是 index where result == 1，意味着 result[pc][pr] == 1，pc 和 pr 中相同位置的项目是 result 等于 1二维数组