ConcaveHullOfPolygons.java

/*
 * Copyright (c) 2022 Martin Davis.
 *
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License 2.0
 * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
 * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
 * and the Eclipse Distribution License is available at
 *
 * http://www.eclipse.org/org/documents/edl-v10.php.
 */
package org.locationtech.jts.algorithm.hull;

import java.util.ArrayDeque;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.MultiPolygon;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.operation.overlayng.CoverageUnion;
import org.locationtech.jts.triangulate.polygon.ConstrainedDelaunayTriangulator;
import org.locationtech.jts.triangulate.tri.Tri;

/**
 * Constructs a concave hull of a set of polygons, respecting 
 * the polygons as constraints.
 * A concave hull is a concave or convex polygon containing all the input polygons,
 * whose vertices are a subset of the vertices in the input.
 * A given set of polygons has a sequence of hulls of increasing concaveness,
 * determined by a numeric target parameter.
 * The computed hull "fills the gap" between the polygons,
 * and does not intersect their interior.
 * <p>
 * The concave hull is constructed by removing the longest outer edges 
 * of the constrained Delaunay Triangulation of the space between the polygons,
 * until the target criterion parameter is reached.
 * <p>
 * The target criteria are:
 * <ul>
 * <li><b>Maximum Edge Length</b> - the length of the longest edge between the polygons is no larger
 * than this value.
 * <li><b>Maximum Edge Length Ratio</b> - determine the Maximum Edge Length 
 * as a fraction of the difference between the longest and shortest edge lengths 
 * between the polygons.  
 * This provides a scale-free parameter.
 * A value of 1 produces the convex hull; a value of 0 produces the original polygons.
 * </ul>
 * Optionally the concave hull can be allowed to contain holes, 
 * via {@link #setHolesAllowed(boolean)}.
 * <p>
 * The hull can be specified as being "tight", via {@link #setTight(boolean)}.
 * This causes the result to follow the outer boundaries of the input polygons. 
 * <p>
 * Instead of the complete hull, the "fill area" between the input polygons 
 * can be computed using {@link #getFill()}.
 * <p>
 * The input polygons must form a valid {@link MultiPolygon}
 * (i.e. they must be non-overlapping and non-edge-adjacent).
 * If needed, a set of possibly-overlapping Polygons 
 * can be converted to a valid MultiPolygon
 * by using {@link Geometry#union()};
 * 
 * @author Martin Davis
 *
 */
public class ConcaveHullOfPolygons {
 
  /**
   * Computes a concave hull of set of polygons
   * using the target criterion of maximum edge length.
   * 
   * @param polygons the input polygons
   * @param maxLength the target maximum edge length
   * @return the concave hull
   */
  public static Geometry concaveHullByLength(Geometry polygons, double maxLength) {
    return concaveHullByLength(polygons, maxLength, false, false);
  }
  
  /**
   * Computes a concave hull of set of polygons
   * using the target criterion of maximum edge length,
   * and allowing control over whether the hull boundary is tight
   * and can contain holes.
   * 
   * @param polygons the input polygons
   * @param maxLength the target maximum edge length
   * @param isTight true if the hull should be tight to the outside of the polygons
   * @param isHolesAllowed true if holes are allowed in the hull polygon
   * @return the concave hull
   */
  public static Geometry concaveHullByLength(Geometry polygons, double maxLength,
      boolean isTight, boolean isHolesAllowed) {
    ConcaveHullOfPolygons hull = new ConcaveHullOfPolygons(polygons);
    hull.setMaximumEdgeLength(maxLength);
    hull.setHolesAllowed(isHolesAllowed);
    hull.setTight(isTight);
    return hull.getHull();
  }
  
  /**
   * Computes a concave hull of set of polygons
   * using the target criterion of maximum edge length ratio.
   * 
   * @param polygons the input polygons
   * @param lengthRatio the target maximum edge length ratio
   * @return the concave hull
   */
  public static Geometry concaveHullByLengthRatio(Geometry polygons, double lengthRatio) {
    return concaveHullByLengthRatio(polygons, lengthRatio, false, false);
  }
  
