shapely for python で、polygonをmesh化後、中心線を算出

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path        import Path
from matplotlib.patches     import PathPatch
from matplotlib.collections import PatchCollection
from matplotlib.lines       import Line2D

# cf. https://shapely.readthedocs.io/en/stable/manual.html
from shapely.geometry       import Polygon
from shapely.geometry       import MultiPolygon


def main():
    polygon =  Polygon(
        shell=[ (10,10),(10,40),(70,40),(70,10),(60,10),
                (60,30),(50,30),(50,10),(10,10) ],
        holes=[ [(20,20),(40,20),(40,30),(20,30),(20,20)] ]
    )

    meshes = divide_by_mesh(polygon,10)  # メッシュ化
    mid_lines = calc_mid_lines( meshes ) # 中心線算出
    
    fig, ax = plt.subplots()

    x_s = []
    y_s = []
    for mid_line in mid_lines:
        ax.add_line( Line2D([mid_line[0],mid_line[2]],
                            [mid_line[1],mid_line[3]] ))
    ax.set_xlim(10, 70)
    ax.set_ylim(10, 40)

    # for mesh in meshes:
    #     plot_polygon(ax, mesh, facecolor='lightblue', edgecolor='red')
    plt.show()


def calc_mid_lines( meshes ):
    mid_lines = []

    for mesh in meshes:
        # 矩形でない場合、pass
        if len( list( mesh.exterior.coords ) ) > 5:
            continue
        (min_x,min_y, max_x,max_y) = mesh.bounds
        mid_x = (min_x+max_x) / 2
        mid_y = (min_y+max_y) / 2

        mid_lines += [(min_x,mid_y,mid_x,mid_y),
                      (mid_x,mid_y,max_x,mid_y),
                      (mid_x,min_y,mid_x,mid_y),
                      (mid_x,mid_y,mid_x,max_y)]

    return mid_lines
    
def divide_by_mesh(polygon,mesh_size):
    
    (min_x,y, max_x,max_y) = polygon.bounds
    print( min_x,y, max_x,max_y )
    x = min_x
    meshes = []
    
    while y < max_y:
        while x < max_x:
            mesh = Polygon(
                [ (x, y+mesh_size),(x+mesh_size,y+mesh_size),
                  (x+mesh_size,y ),(x, y)] )
            
            intersects = polygon.intersection(mesh)
            if intersects.type == "Polygon":
                meshes.append( intersects )
            else:
                for intersect in intersects.geoms:
                    if intersect.type != "Polygon":
                        continue
                    meshes.append( intersect )
                    
            x += mesh_size
        x = min_x
        y += mesh_size
    return meshes
    
# Plots a Polygon to pyplot `ax`
# cf. https://stackoverflow.com/questions/55522395
def plot_polygon(ax, poly, **kwargs):
    path = Path.make_compound_path(
        Path(np.asarray(poly.exterior.coords)[:, :2]),
        *[Path(np.asarray(ring.coords)[:, :2]) for ring in poly.interiors])

    patch = PathPatch(path, **kwargs)
    collection = PatchCollection([patch], **kwargs)
    
    ax.add_collection(collection, autolim=True)
    ax.autoscale_view()
    return collection

if __name__ == '__main__':
    main()