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

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/21339448/
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-18 22:37:24  来源:igfitidea点击:

How to get list of points inside a polygon in python?

pythonalgorithmgispolygonfill

提问by Farsheed

I searched a lot and cant find any practical answer to my question. I have a polygon. For example:

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

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

I want to get a list of all the points inside this border polygon. I heard alot about polygon triangulation techniques or linear/flood/intersection/... filling algorithms. but i cant really come up with an efficient way of implementing this. This poly is small, imagine a polygon with 1 billion points. I am now using PIL draw polygon to fill the poly with red color and loop inside it to find red points. This is a horribly slow technique:

我想得到这个边界多边形内所有点的列表。我听说了很多关于多边形三角测量技术或线性/洪水/交叉/...填充算法。但我真的无法想出一种有效的方法来实现这一点。这个多边形很小,想象一个有 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

I am searching for a Pythonic way of solving this problem. Thank you all.

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

回答by RemcoGerlich

I think drawing the polygon and filling it is a good start, you're going to need something like that anyway and those algorithms are usually fine tuned in C. But don't use a RGB image, use a black/white image, and use numpy.where()to find the pixels where it's 1.

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

According to this question, the mahotaslibrary has a fill_polygonfunction that works with numpy arrays.

根据这个问题,该mahotas库有一个fill_polygon适用于 numpy 数组的函数。

I'm starting the following code from your function (I would subtract the minxand maxxtoo) but note that I can't test it at all, I'm not on my dev machine.

我正在从您的函数中开始以下代码(我也会减去minxmaxx)但请注意,我根本无法测试它,我不在我的开发机器上。

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

回答by Cherif KAOUA

You can use a numpy matrix like a binary image, which can be used with Opencv for example or other image processing libs, Solution 1So a matrix which size is L x H where

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

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

As entry we have your list of tuple(x,y) you gave above which name is polyfor example :

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

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

So we have now a matrix of Size L x H filled with 0, we put now 1 at polygon points positions

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

I think you can simple do that

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

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

we interpret this as a binary (black/white) image which have a contour drawn on it Assume we want to detect this new contour

我们将其解释为一个二值(黑/白)图像,上面绘制了一个轮廓假设我们想要检测这个新的轮廓

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

Note : poly2is containing a list of points of your polygon and all points forming it, i mean all points of each vertices of your polygon which is what you need can find useful !! you can use cv2.CHAIN_APPROX_SIMPLEparameter to get poly2 containing only end points of the polygon lines which is lighter and which was our input :) important:the type of poly2 is numpy array ,its shape is (n,1,2) and not (n,2)

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

Now we draw this contour on this image(matrix) and will fill it too :)

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

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

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

now we have a matrix where there is 1on each points forming and filling the polygon , "thickness=-1" has forced to fill this contour, you can put set thickness = 1 to get only the borders if you want to translate, you can do it by adding parameter offset(xtran,ytrans)

现在我们有一个矩阵,其中形成和填充多边形的每个点上都有1个,“thickness=-1”已强制填充此轮廓,如果要翻译,您可以设置厚度 = 1 以仅获取边框,您可以通过添加参数offset(xtran,ytrans) 来实现

to get the indices of all theses points simply call

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

list_of_points_indices=numpy.nonzero(matrix)


Solution 2

解决方案2

Which is smarter is to directly transform your list of points (poly) to a contour format (poly2) and draw it on the matrix

哪个更聪明是直接将您的点列表(poly)转换为轮廓格式(poly2)并将其绘制在矩阵上

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

and draw it on the Matrix matrix

并将其绘制在 Matrix 矩阵上

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

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

And get the list of this points with :

并使用以下命令获取此点的列表:

list_of_points_indices=numpy.nonzero(matrix)

Play with thicknessto fill or not the polygon , see solution 1 for more details.

使用厚度来填充或不填充多边形,有关详细信息,请参阅解决方案 1。

回答by Ulf Aslak

Building upon RemcoGerlich's answer here's a validated function:

基于 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))]

Example:

例子:

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

enter image description here

在此处输入图片说明

回答by Ishan Tomar

Try this code. poly_coords are the coordinates of your polygon, 'coord' is coordinate of point you want to check if it is inside the polygon.

试试这个代码。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

If you want to make it even faster, take out the part of algo related to polygon, slope and offset and run the rest of code using 'map' function. Something like this:

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

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)

You can remove 'print' commands if you want. This code is made for python 2.7

如果需要,您可以删除“打印”命令。此代码是为 python 2.7 制作的

回答by Stanpol

I suggest to use matplotlib contains_points()

我建议使用 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