Qual è il modo più veloce per verificare se un punto si trova all'interno di un poligono in python

Question

Ho trovato due metodi principali per verificare se un punto appartiene all'interno di un poligono. Uno sta usando il metodo di tracciamento dei raggi usato here, che è la risposta più raccomandata, l'altra sta usando matplotlib path.contains_points (che a me sembra un po 'oscuro). Dovrò controllare molti punti continuamente. Qualcuno sa se qualcuno di questi due è più raccomandabile dell'altro o se ci sono ancora terze opzioni migliori?

UPDATE:

Ho controllato i due metodi e matplotlib sembra molto più veloce.

from time import time 
import numpy as np 
import matplotlib.path as mpltPath 

# regular polygon for testing 
lenpoly = 100 
polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)[:-1]] 

# random points set of points to test 
N = 10000 
points = zip(np.random.random(N),np.random.random(N)) 


# Ray tracing 
def ray_tracing_method(x,y,poly): 

    n = len(poly) 
    inside = False 

    p1x,p1y = poly[0] 
    for i in range(n+1): 
     p2x,p2y = poly[i % n] 
     if y > min(p1y,p2y): 
      if y <= max(p1y,p2y): 
       if x <= max(p1x,p2x): 
        if p1y != p2y: 
         xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x 
        if p1x == p2x or x <= xints: 
         inside = not inside 
     p1x,p1y = p2x,p2y 

    return inside 

start_time = time() 
inside1 = [ray_tracing_method(point[0], point[1], polygon) for point in points] 
print "Ray Tracing Elapsed time: " + str(time()-start_time) 

# Matplotlib mplPath 
start_time = time() 
path = mpltPath.Path(polygon) 
inside2 = path.contains_points(points) 
print "Matplotlib contains_points Elapsed time: " + str(time()-start_time) 

che dà,

Ray Tracing Elapsed time: 0.441395998001 
Matplotlib contains_points Elapsed time: 0.00994491577148 

stessa differenza relativa è stata ottenuta una utilizzando un triangolo invece del 100 lati poligono. Vorrei anche verificare formosa dal momento che sembra un pacchetto solo dedicato a questo tipo di problemi

Answer 1

Si può considerare shapely:

from shapely.geometry import Point 
from shapely.geometry.polygon import Polygon 

point = Point(0.5, 0.5) 
polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]) 
print(polygon.contains(point)) 

Dai metodi che hai citato ho usato solo la seconda, path.contains_points, e funziona bene. In ogni caso, a seconda della precisione necessaria per il test, suggerirei di creare una griglia bool numerosa con tutti i nodi all'interno del poligono per essere True (Falso se non). Se avete intenzione di fare un test per un sacco di punti questo potrebbe essere più veloce (anche se preavviso questo si basa si stanno facendo un test all'interno di un "pixel" tolleranza):

from matplotlib import path 
import matplotlib.pyplot as plt 
import numpy as np 

first = -3 
size = (3-first)/100 
xv,yv = np.meshgrid(np.linspace(-3,3,100),np.linspace(-3,3,100)) 
p = path.Path([(0,0), (0, 1), (1, 1), (1, 0)]) # square with legs length 1 and bottom left corner at the origin 
flags = p.contains_points(np.hstack((xv.flatten()[:,np.newaxis],yv.flatten()[:,np.newaxis]))) 
grid = np.zeros((101,101),dtype='bool') 
grid[((xv.flatten()-first)/size).astype('int'),((yv.flatten()-first)/size).astype('int')] = flags 

xi,yi = np.random.randint(-300,300,100)/100,np.random.randint(-300,300,100)/100 
vflag = grid[((xi-first)/size).astype('int'),((yi-first)/size).astype('int')] 
plt.imshow(grid.T,origin='lower',interpolation='nearest',cmap='binary') 
plt.scatter(((xi-first)/size).astype('int'),((yi-first)/size).astype('int'),c=vflag,cmap='Greens',s=90) 
plt.show() 

, i risultati è questo:

Answer 2

Il test è buono, ma misura solo qualche situazione specifica: abbiamo una poligono con molti vertici, e lunga serie di punti per controllare loro all'interno del poligono.

Inoltre, suppongo che si sta misurando non matplotlib-dentro-poligono-metodo vs ray-metodo, ma matplotlib-in qualche modo ottimizzato-iterazione vs semplice-list-iterazione

Facciamo N confronti indipendenti (N coppie di punti e poligoni)?

# ... your code... 
lenpoly = 100 
polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)[:-1]] 

