Trova gli angoli di un poligono rappresentata da una regione di maschera

Question

BW = poly2mask(x, y, m, n) calcola una regione binario di interesse (ROI) maschera, BW, da un poligono ROI, rappresentata dai vettori x ed y. La dimensione di BW è m-by-n.

poly2mask insiemi pixel in BW che si trovano all'interno del poligono (X, Y) per 1 e imposta i pixel all'esterno del poligono di 0.

Problema: Dato una maschera binaria tale BW di un quadrilatero convesso, quale sarebbe il modo più efficace per determinare i quattro angoli?

Ad esempio,

soluzione migliore finora: Usa edge per trovare le linee di delimitazione, la trasformata di Hough per trovare le 4 linee l'immagine bordo in e poi trovare i punti di intersezione di quelle 4 linee o utilizzare un rivelatore d'angolo sull'immagine del bordo. Sembra complicato, e non posso fare a meno di pensare che ci sia una soluzione più semplice là fuori.

Btw, convhull non restituisce sempre 4 punti (forse qualcuno può suggerire le opzioni qhull per impedirlo): restituisce anche alcuni punti lungo i bordi.

MODIFICA: Amro's answer sembra abbastanza elegante ed efficiente. Ma potrebbero esserci più "angoli" in ogni angolo reale poiché i picchi non sono unici. Potrei raggrupparli in base a θ e fare la media degli "angoli" attorno a un angolo reale, ma il problema principale è l'uso di order(1:10).

È il numero 10 sufficiente per tenere conto di tutti gli angoli o questo esclude un "angolo" in un angolo reale?

Answer 1

Questo è un po 'simile a quello suggerito da @AndyL. Comunque sto usando la firma del confine in coordinate polari invece della tangente.

Nota che comincio estraendo i bordi, ottenendo il limite, quindi convertendolo in firma.Finalmente troviamo i punti sul confine più lontani dal centroide, quei punti costituiscono gli angoli trovati. (In alternativa possiamo anche rilevare i picchi nella firma per gli angoli).

Quello che segue è una completa implementazione:

I = imread('oxyjj.png'); 
if ndims(I)==3 
    I = rgb2gray(I); 
end 
subplot(221), imshow(I), title('org') 

%%# Process Image 
%# edge detection 
BW = edge(I, 'sobel'); 
subplot(222), imshow(BW), title('edge') 

%# dilation-erosion 
se = strel('disk', 2); 
BW = imdilate(BW,se); 
BW = imerode(BW,se); 
subplot(223), imshow(BW), title('dilation-erosion') 

%# fill holes 
BW = imfill(BW, 'holes'); 
subplot(224), imshow(BW), title('fill') 

%# get boundary 
B = bwboundaries(BW, 8, 'noholes'); 
B = B{1}; 

%%# boudary signature 
%# convert boundary from cartesian to ploar coordinates 
objB = bsxfun(@minus, B, mean(B)); 
[theta, rho] = cart2pol(objB(:,2), objB(:,1)); 

%# find corners 
%#corners = find(diff(diff(rho)>0) < 0);  %# find peaks 
[~,order] = sort(rho, 'descend'); 
corners = order(1:10); 

%# plot boundary signature + corners 
figure, plot(theta, rho, '.'), hold on 
plot(theta(corners), rho(corners), 'ro'), hold off 
xlim([-pi pi]), title('Boundary Signature'), xlabel('\theta'), ylabel('\rho') 

%# plot image + corners 
figure, imshow(BW), hold on 
plot(B(corners,2), B(corners,1), 's', 'MarkerSize',10, 'MarkerFaceColor','r') 
hold off, title('Corners') 

---

EDIT: In risposta al commento di Jacob, dovrei spiegare che ho cercato di trovare le vette nella firma usando derivate primo/secondo, ma finì per prendere gli N-points più lontani. 10 era solo un valore ad-hoc, e sarebbe difficile da generalizzare (ho provato a prendere 4 come numero di curve, ma non le copriva tutte). Penso che valga l'idea di raggrupparli per rimuovere i duplicati.

Per quanto vedo io, il problema con il 1 ° approccio è stato che se si traccia rho senza prendere θ in considerazione, si otterrà una forma diversa (non gli stessi picchi), dal momento che la velocità con cui tracciare il confine è diverso e dipende dalla curvatura. Se potessimo capire come normalizzare in questo modo, possiamo ottenere risultati più accurati usando derivati.

