多边形封闭区域分割方法
把封闭区域分割成指定的份数,可以使用泰森多边形进行解决。泰森多边形又叫维诺图(Voronoi Diagram),它是对空间平面的一种剖分,其特点是多边形内的任何点到该多边形的控制点的距离最近,而离相邻多边形内控制点的距离较远,且每个多边形内有且仅有一个控制点(具体如下图所示)。

因此,对于前述的建筑物外轮廓,如果建筑物在某层要分成n部分,我们可以随机在其中找n个点做为控制点,然后就可以使用泰森多边形方法将它分成n部分。这就解决我们的封闭区域分割问题。
我们想实现尽可能平均的分割,那么n个控制点就不能随机选取。对于整个封闭区域,可以看成是若干点的集合,我们可以使用k-means算法计算n个聚类中心,将这个点集分成n个子集。这n个聚类中心实际是每个子集的质心,以它们为维诺图中的控制点,分割出的n个部分将更加平均。
上述算法还有一个问题,即如何随机选取若干样点组成点集,来表示封闭区域。我们可以计算封闭区域的最小外包矩形,这个最小外包矩形可以确定选取样点的横坐标和纵坐标的范围。然后在此范围内随机给出样点的横坐标和纵坐标,这样生成的样点一定在最小外包矩形内,但不一定都在封闭区域内。这里我们使用射线法判断生成的点是否在封闭区域内。
在计算机表达中,封闭区域的外轮廓实际是一个多边形。射线法以判断点开始,向右(或向左)的水平方向作一射线,计算该射线与多边形每条边的交点个数,如果交点个数为奇数,则点位于多边形内,偶数则在多边形外。该算法对于复合多边形也能正确判断(如下图所示)。

综上所述,我们的算法如下:
- 在封闭区域内取m个样点表示整个区域(使用射线法判断随机样点是否在封闭区域内);
- 使用k-means算法对m个样点进行聚成n类,求出每个子类的聚类中心点;
- 将n个聚类中心点做为维诺图的控制点,生成封闭区域的维诺图;
- 根据生成维诺图分割封闭区域,得到每个分割部分外轮廓的点集合。
具体实现:
一、引入相关依赖包
<!-- geo tools --> <dependency> <groupId>org.locationtech.jts</groupId> <artifactId>jts-core</artifactId> <version>1.16.1</version> </dependency> <dependency> <groupId>org.geotools</groupId> <artifactId>gt-geojson</artifactId> <version>21.2</version> </dependency> <!-- k-means --> <dependency> <groupId>ca.pjer</groupId> <artifactId>ekmeans</artifactId> <version>2.0.0</version> </dependency>
二、实现代码
/** * 多边形平分 * * @param polygon 多边形 * @param count 平分份数 * @param random 点集数量(kmeans算法用,越多越精确但速度越慢) */ public static MultiPolygon splitPolygon(Polygon polygon, int count, int random) { GeometryFactory geometryFactory = new GeometryFactory(); // 1. 构建随机点 double minY = 90; double maxY = -90; double minX = 180; double maxX = -180; for (Coordinate point : polygon.getCoordinates()) { // 最小 if (point.getY() < minY) { minY = point.getY(); } if (point.getX() < minX) { minX = point.getX(); } // 最大 if (point.getY() > maxY) { maxY = point.getY(); } if (point.getX() > maxX) { maxX = point.getX(); } } double[][] points = new double[random][2]; int index = 0; do { double x = RandomUtils.nextDouble(minX, maxX); double y = RandomUtils.nextDouble(minY, maxY); if (polygon.contains(geometryFactory.createPoint(new Coordinate(x, y)))) { points[index++] = new double[]{x, y}; } } while (index < random); // 2. 利用EKmeans 获取分组和簇的质心 double[][] centroids = new double[count][2]; EKmeans eKmeans = new EKmeans(centroids, points); eKmeans.setEqual(true); // 替换距离计算方法,使用基于经纬度的距离计算 eKmeans.setDistanceFunction((p1, p2) -> new Coordinate(p1[0], p1[1]).distance(new Coordinate(p1[0], p1[1]))); eKmeans.run(); // 获取分组id和中心聚合 centroids = eKmeans.getCentroids(); // 3. 构建泰森多边形 List<Coordinate> coords = new ArrayList<>(); for (double[] p : centroids) { coords.add(new Coordinate(p[0], p[1])); } VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder(); Envelope clipEnvelpoe = new Envelope(); voronoiDiagramBuilder.setSites(coords); voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe); Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory()); // 4. 利用封闭面切割泰森多边形 Coordinate[] coordinates = polygon.getCoordinates(); Polygon pyg = geometryFactory.createPolygon(coordinates); List<Geometry> geometries = new ArrayList<>(); for (int i = 0; i < geom.getNumGeometries(); i++) { Geometry geometry = geom.getGeometryN(i); geometries.add(pyg.intersection(geometry)); } MultiPolygon multiPolygon = geometryFactory.createMultiPolygon(geometries.toArray(new Polygon[0])); return multiPolygon; } /** * Geometry平分 * * @param geometry Geometry * @param count 平分份数 * @param random 点集数量(kmeans算法用,越多越精确但速度越慢) * @return */ public static List<Geometry> splitGeometry(Geometry geometry, int count, int random) { GeometryFactory geometryFactory = new GeometryFactory(); // 1. 构建随机点 double minY = 90; double maxY = -90; double minX = 180; double maxX = -180; for (Coordinate point : geometry.getCoordinates()) { // 最小 if (point.getY() < minY) { minY = point.getY(); } if (point.getX() < minX) { minX = point.getX(); } // 最大 if (point.getY() > maxY) { maxY = point.getY(); } if (point.getX() > maxX) { maxX = point.getX(); } } double[][] points = new double[random][2]; int index = 0; do { double x = RandomUtils.nextDouble(minX, maxX); double y = RandomUtils.nextDouble(minY, maxY); if (geometry.contains(geometryFactory.createPoint(new Coordinate(x, y)))) { points[index++] = new double[]{x, y}; } } while (index < random); // 2. 利用EKmeans 获取分组和簇的质心 double[][] centroids = new double[count][2]; EKmeans eKmeans = new EKmeans(centroids, points); eKmeans.setEqual(true); // 替换距离计算方法,使用基于经纬度的距离计算 eKmeans.setDistanceFunction((p1, p2) -> Coordinates.getDistance(new CoordinatesPoint(p1[0], p1[1]), new CoordinatesPoint(p2[0], p2[1]))); eKmeans.run(); // 获取分组id和中心聚合 centroids = eKmeans.getCentroids(); // 3. 构建泰森多边形 List<Coordinate> coords = new ArrayList<>(); for (double[] p : centroids) { coords.add(new Coordinate(p[0], p[1])); } VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder(); Envelope clipEnvelpoe = new Envelope(); voronoiDiagramBuilder.setSites(coords); voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe); Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory()); // 4. 利用封闭面切割泰森多边形 List<Geometry> geometries = new ArrayList<>(); for (int i = 0; i < geom.getNumGeometries(); i++) { geometries.add(geometry.intersection(geom.getGeometryN(i))); } return geometries; }
三、切分结果

