如何在python中获取多边形内的点列表？

Posted 2023-02-14

【中文标题】如何在python中获取多边形内的点列表？

【英文标题】：How to get list of points inside a polygon in python?

【发布时间】：2014-02-15 20:19:53

【问题描述】：

我搜索了很多，但找不到任何实用的答案。我有一个多边形。例如：

    [(86, 52), (85, 52), (81, 53), (80, 52), (79, 48), (81, 49), (86, 53),
     (85, 51), (82, 54), (84, 54), (83, 49), (81, 52), (80, 50), (81, 48),
     (85, 50), (86, 54), (85, 54), (80, 48), (79, 50), (85, 49), (80, 51),
     (85, 53), (82, 49), (83, 54), (82, 53), (84, 49), (79, 49)]

我想得到这个边界多边形内所有点的列表。我听说了很多关于多边形三角测量技术或线性/洪水/交叉/...填充算法。但我真的想不出一种有效的方法来实现这一点。这个多边形很小，想象一个有 10 亿个点的多边形。我现在使用 PIL 绘制多边形用红色填充多边形并在其中循环以找到红点。这是一个非常缓慢的技术：

def render(poly, z):
    xs = [i[0] for i in poly]
    ys = [i[1] for i in poly]
    minx, maxx = min(xs), max(xs)
    miny, maxy = min(ys), max(ys)
    X = maxx - minx + 1
    Y = maxy - miny + 1
    newPoly = [(x - minx, y - miny) for (x, y) in polygons]
    i = Image.new("RGB", (X, Y))
    draw = ImageDraw.Draw(i)
    draw.polygon(newPoly, fill="red")
    # i.show()
    tiles = list()
    w, h = i.size
    print w, h
    for x in range(w):
        for y in range(h):
            data = i.getpixel((x, y))
            if data != (0, 0, 0):
                tiles.append((x + minx, y + miny))

    return tiles

我正在寻找解决这个问题的 Pythonic 方法。 谢谢大家。

【问题讨论】：

shapely 有什么可以提供给您的吗？ pypi.python.org/pypi/Shapely

我使用了 shapely，但找不到任何解决此问题的方法。谢谢，我会搜索的。

您可能正在寻找这个：rosettacode.org/wiki/Ray-casting_algorithm.

好吧，链接中的函数给出并指向和一个多边形来检查它是否在里面。这不是我需要的。我可以为所有项目手动创建一个网格和循环。但我正在寻找一种更直接的方法。谢谢顺便说一句。

@Farsheed 您找到问题的答案了吗？我现在正在寻找类似的解决方案，我有 100,000 个点的坐标和几个多边形的坐标。我需要删除这些多边形内的那些点

【参考方案1】：

我建议使用matplotlib contains_points()

from matplotlib.path import Path

tupVerts=[(86, 52), (85, 52), (81, 53), (80, 52), (79, 48), (81, 49), (86, 53),
 (85, 51), (82, 54), (84, 54), (83, 49), (81, 52), (80, 50), (81, 48),
 (85, 50), (86, 54), (85, 54), (80, 48), (79, 50), (85, 49), (80, 51),
 (85, 53), (82, 49), (83, 54), (82, 53), (84, 49), (79, 49)]


x, y = np.meshgrid(np.arange(300), np.arange(300)) # make a canvas with coordinates
x, y = x.flatten(), y.flatten()
points = np.vstack((x,y)).T 

p = Path(tupVerts) # make a polygon
grid = p.contains_points(points)
mask = grid.reshape(300,300) # now you have a mask with points inside a polygon

【讨论】：

你可以使用 np.indices() 来创建一个带有坐标的数组

这比使用 shapely 评估多边形内的每个点要快得多。使用 1024x1024 图像时，使用 shapely 计算每个像素相对于多边形的输入或输出值需要 4 到 6 分钟。使用上面的 matplotlib contains_points() 我可以在平均 5 秒内完成。

