Я нашел два основных метода определения принадлежности точки к многоугольнику. Один использует метод трассировки лучей, используемый здесь , что является наиболее рекомендуемым ответом, другой использует matplotlib path.contains_points(что мне кажется немного непонятным). Мне придется постоянно проверять множество точек. Кто-нибудь знает, является ли один из этих двух более рекомендуемым, чем другой, или есть еще лучшие третьи варианты?
ОБНОВИТЬ:
Я проверил два метода, и matplotlib выглядит намного быстрее.
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))
который дает,
Ray Tracing Elapsed time: 0.441395998001
Matplotlib contains_points Elapsed time: 0.00994491577148
Такая же относительная разница была получена при использовании треугольника вместо многоугольника со 100 сторонами. Я также проверю форму, так как это выглядит как пакет, предназначенный только для такого рода проблем
python
matplotlib
Рубен Перес-Карраско
источник
источник
Ответы:
Вы можете считать стройными :
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, и он отлично работает. В любом случае, в зависимости от точности, которая вам нужна для вашего теста, я бы предложил создать сетку numpy bool со всеми узлами внутри многоугольника как True (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(), результаты следующие:
источник
Если скорость - это то, что вам нужно, а дополнительные зависимости не проблема, возможно, вы найдете его
numbaвесьма полезным (теперь его довольно легко установить на любой платформе). Предложенныйray_tracingвами классический подход можно легко перенестиnumbaс помощьюnumba @jitдекоратора и преобразования многоугольника в массив numpy. Код должен выглядеть так:@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.В коде prebuild src строка:
@cc.export('ray_tracing' , 'b1(f8, f8, f8[:,:])')Используется для объявления имени функции и ее типов ввода-вывода var, логического вывода
b1и двух чисел с плавающей запятойf8и двумерного массива чисел с плавающей запятой вf8[:,:]качестве ввода.Редактировать 4 января 2021 г.
Для моего варианта использования мне нужно проверить, находятся ли несколько точек внутри одного многоугольника - в таком контексте полезно воспользоваться преимуществами параллельных возможностей numba для перебора серии точек. Приведенный выше пример можно изменить на:
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Примечание: предварительная компиляция приведенного выше кода не активирует параллельные возможности numba (параллельный целевой ЦП не поддерживается
pycc/AOTкомпиляцией) см .: https://github.com/numba/numba/issues/3336Контрольная работа:
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))Для
N=1000072-ядерной машины возвращает:%%timeit parallelpointinpolygon(points, polygon) # 480 µs ± 8.19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)источник
Ваш тест хорош, но он измеряет только некоторую конкретную ситуацию: у нас есть один многоугольник с множеством вершин и длинный массив точек для проверки их внутри многоугольника.
Более того, я полагаю, что вы измеряете не matplotlib-inside-polygon-method vs ray-method, а matplotlib-some-optimized-iteration vs simple-list-iteration
Сделаем 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.103765010834Matplotlib по-прежнему намного лучше, но не в 100 раз лучше. Теперь попробуем многоугольник намного проще ...
lenpoly = 5 # ... same codeрезультат:
Ray Tracing Elapsed time: 0.0727779865265 Matplotlib contains_points Elapsed time: 0.105288982391источник
Я просто оставлю это здесь, просто переписав приведенный выше код с помощью 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источник