基于K-means 切割多边形 JAVA实现

共计 8246 个字符，预计需要花费 21 分钟才能阅读完成。

基于 K -means 切割多边形 JAVA 实现
思路初稿详见多边形等分
依赖

geotools
ekmeans

<!–geotools–>
<!– https://mvnrepository.com/artifact/org.locationtech.jts/jts-core –>
<dependency>
<groupId>org.locationtech.jts</groupId>
<artifactId>jts-core</artifactId>
<version>1.16.0</version>
</dependency>
<dependency>
<groupId>junit</groupId>
<artifactId>junit</artifactId>
<version>4.12</version>
<scope>test</scope>
</dependency>

<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-shapefile</artifactId>
<version>${geotools.version}</version>
</dependency>

<!– https://mvnrepository.com/artifact/org.geotools/gt-main –>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-main</artifactId>
<version>${geotools.version}</version>
</dependency>

<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-swing</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geojson</artifactId>
<version>${geotools.version}</version>
</dependency>
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geometry</artifactId>
<version>${geotools.version}</version>
</dependency>
<!– https://mvnrepository.com/artifact/com.vividsolutions/jts –>
<dependency>
<groupId>com.vividsolutions</groupId>
<artifactId>jts</artifactId>
<version>1.13</version>
</dependency>

<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geojson</artifactId>
<version>${geotools.version}</version>
</dependency>
<!–kmeans 实现 –>
<dependency>
<groupId>ca.pjer</groupId>
<artifactId>ekmeans</artifactId>
<version>2.0.0</version>
<scope>compile</scope>
</dependency>
实现过程
结果类 KmeanPolygonResult
package com.huifer.planar.asEt.entity;

import java.util.ArrayList;
import java.util.List;
import lombok.AllArgsConstructor;
import lombok.Data;
import lombok.NoArgsConstructor;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;

/**
* <p>Title : KmeanPolygonResult </p>
* <p>Description : Kmean+ 泰森多边形平分多边形结果 </p>
*
* @author huifer
* @date 2019-01-16
*/
@Data
@NoArgsConstructor
@AllArgsConstructor
public class KmeanPolygonResult {

/**
* 原始平面
*/
private Polygon polygon;
/**
* 随机点集合
*/
private ArrayList<Point> pointList;
/**
* 分组后组 id
*/
private int[] assignments;
/**
* 中心集合
*/
private double[][] centroids;
/**
* xlist
*/
private ArrayList<Double> xlist;
/**
* ylist
*/
private ArrayList<Double> ylist;
/**
* 构造泰森多边形
*/
private List<Geometry> voronoi;
}

k-means 聚合类
package com.huifer.planar.asEt.utils;

import ca.pjer.ekmeans.EKmeans;
import lombok.Data;

/**
* <p>Title : Kmeans </p>
* <p>Description : Kmeans 聚合 </p>
*
* @author huifer
* @date 2019-01-15
*/
@Data
public class Kmeans {

/***
* 聚合总量
*/
private int n;
/***
* 数据集合
*/
private double[][] points;
/**
* 簇族数量
*/
private int k;

/**
* 中心集合
*/
private double[][] centroids;
/**
* 分组后组 id
*/
private int[] assignments;

public Kmeans(double[][] points, int k) {
this.n = points.length;
this.points = points;
this.k = k;
this.centroids = new double[k][2];
init();
}

/**
* 求得分组数据
*/
public void init() {
EKmeans eKmeans = new EKmeans(centroids, points);
eKmeans.setIteration(128);
eKmeans.setEqual(true);
eKmeans.setDistanceFunction(EKmeans.EUCLIDEAN_DISTANCE_FUNCTION);
eKmeans.run();
int[] assignments = eKmeans.getAssignments();
this.assignments = assignments;
double[][] centroids = eKmeans.getCentroids();
this.centroids = centroids;
}

}
泰森多边形构造
package com.huifer.planar.asEt.algo.impl;

import com.huifer.planar.asEt.algo.VoronoiInterface;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;
import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;

/**
* <p>Title : VoronoiInterfaceImpl </p>
* <p>Description : 泰森多边形 & 德劳内三角形 </p>
*
* @author huifer
* @date 2019-01-16
*/
public class VoronoiInterfaceImpl implements VoronoiInterface {

@Override
public List<Geometry> voronoi(double[][] doubles) {
return voronoi(doublesToCoordinate(doubles));
}

@Override
public Collection delaunay(double[][] doubles) {
return delaunay(doublesToCoordinate(doubles));
}

@Override
public List<Geometry> voronoi(ArrayList<Point> points) {
List<Geometry> voronoi = voronoi(pointToCoordinate(points));
return voronoi;
}

@Override
public Collection delaunay(ArrayList<Point> points) {
Collection delaunay = delaunay(pointToCoordinate(points));
return delaunay;
}

private List<Coordinate> doublesToCoordinate(double[][] doubles) {
List<Coordinate> coords = new ArrayList<Coordinate>();
for (int i = 0; i < doubles.length; i++) {
Coordinate coord = new Coordinate(doubles[i][0], doubles[i][1], i);
coords.add(coord);
}
return coords;
}

private List<Coordinate> pointToCoordinate(ArrayList<Point> points) {
List<Coordinate> coords = new ArrayList<Coordinate>();
for (int i = 0; i < points.size(); i++) {
Coordinate coord = new Coordinate(points.get(i).getX(), points.get(i).getY(), i);
coords.add(coord);
}
return coords;
}

private VoronoiDiagramBuilder getVoronoiDiagramBuilder(List<Coordinate> coords) {
VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder();
Envelope clipEnvelpoe = new Envelope();
voronoiDiagramBuilder.setSites(coords);
voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe);
return voronoiDiagramBuilder;
}

