Qual è un modo semplice per calcolare la sovrapposizione tra un'immagine e un poligono?

Question

Ho un poligono chiuso non autointersecante. I suoi vertici vengono salvati in due vettori X e Y. Infine i valori di X e Y sono rilegati tra 0 e 22.

Mi piacerebbe costruire una matrice di dimensione 22x22 e impostare il valore di ciascun bin uguale a vero se parte del poligono si sovrappone a quel bin, altrimenti falso.

Il mio primo pensiero era di generare una griglia di punti definita con [a, b] = meshgrid(1:22) e quindi di utilizzare inpolygon per determinare quali punti della griglia erano nel poligono.

[a b] = meshgrid(1:22); 
inPoly1 = inpolygon(a,b,X,Y); 

Tuttavia questo solo restituisce vero se se il centro del cassone è contenuta nel poligono, cioè restituisce la forma rosso nell'immagine sottostante. Tuttavia quale bisogno è più lungo la linea della forma verde (anche se è ancora una soluzione incompleta).

Per ottenere il blob verde ho eseguito quattro chiamate allo inpolygon. Per ogni confronto ho spostato la griglia dei punti NE, NW, SE o SW di 1/2. Questo è equivalente al test se gli angoli di un raccoglitore si trovano nel poligono.

inPoly2 = inpolygon(a-.5,b-.5,X,Y) | inpolygon(a+.5,b-.5,X,Y) | inpolygon(a-.5,b+5,X,Y) | inpolygon(a+.5,b+.5,X,Y); 

Mentre questo mi forniscono una soluzione parziale non riesce nel caso in cui un vertice è contiene in un bidone, ma nessuno degli angoli bin sono.

C'è un modo più diretto di attaccare questo problema, preferibilmente una soluzione che produce codice più leggibile?

Questo terreno è stato disegnato con:

imagesc(inPoly1 + inPoly2); hold on; 
line(a, b, 'w.'); 
line(X, Y, 'y); 

Answer 1

Un suggerimento è utilizzare la funzione polybool (non disponibile in 2008b o precedenti). Trova l'intersezione di due poligoni e restituisce i vertici risultanti (o un vettore vuoto se non esistono vertici). Per usarlo qui, iteriamo (usando arrayfun) su tutti i quadrati del controllo della griglia per vedere se l'argomento di output su polybool è vuoto (ad esempio non si sovrappone).

N=22; 
sqX = repmat([1:N]',1,N); 
sqX = sqX(:); 
sqY = repmat(1:N,N,1); 
sqY = sqY(:); 

intersects = arrayfun((@(xs,ys) ... 
     (~isempty(polybool('intersection',X,Y,[xs-1 xs-1 xs xs],[ys-1 ys ys ys-1])))),... 
     sqX,sqY); 

intersects = reshape(intersects,22,22); 

Ecco l'immagine risultante:

Codice per il tracciato:

imagesc(.5:1:N-.5,.5:1:N-.5,intersects'); 
hold on; 
plot(X,Y,'w'); 
for x = 1:N 
    plot([0 N],[x x],'-k'); 
    plot([x x],[0 N],'-k'); 
end 
hold off; 

Answer 2

Che ne dite di questo algoritmo pseudocodice:

For each pair of points p1=p(i), p2=p(i+1), i = 1..n-1 
    Find the line passing through p1 and p2 
    Find every tile this line intersects // See note 
    Add intersecting tiles to the list of contained tiles 

Find the red area using the centers of each tile, and add these to the list of contained tiles 

Nota: Questa linea avrà un po 'di sforzo da implementare, ma Penso che ci sia un algoritmo abbastanza semplice e ben noto per questo.

Inoltre, se si utilizzasse .NET, definirei semplicemente un rettangolo corrispondente a ciascun riquadro della griglia e quindi vedere quali intersecano il poligono. Non so se il controllo dell'intersezione sia facile in Matlab, comunque.

Answer 3

prima cosa definire un cerchio a bassa risoluzione per questo esempio

X=11+cos(linspace(0,2*pi,10))*5; 
Y=11+sin(linspace(0,2.01*pi,10))*5; 

Come il vostro esempio si adatta con in una griglia di ~ 22 unità. Quindi, seguendo il tuo comando, dichiariamo un meshgrid e controlliamo se i punti sono nel poligono.

L'unica differenza è che, laddove la soluzione originale ha impiegato uno step, questa griglia può richiedere passi più piccoli. E, infine, per includere qualsiasi cosa entro i "bordi" delle piazze

inPolyFull=unique(round([a(inPoly1) b(inPoly1)]) ,'rows'); 

Il round prende semplicemente la nostra griglia di alta risoluzione e arrotonda i punti in modo appropriato alle loro più vicine equivalenti a bassa risoluzione. Quindi rimuoviamo tutti i duplicati in uno stile vettoriale o in modalità pair-wise chiamando lo unique con il qualificatore 'rows'.E il gioco è fatto

Per visualizzare il risultato,

[aOrig bOrig] = meshgrid(1:22); 
imagesc(1:stepSize:22,1:stepSize:22,inPoly1); hold on; 
plot(X,Y,'y'); 
plot(aOrig,bOrig,'k.'); 
plot(inPolyFull(:,1),inPolyFull(:,2),'w.'); hold off; 

Modifica della stepSize ha l'effetto previsto di migliorare il risultato a costo di velocità e della memoria.

Se è necessario il risultato di essere nello stesso formato come l'inPoly2 nel tuo esempio, è possibile utilizzare

inPoly2=zeros(22); 
inPoly2(inPolyFull(:,1),inPolyFull(:,2))=1 

Speranza che aiuta. Posso pensare ad altri modi per farlo, ma questo sembra il più semplice.

Answer 4

Io suggerirei di usare poly2mask in Image Processing Toolbox, che fa più o meno quello tu vuoi, io penso, e anche più o meno ciò che tu e Salain avete suggerito.

Answer 5

leggero miglioramento

In primo luogo, per semplificare la vostra "soluzione parziale" - quello che stai facendo è solo in cerca agli angoli. Se invece di considerare la griglia di punti 22x22, potresti considerare la griglia degli angoli 23x23 (che sarà compensata dalla griglia più piccola di (-0,5, -0,5). Una volta ottenuto ciò, puoi segnare i punti sulla griglia 22x22 che hanno almeno un angolo nel poligono

soluzione completa:..

Tuttavia, quello che stai veramente cercando è se il poligono interseca con la scatola 1x1 che circonda ogni pixel Questo non significa necessariamente includere uno qualsiasi degli angoli, ma è richiede che il poligono intersechi uno dei quattro lati della casella

Un modo si potrebbe trovare i pixel in cui il poligono interseca con la scatola contenente è con il seguente algoritmo:

For each pair of adjacent points in the polygon, calling them pA and pB: 
    Calculate rounded Y-values: Round(pA.y) and Round(pB.y) 
    For each horizontal pixel edge between these two values: 
     * Solve the simple linear equation to find out at what X-coordinate 
      the line between pA and pB crosses this edge 
     * Round the X-coordinate 
     * Use the rounded X-coordinate to mark the pixels above and below 
      where it crosses the edge 
    Do a similar thing for the other axis 

Così, per esempio, diciamo che stiamo guardando pA = (1, 1) e pB = (2, 3).

• In primo luogo, abbiamo calcolato i valori Y arrotondati: 1 e 3.
• Poi, guardiamo i bordi dei pixel tra questi valori: y = 1.5 e y = 2.5 (bordi di pixel sono semi-offset da pixel
• Per ciascuno di questi, risolviamo l'equazione lineare per trovare dove pA ->pB interseca con la nostra . bordi Questo ci dà: x = 1.25, y = 1.5 e x = 1.75, y = 2.5
• per ciascuno di questi incroci, prendiamo il valore X arrotondato, e usarlo per marcare i pixel entrambi i lati del bordo
  • x = 1.25 è arrotondato a 1.. (per il bordo y = 1.5). Pertanto possiamo contrassegnare il file pixel a (1, 1) e (1, 2) come parte del nostro set.
  • x = 1.75 è arrotondato a 2 (per il bordo y = 2.5). Pertanto, possiamo contrassegnare i pixel a (2, 2) e (2, 3).

Ecco i bordi orizzontali curato. Quindi, diamo un'occhiata a quelle verticali:

• Per prima cosa calcolare gli X-valori arrotondati: 1 e 2
• Poi, guardiamo i bordi dei pixel. Qui ce n'è uno solo: x = 1.5.
• Per questo bordo, troviamo il punto in cui incontra la linea pA ->pB. Questo ci dà x = 1.5, y = 2.
• Per questa intersezione, prendiamo il valore Y arrotondato, e usarlo per marcare pixel lati del bordo:
  • y = 2 viene arrotondato a 2. Possiamo quindi segnare i pixel a (1, 2) e (2, 2).

Fatto!

Bene, una specie di. Questo ti darà i bordi, ma non riempirà il corpo del poligono. Tuttavia, puoi combinare questi con i tuoi precedenti (rossi) risultati per ottenere il set completo.

Answer 6

Beh, immagino di essere in ritardo, anche se a rigor di termini il tempo di taglia è stato fino a domani;). Ma ecco il mio tentativo. Innanzitutto, una funzione che contrassegna le celle che contengono/toccano un punto. Data una griglia strutturato con spaziatura lx, ly, e un insieme di punti di coordinate (Xp, Yp), impostare cellule contenenti:

function cells = mark_cells(lx, ly, Xp, Yp, cells) 

% Find cell numbers to which points belong. 
% Search by subtracting point coordinates from 
% grid coordinates and observing the sign of the result. 
% Points lying on edges/grid points are assumed 
% to belong to all surrounding cells. 

sx=sign(bsxfun(@minus, lx, Xp')); 
sy=sign(bsxfun(@minus, ly, Yp')); 
cx=diff(sx, 1, 2); 
cy=diff(sy, 1, 2); 

% for every point, mark the surrounding cells 
for i=1:size(cy, 1) 
    cells(find(cx(i,:)), find(cy(i,:)))=1; 
end 
end 

Ora, il resto del codice. Per ogni segmento nel poligono (devi attraversare i segmenti uno per uno), intersecare il segmento con le linee della griglia. L'intersezione viene eseguita attentamente, per le linee orizzontali e verticali separatamente, utilizzando le coordinate del punto della griglia fornite per evitare imprecisioni numeriche. Per i punti di intersezione trovati io chiamo mark_cells per marcare le cellule circostanti per 1:

% example grid 
nx=21; 
ny=51; 
lx = linspace(0, 1, nx); 
ly = linspace(0, 1, ny); 
dx=1/(nx-1); 
dy=1/(ny-1); 
cells = zeros(nx-1, ny-1); 

% for every line in the polygon... 
% Xp and Yp contain start-end points of a single segment 
Xp = [0.15 0.61]; 
Yp = [0.1 0.78]; 

% line equation 
slope = diff(Yp)/diff(Xp); 
inter = Yp(1) - (slope*Xp(1)); 

if isinf(slope) 
    % SPECIAL CASE: vertical polygon segments 
    % intersect horizontal grid lines 
    ymax = 1+floor(max(Yp)/dy); 
    ymin = 1+ceil(min(Yp)/dy); 
    x=repmat(Xp(1), 1, ymax-ymin+1); 
    y=ly(ymin:ymax); 
    cells = mark_cells(lx, ly, x, y, cells); 
else 
    % SPECIAL CASE: not horizontal polygon segments 
    if slope ~= 0 
     % intersect horizontal grid lines 
     ymax = 1+floor(max(Yp)/dy); 
     ymin = 1+ceil(min(Yp)/dy); 
     xmax = (ly(ymax)-inter)/slope; 
     xmin = (ly(ymin)-inter)/slope; 
     % interpolate in x... 
     x=linspace(xmin, xmax, ymax-ymin+1); 
     % use exact grid point y-coordinates! 
     y=ly(ymin:ymax); 
     cells = mark_cells(lx, ly, x, y, cells); 
    end 

    % intersect vertical grid lines 
    xmax = 1+floor(max(Xp)/dx); 
    xmin = 1+ceil(min(Xp)/dx); 
    % interpolate in y... 
    ymax = inter+slope*lx(xmax); 
    ymin = inter+slope*lx(xmin); 
    % use exact grid point x-coordinates! 
    x=lx(xmin:xmax); 
    y=linspace(ymin, ymax, xmax-xmin+1); 
    cells = mark_cells(lx, ly, x, y, cells); 
end 

uscita per l'esempio rete/segmento:

Passeggiando per tutti i segmenti del poligono ti dà il poligono 'Halo'.Infine, l'interno del poligono viene ottenuto utilizzando la funzione standard inpolygon. Fammi sapere se hai bisogno di maggiori dettagli sul codice.