这也比在 Shapely 中使用 MultiPoint(points).intersection(polygon) 快得多。

【参考方案2】：

根据 RemcoGerlich 的回答，这里是一个经过验证的函数：

import numpy as np
import mahotas

def render(poly):
    """Return polygon as grid of points inside polygon.

    Input : poly (list of lists)
    Output : output (list of lists)
    """
    xs, ys = zip(*poly)
    minx, maxx = min(xs), max(xs)
    miny, maxy = min(ys), max(ys)

    newPoly = [(int(x - minx), int(y - miny)) for (x, y) in poly]

    X = maxx - minx + 1
    Y = maxy - miny + 1

    grid = np.zeros((X, Y), dtype=np.int8)
    mahotas.polygon.fill_polygon(newPoly, grid)

    return [(x + minx, y + miny) for (x, y) in zip(*np.nonzero(grid))]

示例：

poly = [
    [0, 0],
    [0, 10],
    [10, 10],
    [10, 0]
]

plt.figure(None, (5, 5))
x, y = zip(*render(poly))
plt.scatter(x, y)
x, y = zip(*poly)
plt.plot(x, y, c="r")
plt.show()

【参考方案3】：

我认为绘制多边形并填充它是一个好的开始，无论如何你都需要这样的东西，而且这些算法通常在 C 中进行了微调。但不要使用 RGB 图像，使用黑白图像图像，并使用numpy.where() 找到它为 1 的像素。

根据this question，mahotas 库有一个fill_polygon 函数，可用于 numpy 数组。

我从您的函数开始以下代码（我也会减去 minx 和 maxx）但请注意，我根本无法测试它，我不在我的开发机器上。

import numpy as np
import mahotas

def render(poly): # removed parameter 'z'
    xs = [i[0] for i in poly]
    ys = [i[1] for i in poly]
    minx, maxx = min(xs), max(xs)
    miny, maxy = min(ys), max(ys)
    X = maxx - minx + 1
    Y = maxy - miny + 1
    newPoly = [(x - minx, y - miny) for (x, y) in poly]           

    grid = np.zeros((X, Y), dtype=np.int8)
    mahotas.polygon.fill_polygon(newPoly, grid)

    return [(x + minx, y + miny) for (x, y) in np.where(grid)]

【讨论】：

我正在测试它。第一点是它应该是“从 mahotas 导入多边形”。让我再测试一下。

好的，您的代码比 PIL 填充方法慢 38%。使用 1.3M 多边形进行测试。 PIL 需要 5.6 秒，np+mahotas 需要 6.7 秒。

哦，酷。不过，这相差不到 38% :-)。你能用 PIL 填充多边形并使用 np.where 来查找点吗？我认为 PIL 的图像也是 numpy 数组？

我想知道有多少时间花在减去 minx 和 miny 并稍后将其添加回来。但是，如果您的多边形具有任意索引，则无法真正想到避免它的方法。

我现在的问题是如何将 np.where 与 img 一起使用。我们可以使用 np.aaray(img) 将图像转换为数组。但是如何在哪里使用以及如何将其转换回正常的瓷砖列表？

【参考方案4】：

您可以使用像二进制图像这样的 numpy 矩阵，例如可以与 Opencv 或其他图像处理库一起使用， 解决方案 1 所以一个大小为 L x H 的矩阵，其中

L=max(x) - min (x)
H=max(y) - min (y)

作为条目，我们有您在上面给出的元组 (x,y) 列表，其中名称是 poly 例如：

import numpy as np
matrix =np.zeros((L,H),dtype=np.int32) # you can use np.uint8 if unsigned x ,y

所以我们现在有一个大小为 L x H 的矩阵，用 0 填充，我们现在将 1 放在多边形点的位置

我认为你可以简单地做到这一点

matrix[poly]=1  # which will put 1 at each (x,y) of the list **poly**

