java jts 针对shp含洞多边形进行三角剖分切分成可行区域

前言

java jts 提供了Delaunay三角剖分的相关方法,但是该方法不考虑含洞的多边形的。虽然 jts 的 ConformingDelaunayTriangulationBuilder 类可以通过提供线段约束的方式防止切割到洞内，但是仅支持最多99条线段，虽然可以通过重写破除99条线段的约束，但是线段多了性能会变得很差。且在网上没找到相关的解决方案。

本文主要是基于 jts 实现对含洞多边形进行三角剖分，并且可以不切割到洞内的方法。

三角剖分是生成寻路用的 navmesh 的步骤之一

原生工具切割效果展示

首先是初始shp文件，就是一系列的多边形障碍物，三角剖分时不能切割到障碍物内部。

将上面的shp取反后就是一个含洞的多边形，洞内就是上面的多边形障碍物。

切分后的效果如下图所示

放大局部后如下图所示，可以看出切割的时候是不考虑多边形障碍物的，会切割到多边形内部。

实现代码

pom.xml依赖

 <dependency><groupId>org.geotools</groupId><artifactId>gt-geotiff</artifactId><version>21.1</version>
</dependency>
<dependency><groupId>org.geotools</groupId><artifactId>gt-shapefile</artifactId><version>22.1</version>
</dependency>
<dependency><groupId>cn.hutool</groupId><artifactId>hutool-all</artifactId><version>5.7.19</version>
</dependency>

其他工具

QGIS：用于打开 shp 文件可视化看效果

代码

代码目前测试暂未发现还有什么问题，但不确保一点问题都没。

