Внедряване на ъглов детектор на Харис

Внедрявам ъгълов детектор на Харис за образователни цели, но съм заседнал в частта за отговор на Харис. По принцип това, което правя, е:

  1. Изчислете градиентите на интензитета на изображението в посока x и y
  2. Изход за замъгляване на (1)
  3. Изчислете реакцията на Харис върху резултата от (2)
  4. Потискане на не-максимумите в изхода на (3) в 3x3-съседство и прагов изход

1 и 2 изглежда работят добре; обаче получавам много малки стойности като реакция на Харис и нито една точка не достига прага. Входът е стандартна фотография на открито.

[...]
[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

За примерната снимка максимумът в крайна сметка е 6.4163e-018, което изглежда твърде ниско.


person Etan    schedule 05.10.2010    source източник


Отговори (5)


Ъгъл в откриването на ъгли на Харис се определя като „пиксел с най-висока стойност в даден регион“ (обикновено 3X3 или 5x5), така че коментарът ви за това, че няма точка за достигане на „праг“, ми изглежда странен. Просто съберете всички пиксели, които имат по-висока стойност от всички други пиксели в 5x5 квартала около тях.

Освен това: не съм 100% сигурен, но мисля, че трябва да имате:

K(j,i) = det(H) - lambda*(trace(H)^2) Където ламбда е положителна константа, която работи във вашия случай (а предложената от Харис стойност е 0,04).

По принцип единственият разумен момент за филтриране на вашия вход е преди тази точка:

[Ix, Iy] = intensityGradients(img);

Филтрирането на Ix2, Iy2 и Ixy няма много смисъл за мен.

Освен това мисля, че вашият примерен код тук е грешен (има ли функция harrisResponse две или три входни променливи?):

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

function K = harrisResponse(Ix, Iy)
person jilles de wit    schedule 05.10.2010
comment
Върнах се да не филтрирам повече Ix2 и т.н., следователно имаше оставена грешка в копието на stackoverflow. - person Etan; 05.10.2010
comment
Проблемът беше, че не сумирах всички пиксели в квадрата 3x3, за да разбера Ix2 и т.н.; вместо това току-що използвах съответния пиксел. След промяна на H по начин, който сумира всички Ix2, Ixy и Iy2 за всичките 9 пиксела, изглежда много добре. - person Etan; 05.10.2010
comment
det(H)/trace(H) е често използвано приближение в случай, че няма да имате ламбда. - person Etan; 05.10.2010
comment
Не знаех за последния трик. хубаво. Изглежда сам си решил проблема, браво! (и старият трик все още работи: самото обяснение на проблема на някого ви помага да стигнете до решението) това ли е работещият код? - person jilles de wit; 05.10.2010

Решението, което внедрих с 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)
person Walid Bousseta    schedule 07.05.2017
comment
какво е pc и pr какво дава? - person acoustic python; 16.11.2020
comment
@acousticpython pc и pr са индекс, където резултат == 1, означава, че резултат [pc][pr] == 1, елементи в pc и pr на една и съща позиция са, където резултатът е равен на единица в 2d масива - person Walid Bousseta; 16.11.2020

По принцип откриването на ъгъла на Харис ще има 5 стъпки:

  1. Градиентно изчисление
  2. Гаусово изглаждане
  3. Изчисление на Харис
  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);
person Megha Billure    schedule 21.04.2016

Предложената реализация е ужасно неефективна. Нека започнем след изчисляването на градиенти (които също могат да бъдат оптимизирани):

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 мрази цикли.

person Piotr    schedule 17.09.2013
comment
Но защо да изчисляваме C = (Ix.*Iy).^4 вместо просто C = (Ix.*Iy)? - person Hilder Vitor Lima Pereira; 31.03.2015

Има функция за това в Computer Vision System Toolbox, наречена detectHarrisFeatures.

person Dima    schedule 10.04.2014