@Override
public List<Geometry> voronoi(List<Coordinate> coords) {
VoronoiDiagramBuilder voronoiDiagramBuilder = getVoronoiDiagramBuilder(
coords);
Geometry geom = voronoiDiagramBuilder
.getDiagram(JTSFactoryFinder.getGeometryFactory());

List<Geometry> voronoi = new ArrayList<>();
int numGeometries = geom.getNumGeometries();
for (int i = 0; i < numGeometries; i++) {
Geometry geometryN = geom.getGeometryN(i);
voronoi.add(geometryN);
}
return voronoi;
}

@Override
public Collection delaunay(List<Coordinate> coords) {
VoronoiDiagramBuilder voronoiDiagramBuilder = getVoronoiDiagramBuilder(coords);
QuadEdgeSubdivision subdivision = voronoiDiagramBuilder.getSubdivision();
Collection delaunay = subdivision.getEdges();
return delaunay;
}
}

Kmeans 切割面类
package com.huifer.planar.asEt.algo.impl;

import com.huifer.planar.asEt.algo.KmeanPolygonSplitInterface;
import com.huifer.planar.asEt.algo.VoronoiInterface;
import com.huifer.planar.asEt.entity.KmeanPolygonResult;
import com.huifer.planar.asEt.utils.Kmeans;
import com.huifer.planar.asEt.utils.shptools.overlay.Operation;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.io.ParseException;

/**
* <p>Title : KmeanPolygonSplitCore </p>
* <p>Description : Kmean split polygon core </p>
*
* @author huifer
* @date 2019-01-16
*/
public class KmeanPolygonSplitCore implements KmeanPolygonSplitInterface {

private static double random(double max, double min) {
double d = (Math.random() * (max – min) + min);
return d;
}

@Override
public KmeanPolygonResult splitPolygon(String wkt, int setp, int k) throws ParseException {
Operation op = new Operation();
KmeanPolygonResult result = new KmeanPolygonResult();
Polygon polygon = op.createPolygonByWKT(wkt);
Coordinate[] coordinates = polygon.getCoordinates();
ArrayList<Double> xList = new ArrayList<>();
ArrayList<Double> yList = new ArrayList<>();

Arrays.stream(coordinates).forEach(
s -> {
xList.add(s.x);
yList.add(s.y);
}
);
// xy 最大最小值
Double xMax = xList.stream().reduce(Double::max).get();
Double xMin = xList.stream().reduce(Double::min).get();
Double yMax = yList.stream().reduce(Double::max).get();
Double yMin = yList.stream().reduce(Double::min).get();

// 当前点数量
int pointCount = 0;
ArrayList<Point> pointArrayList = new ArrayList<>();

for (int i = 0; i < Integer.MAX_VALUE; i++) {
// 最大最小值随机
if (pointCount <= setp) {
double rx = random(xMax, xMin);
double ry = random(yMax, yMin);
Point nowPoint = op.createPointByWkt(“POINT(” + rx + ” ” + ry + “)”);
boolean contains = polygon.contains(nowPoint);
if (contains) {
pointArrayList.add(nowPoint);
pointCount++;
}
} else {
break;
}
}
// k-means 数据 构造
double[][] kmData = new double[pointArrayList.size()][2];
for (int i = 0; i < pointArrayList.size(); i++) {
Point point = pointArrayList.get(i);
double[] oneData = new double[2];
oneData[0] = point.getX();
oneData[1] = point.getY();
kmData[i] = oneData;
}

// k-means 结果

Kmeans kmeans = new Kmeans(kmData, k);

// 构造泰森多边形
VoronoiInterface vo = new VoronoiInterfaceImpl();
List<Geometry> voronoi = vo.voronoi(kmeans.getCentroids());

result.setPolygon(polygon);
result.setPointList(pointArrayList);
result.setAssignments(kmeans.getAssignments());
result.setCentroids(kmeans.getCentroids());
result.setXlist(xList);
result.setYlist(yList);
result.setVoronoi(voronoi);
return result;
}

}

结果展示

注
本文代码及可视化代码均放在 gitee 码云上 欢迎 star & fork