La implementación de un detector de esquina Harris

Question

Me estoy poniendo en práctica un detector de esquina Harris con fines educativos, pero estoy atascado en la parte respuesta Harris. Básicamente, lo que estoy haciendo, es la siguiente:

1. Calcular gradientes de intensidad de la imagen en x e y-dirección
2. salida de la falta de definición de (1)
3. Compute Harris respuesta sobre la producción de (2)
4. Suppress no maximas en la producción de (3) en una salida de 3x3-barrio y umbral

1 y 2 parece que funcionan bien; Sin embargo, me da valores muy pequeños como la respuesta de Harris, y ningún momento alcanzar el umbral. Entrada es una fotografía al aire libre estándar.

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

En la imagen de ejemplo, extremos max siendo 6.4163e-018 que parece demasiado bajo.

Answer 1

Un rincón en la detección de esquina Harris se define como "el pixel más alto valor en una región" (por lo general 3X3 o 5x5) para que su comentario sobre ningún punto de llegar a un "umbral" me parece extraño. Para ello, consigue todos los píxeles que tienen un valor más alto que todos los demás píxeles en el barrio 5x5 alrededor de ellos.

Aparte de eso: No estoy 100% seguro, pero creo que debe tener:

K(j,i) = det(H) - lambda*(trace(H)^2) Donde lambda es una constante positiva que las obras en su caso (y Harris sugirieron valor es 0,04).

En general el momento sólo tiene sentido para filtrar la entrada es antes de este punto:

[Ix, Iy] = intensityGradients(img);

Filtrado Ix2, Iy2 y Ixy no tiene mucho sentido para mí.

Además, creo que el código de ejemplo es incorrecto aquí (Cómo funciona la harrisResponse tener dos o tres variables de entrada?):

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

function K = harrisResponse(Ix, Iy)

Answer 2

La solución que he implementado con el pitón, que funciona para mí espero que encuentre lo que busca

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)

Answer 3

aplicación propuesta es terriblemente ineficiente. Vamos a comenzar después de calcular los gradientes (que pueden ser optimizados también):

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(:))

No se requiere bucles, porque Matlab odia bucles.

Answer 4

Básicamente, la detección esquina Harris tendrá 5 pasos:

1. Degradado cálculo
2. Gaussian suavizado
3. Harris medida cálculo
4. No-máxima supresión
5. Umbral

Si va a implementar en MATLAB, que será fácil de entender el algoritmo y obtener los resultados.

El siguiente código de MATLAB puede ayudarle a resolver sus dudas:

% 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);

Answer 5

Hay una función para que en el sistema de la caja de herramientas de visión de computadora llamado detectHarrisFeatures.