多边形封闭区域分割方法
把封闭区域分割成指定的份数，可以使用泰森多边形进行解决。泰森多边形又叫维诺图（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;
}
三、切分结果
