将多边形细分为较小多边形的算法

Posted 2023-02-19

【中文标题】将多边形细分为较小多边形的算法

【英文标题】：Algorithm to subdivide a polygon in smaller polygons

【发布时间】：2012-01-19 11:40:36

【问题描述】：

我有一个由平面上的连续边组成的多边形，并且想将其细分为三角形或矩形的子多边形。 我在哪里可以找到一个算法来做到这一点？ 谢谢！

【问题讨论】：

您能否发布一张图片，并附上您正在尝试做的事情的示例？在绘制图片的过程中，您很有可能会提出一种可接受的算法。

抱歉不清楚。基本上我正在开发一个生成街道和地块的应用程序，地块是街道交叉口内的飞机。我想将地块细分为较小的地块，以便将房屋放在较小的地块上。

@LaurentCrivello：正如missigno 提到的，这是一个已知问题。如果您正在寻找示例代码，here's Haskell 中的一个简短实现（查看triangulatePolygon）

你可能会觉得this Soulwire blog post很有趣。

【参考方案1】：

经过多次搜索偶然发现。

感谢@Aditya Chhabra 的提交，它工作得很好，但get_squares_from_rect 由于迭代剪辑，对于小边长来说非常慢。

如果我们将所有 LineStrings 合并到一个集合中，然后一步进行剪辑和多边形化，我们可以立即做到这一点，我在 this question 中找到了这一点。

以前边长 0.0001 (EPSG:4326) 需要 > 1 分钟，现在不需要时间。

from shapely.ops import unary_union, polygonize, linemerge
from shapely.geometry import LineString
import numpy as np


def get_squares_from_rect_faster(RectangularPolygon, side_length=0.0025):
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)

    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]

    lines = horizontal_splitters + vertical_splitters
    lines.append(RectangularPolygon.boundary)
    lines = unary_union(lines)
    lines = linemerge(lines)

    return list(polygonize(lines))

【参考方案2】：

我自己也在寻找这个问题的答案，但找不到。尝试将几块拼接在一起，结果如下。

这不一定是最理想的例程，但它为我完成了工作。如果您想提高性能，请尝试使用代码。

算法的简要说明：

使用原始几何本身的边界、凸包的边界和最小旋转矩形，导出所有可能的矩形。

将所有矩形分成指定边长的小方块。

删除重复项使用四舍五入的质心。 (r: 四舍五入参数)

保留那些“在”几何体中的方格，或者那些“相交”几何体的方格，取决于哪个更接近所需方格的总数。

已编辑

新解决方案

#### Python script for dividing any shapely polygon into smaller equal sized polygons

import numpy as np
from shapely.ops import split
import geopandas
from shapely.geometry import MultiPolygon, Polygon

def rhombus(square):
    """
    Naively transform the square into a Rhombus at a 45 degree angle
    """
    coords = square.boundary.coords.xy
    xx = list(coords[0])
    yy = list(coords[1])
    radians = 1
    points = list(zip(xx, yy))
    Rhombus = Polygon(
        [
            points[0],
            points[1],
            points[3],
            ((2 * points[3][0]) - points[2][0], (2 * points[3][1]) - points[2][1]),
            points[4],
        ]
    )
    return Rhombus


def get_squares_from_rect(RectangularPolygon, side_length=0.0025):
    """
    Divide a Rectangle (Shapely Polygon) into squares of equal area.

    `side_length` : required side of square

    """
    rect_coords = np.array(RectangularPolygon.boundary.coords.xy)
    y_list = rect_coords[1]
    x_list = rect_coords[0]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells = int(np.round(width / side_length))
    ycells = int(np.round(height / side_length))

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = RectangularPolygon
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)

    return square_polygons


def split_polygon(G, side_length=0.025, shape="square", thresh=0.9):
    """
    Using a rectangular envelope around `G`, creates a mesh of squares of required length.
    
    Removes non-intersecting polygons. 
            

    Args:
    
    - `thresh` : Range - [0,1]

        This controls - the number of smaller polygons at the boundaries.
        
        A thresh == 1 will only create (or retain) smaller polygons that are 
        completely enclosed (area of intersection=area of smaller polygon) 
        by the original Geometry - `G`.
        
        A thresh == 0 will create (or retain) smaller polygons that 
        have a non-zero intersection (area of intersection>0) with the
        original geometry - `G` 

    - `side_length` : Range - (0,infinity)
        side_length must be such that the resultant geometries are smaller 
        than the original geometry - `G`, for a useful result.

        side_length should be >0 (non-zero positive)

    - `shape` : square/rhombus
        Desired shape of subset geometries. 


    """
    assert side_length>0, "side_length must be a float>0"
    Rectangle    = G.envelope
    squares      = get_squares_from_rect(Rectangle, side_length=side_length)
    SquareGeoDF  = geopandas.GeoDataFrame(squares).rename(columns=0: "geometry")
    Geoms        = SquareGeoDF[SquareGeoDF.intersects(G)].geometry.values
    if shape == "rhombus":
        Geoms = [rhombus(g) for g in Geoms]
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    elif shape == "square":
        geoms = [g for g in Geoms if ((g.intersection(G)).area / g.area) >= thresh]
    return geoms

# Reading geometric data

geo_filepath = "/data/geojson/pc_14.geojson"
GeoDF        = geopandas.read_file(geo_filepath)

# Selecting random shapely-geometry

G = np.random.choice(GeoDF.geometry.values)


squares   = split_polygon(G,shape='square',thresh=0.5,side_length=0.025)
rhombuses = split_polygon(G,shape='rhombus',thresh=0.5,side_length=0.025)

以前的解决方案：

