Jaki jest najszybszy sposób sprawdzenia, czy punkt znajduje się wewnątrz wielokąta w Pythonie

87

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
Ponieważ implementacją matplotlib jest C ++, prawdopodobnie można oczekiwać, że będzie szybsza. Biorąc pod uwagę, że matplotlib jest bardzo szeroko stosowany i ponieważ jest to bardzo podstawowa funkcja, prawdopodobnie można bezpiecznie założyć, że działa poprawnie (nawet jeśli może się wydawać „niejasny”). I na koniec: dlaczego nie po prostu go przetestować?
— sebastian
Zaktualizowałem pytanie testem, tak jak przewidziałeś, matplotlib jest znacznie szybsze. Martwiłem się, ponieważ matplotlib nie jest najbardziej znaną odpowiedzią w różnych miejscach, w których szukałem, i chciałem wiedzieć, czy coś przeoczyłem (lub lepszy pakiet). Również matplotlib wyglądał na dużego gościa jak na tak proste pytanie.
— Ruben Perez-Carrasco

Odpowiedzi:

106

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:

punkt wewnątrz wielokąta z tolerancją pikseli

— armatita
źródło
1
Dzięki za to, na razie pozostanę przy matplotlib, ponieważ wydaje się znacznie szybszy niż niestandardowy ray tracing. Niemniej jednak bardzo podoba mi się odpowiedź na dyskretyzację przestrzeni, być może będę jej potrzebować w przyszłości. Sprawdzę też zgrabnie, bo wygląda to na pakiet poświęcony tego typu problemom
— Ruben Perez-Carrasco
20

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 inside

Pierwsze 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 ms

Który po kompilacji zmniejszy się do:

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

Jeś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 ms

W kodzie numba użyłem: „b1 (f8, f8, f8 [:,:])”

Aby się skompilować z nopython=True, każda zmienna musi zostać zadeklarowana przed for 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 zmiennoprzecinkowych f8oraz dwuwymiarowej tablicy wartości zmiennoprzecinkowych f8[:,:]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 D

Uwaga: 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/3336

Test:


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)
— epifanio
źródło
11

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.103765010834

Matplotlib jest nadal znacznie lepszy, ale nie 100 razy lepszy. Teraz spróbujmy znacznie prostszego wielokąta ...

lenpoly = 5
# ... same code

wynik:

Ray Tracing Elapsed time: 0.0727779865265
Matplotlib contains_points Elapsed time: 0.105288982391
— Даниил Тарарухин
źródło
6

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 inside

Opakowano 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
— user3274748
źródło
jak mogę zwrócić tylko prawdę lub fałsz dla jednego poli i jednego x, y?
— Jasar Orion
Użyłbym rozwiązania @epifanio, jeśli robisz tylko jeden poli. Rozwiązanie NumPy jest lepsze do obliczeń w większych partiach.
— Can Hicabi Tartanoglu
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.