  /**
   * Computes a concave hull of set of polygons
   * using the target criterion of maximum edge length ratio,
   * and allowing control over whether the hull boundary is tight
   * and can contain holes.
   * 
   * @param polygons the input polygons
   * @param lengthRatio the target maximum edge length ratio
   * @param isTight true if the hull should be tight to the outside of the polygons
   * @param isHolesAllowed true if holes are allowed in the hull polygon
   * @return the concave hull
   */
  public static Geometry concaveHullByLengthRatio(Geometry polygons, double lengthRatio,
      boolean isTight, boolean isHolesAllowed) {
    ConcaveHullOfPolygons hull = new ConcaveHullOfPolygons(polygons);
    hull.setMaximumEdgeLengthRatio(lengthRatio);
    hull.setHolesAllowed(isHolesAllowed);
    hull.setTight(isTight);
    return hull.getHull();
  }
  
  /**
   * Computes a concave fill area between a set of polygons,
   * using the target criterion of maximum edge length.
   * 
   * @param polygons the input polygons
   * @param maxLength the target maximum edge length
   * @return the concave fill
   */
  public static Geometry concaveFillByLength(Geometry polygons, double maxLength) {
    ConcaveHullOfPolygons hull = new ConcaveHullOfPolygons(polygons);
    hull.setMaximumEdgeLength(maxLength);
    return hull.getFill();
  }
  
  /**
   * Computes a concave fill area between a set of polygons,
   * using the target criterion of maximum edge length ratio.
   * 
   * @param polygons the input polygons
   * @param lengthRatio the target maximum edge length ratio
   * @return the concave fill
   */
  public static Geometry concaveFillByLengthRatio(Geometry polygons, double lengthRatio) {
    ConcaveHullOfPolygons hull = new ConcaveHullOfPolygons(polygons);
    hull.setMaximumEdgeLengthRatio(lengthRatio);
    return hull.getFill();
  }
  
  private static final int FRAME_EXPAND_FACTOR = 4;
  private static final int NOT_SPECIFIED = -1;
  private static final int NOT_FOUND = -1;

  private Geometry inputPolygons;
  private double maxEdgeLength = 0.0;
  private double maxEdgeLengthRatio = NOT_SPECIFIED;
  private boolean isHolesAllowed = false;
  private boolean isTight = false;
  
  private GeometryFactory geomFactory;
  private LinearRing[] polygonRings;
  
  private Set<Tri> hullTris;
  private ArrayDeque<Tri> borderTriQue;
  /**
   * Records the edge index of the longest border edge for border tris,
   * so it can be tested for length and possible removal.
   */
  private Map<Tri, Integer> borderEdgeMap = new HashMap<Tri, Integer>();
  
  /**
   * Creates a new instance for a given geometry.
   * 
   * @param geom the input geometry
   */
  public ConcaveHullOfPolygons(Geometry polygons) {
    if (! (polygons instanceof Polygon || polygons instanceof MultiPolygon)) {
      throw new IllegalArgumentException("Input must be polygonal");
    }
    this.inputPolygons = polygons;
    geomFactory = inputPolygons.getFactory();
  }

  /**
   * Sets the target maximum edge length for the concave hull.
   * The length value must be zero or greater.
   * <ul>
   * <li>The value 0.0 produces the input polygons.
   * <li>Larger values produce less concave results.
   * Above a certain large value the result is the convex hull of the input.
   * <p>
   * The edge length ratio provides a scale-free parameter which
   * is intended to produce similar concave results for a variety of inputs.
   * 
   * @param edgeLength a non-negative length
   */
  public void setMaximumEdgeLength(double edgeLength) {
    if (edgeLength < 0)
      throw new IllegalArgumentException("Edge length must be non-negative");
    this.maxEdgeLength = edgeLength;
    maxEdgeLengthRatio = NOT_SPECIFIED;
  }
  
