【问题标题】:Implementing a Harris corner detector实现 Harris 角点检测器
【发布时间】:2011-04-21 04:55:55
【问题描述】:

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

  1. 计算 x 和 y 方向的图像强度梯度
  2. (1) 的模糊输出
  3. 根据 (2) 的输出计算 Harris 响应
  4. 在 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,这似乎太低了。

【问题讨论】:

    标签: algorithm matlab computer-vision feature-detection corner-detection


    【解决方案1】:

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

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

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

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

    [Ix, Iy] = intensityGradients(img);

    过滤 Ix2Iy2Ixy 对我来说没有多大意义。

    此外,我认为您的示例代码在这里是错误的(函数harrisResponse 有两个或三个输入变量吗?):

    H = harrisResponse(Ix2, Ixy, Iy2);
    [...]
    
    function K = harrisResponse(Ix, Iy)
    

    【讨论】:

    • 我已恢复为不再过滤 Ix2 等,因此在 stackoverflow 上的副本中留下了一些错误。
    • 问题是我没有将 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 个步骤:

      1. 梯度计算
      2. 高斯平滑
      3. Harris 度量计算
      4. 非最大抑制
      5. 阈值

      如果你是在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二维数组
        猜你喜欢
        • 1970-01-01
        • 2023-03-18
        • 1970-01-01
        • 2012-10-24
        • 2018-10-22
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2013-08-17
        相关资源
        最近更新 更多