我们将此解释为二值（黑/白）图像，其上绘制有轮廓 假设我们要检测这个新轮廓

import cv2 # opencv import
ContoursListe,hierarchy = cv2.findContours(self.thresh,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_NONE)
poly2=ContoursListe[0] # we take the first only contour

注意：poly2 包含多边形的点列表以及构成它的所有点，我的意思是多边形的每个顶点的所有点，您需要这些点会发现有用 ！！您可以使用 cv2.CHAIN_APPROX_SIMPLE 参数来获取 poly2 仅包含多边形线的端点，这是我们的输入：） 重要：poly2的类型是numpy数组，它的形状是(n,1,2)而不是(n,2)

现在我们在这个图像（矩阵）上绘制这个轮廓并且也将填充它:)

cv2.drawContours(matrix,[poly2],-1,(1),thickness=-1) thickness=-1

现在我们有一个矩阵，其中每个点都有 1 形成并填充多边形，“thickness=-1”已强制填充此轮廓，您可以设置厚度 = 1 得到只有边界 如果要翻译，可以通过添加参数offset(xtran,ytrans)

来实现

要获得所有这些点的索引，只需调用

list_of_points_indices=numpy.nonzero(matrix)

---

解决方案 2

更聪明的是直接将你的点列表（poly）转换为轮廓格式（poly2）并在矩阵上绘制

poly2=poly.reshape(-1,1,2).astype(np.int32)

并在Matrix矩阵上绘制

matrix =np.zeros((L,H),dtype=np.int32)

cv2.drawContours(matrix,[poly2],-1,(1),thickness=-1)

并通过以下方式获取这些点的列表：

list_of_points_indices=numpy.nonzero(matrix)

使用 thickness 填充或不填充多边形，请参阅解决方案 1 了解更多详情。

【参考方案5】：

试试这个代码。 poly_coords 是多边形的坐标，'coord' 是要检查它是否在多边形内的点的坐标。

def testP(coord, poly_coords):
    """
    The coordinates should be in the form of list of x and y
    """
    test1 = n.array(poly_coords)
    test2 = n.vstack((poly_coords[1:], poly_coords[:1]))
    test  = test2-test1
    m = test[:,1]/test[:,0]
    c = test1[:,1]-m*test1[:,0]
    xval = (coord[1]-c)/m
    print 'xVal:\t'; print xval
    print (test1[:,0]-xval)*(test2[:,0]-xval)
    check = n.where((xval>=coord[0])&((test1[:,0]-xval)*(test2[:,0]-xval)<0))[0]
    print check
    print len(check)
    if len(check)%2==0:
        return False
    else:
        return True

如果您想让它更快，请取出与多边形、坡度和偏移量相关的算法部分，然后使用“地图”函数运行其余代码。像这样的：

test1 = n.array( your polygon)
test2 = n.vstack((test1[1:], test1[:1]))
test  = test2-test1
m = test[:,1]/test[:,0]
c = test1[:,1]-m*test1[:,0]

def testP(coord):
    """
    The coordinates should be in the form of list of x and y
    """
    global test, test1, test2, m,c
    xval = (coord[1]-c)/m
    check = n.where((xval>=coord[0])&((test1[:,0]-xval)*(test2[:,0]-xval)<0))[0]
    if len(check)%2==0:
        return False
    else:
        return True
coords = n.array(( your coords in x,y ))
map (testP, coords)

您可以根据需要删除“打印”命令。此代码是为 python 2.7 编写的

【讨论】：

您能解释一下your polygon 和your coords in x,y 吗？似乎your polygon 的意思类似于poly = [[0, 0], [0, 10], [10, 10], [10, 0]]。但我不知道your coords in x,y到底是什么以及怎么样！

对于多边形，你是绝对正确的。 “x,y 坐标”是指要检查是否在多边形内的点坐标。坐标应该是一个 numpy 数组，例如 (5,6)，其中点在 x 轴上的位置为 5，y 轴为 6