多边形等分

多边形等分思路

前提条件

• 封闭面，不可以有空洞

解题思路

• 封闭面中随机构造点
• 利用 k-means 分组（在此方法中设置分组数量 簇数）
• 计算每一个簇的质心
• 利用质心绘制 voronoi 泰森多边形
• 利用封闭面切割泰森多边形

注

本文中判断点是否在面中利用了 这篇文章

实现

#! /usr/bin/env python
# -*- coding: utf-8 -*-
import shapefile as shp
import matplotlib.pyplot as plt
import random
import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial import Voronoi, voronoi_plot_2d


def isPoiWithinPoly(poi, poly):
    """ 判断是否在面内 :param poi: [x,y] :param poly: [ [ [x1,y1],[x2,y2],[] ] ] :return: boolean """
    sinsc = 0
    for epoly in poly:
        for i in range(len(epoly) - 1):
            s_poi = epoly[i]
            e_poi = epoly[i + 1]
            if isRayIntersectsSegment(poi, s_poi, e_poi):
                sinsc += 1

    return True if sinsc % 2 == 1 else False


def isRayIntersectsSegment(poi, s_poi, e_poi):
    if s_poi[1] == e_poi[1]:
        return False
    if s_poi[1] > poi[1] and e_poi[1] > poi[1]:
        return False
    if s_poi[1] < poi[1] and e_poi[1] < poi[1]:
        return False
    if s_poi[1] == poi[1] and e_poi[1] > poi[1]:
        return False
    if e_poi[1] == poi[1] and s_poi[1] > poi[1]:
        return False
    if s_poi[0] < poi[0] and e_poi[1] < poi[1]:
        return False
    xseg = e_poi[0] - (e_poi[0] - s_poi[0]) * (e_poi[1] - poi[1]) / (e_poi[1] - s_poi[1])  # 求交
    if xseg < poi[0]:
        return False
    return True


# 开始
sf = shp.Reader("浙江省.shp")

plt.figure()

x = None
y = None

# 将SHP中的图形展示到Plots中
for shape in sf.shapeRecords():
    for i in range(len(shape.shape.parts)):
        i_start = shape.shape.parts[i]
        if i == len(shape.shape.parts) - 1:
            i_end = len(shape.shape.points)
        else:
            i_end = shape.shape.parts[i + 1]
        x = [i[0] for i in shape.shape.points[i_start:i_end]]
        y = [i[1] for i in shape.shape.points[i_start:i_end]]
        plt.plot(x, y)

# 获取当前面的最大最小XY
minX = min(x)
maxX = max(x)
minY = min(y)
maxY = max(y)

# 随机生成一定数量的散点 在这个平面内
# step 总共要多少个点
step = 10000
# count 记录当前点数量
count = 1

point_random = []
for i in range(10000):
    rx = random.uniform(minX, maxX)
    ry = random.uniform(minY, maxY)
    if isPoiWithinPoly([rx, ry], [list(zip(x, y))]) and count < step:
        count += 1
        point_random.append([rx, ry])
        # plt.scatter(rx, ry, s=20, c='aqua', alpha=0.5)

print("当前共有", count, "点")
print("随机点", point_random)
# kmean 簇族 分组量
group_size = 10
clf = KMeans(n_clusters=group_size)
cluster_group = clf.fit_predict(point_random)
# 遍历构造新的随机点 + 分组 [x,y,group_id]
point_kmean = []
# 设置颜色集合
cnames = ['black', 'blue', 'red', 'khaki', 'aliceblue', 'antiquewhite', 'aqua', 'aquamarine', 'azure', 'beige',
          'bisque', 'blanchedalmond', 'blueviolet', 'brown', 'burlywood', 'cadetblue', 'chartreuse', 'chocolate',
          'coral', 'cornflowerblue', 'cornsilk', 'crimson', 'cyan', 'darkblue', 'darkcyan', 'darkgoldenrod', 'darkgray',
          'darkgreen', 'darkkhaki', 'darkmagenta', 'darkolivegreen', 'darkorange', 'darkorchid', 'darkred',
          'darksalmon', 'darkseagreen', 'darkslateblue', 'darkslategray', 'darkturquoise', 'darkviolet', 'deeppink',
          'deepskyblue', 'dimgray', 'dodgerblue', 'firebrick', 'floralwhite', 'forestgreen', 'fuchsia', 'gainsboro',
          'ghostwhite', 'gold', 'goldenrod', 'gray', 'green', 'greenyellow', 'honeydew', 'hotpink', 'indianred',
          'indigo', 'ivory', 'khaki', 'lavender', 'lavenderblush', 'lawngreen', 'lemonchiffon', 'lightblue',
          'lightcoral', 'lightcyan', 'lightgoldenrodyellow', 'lightgreen', 'lightgray', 'lightpink', 'lightsalmon',
          'lightseagreen', 'lightskyblue', 'lightslategray', 'lightsteelblue', 'lightyellow', 'lime', 'limegreen',
          'linen', 'magenta', 'maroon', 'mediumaquamarine', 'mediumblue', 'mediumorchid', 'mediumpurple',
          'mediumseagreen', 'mediumslateblue', 'mediumspringgreen', 'mediumturquoise', 'mediumvioletred',
          'midnightblue', 'mintcream', 'mistyrose', 'moccasin', 'navajowhite', 'navy', 'oldlace', 'olive', 'olivedrab',
          'orange', 'orangered', 'orchid', 'palegoldenrod', 'palegreen', 'palevioletred', 'papayawhip', 'peachpuff',
          'peru', 'pink', 'plum', 'powderblue', 'purple', 'red', 'rosybrown', 'royalblue', 'saddlebrown', 'salmon',
          'sandybrown', 'seagreen', 'seashell', 'sienna', 'silver', 'skyblue', 'slateblue', 'slategray', 'snow',
          'springgreen', 'steelblue', 'tan', 'teal', 'thistle', 'tomato', 'turquoise', 'violet', 'wheat', 'white',
          'whitesmoke', 'yellow', 'yellowgreen']

for point, gp_id in zip(point_random, cluster_group):
    point_kmean.append([point, gp_id])
    # 放到 plt 中展示
    plt.scatter(point[0], point[1], s=5, c=cnames[gp_id], alpha=1)

# 通过 kmean 质心绘制泰森多边形
cluster_center = clf.cluster_centers_
print("质心")
print(cluster_center)

for i in cluster_center:
    plt.scatter(i[0], i[1], s=25, c='green', alpha=1)

# 泰森多边形加载到plot上
points = np.array(cluster_center)
vor = Voronoi(points)

voronoi_plot_2d(vor)
for region in vor.regions:
    if not -1 in region:
        polygon = [vor.vertices[i] for i in region]
        plt.fill(*zip(*polygon))

plt.show()