Znalazłem dwie główne metody sprawdzania, czy punkt należy do wielokąta. Jeden z nich wykorzystuje zastosowaną tutaj metodę śledzenia promieni , który jest najbardziej zalecaną odpowiedzią, drugi to matplotlib path.contains_points(co wydaje mi się nieco niejasne). Będę musiał ciągle sprawdzać wiele punktów. Czy ktoś wie, czy którakolwiek z tych dwóch jest bardziej godna polecenia niż druga, czy też istnieją jeszcze lepsze opcje trzecie?
AKTUALIZACJA:
Sprawdziłem dwie metody i matplotlib wygląda znacznie szybciej.
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 = np.random.rand(N,2)
# 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))
co daje,
Ray Tracing Elapsed time: 0.441395998001
Matplotlib contains_points Elapsed time: 0.00994491577148
Tę samą względną różnicę uzyskano przy użyciu trójkąta zamiast wielokąta o 100 bokach. Sprawdzę też zgrabnie, bo wygląda na pakiet właśnie poświęcony tego typu problemom
python
matplotlib
Ruben Perez-Carrasco
źródło
źródło
Odpowiedzi:
Możesz rozważyć zgrabną :
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))Z metod, o których wspomniałeś, użyłem tylko drugiej
path.contains_pointsi działa dobrze. W każdym razie, w zależności od precyzji potrzebnej do testu, sugerowałbym utworzenie numpy bool grid ze wszystkimi węzłami wewnątrz wielokąta, aby były True (False, jeśli nie). Jeśli masz zamiar wykonać test dla wielu punktów, może to być szybsze ( chociaż zwróć uwagę, że wymaga to wykonywania testu z tolerancją „piksela” ):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()wyniki są następujące:
źródło
Jeśli potrzebujesz szybkości, a dodatkowe zależności nie stanowią problemu, być może uznasz to za
numbacałkiem przydatne (teraz jest to dość łatwe do zainstalowania na dowolnej platformie).ray_tracingZaproponowane przez ciebie klasyczne podejście można łatwo przenieśćnumbaza pomocąnumba @jitdekoratora i rzutowania wielokąta na tablicę numpy. Kod powinien wyglądać następująco:@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 insidePierwsze wykonanie zajmie trochę więcej czasu niż jakiekolwiek następne wywołanie:
%%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 msKtóry po kompilacji zmniejszy się do:
CPU times: user 18.7 ms, sys: 320 µs, total: 19.1 ms Wall time: 18.4 msJeśli potrzebujesz szybkości przy pierwszym wywołaniu funkcji, możesz wstępnie skompilować kod w module przy użyciu
pycc. Przechowuj funkcję w src.py, na przykład: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()Zbuduj go
python src.pyi uruchom: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 msW kodzie numba użyłem: „b1 (f8, f8, f8 [:,:])”
Aby się skompilować z
nopython=True, każda zmienna musi zostać zadeklarowana przedfor loop.W prekompilowanym kodzie źródłowym wiersz:
@cc.export('ray_tracing' , 'b1(f8, f8, f8[:,:])')Służy do zadeklarowania nazwy funkcji i jej typów zmiennoprzecinkowych we / wy, wyjścia logicznego
b1i dwóch wartości zmiennoprzecinkowychf8oraz dwuwymiarowej tablicy wartości zmiennoprzecinkowychf8[:,:]jako danych wejściowych.Edytuj 4 stycznia 2021 r
W moim przypadku muszę sprawdzić, czy wiele punktów znajduje się w jednym wielokącie - w takim kontekście przydatne jest wykorzystanie możliwości równoległych numba do zapętlenia serii punktów. Powyższy przykład można zmienić na:
from numba import jit, njit import numba import numpy as np @jit(nopython=True) def pointinpolygon(x,y,poly): n = len(poly) inside = False p2x = 0.0 p2y = 0.0 xints = 0.0 p1x,p1y = poly[0] for i in numba.prange(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 @njit(parallel=True) def parallelpointinpolygon(points, polygon): D = np.empty(len(points), dtype=numba.boolean) for i in numba.prange(1, len(D)): D[i] = pointinpolygon(points[i,0], points[i,1], polygon) return DUwaga: prekompilacja powyższego kodu nie włączy równoległych możliwości numba (równoległy procesor CPU nie jest obsługiwany przez
pycc/AOTkompilację) patrz: https://github.com/numba/numba/issues/3336Test:
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]] polygon = np.array(polygon) N = 10000 points = np.random.uniform(-1.5, 1.5, size=(N, 2))W przypadku maszyny
N=10000z 72 rdzeniami zwraca:%%timeit parallelpointinpolygon(points, polygon) # 480 µs ± 8.19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)źródło
Twój test jest dobry, ale mierzy tylko określoną sytuację: mamy jeden wielokąt z wieloma wierzchołkami i długą tablicą punktów, aby sprawdzić je w wielokącie.
Co więcej, przypuszczam, że mierzysz nie matplotlib-inside-polygon-method vs ray-method, ale matplotlib-w jakiś sposób zoptymalizowany-iteracja vs prosta-iteracja-listy
Zróbmy N niezależnych porównań (N par punktów i wielokątów)?
# ... 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)Wynik:
Ray Tracing Elapsed time: 0.548588991165 Matplotlib contains_points Elapsed time: 0.103765010834Matplotlib jest nadal znacznie lepszy, ale nie 100 razy lepszy. Teraz spróbujmy znacznie prostszego wielokąta ...
lenpoly = 5 # ... same codewynik:
Ray Tracing Elapsed time: 0.0727779865265 Matplotlib contains_points Elapsed time: 0.105288982391źródło
Zostawię to tutaj, po prostu przepisałem powyższy kod za pomocą numpy, może ktoś uzna to za przydatne:
def ray_tracing_numpy(x,y,poly): n = len(poly) inside = np.zeros(len(x),np.bool_) p2x = 0.0 p2y = 0.0 xints = 0.0 p1x,p1y = poly[0] for i in range(n+1): p2x,p2y = poly[i % n] idx = np.nonzero((y > min(p1y,p2y)) & (y <= max(p1y,p2y)) & (x <= max(p1x,p2x)))[0] if p1y != p2y: xints = (y[idx]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x if p1x == p2x: inside[idx] = ~inside[idx] else: idxx = idx[x[idx] <= xints] inside[idxx] = ~inside[idxx] p1x,p1y = p2x,p2y return insideOpakowano ray_tracing w
def ray_tracing_mult(x,y,poly): return [ray_tracing(xi, yi, poly[:-1,:]) for xi,yi in zip(x,y)]Testowany na 100000 punktów, wyniki:
ray_tracing_mult 0:00:00.850656 ray_tracing_numpy 0:00:00.003769źródło