在python中检查点是否在多边形内的最快方法是什么

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/36399381/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-19 17:50:05  来源:igfitidea点击:

What's the fastest way of checking if a point is inside a polygon in python

pythonmatplotlib

提问by Ruben Perez-Carrasco

I found two main methods to look if a point belongs inside a polygon. One is using the ray tracing method used here, which is the most recommended answer, the other is using matplotlib path.contains_points(which seems a bit obscure to me). I will have to check lots of points continuously. Does anybody know if any of these two is more recommendable than the other or if there are even better third options?

我找到了两种主要方法来查看一个点是否属于多边形内。一种是使用这里使用的光线追踪方法,这是最推荐的答案,另一种是使用 matplotlib path.contains_points(对我来说似乎有点晦涩)。我将不得不连续检查很多点。有谁知道这两个中的任何一个是否比另一个更值得推荐,或者是否有更好的第三种选择?

UPDATE:

更新:

I checked the two methods and matplotlib looks much faster.

我检查了这两种方法,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 = 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)

which gives,

这使,

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

Same relative difference was obtained one using a triangle instead of the 100 sides polygon. I will also check shapely since it looks a package just devoted to these kind of problems

使用三角形而不是 100 边多边形获得了相同的相对差异。我也会好好检查一下,因为它看起来是一个专门解决这类问题的包

回答by armatita

You can consider 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))

From the methods you've mentioned I've only used the second, path.contains_points, and it works fine. In any case depending on the precision you need for your test I would suggest creating a numpy bool grid with all nodes inside the polygon to be True (False if not). If you are going to make a test for a lot of points this might be faster (although notice this relies you are making a test within a "pixel" tolerance):

从你提到的方法中,我只使用了第二种,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()

, the results is this:

,结果是这样的:

point inside polygon within pixel tolerance

像素容差内的多边形内点

回答by epifanio

If speed is what you need and extra dependencies are not a problem, you maybe find numbaquite useful (now it is pretty easy to install, on any platform). The classic ray_tracingapproach you proposed can be easily ported to numbaby using numba @jitdecorator and casting the polygon to a numpy array. The code should look like:

如果速度是您所需要的并且额外的依赖不是问题,那么您可能会发现它numba非常有用(现在它很容易在任何平台上安装)。通过使用装饰器并将多边形转换为 numpy 数组,ray_tracing可以轻松移植到您提出的经典方法。代码应如下所示: numbanumba @jit

@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

The first execution will take a little longer than any subsequent call:

第一次执行将比任何后续调用花费更长的时间:

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

Which, after compilation will decrease to:

其中,编译后将减少为:

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

If you need speed at the first call of the function you can then pre-compile the code in a module using pycc. Store the function in a src.py like:

如果您在第一次调用函数时需要速度,则可以使用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()

Build it with python src.pyand run:

使用以下命令构建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

In the numba code I used: 'b1(f8, f8, f8[:,:])'

在我使用的 numba 代码中:'b1(f8, f8, f8[:,:])'

In order to compile with nopython=True, each var needs to be declared before the for loop.

为了用 编译nopython=True,每个 var 需要在for loop.

In the prebuild src code the line:

在预构建 src 代码中,该行:

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

Is used to declare the function name and its I/O var types, a boolean output b1and two floats f8and a two-dimensional array of floats f8[:,:]as input.

用于声明函数名称及其 I/O 变量类型、一个布尔输出b1和两个浮点数f8以及一个二维浮点数数组f8[:,:]作为输入。

回答by Даниил Тарарухин

Your test is good, but it measures only some specific situation: we have one polygon with many vertices, and long array of points to check them within polygon.

你的测试很好,但它只测量了一些特定的情况:我们有一个有很多顶点的多边形,还有一长串点来检查它们在多边形内。

Moreover, I suppose that you're measuring not matplotlib-inside-polygon-method vs ray-method, but matplotlib-somehow-optimized-iteration vs simple-list-iteration

此外,我想你测量的不是 matplotlib-inside-polygon-method 与 ray-method,而是 matplotlib-somehow-optimized-iteration vs simple-list-iteration

Let's make N independent comparisons (N pairs of point and polygon)?

让我们进行 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)

Result:

结果:

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

Matplotlib is still much better, but not 100 times better. Now let's try much simpler polygon...

Matplotlib 仍然要好得多,但不是好 100 倍。现在让我们尝试更简单的多边形...

lenpoly = 5
# ... same code

result:

结果:

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

回答by user3274748

I will just leave it here, just rewrote the code above using numpy, maybe somebody finds it useful:

我就留在这里,用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    

Wrapped ray_tracing into

包裹 ray_tracing 到

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

Tested on 100000 points, results:

在 100000 点上测试,结果:

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