Come posso determinare se un punto 2D si trova all'interno di un poligono?
https://stackoverflow.com/questions/217578
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Sto provando a creare un punto 2D veloce all'interno dell'algoritmo poligonale, da utilizzare in hit-test (ad es. Polygon.contains(p:Point)). Suggerimenti per tecniche efficaci sarebbero apprezzati.
Soluzione
Per la grafica, preferirei non preferire numeri interi. Molti sistemi usano numeri interi per la pittura dell'interfaccia utente (dopo tutto i pixel sono ints), ma macOS ad esempio usa float per tutto. macOS conosce solo i punti e un punto può essere tradotto in un pixel, ma a seconda della risoluzione del monitor, potrebbe tradursi in qualcos'altro. Sugli schermi retina mezzo punto (0,5 / 0,5) è pixel. Tuttavia, non ho mai notato che le interfacce utente di macOS sono significativamente più lente di altre interfacce utente. Dopo che tutte le API 3D (OpenGL o Direct3D) funzionano anche con float e librerie grafiche moderne sfruttano molto spesso l'accelerazione GPU.
Ora hai detto che la velocità è la tua principale preoccupazione, ok, andiamo per la velocità. Prima di eseguire qualsiasi algoritmo sofisticato, eseguire innanzitutto un semplice test. Crea un riquadro di selezione allineato all'asse attorno al poligono. Questo è molto semplice, veloce e può già essere sicuro per molti calcoli. Come funziona? Scorri su tutti i punti del poligono e trova i valori min / max di X e Y.
es. hai i punti (9/1), (4/3), (2/7), (8/2), (3/6). Ciò significa che Xmin è 2, Xmax è 9, Ymin è 1 e Ymax è 7. Un punto esterno al rettangolo con i due bordi (2/1) e (9/7) non può essere all'interno del poligono.
// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
Questo è il primo test da eseguire per qualsiasi punto. Come puoi vedere, questo test è ultra veloce ma è anche molto grossolano. Per gestire i punti che si trovano all'interno del rettangolo di delimitazione, abbiamo bisogno di un algoritmo più sofisticato. Ci sono un paio di modi in cui questo può essere calcolato. Quale metodo funziona dipende anche dal fatto se il poligono può avere buchi o sarà sempre solido. Ecco alcuni esempi di quelli solidi (uno convesso, uno concavo):
Ed eccone uno con un buco:
Quello verde ha un buco nel mezzo!
L'algoritmo più semplice, in grado di gestire tutti e tre i casi sopra ed è ancora piuttosto veloce, si chiama ray casting . L'idea dell'algoritmo è piuttosto semplice: disegna un raggio virtuale da qualsiasi punto esterno al poligono fino al punto e conta la frequenza con cui colpisce un lato del poligono. Se il numero di colpi è pari, è al di fuori del poligono, se è dispari, è dentro.
L'algoritmo numero di avvolgimento sarebbe un'alternativa, è più preciso per i punti molto vicini a una linea poligonale ma è anche molto più lento. Il lancio del raggio potrebbe non riuscire per punti troppo vicini a un lato del poligono a causa di problemi di precisione e arrotondamento in virgola mobile limitati, ma in realtà non è certo un problema, come se un punto si trovi così vicino a un lato, spesso visivamente non è nemmeno possibile per un spettatore per riconoscere se è già dentro o ancora fuori.
Hai ancora il riquadro di selezione sopra, ricordi? Scegli un punto fuori dal riquadro di delimitazione e usalo come punto di partenza per il tuo raggio. Per esempio. il punto (Xmin - e/p.y) è sicuramente fuori dal poligono.
Ma cos'è e? Bene, v1x1/v1y1 (in realtà epsilon) dà al riquadro di delimitazione un riempimento . Come ho detto, il ray tracing fallisce se iniziamo troppo vicino a una linea poligonale. Poiché il rettangolo di selezione potrebbe eguagliare il poligono (se il poligono è un rettangolo allineato all'asse, il rettangolo di selezione è uguale al poligono stesso!), Abbiamo bisogno di un po 'di riempimento per renderlo sicuro, tutto qui. Quanto dovresti scegliere v1x2/v1y2? Non troppo grosso. Dipende dalla scala del sistema di coordinate utilizzata per il disegno. Se la larghezza del passo del pixel è 1.0, scegli semplicemente 1.0 (eppure 0.1 avrebbe funzionato pure)
Ora che abbiamo il raggio con le sue coordinate di inizio e fine, il problema si sposta da " è il punto all'interno del poligono " a " con quale frequenza il raggio interseca un lato poligonale " ;. Pertanto non possiamo semplicemente lavorare con i punti poligonali come prima, ora abbiamo bisogno dei lati effettivi. UNil lato è sempre definito da due punti.
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
Devi testare il raggio su tutti i lati. Considera il raggio come un vettore e ogni lato è un vettore. Il raggio deve colpire ogni lato esattamente una volta o mai. Non può colpire due volte lo stesso lato. Due linee nello spazio 2D si intersecano sempre esattamente una volta, a meno che non siano parallele, nel qual caso non si intersecano mai. Tuttavia, poiché i vettori hanno una lunghezza limitata, due vettori potrebbero non essere paralleli e non si intersecano mai perché sono troppo corti per incontrarsi mai.
// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
// Test if current side intersects with ray.
// If yes, intersections++;
}
if ((intersections & 1) == 1) {
// Inside of polygon
} else {
// Outside of polygon
}
Finora tutto bene, ma come si verifica se due vettori si intersecano? Ecco un po 'di codice C (non testato), che dovrebbe fare il trucco:
#define NO 0
#define YES 1
#define COLLINEAR 2
int areIntersecting(
float v1x1, float v1y1, float v1x2, float v1y2,
float v2x1, float v2y1, float v2x2, float v2y2
) {
float d1, d2;
float a1, a2, b1, b2, c1, c2;
// Convert vector 1 to a line (line 1) of infinite length.
// We want the line in linear equation standard form: A*x + B*y + C = 0
// See: http://en.wikipedia.org/wiki/Linear_equation
a1 = v1y2 - v1y1;
b1 = v1x1 - v1x2;
c1 = (v1x2 * v1y1) - (v1x1 * v1y2);
// Every point (x,y), that solves the equation above, is on the line,
// every point that does not solve it, is not. The equation will have a
// positive result if it is on one side of the line and a negative one
// if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
// 2 into the equation above.
d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
d2 = (a1 * v2x2) + (b1 * v2y2) + c1;
// If d1 and d2 both have the same sign, they are both on the same side
// of our line 1 and in that case no intersection is possible. Careful,
// 0 is a special case, that's why we don't test ">=" and "<=",
// but "<" and ">".
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// The fact that vector 2 intersected the infinite line 1 above doesn't
// mean it also intersects the vector 1. Vector 1 is only a subset of that
// infinite line 1, so it may have intersected that line before the vector
// started or after it ended. To know for sure, we have to repeat the
// the same test the other way round. We start by calculating the
// infinite line 2 in linear equation standard form.
a2 = v2y2 - v2y1;
b2 = v2x1 - v2x2;
c2 = (v2x2 * v2y1) - (v2x1 * v2y2);
// Calculate d1 and d2 again, this time using points of vector 1.
d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
d2 = (a2 * v1x2) + (b2 * v1y2) + c2;
// Again, if both have the same sign (and neither one is 0),
// no intersection is possible.
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// If we get here, only two possibilities are left. Either the two
// vectors intersect in exactly one point or they are collinear, which
// means they intersect in any number of points from zero to infinite.
if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;
// If they are not collinear, they must intersect in exactly one point.
return YES;
}
I valori di input sono i due endpoint del vettore 1 (v2x1/v2y1 e v2x2/v2y2) e del vettore 2 (YES e NO). Quindi hai 2 vettori, 4 punti, 8 coordinate. == 0.0f e < epsilon sono chiari. <=> aumenta le intersezioni, <=> non fa nulla.
Che dire di COLLINEAR? Significa che entrambi i vettori si trovano sulla stessa linea infinita, a seconda della posizione e della lunghezza, non si intersecano affatto o si intersecano in un numero infinito di punti. Non sono assolutamente sicuro di come gestire questo caso, non lo considererei come intersezione in entrambi i casi. Bene, questo caso è piuttosto raro nella pratica comunque a causa di errori di arrotondamento in virgola mobile; il codice migliore probabilmente non testerebbe per <=> ma invece per qualcosa come <=>, dove epsilon è un numero piuttosto piccolo.
Se devi testare un numero maggiore di punti, puoi certamente accelerare un po 'il tutto mantenendo in memoria le forme standard di equazione lineare dei lati del poligono, quindi non devi ricalcolarle ogni volta. Ciò consentirà di risparmiare due moltiplicazioni in virgola mobile e tre sottrazioni in virgola mobile su ogni test in cambio della memorizzazione di tre valori in virgola mobile per lato poligonale in memoria. È un tipico compromesso tra memoria e tempo di calcolo.
Ultimo ma non meno importante: se è possibile utilizzare l'hardware 3D per risolvere il problema, esiste un'alternativa interessante. Lascia che la GPU faccia tutto il lavoro per te. Crea una superficie pittorica fuori dallo schermo. Riempilo completamente con il colore nero. Ora lascia che OpenGL o Direct3D dipingano il tuo poligono (o anche tutti i tuoi poligoni se vuoi solo verificare se il punto si trova all'interno di uno di essi, ma non ti interessa quale) e riempi il poligono con un altro colore, ad es bianca. Per verificare se un punto si trova all'interno del poligono, ottenere il colore di questo punto dalla superficie del disegno. Questo è solo un recupero di memoria O (1).
Naturalmente questo metodo è utilizzabile solo se la superficie del disegno non deve essere enorme. Se non si adatta alla memoria della GPU, questo metodo è più lento rispetto alla CPU. Se dovrebbe essere enorme e la tua GPU supporta shader moderni, puoi comunque utilizzare la GPU implementando il ray casting mostrato sopra come shader GPU, il che è assolutamente possibile. Per un numero maggiore di poligoni o un numero elevato di punti da testare, questo pagherà, considerando che alcune GPU saranno in grado di testare da 64 a 256 punti in parallelo. Si noti tuttavia che il trasferimento di dati dalla CPU alla GPU e viceversa è sempre costoso, quindi per testare solo un paio di punti contro un paio di poligoni semplici, in cui i punti o i poligoni sono dinamici e cambieranno frequentemente, un approccio GPU raramente pagherà off.
Altri suggerimenti
Penso che il seguente pezzo di codice sia la migliore soluzione (presa da qui ):
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
Argomenti
- annulla : numero di vertici nel poligono. Se ripetere il primo vertice alla fine è stato discusso nell'articolo di cui sopra.
- vertx, verty : array contenenti le coordinate xey dei vertici del poligono.
- testx, testy : coordinate X e y del punto di test.
È allo stesso tempo corto ed efficiente e funziona sia per poligoni convessi che concavi. Come suggerito prima, dovresti prima controllare il rettangolo di delimitazione e trattare i fori poligonali separatamente.
L'idea alla base è piuttosto semplice. L'autore lo descrive come segue:
Eseguo un raggio semi-infinito in orizzontale (aumentando x, fisso y) fuori dal punto di prova e conto quanti spigoli attraversa. Ad ogni incrocio, il raggio cambia tra interno ed esterno. Questo è chiamato teorema della curva giordana.
La variabile c passa da 0 a 1 e da 1 a 0 ogni volta che il raggio orizzontale attraversa un bordo. Quindi, in sostanza, sta monitorando se il numero di spigoli incrociati è pari o dispari. 0 significa pari e 1 significa dispari.
Ecco una versione C # della risposta data da nirg , che proviene da questo professore RPI . Si noti che l'uso del codice da quella fonte RPI richiede l'attribuzione.
Nella parte superiore è stato aggiunto un segno di spunta. Tuttavia, come sottolinea James Brown, il codice principale è quasi veloce quanto il controllo del riquadro di selezione stesso, quindi il controllo del riquadro di selezione può effettivamente rallentare l'operazione complessiva, nel caso in cui la maggior parte dei punti che stai controllando si trovino all'interno del riquadro di selezione . Quindi potresti lasciare il riquadro di selezione check out o un'alternativa sarebbe quella di precompilare i box di delimitazione dei tuoi poligoni se non cambiano forma troppo spesso.
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
minX = Math.Min( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
{
return false;
}
// http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
{
if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
Ecco una variante JavaScript della risposta di M. Katz basata sull'approccio di Nirg:
function pointIsInPoly(p, polygon) {
var isInside = false;
var minX = polygon[0].x, maxX = polygon[0].x;
var minY = polygon[0].y, maxY = polygon[0].y;
for (var n = 1; n < polygon.length; n++) {
var q = polygon[n];
minX = Math.min(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
Calcola la somma orientata degli angoli tra il punto p e ciascuno degli apici del poligono. Se l'angolo orientato totale è di 360 gradi, il punto si trova all'interno. Se il totale è 0, il punto è esterno.
Mi piace meglio questo metodo perché è più robusto e meno dipendente dalla precisione numerica.
I metodi che calcolano l'uniformità del numero di intersezioni sono limitati perché puoi "colpire" un apice durante il calcolo del numero di intersezioni.
EDIT: A proposito, questo metodo funziona con poligoni concavi e convessi.
EDIT: di recente ho trovato un intero articolo di Wikipedia sull'argomento.
L'articolo Eric Haines citato da bobobobo è davvero eccellente. Particolarmente interessanti sono le tabelle che confrontano le prestazioni degli algoritmi; il metodo di somma degli angoli è davvero pessimo rispetto agli altri. È anche interessante notare che le ottimizzazioni come l'uso di una griglia di ricerca per suddividere ulteriormente il poligono in & Quot; in & Quot; e " out " settori possono rendere il test incredibilmente veloce anche su poligoni con > 1000 lati.
Comunque, sono i primi giorni ma il mio voto va al " incroci " il metodo, che è praticamente quello che descrive Mecki penso. Tuttavia l'ho trovato più succinta descritto e codificato da David Bourke . Adoro il fatto che non sia necessaria una vera trigonometria, che funzioni per convesso e concavo, e si comporta in modo ragionevole all'aumentare del numero di lati.
A proposito, ecco una delle tabelle delle prestazioni dall'articolo di Eric Haines per interesse, test su poligoni casuali.
number of edges per polygon
3 4 10 100 1000
MacMartin 2.9 3.2 5.9 50.6 485
Crossings 3.1 3.4 6.8 60.0 624
Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787
Triangle Fan 1.2 2.1 7.3 85.4 865
Barycentric 2.1 3.8 13.8 160.7 1665
Angle Summation 56.2 70.4 153.6 1403.8 14693
Grid (100x100) 1.5 1.5 1.6 2.1 9.8
Grid (20x20) 1.7 1.7 1.9 5.7 42.2
Bins (100) 1.8 1.9 2.7 15.1 117
Bins (20) 2.1 2.2 3.7 26.3 278
Questa domanda è così interessante. Ho un'altra idea realizzabile diversa dalle altre risposte di questo post.
L'idea è di usare la somma degli angoli per decidere se il bersaglio è dentro o fuori. Se il bersaglio si trova all'interno di un'area, la somma della forma angolare del bersaglio e ogni due punti del bordo saranno 360. Se il bersaglio è esterno, la somma non sarà 360. L'angolo ha direzione. Se l'angolo va indietro, l'angolo è negativo. È come calcolare il numero di avvolgimento .
Fai riferimento a questa immagine per ottenere una comprensione di base dell'idea:
Il mio algoritmo presuppone che la direzione positiva sia in senso orario. Ecco un potenziale input:
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
Quello che segue è il codice Python che implementa l'idea:
def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
a = border[i]
b = border[i + 1]
# calculate distance of vector
A = getDistance(a[0], a[1], b[0], b[1]);
B = getDistance(target[0], target[1], a[0], a[1])
C = getDistance(target[0], target[1], b[0], b[1])
# calculate direction of vector
ta_x = a[0] - target[0]
ta_y = a[1] - target[1]
tb_x = b[0] - target[0]
tb_y = b[1] - target[1]
cross = tb_y * ta_x - tb_x * ta_y
clockwise = cross < 0
# calculate sum of angles
if(clockwise):
degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
else:
degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
if(abs(round(degree) - 360) <= 3):
return True
return False
Versione rapida della risposta di nirg :
extension CGPoint {
func isInsidePolygon(vertices: [CGPoint]) -> Bool {
guard !vertices.isEmpty else { return false }
var j = vertices.last!, c = false
for i in vertices {
let a = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
Ho lavorato un po 'su questo argomento quando ero un ricercatore in Michael Stonebraker - tu sapete, il professore che ha ideato Ingres , PostgreSQL , ecc.
Ci siamo resi conto che il modo più veloce era fare prima un rettangolo di selezione perché è SUPER veloce. Se è fuori dal riquadro di delimitazione, è fuori. Altrimenti, fai il lavoro più duro ...
Se vuoi un ottimo algoritmo, guarda il codice sorgente PostgreSQL del progetto open source per il geo work ...
Voglio sottolineare, non abbiamo mai avuto alcuna comprensione della destra o della mano sinistra (anche espressa come " dentro " vs " outside " problema .. .
UPDATE
Il link di BKB ha fornito un buon numero di algoritmi ragionevoli. Stavo lavorando su problemi di scienze della terra e quindi avevo bisogno di una soluzione che funzioni in latitudine / longitudine e che abbia il problema peculiare della mano: l'area all'interno dell'area più piccola o più grande? La risposta è che & Quot; direzione & Quot; delle cose importanti: è per mancini o destrorsi e in questo modo puoi indicare entrambe le aree come " dentro " qualsiasi dato poligono. Come tale, il mio lavoro ha utilizzato la soluzione tre elencata in quella pagina.
Inoltre, il mio lavoro utilizzava funzioni separate per " sulla riga " test.
... Dal momento che qualcuno ha chiesto: abbiamo capito che i test del riquadro di selezione erano i migliori quando il numero di vertici è andato oltre un certo numero - fai un test molto veloce prima di fare il test più lungo se necessario ... Viene creato un riquadro di selezione da semplicemente prendendo la più grande x, la più piccola x, la più grande y e la più piccola y e mettendole insieme per fare quattro punti di una scatola ...
Un altro suggerimento per quelli che seguono: abbiamo fatto tutto il nostro più sofisticato e " light-dimming " calcolare in uno spazio della griglia tutto in punti positivi su un piano e quindi riproiettare nuovamente in " real " longitudine / latitudine, evitando così possibili errori di avvolgimento quando una linea incrociata 180 di longitudine e durante la manipolazione di regioni polari. Ha funzionato alla grande!
Mi piace molto la soluzione pubblicata da Nirg e modificata da bobobobo. L'ho appena reso javascript e un po 'più leggibile per il mio uso:
function insidePoly(poly, pointx, pointy) {
var i, j;
var inside = false;
for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
}
return inside;
}
La risposta di David Segond è praticamente la risposta generale standard e quella di Richard T è l'ottimizzazione più comune, sebbene ce ne siano altre. Altre forti ottimizzazioni si basano su soluzioni meno generali. Ad esempio, se stai per controllare lo stesso poligono con molti punti, la triangolazione del poligono può accelerare enormemente le cose in quanto vi sono numerosi algoritmi di ricerca TIN molto veloci. Un altro è se il poligono e i punti si trovano su un piano limitato a bassa risoluzione, ad esempio una visualizzazione a schermo, è possibile dipingere il poligono su un buffer di visualizzazione mappato in memoria in un determinato colore e controllare il colore di un dato pixel per vedere se si trova nei poligoni.
Come molte ottimizzazioni, queste si basano su casi specifici piuttosto che generali e producono vantaggi in base al tempo ammortizzato anziché al singolo utilizzo.
Lavorando in questo campo, ho trovato che la geometria di calcolo di Joeseph O'Rourkes in C 'ISBN 0-521-44034-3 era di grande aiuto.
La banale soluzione sarebbe quella di dividere il poligono in triangoli e premere test i triangoli come spiegato qui
Se il tuo poligono è CONVESSO , potrebbe esserci un approccio migliore. Guarda il poligono come una raccolta di linee infinite. Ogni linea che divide lo spazio in due. per ogni punto è facile dire se è da un lato o dall'altro lato della linea. Se un punto si trova sullo stesso lato di tutte le linee, allora si trova all'interno del poligono.
Mi rendo conto che questo è vecchio, ma ecco un algoritmo di ray casting implementato in Cocoa, nel caso qualcuno fosse interessato. Non sono sicuro che sia il modo più efficiente di fare le cose, ma potrebbe aiutare qualcuno.
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
BOOL result;
float aggregateX = 0; //I use these to calculate the centroid of the shape
float aggregateY = 0;
NSPoint firstPoint[1];
[currentPath elementAtIndex:0 associatedPoints:firstPoint];
float olderX = firstPoint[0].x;
float olderY = firstPoint[0].y;
NSPoint interPoint;
int noOfIntersections = 0;
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
aggregateX += points[0].x;
aggregateY += points[0].y;
}
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
//line equations in Ax + By = C form
float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;
float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);
float _A_BAR = olderY - points[0].y;
float _B_BAR = points[0].x - olderX;
float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);
float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
if (det != 0) {
//intersection points with the edges
float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
if (olderX <= points[0].x) {
//doesn't matter in which direction the ray goes, so I send it right-ward.
if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {
noOfIntersections++;
}
} else {
if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
noOfIntersections++;
}
}
}
olderX = points[0].x;
olderY = points[0].y;
}
if (noOfIntersections % 2 == 0) {
result = FALSE;
} else {
result = TRUE;
}
return result;
}
Versione obj-C della risposta di nirg con metodo di esempio per i punti di prova. La risposta di Nirg ha funzionato bene per me.
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
NSUInteger nvert = [vertices count];
NSInteger i, j, c = 0;
CGPoint verti, vertj;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
c = !c;
}
return (c ? YES : NO);
}
- (void)testPoint {
NSArray *polygonVertices = [NSArray arrayWithObjects:
[NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
[NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
[NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
[NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
nil
];
CGPoint tappedPoint = CGPointMake(23.0, 70.0);
if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
NSLog(@"YES");
} else {
NSLog(@"NO");
}
}
La versione C # della risposta di nirg è qui: condividerò semplicemente il codice. Potrebbe far risparmiare tempo a qualcuno.
public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
bool result = false;
int j = polygon.Count() - 1;
for (int i = 0; i < polygon.Count(); i++) {
if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
result = !result;
}
}
j = i;
}
return result;
}
Non c'è niente di più bello di una definizione induttiva di un problema. Per completezza qui hai una versione in prologo che potrebbe anche chiarire i pensieri dietro ray casting :
Basato sulla simulazione dell'algoritmo di semplicità in http: // www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
Alcuni predicati di supporto:
exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).
inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).
get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
L'equazione di una linea con 2 punti A e B (Linea (A, B)) è:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
È importante che la direzione di rotazione per la linea sia impostato in senso orario per i confini e in senso antiorario per i buchi. Verificheremo se il punto (X, Y), ovvero il punto testato è a sinistra mezzo piano della nostra linea (è una questione di gusti, potrebbe anche essere il lato destro, ma anche la direzione delle linee di confine deve essere cambiata quel caso), questo è per proiettare il raggio dal punto a destra (o sinistra) e riconoscere l'intersezione con la linea. Abbiamo scelto di proiettare il raggio in direzione orizzontale (di nuovo è una questione di gusti, potrebbe anche essere fatto in verticale con restrizioni simili), quindi abbiamo:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Ora dobbiamo sapere se il punto è sul lato sinistro (o destro) di solo il segmento di linea, non l'intero piano, quindi è necessario limitare la ricerca solo a questo segmento, ma da allora è facile essere all'interno del segmento solo un punto nella linea può essere più alto di Y nell'asse verticale. Dal momento che questa è una restrizione più forte deve essere il primo a controllare, quindi prendiamo prima solo quelle righe soddisfare questo requisito e quindi verificarne il possesso. Dal Giordano Teorema della curva qualsiasi raggio proiettato su un poligono deve intersecarsi in un numero pari di righe. Quindi abbiamo finito, lanceremo il raggio su a destra e poi ogni volta che interseca una linea, commuta il suo stato. Tuttavia, nella nostra implementazione, siamo in grado di verificare la lunghezza di il pacchetto di soluzioni che soddisfa le restrizioni indicate e decide il amicizia su di esso. per ogni riga del poligono è necessario eseguire questo.
is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA));
is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).
in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines).
traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).
% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
Versione Java:
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
.Net port:
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
int i, j, c = 0;
vertx.Add(1);
vertx.Add(2);
vertx.Add(1);
vertx.Add(4);
vertx.Add(4);
vertx.Add(1);
verty.Add(1);
verty.Add(2);
verty.Add(4);
verty.Add(4);
verty.Add(1);
verty.Add(1);
int nvert = 6; //Vértices del poligono
double testx = 2;
double testy = 5;
for (i = 0, j = nvert - 1; i < nvert; j = i++)
{
if (((verty[i] > testy) != (verty[j] > testy)) &&
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
c = 1;
}
}
VERSIONE VBA:
Nota: Ricorda che se il tuo poligono è un'area all'interno di una mappa, Latitudine / Longitudine sono valori Y / X rispetto a X / Y (Latitudine = Y, Longitudine = X) a causa di ciò che comprendo sono implicazioni storiche da lontano quando Longitudine non era una misura.
MODULO DI CLASSE: CPoint
Private pXValue As Double
Private pYValue As Double
'''''X Value Property'''''
Public Property Get X() As Double
X = pXValue
End Property
Public Property Let X(Value As Double)
pXValue = Value
End Property
'''''Y Value Property'''''
Public Property Get Y() As Double
Y = pYValue
End Property
Public Property Let Y(Value As Double)
pYValue = Value
End Property
MODULO:
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean
Dim i As Integer
Dim j As Integer
Dim q As Object
Dim minX As Double
Dim maxX As Double
Dim minY As Double
Dim maxY As Double
minX = polygon(0).X
maxX = polygon(0).X
minY = polygon(0).Y
maxY = polygon(0).Y
For i = 1 To UBound(polygon)
Set q = polygon(i)
minX = vbMin(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
isPointInPolygon = False
Exit Function
End If
' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
isPointInPolygon = False
i = 0
j = UBound(polygon)
Do While i < UBound(polygon) + 1
If (polygon(i).Y > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
ElseIf (polygon(i).Y < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
End If
j = i
i = i + 1
Loop
End Function
Function vbMax(n1, n2) As Double
vbMax = IIf(n1 > n2, n1, n2)
End Function
Function vbMin(n1, n2) As Double
vbMin = IIf(n1 > n2, n2, n1)
End Function
Sub TestPointInPolygon()
Dim i As Integer
Dim InPolygon As Boolean
' MARKER Object
Dim p As CPoint
Set p = New CPoint
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
Ho realizzato un'implementazione Python di nirg's c ++ codice :
Ingressi
- bounding_points: nodi che compongono il poligono.
bounding_box_positions: punti candidati da filtrare. (Nella mia implementazione creata dal rettangolo di selezione.
(Gli input sono elenchi di tuple nel formato:
[(xcord, ycord), ...])
I ritorni
- Tutti i punti all'interno del poligono.
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
Ancora una volta, l'idea è presa da qui
Ecco un punto nel test poligonale in C che non utilizza il ray-casting. E può funzionare per aree sovrapposte (autointersezioni), vedere l'argomento use_holes.
/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);
/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
const bool use_holes)
{
/* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
float angletot = 0.0;
float fp1[2], fp2[2];
unsigned int i;
const float *p1, *p2;
p1 = verts[nr - 1];
/* first vector */
fp1[0] = p1[0] - pt[0];
fp1[1] = p1[1] - pt[1];
for (i = 0; i < nr; i++) {
p2 = verts[i];
/* second vector */
fp2[0] = p2[0] - pt[0];
fp2[1] = p2[1] - pt[1];
/* dot and angle and cross */
angletot += angle_signed_v2v2(fp1, fp2);
/* circulate */
copy_v2_v2(fp1, fp2);
p1 = p2;
}
angletot = fabsf(angletot);
if (use_holes) {
const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
angletot -= nested * (float)(M_PI * 2.0);
return (angletot > 4.0f) != ((int)nested % 2);
}
else {
return (angletot > 4.0f);
}
}
/* math lib */
static float dot_v2v2(const float a[2], const float b[2])
{
return a[0] * b[0] + a[1] * b[1];
}
static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
return atan2f(perp_dot, dot_v2v2(v1, v2));
}
static void copy_v2_v2(float r[2], const float a[2])
{
r[0] = a[0];
r[1] = a[1];
}
Nota: questo è uno dei metodi meno ottimali poiché include molte chiamate a atan2f, ma potrebbe essere interessante per gli sviluppatori che leggono questo thread (nei miei test è ~ 23x più lento rispetto all'uso del metodo dell'intersezione di linea ).
Per rilevare hit su Polygon dobbiamo testare due cose:
- Se Punto è all'interno dell'area del poligono. (può essere realizzato con Ray-Casting Algorithm)
- Se Point si trova sul bordo del poligono (può essere realizzato dallo stesso algoritmo utilizzato per il rilevamento del punto sulla polilinea (linea)).
Per gestire i seguenti casi speciali in Algoritmo di ray casting :
- Il raggio si sovrappone a uno dei lati del poligono.
- Il punto si trova all'interno del poligono e il raggio passa attraverso un vertice del poligono.
- Il punto si trova all'esterno del poligono e il raggio tocca appena uno degli angoli del poligono.
Controlla Determinare se un punto si trova all'interno di un poligono complesso . L'articolo fornisce un modo semplice per risolverli in modo che non ci sia alcun trattamento speciale richiesto per i casi di cui sopra.
Puoi farlo controllando se l'area formata collegando il punto desiderato ai vertici del tuo poligono corrisponde all'area del poligono stesso.
Oppure potresti verificare se la somma degli angoli interni dal tuo punto a ciascuna coppia di due vertici poligonali consecutivi al tuo punto di controllo si somma a 360, ma ho la sensazione che la prima opzione sia più veloce perché non comporta divisioni né calcoli dell'inverso delle funzioni trigonometriche.
Non so cosa succede se il tuo poligono ha un buco al suo interno, ma mi sembra che l'idea principale possa essere adattata a questa situazione
Puoi anche pubblicare la domanda in una comunità matematica. Scommetto che hanno un milione di modi per farlo
Se stai cercando una libreria java-script c'è un'estensione javascript google maps v3 per la classe Polygon per rilevare se un punto risiede o meno al suo interno.
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
Quando si utilizza qt (Qt 4.3+) , si può usare la funzione di QPolygon contienePoint
La risposta dipende dal fatto che tu abbia i poligoni semplici o complessi. I poligoni semplici non devono avere intersezioni di segmenti di linea. Quindi possono avere i buchi ma le linee non possono incrociarsi. Le regioni complesse possono avere le intersezioni di linea, quindi possono avere le regioni sovrapposte o le regioni che si toccano solo per un singolo punto.
Per i poligoni semplici l'algoritmo migliore è l'algoritmo Ray casting (numero di incrocio). Per i poligoni complessi, questo algoritmo non rileva i punti all'interno delle regioni sovrapposte. Quindi per i poligoni complessi devi usare l'algoritmo del numero di Winding.
Ecco un eccellente articolo con l'implementazione in C di entrambi gli algoritmi. Li ho provati e funzionano bene.
Versione scala della soluzione di nirg (presuppone che il controllo preliminare del rettangolo di delimitazione sia eseguito separatamente):
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
Ecco la versione golang della risposta @nirg (ispirata al codice C # di @@ m-katz)
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
Nessuno sorpreso ne ha parlato prima, ma per i pragmatici che richiedono un database: MongoDB ha un eccellente supporto per le query Geo, incluso questo.
Quello che stai cercando è:
db.neighbors.findOne ({geometry: {$ geoIntersects: {$ geometry: { digitare: " Punto " ;, coordinate: [" longitudine " ;, " latitudine " ]}}} })
Neighborhoods è la raccolta che memorizza uno o più poligoni nel formato GeoJson standard. Se la query restituisce null, non viene intersecata, altrimenti lo è.
Molto ben documentato qui: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
Le prestazioni per oltre 6.000 punti classificati in una griglia a 330 poligoni irregolari sono state meno di un minuto senza ottimizzazione e incluso il tempo di aggiornamento dei documenti con il rispettivo poligono.
