Implementación de un detector de esquina Harris

Question

Estoy implementando un detector de esquina Harris con fines educativos, pero estoy atascado en la parte de respuesta de harris. Básicamente, lo que estoy haciendo, es:

1. Compute gradientes de intensidad de la imagen en x e y-dirección
salida
2. la falta de definición de (1)
respuesta
3. Compute Harris sobre la producción de (2)
4. Suprimir no máximas en la salida de (3) en una vecindad de 3x3 y salida de umbral

1 y 2 parecen funcionar bien; sin embargo, obtengo valores muy pequeños como la respuesta de Harris, y ningún punto alcanza el umbral. La entrada es una fotografía estándar al aire libre.

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

Para la imagen de muestra, max termina siendo 6.4163e-018 que parece demasiado bajo.

Answer 1

Una esquina en la detección de esquina de Harris se define como "el píxel de más alto valor en una región" (generalmente 3X3 o 5x5) por lo que su comentario sobre ningún punto que alcanza un "umbral" me parece extraño. Simplemente recolecte todos los píxeles que tengan un valor más alto que todos los demás píxeles en el vecindario 5x5 que los rodea.

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 funciona en su caso (y Harris sugirió valor es 0,04).

En general, el único momento sensata 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í (qué función harrisResponse tener dos o tres variables de entrada?):

aplicación

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

function K = harrisResponse(Ix, Iy) 

Answer 2

propuesto es terriblemente ineficiente. Lets' empezar 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 hay bucles necesarios, porque Matlab odia bucles.

Answer 3

Hay una función para eso en Computer Vision System Toolbox llamada detectHarrisFeatures.

Answer 4

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

1. de cálculo de gradiente
2. gaussiana suavizado
3. Harris medida cálculo
4. supresión no máximo
5. Umbral

Si estás implementando en MATLAB, será fácil 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

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)