import numpy as np
import geopandas
from shapely.ops import split
from shapely.geometry import MultiPolygon, Polygon, Point, MultiPoint

def get_rect_from_geom(G, r=2):
    """
    Get rectangles from a geometry.
    r = rounding factor.
    small r ==> more rounding off ==> more rectangles
    """
    coordinate_arrays = G.exterior.coords.xy
    coordinates = list(
        zip(
            [np.round(c, r) for c in coordinate_arrays[0]],
            [np.round(c, r) for c in coordinate_arrays[1]],
        )
    )
    Rectangles = []
    for c1 in coordinates:
        Coords1 = [a for a in coordinates if a != c1]
        for c2 in Coords1:
            Coords2 = [b for b in Coords1 if b != c2]
            x1, y1 = c1[0], c1[1]
            x2, y2 = c2[0], c2[1]
            K1 = [k for k in Coords2 if k == (x1, y2)]
            K2 = [k for k in Coords2 if k == (x2, y1)]
            if (len(K1) > 0) & (len(K2) > 0):
                rect = [list(c1), list(K1[0]), list(c2), list(K2[0])]
                Rectangles.append(rect)
    return Rectangles


def get_squares_from_rect(rect, side_length=0.0025):
    """
    Divide a rectangle into equal area squares
    side_length = required side of square
    """
    y_list = [r[1] for r in rect]
    x_list = [r[0] for r in rect]
    y1 = min(y_list)
    y2 = max(y_list)
    x1 = min(x_list)
    x2 = max(x_list)
    width = x2 - x1
    height = y2 - y1

    xcells, ycells = int(np.round(width / side_length)), int(
        np.round(height / side_length)
    )

    yindices = np.linspace(y1, y2, ycells + 1)
    xindices = np.linspace(x1, x2, xcells + 1)
    horizontal_splitters = [
        LineString([(x, yindices[0]), (x, yindices[-1])]) for x in xindices
    ]
    vertical_splitters = [
        LineString([(xindices[0], y), (xindices[-1], y)]) for y in yindices
    ]
    result = Polygon(rect)
    for splitter in vertical_splitters:
        result = MultiPolygon(split(result, splitter))
    for splitter in horizontal_splitters:
        result = MultiPolygon(split(result, splitter))
    square_polygons = list(result)
    return [np.stack(SQPOLY.exterior.coords.xy, axis=1) for SQPOLY in square_polygons]


def round_centroid(g, r=10):
    """
    Get Centroids.
    Round off centroid coordinates to `r` decimal points.
    """
    C = g.centroid.coords.xy
    return (np.round(C[0][0], r), np.round(C[1][0], r))


def subdivide_polygon(g, side_length=0.0025, r=10):
    """
    1. Create all possible rectangles coordinates from the geometry, its minimum rotated rectangle, and its convex hull.
    2. Divide all rectangles into smaller squares.

    small r ==> more rounding off ==> fewer overlapping squares. (these are dropped as duplicates)
    large r ==> may lead to a few overlapping squares.
    """
    #   Number of squares required.
    num_squares_reqd = g.area // (side_length ** 2)

    # Some of these combinations can be dropped to improve performance.
    Rectangles = []
    Rectangles.extend(get_rect_from_geom(g))
    Rectangles.extend(get_rect_from_geom(g.minimum_rotated_rectangle))
    Rectangles.extend(get_rect_from_geom(g.convex_hull))

    Squares = []
    for r in range(len(Rectangles)):
        rect = Rectangles[r]
        Squares.extend(get_squares_from_rect(rect, side_length=side_length))
    SquarePolygons = [Polygon(square) for square in Squares]
    GDF = geopandas.GeoDataFrame(SquarePolygons).rename(columns=0: "geometry")
    GDF.loc[:, "centroid"] = GDF.geometry.apply(round_centroid, r=r)
    GDF = GDF.drop_duplicates(subset=["centroid"])

    wgeoms = GDF[GDF.within(g)].geometry.values
    igeoms = GDF[GDF.intersects(g)].geometry.values

    w = abs(num_squares_reqd - len(wgeoms))
    i = abs(num_squares_reqd - len(igeoms))

    print(w, i)
    if w <= i:
        return wgeoms
    else:
        return igeoms



geoms = subdivide(g)

【讨论】：

这看起来真不错！不过我有一个问题：当我将它用于我的数据 (POLYGON ((122.395 32.221, 122.395 32.221, 122.34 35.444, 193.444 35.000 ) 时，我收到错误 AttributeError: 'GeoSeries' object has no attribute 'coords'。您输入为g 的几何图形是什么样的？

@SergedeGossondeVarennes，根据我的评论，您可能使用的是 Geoseries（geopandas 库的一个类），而不是 shapely 库中的几何图形。上述函数中的 g 采用单个 Shapely Polygon 对象。如果你能分享一些你的代码，我可以说更多。

我发现并纠正了它。它对某些几何图形非常有效。但对于某些人来说，它没有，我实际上为此创建了一个问题。 ***.com/q/67965664/5363686 。在某些情况下，get_rect_from_geom(G, r=1) 返回一个空字段。你经历过吗？顺便说一句....真的很棒的解决方案。

@SergedeGossondeVarennes - 我已经编辑了我的答案。这要简单得多，而且速度要快得多。希望这会有所帮助。

它就像一个魅力！各种形状！ @Aditya Chhabra，纯粹的天才！

【参考方案3】：

在computational geometry中，你要解决的问题叫做triangulation。

有一些算法可以解决这个问题，给出具有不同属性的三角剖分。您需要决定哪一个最合适。

【讨论】：

谢谢。然而，以三角形结尾并不是我的最终目标，因为矩形更符合我的定义。不过我还是看看吧，谢谢！