M = 10000 
start_time = time() 
# Ray tracing 
for i in range(M): 
    x,y = np.random.random(), np.random.random() 
    inside1 = ray_tracing_method(x,y, polygon) 
print "Ray Tracing Elapsed time: " + str(time()-start_time) 

# Matplotlib mplPath 
start_time = time() 
for i in range(M): 
    x,y = np.random.random(), np.random.random() 
    inside2 = path.contains_points([[x,y]]) 
print "Matplotlib contains_points Elapsed time: " + str(time()-start_time) 

Risultato:

Ray Tracing Elapsed time: 0.548588991165 
Matplotlib contains_points Elapsed time: 0.103765010834 

Matplotlib è ancora molto meglio, ma non al 100 volte meglio. Ora proviamo molto più semplice poligono ...

lenpoly = 5 
# ... same code 

risultato:

Ray Tracing Elapsed time: 0.0727779865265 
Matplotlib contains_points Elapsed time: 0.105288982391 

Answer 3

Se la velocità è quello che ti serve e le dipendenze in più non sono un problema, magari trovare numba molto utile (ora è abbastanza facile da installare, su qualsiasi piattaforma).Il classico approccio ray_tracing che hai proposto può essere facilmente portato su numba usando il decoratore numba @jit e convertendo il poligono in una matrice numpy. Il codice dovrebbe essere simile:

@jit(nopython=True) 
def ray_tracing(x,y,poly): 
    n = len(poly) 
    inside = False 
    p2x = 0.0 
    p2y = 0.0 
    xints = 0.0 
    p1x,p1y = poly[0] 
    for i in range(n+1): 
     p2x,p2y = poly[i % n] 
     if y > min(p1y,p2y): 
      if y <= max(p1y,p2y): 
       if x <= max(p1x,p2x): 
        if p1y != p2y: 
        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x 
        if p1x == p2x or x <= xints: 
         inside = not inside 
     p1x,p1y = p2x,p2y 

    return inside 

La prima esecuzione avrà un po 'di più rispetto a qualsiasi chiamata successiva:

%%time 
polygon=np.array(polygon) 
inside1 = [numba_ray_tracing_method(point[0], point[1], polygon) for 
point in points] 

CPU times: user 129 ms, sys: 4.08 ms, total: 133 ms 
Wall time: 132 ms 

Il che, dopo la compilazione diminuirà a:

CPU times: user 18.7 ms, sys: 320 µs, total: 19.1 ms 
Wall time: 18.4 ms 

Se bisogno di velocità al primo richiamo della funzione è possibile quindi precompilare il codice in un modulo utilizzando pycc. Conservare la funzione in uno src.py come:

from numba import jit 
from numba.pycc import CC 
cc = CC('nbspatial') 


@cc.export('ray_tracing', 'b1(f8, f8, f8[:,:])') 
@jit(nopython=True) 
def ray_tracing(x,y,poly): 
    n = len(poly) 
    inside = False 
    p2x = 0.0 
    p2y = 0.0 
    xints = 0.0 
    p1x,p1y = poly[0] 
    for i in range(n+1): 
     p2x,p2y = poly[i % n] 
     if y > min(p1y,p2y): 
      if y <= max(p1y,p2y): 
       if x <= max(p1x,p2x): 
        if p1y != p2y: 
         xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x 
        if p1x == p2x or x <= xints: 
         inside = not inside 
     p1x,p1y = p2x,p2y 

    return inside 


if __name__ == "__main__": 
    cc.compile() 

costruire con python src.py ed eseguire:

import nbspatial 

import numpy as np 
lenpoly = 100 
polygon = [[np.sin(x)+0.5,np.cos(x)+0.5] for x in 
np.linspace(0,2*np.pi,lenpoly)[:-1]] 

# random points set of points to test 
N = 10000 
# making a list instead of a generator to help debug 
points = zip(np.random.random(N),np.random.random(N)) 

polygon = np.array(polygon) 

%%time 
result = [nbspatial.ray_tracing(point[0], point[1], polygon) for point in points] 

CPU times: user 20.7 ms, sys: 64 µs, total: 20.8 ms 
Wall time: 19.9 ms 

Nel codice Numba ho usato: 'b1 (F8, f8, f8 [:, :]) '

Per compilare con nopython=True, ogni var deve essere dichiarato prima dello for loop.

Nel codice prebuild src linea:

@cc.export('ray_tracing', 'b1(f8, f8, f8[:,:])') 

Si usa per dichiarare il nome della funzione e suoi tipi I/O var, un'uscita booleana b1 e due galleggianti f8 e una matrice bidimensionale di carri f8[:,:] come input.