  /**
   * Sets the target maximum edge length ratio for the concave hull.
   * The edge length ratio is a fraction of the difference
   * between the longest and shortest edge lengths 
   * in the Delaunay Triangulation of the area between the input polygons.
   * (Roughly speaking, it is a fraction of the difference between
   * the shortest and longest distances between the input polygons.)
   * It is a value in the range 0 to 1. 
   * <ul>
   * <li>The value 0.0 produces the original input polygons.
   * <li>The value 1.0 produces the convex hull.
   * <ul> 
   * 
   * @param edgeLengthRatio a length factor value between 0 and 1
   */
  public void setMaximumEdgeLengthRatio(double edgeLengthRatio) {
    if (edgeLengthRatio < 0 || edgeLengthRatio > 1)
      throw new IllegalArgumentException("Edge length ratio must be in range [0,1]");
    this.maxEdgeLengthRatio = edgeLengthRatio;
  }
  
  /**
   * Sets whether holes are allowed in the concave hull polygon.
   * 
   * @param isHolesAllowed true if holes are allowed in the result
   */
  public void setHolesAllowed(boolean isHolesAllowed) {
    this.isHolesAllowed = isHolesAllowed;
  }
  
  /**
   * Sets whether the boundary of the hull polygon is kept
   * tight to the outer edges of the input polygons.
   * 
   * @param isTight true if the boundary is kept tight
   */
  public void setTight(boolean isTight) {
    this.isTight = isTight;
  }
  
  /**
   * Gets the computed concave hull.
   * 
   * @return the concave hull
   */
  public Geometry getHull() {
    if (inputPolygons.isEmpty()) {
      return createEmptyHull();
    }
    buildHullTris();
    Geometry hull = createHullGeometry(hullTris, true);
    return hull;
  }

  /**
   * Gets the concave fill, which is the area between the input polygons, 
   * subject to the concaveness control parameter.
   * 
   * @return the concave fill
   */
  public Geometry getFill() {
    isTight = true;
    if (inputPolygons.isEmpty()) {
      return createEmptyHull();
    }
    buildHullTris();
    Geometry fill = createHullGeometry(hullTris, false);
    return fill;
  }
  
  private Geometry createEmptyHull() {
    return geomFactory.createPolygon();
  }
  
  private void buildHullTris() {
    polygonRings = extractShellRings(inputPolygons);
    Polygon frame = createFrame(inputPolygons.getEnvelopeInternal(), polygonRings, geomFactory);
    ConstrainedDelaunayTriangulator cdt = new ConstrainedDelaunayTriangulator(frame);
    List<Tri> tris = cdt.getTriangles();
    //System.out.println(tris);
    
    Coordinate[] framePts = frame.getExteriorRing().getCoordinates();
    if (maxEdgeLengthRatio >= 0) {
      maxEdgeLength = computeTargetEdgeLength(tris, framePts, maxEdgeLengthRatio);
    }
    
    hullTris = removeFrameCornerTris(tris, framePts);
    
    removeBorderTris();
    if (isHolesAllowed) removeHoleTris();
  }

  private static double computeTargetEdgeLength(List<Tri> triList, 
      Coordinate[] frameCorners,
      double edgeLengthRatio) {
    if (edgeLengthRatio == 0) return 0;
    double maxEdgeLen = -1;
    double minEdgeLen = -1;
    for (Tri tri : triList) {
      //-- don't include frame triangles
      if (isFrameTri(tri, frameCorners))
        continue;
      
      for (int i = 0; i < 3; i++) {
        //-- constraint edges are not used to determine ratio
        if (! tri.hasAdjacent(i))
          continue;
        
        double len = tri.getLength(i);
        if (len > maxEdgeLen) 
          maxEdgeLen = len;
        if (minEdgeLen < 0 || len < minEdgeLen)
          minEdgeLen = len;
      }
    }
    //-- if ratio = 1 ensure all edges are included
    if (edgeLengthRatio == 1) 
      return 2 * maxEdgeLen;
    
    return edgeLengthRatio * (maxEdgeLen - minEdgeLen) + minEdgeLen;
  }

  private static boolean isFrameTri(Tri tri, Coordinate[] frameCorners) {
    int index = vertexIndex(tri, frameCorners);
    boolean isFrameTri = index >= 0;
    return isFrameTri;
  }
  
