把封闭区域分割成指定的份数,可以使用泰森多边形进行解决。泰森多边形又叫维诺图(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;
}

三、切分结果