Answer 2

Mi piace risolvere questo problema lavorando con un limite, perché riduce questo da un problema 2D a un problema 1D.

Utilizzare bwtraceboundary() dal toolkit di elaborazione immagini per estrarre un elenco di punti sul limite. Quindi converti il ​​limite in una serie di vettori tangenti (ci sono diversi modi per farlo, un modo sarebbe quello di eseguire il sottotitolo del punto lungo il confine dal punto i+delta.) Una volta che hai un elenco di vettori, prendere il prodotto punto di vettori adiacenti. I quattro punti con i prodotti con i punti più piccoli sono i tuoi angoli!

Se si desidera che l'algoritmo funzioni su poligoni con un numero arbitrario di vertici, cercare semplicemente i prodotti punto che sono un certo numero di deviazioni standard al di sotto del punto mediano.

Answer 3

Ho deciso di utilizzare uno Harris corner detector (qui è un more formal description) per ottenere gli angoli. Ciò può essere implementato come segue:

%% Constants 
Window = 3; 
Sigma = 2; 
K = 0.05; 
nCorners = 4; 

%% Derivative masks 
dx = [-1 0 1; -1 0 1; -1 0 1]; 
dy = dx'; %SO code color fix ' 

%% Find the image gradient 
% Mask is the binary image of the quadrilateral 
Ix = conv2(double(Mask),dx,'same'); 
Iy = conv2(double(Mask),dy,'same'); 

%% Use a gaussian windowing function and compute the rest 
Gaussian = fspecial('gaussian',Window,Sigma); 
Ix2 = conv2(Ix.^2, Gaussian, 'same'); 
Iy2 = conv2(Iy.^2, Gaussian, 'same'); 
Ixy = conv2(Ix.*Iy, Gaussian, 'same');  

%% Find the corners 
CornerStrength = (Ix2.*Iy2 - Ixy.^2) - K*(Ix2 + Iy2).^2; 
[val ind] = sort(CornerStrength(:),'descend');  
[Ci Cj] = ind2sub(size(CornerStrength),ind(1:nCorners)); 

%% Display 
imshow(Mask,[]); 
hold on; 
plot(Cj,Ci,'r*'); 

Qui, il problema con molteplici angoli grazie alla funzione di windowing gaussiana che leviga il cambiamento di intensità. Di seguito, c'è una versione ingrandita di un angolo con la mappa di colori hot.

Answer 4

Se avete la Image Processing Toolbox, c'è una funzione chiamata cornermetric che può implementare un rilevatore di angolo Harris o il metodo autovalore minimo Shi e Tomasi di. Questa funzione è presente dalla versione 6.2 di Image Processing Toolbox (versione MATLAB R2008b).

Utilizzando questa funzione, ho trovato un approccio leggermente diverso dalle altre risposte. La soluzione di seguito si basa sull'idea che un'area circolare centrata su ciascun punto d'angolo "vero" si sovrapporrà al poligono di una quantità inferiore rispetto a un'area circolare centrata su un punto d'angolo errato che si trova effettivamente sul bordo. Questa soluzione può anche gestire casi in cui sono rilevati più punti nello stesso angolo ...

Il primo passo è quello di caricare i dati:

rawImage = imread('oxyjj.png'); 
rawImage = rgb2gray(rawImage(7:473, 9:688, :)); % Remove the gray border 
subplot(2, 2, 1); 
imshow(rawImage); 
title('Raw image'); 

Avanti, calcolare l'angolo metriche utilizzando cornermetric. Si noti che sto mascherando la metrica dell'angolo con il poligono originale, in modo che cerchiamo i punti d'angolo che sono all'interno del poligono (vale a dire cercando di trovare i pixel angolari del poligono). imregionalmax viene quindi utilizzato per trovare i massimi locali. Dato che puoi avere cluster di più di 1 pixel con la stessa metrica angolare, aggiungo rumore ai massimi e ricalcolo in modo da ottenere solo 1 pixel in ciascuna regione massima.Ogni regione massima viene quindi all'etichetta bwlabel:

cornerImage = cornermetric(rawImage).*(rawImage > 0); 
maxImage = imregionalmax(cornerImage); 
noise = rand(nnz(maxImage), 1); 
cornerImage(maxImage) = cornerImage(maxImage)+noise; 
maxImage = imregionalmax(cornerImage); 
labeledImage = bwlabel(maxImage); 

