Какой самый быстрый способ проверить, находится ли точка внутри многоугольника в python

Ответ 1

Вы можете рассмотреть стройность:

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))

Из методов, о которых вы упоминали, я использовал только вторую, path.contains_points, и она отлично работает. В любом случае, в зависимости от точности, необходимой для вашего теста, я бы предложил создать сетку с множеством точек с нулевым значением, чтобы все узлы внутри полигона были истинными (False, если нет). Если вы собираетесь провести тест на множество точек, это может быть быстрее (хотя обратите внимание, что это полагает, что вы проводите тест в пределах "пиксельного" допуска):

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()

результаты таковы:

Ответ 2

Если вам нужна скорость и дополнительные зависимости не являются проблемой, вы можете найти numba весьма полезным (теперь его довольно просто установить на любой платформе). Предложенный вами классический подход ray_tracing можно легко перенести на numba, используя декоратор numba @jit и приведя многоangularьник к массиву пустышек. Код должен выглядеть примерно так:

@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

Первое выполнение займет немного больше времени, чем любой последующий вызов:

%%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

Который после компиляции уменьшится до:

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

Если вам нужна скорость при первом вызове функции, вы можете предварительно скомпилировать код в модуле, используя pycc. Сохраните функцию в src.py как:

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()

Создайте его с помощью python src.py и запустите:

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

В коде numba я использовал:   'b1 (f8, f8, f8 [:,:])'

Для компиляции с nopython=True каждый var должен быть объявлен перед for loop.

В предварительно собранном коде src строка:

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

Используется для объявления имени функции и ее типов переменных ввода/вывода, логического выхода b1 и двух чисел с плавающей точкой f8 и двумерного массива чисел с плавающей точкой f8[:,:] в качестве входных данных.

Ответ 3

Ваш тест хорош, но он измеряет только определенную ситуацию: у нас есть один многоугольник со множеством вершин и длинный массив точек для проверки их внутри полигона.

Более того, я полагаю, что вы измеряете метод matplotlib-внутри-polygon-метода vs ray, но matplotlib-somehow-optimized-iteration vs simple-list-итерация

Пусть N независимых сравнений (N пар точки и многоугольника)?

# ... 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)

Результат:

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

Matplotlib все еще намного лучше, но не в 100 раз лучше. Теперь попробуем гораздо более простой полигон...

lenpoly = 5
# ... same code

результат:

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

Ответ 4

Я просто оставлю это здесь, просто переписал приведенный выше код, используя numpy, может быть, кто-то найдет это полезным:

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    

Обернутый ray_tracing в

def ray_tracing_mult(x,y,poly):
    return [ray_tracing(xi, yi, poly[:-1,:]) for xi,yi in zip(x,y)]

Проверено на 100000 баллов, результаты:

ray_tracing_mult 0:00:00.850656
ray_tracing_numpy 0:00:00.003769