  private Set<Tri> removeFrameCornerTris(List<Tri> tris, Coordinate[] frameCorners) {
    Set<Tri> hullTris = new HashSet<Tri>();
    borderTriQue = new ArrayDeque<Tri>();
    for (Tri tri : tris) {
      int index = vertexIndex(tri, frameCorners);
      boolean isFrameTri = index != NOT_FOUND;
      if (isFrameTri) {
        /**
         * Frame tris are adjacent to at most one border tri,
         * which is opposite the frame corner vertex.
         * Or, the opposite tri may be another frame tri,
         * which is not added as a border tri.
         */
        int oppIndex = Tri.oppEdge(index);
        Tri oppTri = tri.getAdjacent(oppIndex);
        boolean isBorderTri = oppTri != null && ! isFrameTri(oppTri, frameCorners);
        if (isBorderTri) {
          addBorderTri(tri, oppIndex);
        }
        //-- remove the frame tri
        tri.remove();
      }
      else {
        hullTris.add(tri);
        //System.out.println(tri);
      }
    }
    return hullTris;
  }

  /**
   * Get the tri vertex index of some point in a list, 
   * or -1 if none are vertices.
   * 
   * @param tri the tri to test for containing a point
   * @param pts the points to test
   * @return the vertex index of a point, or -1
   */
  private static int vertexIndex(Tri tri, Coordinate[] pts) {
    for (Coordinate p : pts) {
      int index = tri.getIndex(p);
      if (index >= 0) 
        return index;
    }
    return NOT_FOUND;
  }
  
  private void removeBorderTris() {
    while (! borderTriQue.isEmpty()) {
      Tri tri = borderTriQue.pop();
      //-- tri might have been removed already
      if (! hullTris.contains(tri)) {
        continue;
      }
      if (isRemovable(tri)) {
        addBorderTris(tri);
        removeBorderTri(tri);
        //System.out.println(tri);
      }
    }
  }

  private void removeHoleTris() {
    while (true) {
      Tri holeTri = findHoleSeedTri(hullTris);
      if (holeTri == null)
        return;
      addBorderTris(holeTri);
      removeBorderTri(holeTri);
      removeBorderTris();
    }
  }
  
  private Tri findHoleSeedTri(Set<Tri> tris) {
    for (Tri tri : tris) {
      if (isHoleSeedTri(tri))
        return tri;
    }
    return null;
  }

  private boolean isHoleSeedTri(Tri tri) {
    if (isBorderTri(tri))
      return false;
    for (int i = 0; i < 3; i++) {
      if (tri.hasAdjacent(i)
          && tri.getLength(i) > maxEdgeLength)
         return true;
    }
    return false;
  }

  private boolean isBorderTri(Tri tri) {
    for (int i = 0; i < 3; i++) {
      if (! tri.hasAdjacent(i))
          return true;
    }
    return false;
  }

  private boolean isRemovable(Tri tri) {
    //-- remove non-bridging tris if keeping hull boundary tight
    if (isTight && isTouchingSinglePolygon(tri))
      return true;
    
    //-- check if outside edge is longer than threshold
    if (borderEdgeMap.containsKey(tri)) {
      int borderEdgeIndex = borderEdgeMap.get(tri);
      double edgeLen = tri.getLength(borderEdgeIndex);
      if (edgeLen > maxEdgeLength)
        return true;
    }
    return false;
  }

  /**
   * Tests whether a triangle touches a single polygon at all vertices.
   * If so, it is a candidate for removal if the hull polygon
   * is being kept tight to the outer boundary of the input polygons.
   * Tris which touch more than one polygon are called "bridging".
   * 
   * @param tri
   * @return true if the tri touches a single polygon
   */
  private boolean isTouchingSinglePolygon(Tri tri) {
    Envelope envTri = envelope(tri);
    for (LinearRing ring : polygonRings) {
      //-- optimization heuristic: a touching tri must be in ring envelope
      if (ring.getEnvelopeInternal().intersects(envTri)) {
        if (hasAllVertices(ring, tri))
          return true;
      }
    }
    return false;
  }