Le regioni marcate vengono quindi dilatate (usando imdilate) con un elemento discoidale strutturazione (creati utilizzando strel):

diskSize = 5; 
dilatedImage = imdilate(labeledImage, strel('disk', diskSize)); 
subplot(2, 2, 2); 
imshow(dilatedImage); 
title('Dilated corner points'); 

Ora che l'etichetta le regioni d'angolo sono state dilatate, si sovrappongono parzialmente al poligono originale. Le regioni su un lato del poligono avranno una sovrapposizione di circa il 50%, mentre le regioni che si trovano su un angolo avranno una sovrapposizione di circa il 25%. La funzione regionprops può essere utilizzata per trovare le aree di sovrapposizione per ogni regione etichettata, e le 4 regioni che hanno la minor quantità di sovrapposizione può quindi essere considerato come i veri angoli:

maskImage = dilatedImage.*(rawImage > 0);  % Overlap with the polygon 
stats = regionprops(maskImage, 'Area');   % Compute the areas 
[sortedValues, index] = sort([stats.Area]);  % Sort in ascending order 
cornerLabels = index(1:4);      % The 4 smallest region labels 
maskImage = ismember(maskImage, cornerLabels); % Mask of the 4 smallest regions 
subplot(2, 2, 3); 
imshow(maskImage); 
title('Regions of minimal overlap'); 

E ora siamo in grado di ottenere il coordinate pixel degli angoli utilizzando find e ismember:

[r, c] = find(ismember(labeledImage, cornerLabels)); 
subplot(2, 2, 4); 
imshow(rawImage); 
hold on; 
plot(c, r, 'r+', 'MarkerSize', 16, 'LineWidth', 2); 
title('Corner points'); 

ed ecco una prova con una regione a forma di diamante:

0.123.

Answer 5

Ecco un esempio utilizzando Ruby e HornetsEye. Fondamentalmente il programma crea un istogramma dell'orientamento del gradiente Sobel quantizzato per trovare gli orientamenti dominanti. Se vengono rilevati quattro orientamenti dominanti, le linee vengono adattate e si presume che le intersezioni tra linee adiacenti siano gli angoli del rettangolo proiettato.

#!/usr/bin/env ruby 
require 'hornetseye' 
include Hornetseye 
Q = 36 
img = MultiArray.load_ubyte 'http://imgur.com/oxyjj.png' 
dx, dy = 8, 6 
box = [ dx ... 688, dy ... 473 ] 
crop = img[ *box ] 
crop.show 
s0, s1 = crop.sobel(0), crop.sobel(1) 
mag = Math.sqrt s0 ** 2 + s1 ** 2 
mag.normalise.show 
arg = Math.atan2 s1, s0 
msk = mag >= 500 
arg_q = ((arg.mask(msk)/Math::PI + 1) * Q/2).to_int % Q 
hist = arg_q.hist_weighted Q, mag.mask(msk) 
segments = (hist >= hist.max/4).components 
lines = arg_q.map segments 
lines.unmask(msk).normalise.show 
if segments.max == 4 
    pos = MultiArray.scomplex *crop.shape 
    pos.real = MultiArray.int(*crop.shape).indgen! % crop.shape[0] 
    pos.imag = MultiArray.int(*crop.shape).indgen!/crop.shape[0] 
    weights = lines.hist(5).major 1.0 
    centre = lines.hist_weighted(5, pos.mask(msk))/weights 
    vector = pos.mask(msk) - lines.map(centre) 
    orientation = lines.hist_weighted(5, vector ** 2) ** 0.5 
    corner = Sequence[ *(0 ... 4).collect do |i| 
    i1, i2 = i + 1, (i + 1) % 4 + 1 
    l1, a1, l2, a2 = centre[i1], orientation[i1], centre[i2], orientation[i2] 
    (l1 * a1.conj * a2 - l2 * a1 * a2.conj - 
     l1.conj * a1 * a2 + l2.conj * a1 * a2)/
     (a1.conj * a2 - a1 * a2.conj) 
    end ] 
    result = MultiArray.ubytergb(*img.shape).fill! 128 
    result[ *box ] = crop 
    corner.to_a.each do |c| 
    result[ c.real.to_i + dx - 1 .. c.real.to_i + dx + 1, 
      c.imag.to_i + dy - 1 .. c.imag.to_i + dy + 1 ] = RGB 255, 0, 0 
    end 
    result.show 
end