如何确定哪些点位于多边形内部哪些不是(大量点)？

Question

我有一大堆数据点(100,000+)存储在一个二维numpy数组中(第一列:x坐标,第二列:y坐标).我还有几个1维数组存储每个数据点的附加信息.我现在想要从这些一维数组的子集创建图,其中只包含给定多边形中的点.

我提出了以下解决方案,既不优雅也不快:

#XY is the 2D array.
#A is one of the 1D arrays.
#poly is a matplotlib.patches.Polygon

mask = np.array([bool(poly.get_path().contains_point(i)) for i in XY])

matplotlib.pylab.hist(A[mask], 100)
matplotlib.pylab.show()

你能帮我改进一下这段代码吗？我尝试使用np.vectorize而不是列表理解,但无法让它工作.

Answer 1

使用matplotlib.nxutils.points_inside_poly,它实现了一个非常有效的测试.

matplotlib常见问题解答中这个有40年历史的算法的例子和进一步说明.

更新:请注意,points_inside_poly自matplotlib版本1.2.0起不推荐使用.请改用matplotlib.path.Path.contains_points.

Answer 2

我担心我不熟悉您正在使用的库,但我认为我对您可以使用的算法有一个合理的想法,我将直接介绍如何用vanilla python实现它然后我就是确保你可以改进它并使用这些库实现它.此外,我并不是说这是实现这一目标的最佳方式,但我希望能够快速得到我的回复,所以这里也是如此.

现在,这个想法来自于在算法中使用两个向量的交叉积来找到一组点的凸集,例如格雷厄姆的扫描.假设我们有两个点p1和p2,它们分别定义了从原点(0,0)到(x1,y1)和(x2,y2)的点向量p1和p2.p1 × p2的叉积给出第三矢量p3,其垂直于p1和p2,并且具有由矢量限定的平行四边形的面积给出的幅度.

一个非常有用的结果是矩阵的行列式

/ x1, x2 \
\ y1, y2 /

...是x1*y2-x2*y1给出向量p3的大小,符号表示p3是"从平面"出来还是"进入"它.这里的关键点是,如果这个幅度是正的那么p2是p1的 "左边",如果它是负的,那么p2是p1的 "右边" .

希望这个ascii艺术示例将有所帮助:

    . p2(4, 5)
   /
  /
 /
/_ _ _ _ _. p1(5, 0)

x1*y2 - x2*y1 = 5*4 - 0*5 = 20,因此p2是p1的 "左侧"

最后说明为什么这对我们有用!如果我们有一个多边形的顶点列表和图形中的一组其他点,那么对于多边形的每个边,我们可以得到该边的矢量.我们还可以获得将起始顶点连接到图中所有其他点的向量,并通过测试这些向量是否位于边的左侧或右侧,我们可以消除每个边的一些点.所有未在过程结束时删除的是多边形内的那些点.无论如何,对一些代码,以使一些更多的意义!

如果您以逆时针方向绘制多边形顶点,请按照您访问它们的顺序获取多边形顶点列表,例如某些五边形可能是:

poly = [(1, 1), (4, 2), (5, 5), (3, 8), (0, 4)]

获取一个包含图中所有其他点的集合,我们将逐渐从该集合中删除无效点,直到在过程结束时留下的那些点正好是多边形内的那些点.

points = set(['(3, 0), (10, -2), (3,3), ...])

代码本身的主要部分实际上非常紧凑,我花了多长时间来写它是如何工作的.to_right采用两个代表向量的元组,True如果v2位于右侧则返回v1.然后循环遍历多边形的所有边缘,如果它们位于任何边缘的右侧,则从工作集中移除点.

def to_right(v1, v2):
    return (v1[0]*v2[1] - v1[1]*v2[0]) < 0

for i in range(len(poly)):
    v1 = poly[i-1]
    v2 = poly[i]
    for p in points:
        if(to_right(v2-v1, p-v1)):
            points.remove(p)

编辑:为了澄清,如果它们在右侧而不是左侧,则它们被移除的事实与指定多边形顶点的顺序相关联.如果它们按顺时针顺序排列,则可能需要消除左侧点.目前我对这个问题没有特别好的解决方案.

无论如何,希望我对这些东西是正确的,即使不是OP,它对某些人也有帮助.该算法的渐近复杂度为O(mn),其中n是图中点的数量,m是多边形顶点的数量,因为在最坏的情况下,所有点都位于多边形内部,我们必须检查每个点对于每一个边缘,没有被删除.