public class ShpUtils {public static final GeometryFactory GEOMETRY_FACTORY = JTSFactoryFinder.getGeometryFactory();public static void main(String[] args) throws Exception {List<Polygon> resultTriangles = genTriangle("D:\\shp\\test.shp", false, Lists.newArrayList());//将剖分的结果写入shp文件，可以用于看剖分效果writeMultiPolygonShp(GEOMETRY_FACTORY.createMultiPolygon(resultTriangles.toArray(new Polygon[resultTriangles.size()])), "D:\\shp\\test_triangle.shp");}/*** 三角剖分** @param shpPath     shp文件地址* @param isNegation  是否取反后再剖分* @param coordinates 外边框(地图的边界范围)，可以为空* @throws Exception*/public static List<Polygon> genTriangle(String shpPath, boolean isNegation, List<Coordinate> coordinates) throws Exception {SimpleFeatureSource featureSource = getFeatureSourceByPath(shpPath);List<MultiPolygon> multiPolygons = getPolygons(featureSource);Geometry triangles;if (!isNegation) {if (CollectionUtil.isEmpty(coordinates)) {//如果没传入外边框,则不设置外边框coordinates = Lists.newArrayList();//取shp外接矩形的外边框的话剖分结果应该是和取反后再三角剖分的效果是一样的,取反默认就是shp文件外接矩形对应的外边框
//            ReferencedEnvelope bounds = featureSource.getBounds();
//            List<Coordinate> coordinates = boundToCoordinate(bounds);}triangles = genTriangle(coordinates, (List) multiPolygons);} else {//取反后三角剖分Geometry negationShp = negationShp(featureSource, coordinates);triangles = genTriangle(Lists.newArrayList(), Lists.newArrayList(negationShp));}//最终剖分得到的三角形List<Polygon> resultTriangles = Lists.newArrayList();List<Polygon> polygons = multiPolygonToPolygon(multiPolygons);for (int i = 0; i < triangles.getNumGeometries(); i++) {//剖分后的三角形Polygon triangle = (Polygon) triangles.getGeometryN(i);List<Polygon> disposePolygons = disposeTriangle(triangle, polygons, 0);resultTriangles.addAll(disposePolygons);//注释掉上面两行使用下面这行即为原生的剖分效果
//            resultTriangles.add(triangle);}return resultTriangles;}/*** 对剖分的三角形存在下面两中情况的进行处理* 1.剖分的三角形可能在障碍物内部* 2.剖分的三角形可能和障碍物有交集，但没被完全包含** @param triangle* @param polygons* @param start* @return*/private static List<Polygon> disposeTriangle(Polygon triangle, List<Polygon> polygons, int start) {List<Polygon> resultTriangles = Lists.newArrayList();for (int i = start; i < polygons.size(); i++) {Polygon polygon = polygons.get(i);if (polygon.contains(triangle)) {//当剖分的三角形完全在障碍物多边形内部，则忽略掉该三角形,在障碍物内部的三角形是不可能为可行区域的return resultTriangles;} else if (polygon.intersects(triangle)) {//当剖分的三角形和障碍物多边形有交集,则相交部分必为不可行区域//取出当前三角形相对于障碍物多边形没有交集的部分,没有交集的部分可能为可行区域，需要对没有交集的部分进行处理Geometry difference = triangle.difference(polygon);if (difference instanceof MultiPolygon) {//障碍物多边形可能将三角形截断成多个多边形，需要对截断后的多个多边形分别进行处理for (int j = 0; j < difference.getNumGeometries(); j++) {//只需要对之后的障碍物多边形校验是否相交resultTriangles.addAll(disposeTriangle((Polygon) difference.getGeometryN(j), polygons, i + 1));}} else if (difference instanceof Polygon) {//如果没有交集的部分依然只是一个多边形,则检测该多边形与后续障碍物多边形是否有交集triangle = (Polygon) difference;} else if (difference instanceof GeometryCollection) {//同 MultiPolygonfor (int k = 0; k < difference.getNumGeometries(); k++) {Geometry geometryN = difference.getGeometryN(k);if (geometryN instanceof Polygon) {resultTriangles.addAll(disposeTriangle((Polygon) geometryN, polygons, i + 1));} else {//TODO 有可能得到的结果是非多边形，目前测试可能得到一条线段，线段应该不需要处理，不确定是否有其他情况System.out.println(MessageFormat.format("不可解析图形#1,{0}", geometryN));}}} else {//TODO 目前测试可能得到一条线段，不予处理System.out.println(MessageFormat.format("不可解析图形#2,{0}", difference));}}}//确保返回的几何图形都为三角形if (triangle.getCoordinates().length == 4) {resultTriangles.add(triangle);} else {//对于非三角形，还需要进行一次原生工具的三角剖分,剖分成多个三角形Geometry geometry = genTriangle(Lists.newArrayList(), Lists.newArrayList(triangle));for (int i = 0; i < geometry.getNumGeometries(); i++) {resultTriangles.add((Polygon) geometry.getGeometryN(i));}}//下方代码是判断是否为凸多边形，返回的结果集集合图形不一定是三角形，但不会是凹多边形// navmesh寻路只要不是凹多边形就可以使用
//        if (isConvex(triangle)) {
//            resultTriangles.add(triangle);
//        } else {
//            //如果是被障碍物多边形截断后的形状，有可能是凹多边形，对于凹多边形，还需要进行一次原生工具的三角剖分
//            Geometry geometry = genTriangle(Lists.newArrayList(), Lists.newArrayList(triangle));
//            for (int i = 0; i < geometry.getNumGeometries(); i++) {
//                resultTriangles.add((Polygon) geometry.getGeometryN(i));
//            }
//        }return resultTriangles;}/*** 多边形是否为凸多边形* 本方法代码是chatgpt生成的，目前测试好像可用** @param polygon 多边形* @return*/public static boolean isConvex(Polygon polygon) {Coordinate[] coords = polygon.getExteriorRing().getCoordinates();if (coords.length < 4) {//之所以是四个是因为第一个点和第四个点是一个样的，所以四个点则为三角形return false; // 少于4个顶点的多边形不可能是凸多边形}boolean direction = getTurnDirection(coords[0], coords[1], coords[2]);for (int i = 1; i < coords.length - 2; i++) {boolean currentDirection = getTurnDirection(coords[i], coords[i + 1], coords[i + 2]);if (currentDirection != direction) {return false;}}return true;}/*** chatgpt 生成代码** @param a* @param b* @param c* @return*/private static boolean getTurnDirection(Coordinate a, Coordinate b, Coordinate c) {double crossProduct = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);return crossProduct > 0; // 返回true表示逆时针转，false表示顺时针转}/*** 将 MultiPolygon 转成 Polygon 集合** @param multiPolygons* @return*/public static List<Polygon> multiPolygonToPolygon(List<MultiPolygon> multiPolygons) {List<Polygon> polygonList = Lists.newArrayList();for (MultiPolygon polygon : multiPolygons) {for (int i = 0; i < polygon.getNumGeometries(); i++) {polygonList.add((Polygon) polygon.getGeometryN(i));}}return polygonList;}/*** 使用jts 原生工具进行三角剖分** @param bounds   外边框的四个坐标点(地图边界范围),可以取所有多边形的最大外边框，也可以自己设置* @param polygons 所有障碍物多边形* @return*/public static Geometry genTriangle(List<Coordinate> bounds, List<Geometry> polygons) {//获取外边框及所有障碍物多边形点位的集合Set<Coordinate> coordinates = Sets.newHashSet();for (Geometry polygon : polygons) {coordinates.addAll(Arrays.asList(polygon.getCoordinates()));}coordinates.addAll(bounds);//对所有点调用 jts 原生工具进行三角剖分DelaunayTriangulationBuilder builder = new DelaunayTriangulationBuilder();builder.setSites(coordinates);return builder.getTriangles(GEOMETRY_FACTORY);}/*** 获取外边框的四个坐标点** @param bounds* @return*/public static List<Coordinate> boundToCoordinate(ReferencedEnvelope bounds) {List<Coordinate> coordinates = Lists.newArrayList();coordinates.add(new Coordinate(bounds.getMinX(), bounds.getMinY()));coordinates.add(new Coordinate(bounds.getMinX(), bounds.getMaxY()));coordinates.add(new Coordinate(bounds.getMaxX(), bounds.getMaxY()));coordinates.add(new Coordinate(bounds.getMaxX(), bounds.getMinY()));return coordinates;}/*** 加载shp文件** @param shpPath shp文件地址* @return* @throws Exception*/public static SimpleFeatureSource getFeatureSourceByPath(String shpPath) throws Exception {ShapefileDataStore store = new ShapefileDataStore(new URL("file:/" + shpPath));String typeName = store.getTypeNames()[0];return store.getFeatureSource(typeName);}/*** 获取shp文件内的所有多边形** @param featureSource* @return* @throws IOException*/private static List<MultiPolygon> getPolygons(SimpleFeatureSource featureSource) throws IOException {List<MultiPolygon> polygons = Lists.newArrayList();SimpleFeatureIterator features = featureSource.getFeatures().features();while (features.hasNext()) {SimpleFeature feature = features.next();Geometry geometry = (Geometry) feature.getDefaultGeometry();if (geometry instanceof MultiPolygon) {MultiPolygon multiPolygon = (MultiPolygon) geometry;polygons.add(multiPolygon);}}return polygons;}/*** 取反* 比如如果一个shp 含有多个多边形障碍物，可以用本方法生成一个多边形含有多个洞,洞表示原来的障碍物多边形** @param featureSource* @param coordinates   外边框(地图边界范围)* @return 取反后的几何图形* @throws IOException*/public static Geometry negationShp(SimpleFeatureSource featureSource, List<Coordinate> coordinates) throws IOException {if (CollectionUtil.isEmpty(coordinates)) {//如果没有自定义外边框，则取shp文件的外接矩形作为外边框ReferencedEnvelope bounds = featureSource.getBounds();coordinates = boundToCoordinate(bounds);}//考虑是否判断最后一个点和第一个点一样时不添加新点coordinates.add(coordinates.get(0));Geometry geometryNegation = GEOMETRY_FACTORY.createPolygon(coordinates.toArray(new Coordinate[0]));SimpleFeatureIterator features = featureSource.getFeatures().features();while (features.hasNext()) {SimpleFeature feature = features.next();Geometry geometry = (Geometry) feature.getDefaultGeometry();geometryNegation = geometryNegation.difference(geometry);}return geometryNegation;}/*** 将几何图形写入shp文件** @param geometry  几何图形(可以是一个集合GeometryCollection)* @param writePath 写入的路径位置* @throws Exception*/public static void writeMultiPolygonShp(Geometry geometry, String writePath) throws Exception {ShapefileDataStore ds = new ShapefileDataStore(new URL("file:/" + writePath));//设置编码ds.setCharset(Charset.forName("UTF-8"));//定义图形信息和属性信息SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();tb.setCRS(DefaultGeographicCRS.WGS84);tb.setName("shapefile");tb.add("the_geom", MultiPolygon.class);tb.add("num", Integer.class);ds.createSchema(tb.buildFeatureType());//设置WriterFeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);//写入文件信息for (int i = 0; i < geometry.getNumGeometries(); i++) {SimpleFeature feature = writer.next();Geometry triangle = geometry.getGeometryN(i);feature.setAttribute("the_geom", triangle);feature.setAttribute("num", i);}writer.write();writer.close();ds.dispose();}
}

效果

上面代码切割的效果如下图所示

上面代码有个取反后再进行三角剖分的代码，效果如下

与原生工具对比如下,不会切割到障碍物内部了。

扩展

求共边

上面生成的三角形数量还是比较多的, navmesh 不要求一定要是三角形，可以是凸多边形，不可以是凹多边形，所以可以基于上面的三角形集合根据共边合并成凸多边形，从而减少寻路的三角形数量。从而提升寻路性能。

耳切法

网上是有说耳切法进行三角剖分是可以切割有洞的多边形的，有兴趣的可以自己试试

寻路

GitHub - jzyong/GameAI4j: Game AI for java.NavMesh、A*、BehaviorTree、FSM

上面地址提供了如何基于 navmesh 寻路的相关代码

核心类

PolygonMeshWindow 多边形寻路的可视化 ui 界面启动类, loadMap 方法会加载 navmesh 文件初始化PolygonNavMesh
PolygonNavMesh 多边形寻路的类,将多边形转为改结构即可用 findPath 进行寻路

TriangleNavMeshWindow 三角形寻路的可视化 ui 界面启动类
TriangleNavMesh 用于三角形寻路的类

三角形寻路没用过，可以考虑直接用多边形寻路。

navmesh文件结构介绍

只需要将本文生成的三角形转换成上面库可以解析的navmesh文件，即可用上面的库进行寻路

navmesh的结构如下图所示

主要就是介绍下 pathTriangles 和 pathVertices

pathVertices s是所有多边形的点集数组,如下图，表示的是多边形所有顶点。

pathTriangles 多边形的所有标识，存的值表示的是 pathVertices 中的数组下标对应的点。具有以下特性

1.每三个点表示一个三角形，比如第1-3个值 0,1,2 表示 pathVertices[0]、pathVertices[1]、pathVertices[2] 组成一个三角形，第4-6个指 0、2、3 又是一个三角形。
2.当连续的三角形有两个下标是一样的，表示是共边三角形，可以组成一个凸多边形。所以生成navmesh的时候需要注意不然让连续的三角形有两个下标一样的却不能合成凸多边形。比如前面 0、1、2 和0、2、3有相同的点0、2,所以就可以组成一个0、1、2、3的四边形，以此类推可以合并成五边形、六边形。

另一种思路

上面的方法是将 shp 转换的三角形或者多边形转成 navmesh 文件，但是这样得按照navmesh的格式规则，其实可以自己去看 P olygonNavMesh 类的代码，直接将三角形和多边形转成 PolygonNavMesh 需要的数据结构。这个也不会很难，无非就修改一些 navmesh的解析逻辑