La implementación de un detector de esquina Harris
-
27-09-2019 - |
Pregunta
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:
- Calcular gradientes de intensidad de la imagen en x e y-dirección
- salida de la falta de definición de (1)
- Compute Harris respuesta sobre la producción de (2)
- 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.
Solución
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)
Otros consejos
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)
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.
Básicamente, la detección esquina Harris tendrá 5 pasos:
- Degradado cálculo
- Gaussian suavizado
- Harris medida cálculo
- No-máxima supresión
- 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);
Hay una función para que en el sistema de la caja de herramientas de visión de computadora llamado detectHarrisFeatures
.