  private void addBorderTris(Tri tri) {
    addBorderTri(tri, 0);
    addBorderTri(tri, 1);
    addBorderTri(tri, 2);
  }
  
  /**
   * Adds an adjacent tri to the current border.
   * The adjacent edge is recorded as the border edge for the tri.
   * Note that only edges adjacent to another tri can become border edges.
   * Since constraint-adjacent edges do not have an adjacent tri,
   * they can never be on the border and thus will not be removed
   * due to being shorter than the length threshold.
   * The tri containing them may still be removed via another edge, however. 
   * 
   * @param tri the tri adjacent to the tri to be added to the border
   * @param index the index of the adjacent tri
   */
  private void addBorderTri(Tri tri, int index) {
    Tri adj = tri.getAdjacent( index );
    if (adj == null) 
      return;
    borderTriQue.add(adj);
    int borderEdgeIndex = adj.getIndex(tri);
    borderEdgeMap.put(adj, borderEdgeIndex);
  }

  private void removeBorderTri(Tri tri) {
    tri.remove();
    hullTris.remove(tri);
    borderEdgeMap.remove(tri);
  }
  
  private static boolean hasAllVertices(LinearRing ring, Tri tri) {
    for (int i = 0; i < 3; i++) {
      Coordinate v = tri.getCoordinate(i);
      if (! hasVertex(ring, v)) {
        return false;
      }
    }
    return true;
  }
  
  private static boolean hasVertex(LinearRing ring, Coordinate v) {
    for(int i = 1; i < ring.getNumPoints(); i++) {
      if (v.equals2D(ring.getCoordinateN(i))) {
        return true;
      }
    }
    return false;
  }

  private static Envelope envelope(Tri tri) {
    Envelope env = new Envelope(tri.getCoordinate(0), tri.getCoordinate(1));
    env.expandToInclude(tri.getCoordinate(2));
    return env;
  }

  private Geometry createHullGeometry(Set<Tri> hullTris, boolean isIncludeInput) {
    if (! isIncludeInput && hullTris.isEmpty())
      return createEmptyHull();
    
    //-- union triangulation
    Geometry triCoverage = Tri.toGeometry(hullTris, geomFactory);
    //System.out.println(triCoverage);
    Geometry fillGeometry = CoverageUnion.union(triCoverage);
    
    if (! isIncludeInput) {
      return fillGeometry;
    }
    if (fillGeometry.isEmpty()) {
      return inputPolygons.copy();
    }
    //-- union with input polygons
    Geometry[] geoms = new Geometry[] { fillGeometry, inputPolygons };
    GeometryCollection geomColl = geomFactory.createGeometryCollection(geoms);
    Geometry hull = CoverageUnion.union(geomColl);
    return hull;
  }
  
  /**
   * Creates a rectangular "frame" around the input polygons,
   * with the input polygons as holes in it.
   * The frame is large enough that the constrained Delaunay triangulation
   * of it should contain the convex hull of the input as edges.
   * The frame corner triangles can be removed to produce a 
   * triangulation of the space around and between the input polygons.
   * 
   * @param polygonsEnv
   * @param polygonRings
   * @param geomFactory 
   * @return the frame polygon
   */
  private static Polygon createFrame(Envelope polygonsEnv, LinearRing[] polygonRings, GeometryFactory geomFactory) {
    double diam = polygonsEnv.getDiameter();
    Envelope envFrame = polygonsEnv.copy();
    envFrame.expandBy(FRAME_EXPAND_FACTOR * diam);
    Polygon frameOuter = (Polygon) geomFactory.toGeometry(envFrame);
    LinearRing shell = (LinearRing) frameOuter.getExteriorRing().copy();
    Polygon frame = geomFactory.createPolygon(shell, polygonRings);
    return frame;
  }

  private static LinearRing[] extractShellRings(Geometry polygons) {
    LinearRing[] rings = new LinearRing[polygons.getNumGeometries()];
    for (int i = 0; i < polygons.getNumGeometries(); i++) {
      Polygon consPoly = (Polygon) polygons.getGeometryN(i);
      rings[i] = (LinearRing) consPoly.getExteriorRing().copy();
    }
    return rings;
  }
}