ConcaveHull.java

/*
 * Copyright (c) 2021 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.ArrayList;
import java.util.List;
import java.util.PriorityQueue;

import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;

/**
 * Constructs a concave hull of a set of points.
 * A concave hull is a concave or convex polygon containing all the input points,
 * whose vertices are a subset of the vertices in the input.
 * A given set of points has a sequence of hulls of increasing concaveness,
 * determined by a numeric target parameter.
 * <p>
 * The hull is constructed by removing border triangles 
 * of the Delaunay Triangulation of the points,
 * as long as their "size" is larger than the target criterion.
 * <p>
 * The target criteria are:
 * <ul>
 * <li><b>Maximum Edge Length</b> - the length of the longest edge of the hull is no larger
 * than this value.
 * <li><b>Maximum Edge Length Ratio</b> - determines the Maximum Edge Length 
 * by a fraction of the difference between the longest and shortest edge lengths 
 * in the Delaunay Triangulation.  
 * This normalizes the <b>Maximum Edge Length</b> to be scale-free.
 * A value of 1 produces the convex hull; a value of 0 produces maximum concaveness.
 * <li><b>Alpha</b> - produces Alpha-shapes, 
 * by removing border triangles with a circumradius greater than alpha. 
 * Large values produce the convex hull; a value of 0 produces maximum concaveness. 
 * </ul>
 * The preferred criterion is the <b>Maximum Edge Length Ratio</b>, since it is 
 * scale-free and local (so that no assumption needs to be made about the 
 * total amount of concaveness present).
 * <p>
 * Other length criteria can be used by setting the Maximum Edge Length directly.
 * For example, use a length relative  to the longest edge length
 * in the Minimum Spanning Tree of the point set.
 * Or, use a length derived from the {@link #uniformGridEdgeLength(Geometry)} value.
 * <p>
 * The computed hull is always a single connected {@link Polygon}
 * (unless it is degenerate, in which case it will be a {@link Point} or a {@link LineString}).
 * This constraint may cause the concave hull to fail to meet the target criterion.
 * <p>
 * Optionally the concave hull can be allowed to contain holes by calling {@link #setHolesAllowed(boolean)}.
 * 
 * @author Martin Davis
 *
 */
public class ConcaveHull 
{
  /**
   * Computes the approximate edge length of
   * a uniform square grid having the same number of
   * points as a geometry and the same area as its convex hull.
   * This value can be used to determine a suitable length threshold value
   * for computing a concave hull.  
   * A value from 2 to 4 times the uniform grid length 
   * seems to produce reasonable results.
   *  
   * @param geom a geometry
   * @return the approximate uniform grid length
   */
  public static double uniformGridEdgeLength(Geometry geom) {
    double areaCH = geom.convexHull().getArea();
    int numPts = geom.getNumPoints();
    return Math.sqrt(areaCH / numPts);
  }
  
  /**
   * Computes a concave hull of the vertices in a geometry
   * using the target criterion of maximum edge length.
   * 
   * @param geom the input geometry
   * @param maxLength the target maximum edge length
   * @return the concave hull
   */
  public static Geometry concaveHullByLength(Geometry geom, double maxLength) {
    return concaveHullByLength(geom, maxLength, false);
  }
  
  /**
   * Computes a concave hull of the vertices in a geometry
   * using the target criterion of maximum edge length,
   * and optionally allowing holes.
   * 
   * @param geom the input geometry
   * @param maxLength the target maximum edge length
   * @param isHolesAllowed whether holes are allowed in the result
   * @return the concave hull
   */
  public static Geometry concaveHullByLength(Geometry geom, double maxLength, boolean isHolesAllowed) {
    ConcaveHull hull = new ConcaveHull(geom);
    hull.setMaximumEdgeLength(maxLength);
    hull.setHolesAllowed(isHolesAllowed);
    return hull.getHull();
  }
  
  /**
   * Computes a concave hull of the vertices in a geometry
   * using the target criterion of maximum edge length ratio.
   * The edge length ratio is a fraction of the length difference
   * between the longest and shortest edges 
   * in the Delaunay Triangulation of the input points. 
   * 
   * @param geom the input geometry
   * @param lengthRatio the target edge length factor
   * @return the concave hull
   */
  public static Geometry concaveHullByLengthRatio(Geometry geom, double lengthRatio) {
    return concaveHullByLengthRatio(geom, lengthRatio, false);
  }
  
  /**
   * Computes a concave hull of the vertices in a geometry
   * using the target criterion of maximum edge length factor,
   * and optionally allowing holes.
   * The edge length factor is a fraction of the length difference
   * between the longest and shortest edges 
   * in the Delaunay Triangulation of the input points. 
   * 
   * @param geom the input geometry
   * @param maxLength the target maximum edge length
   * @param isHolesAllowed whether holes are allowed in the result
   * @return the concave hull
   */
  public static Geometry concaveHullByLengthRatio(Geometry geom, double lengthRatio, boolean isHolesAllowed) {
    ConcaveHull hull = new ConcaveHull(geom);
    hull.setMaximumEdgeLengthRatio(lengthRatio);
    hull.setHolesAllowed(isHolesAllowed);
    return hull.getHull();
  }
  
  /**
   * Computes the alpha shape of a geometry as a polygon.
   * The alpha parameter is the radius of the eroding disc.
   * 
   * @param geom the input geometry
   * @param alpha the radius of the eroding disc
   * @param isHolesAllowed whether holes are allowed in the result
   * @return the alpha shape polygon
   */
  public static Geometry alphaShape(Geometry geom, double alpha, boolean isHolesAllowed) {
    ConcaveHull hull = new ConcaveHull(geom);
    hull.setAlpha(alpha);
    hull.setHolesAllowed(isHolesAllowed);
    return hull.getHull();
  }
  
  private static int PARAM_EDGE_LENGTH = 1;
  private static int PARAM_ALPHA = 2;
  
  private Geometry inputGeometry;
  private double maxEdgeLengthRatio = -1;
  private double alpha = -1;
  private boolean isHolesAllowed = false;
  
  private int criteriaType = PARAM_EDGE_LENGTH;
  private double maxSizeInHull = 0.0;
  private GeometryFactory geomFactory;


  /**
   * Creates a new instance for a given geometry.
   * 
   * @param geom the input geometry
   */
  public ConcaveHull(Geometry geom) {
    this.inputGeometry = geom;
    this.geomFactory = geom.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 concave hull of smallest area
   * that is still connected.
   * <li>Larger values produce less concave results.
   * A value equal or greater than the longest Delaunay Triangulation edge length
   * produces the convex hull.
   * </ul>
   * The {@link #uniformGridEdgeLength(Geometry)} value may be used as
   * the basis for estimating an appropriate target maximum edge length.
   * 
   * @param edgeLength a non-negative length
   * 
   * @see #uniformGridEdgeLength(Geometry)
   */
  public void setMaximumEdgeLength(double edgeLength) {
    if (edgeLength < 0)
      throw new IllegalArgumentException("Edge length must be non-negative");
    this.maxSizeInHull = edgeLength;
    maxEdgeLengthRatio = -1;
    criteriaType = PARAM_EDGE_LENGTH;
  }
  
  /**
   * 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 input points.
   * It is a value in the range 0 to 1. 
   * <ul>
   * <li>The value 0.0 produces a concave hull of minimum area
   * that is still connected.
   * <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;
    criteriaType = PARAM_EDGE_LENGTH;
  }
  
  /**
   * Sets the alpha parameter to compute an alpha shape of the input.
   * Alpha is the radius of the eroding disc.
   * Border triangles with circumradius greater than alpha are removed.
   * 
   * @param alpha the alpha radius
   */
  public void setAlpha(double alpha) {
    this.alpha = alpha;
    maxSizeInHull = alpha;
    criteriaType = PARAM_ALPHA;
  }
  
  /**
   * 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;
  }
  
  /**
   * Gets the computed concave hull.
   * 
   * @return the concave hull
   */
  public Geometry getHull() {
    if (inputGeometry.isEmpty()) {
      return geomFactory.createPolygon();
    }
    List<HullTri> triList = HullTriangulation.createDelaunayTriangulation(inputGeometry);
    setSize(triList);
    
    if (maxEdgeLengthRatio >= 0) {
      maxSizeInHull = computeTargetEdgeLength(triList, maxEdgeLengthRatio);
    }
    if (triList.isEmpty())
      return inputGeometry.convexHull();
    
    computeHull(triList);    

    Geometry hull = toGeometry(triList, geomFactory);
    return hull;
  }

  private void setSize(List<HullTri> triList) {
    for (HullTri tri : triList) {
      if (criteriaType == PARAM_EDGE_LENGTH) {
        tri.setSizeToLongestEdge();
      }
      else {
        tri.setSizeToCircumradius();
      }
    }
  }

  private static double computeTargetEdgeLength(List<HullTri> triList, 
      double edgeLengthRatio) {
    if (edgeLengthRatio == 0) return 0;
    double maxEdgeLen = -1;
    double minEdgeLen = -1;
    for (HullTri tri : triList) {
      for (int i = 0; i < 3; i++) {
        double len = tri.getCoordinate(i).distance(tri.getCoordinate(HullTri.next(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;
  }
  
  /**
   * Computes the concave hull using edge length as the target criterion.
   * The erosion is done in two phases: first the border, then any
   * internal holes (if required).
   * This allows an fast connection check to be used
   * when eroding holes,
   * which makes this much more efficient than the area-based algorithm.
   * 
   * @param triList
   */
  private void computeHull(List<HullTri> triList) {
    computeHullBorder(triList);
    if (isHolesAllowed) {
      computeHullHoles(triList);
    }
  }
  
  private void computeHullBorder(List<HullTri> triList) {
    PriorityQueue<HullTri> queue = createBorderQueue(triList);
    // process tris in order of decreasing size (edge length or circumradius)
    while (! queue.isEmpty()) {
      HullTri tri = queue.poll();
      
      if (isInHull(tri)) 
        break;
      
      if (isRemovableBorder(tri)) {
        //-- the non-null adjacents are now on the border
        HullTri adj0 = (HullTri) tri.getAdjacent(0);
        HullTri adj1 = (HullTri) tri.getAdjacent(1);
        HullTri adj2 = (HullTri) tri.getAdjacent(2);
        
        tri.remove(triList);
        
        //-- add border adjacents to queue
        addBorderTri(adj0, queue);
        addBorderTri(adj1, queue);
        addBorderTri(adj2, queue);
      }
    }
  }
  
  private PriorityQueue<HullTri> createBorderQueue(List<HullTri> triList) {
    PriorityQueue<HullTri> queue = new PriorityQueue<HullTri>();
    for (HullTri tri : triList) {
      addBorderTri(tri, queue);
    }
    return queue;
  }

  /**
   * Adds a Tri to the queue.
   * Only add tris with a single border edge,
   * since otherwise that would risk isolating a vertex if
   * the tri ends up being eroded from the hull.
   * Sets the tri size according to the threshold parameter being used.
   * 
   * @param tri the Tri to add
   * @param queue the priority queue to add to
   */
  private void addBorderTri(HullTri tri, PriorityQueue<HullTri> queue) {
    if (tri == null) return;
    if (tri.numAdjacent() != 2) return;
    setSize(tri);
    queue.add(tri);
  }

  private void setSize(HullTri tri) {
    if (criteriaType == PARAM_EDGE_LENGTH)
      tri.setSizeToBoundary();
    else 
      tri.setSizeToCircumradius();
  }
  
  /**
   * Tests if a tri is included in the hull.
   * Tris with size less than the maximum are included in the hull.
   * 
   * @param tri the tri to test
   * @return true if the tri is included in the hull
   */
  private boolean isInHull(HullTri tri) {
    return tri.getSize() < maxSizeInHull;
  }
  
  private void computeHullHoles(List<HullTri> triList) {
    List<HullTri> candidateHoles = findCandidateHoles(triList, maxSizeInHull);
    // remove tris in order of decreasing size (edge length)
    for (HullTri tri : candidateHoles) {
      if (tri.isRemoved() 
          || tri.isBorder() 
          || tri.hasBoundaryTouch())
        continue;
      removeHole(triList, tri);
    }
  }
  
  /**
   * Finds tris which may be the start of holes.
   * Only tris which have a long enough edge and which do not touch the current hull
   * boundary are included.
   * This avoids the risk of disconnecting the result polygon.
   * The list is sorted in decreasing order of size.
   * 
   * @param triList
   * @param maxSizeInHull maximum tri size which is not in a hole
   * @return
   */
  private static List<HullTri> findCandidateHoles(List<HullTri> triList, double maxSizeInHull) {
    List<HullTri> candidates = new ArrayList<HullTri>();
    for (HullTri tri : triList) {
      //-- tris below the size threshold are in the hull, so NOT in a hole
      if (tri.getSize() < maxSizeInHull) continue;
      
      boolean isTouchingBoundary = tri.isBorder() || tri.hasBoundaryTouch();
      if (! isTouchingBoundary) {
        candidates.add(tri);
      }
    }
    // sort by HullTri comparator - larger sizes first
    candidates.sort(null);
    return candidates;
  }
  
  /**
   * Erodes a hole starting at a given triangle, 
   * and eroding all adjacent triangles with boundary edge length above target.
   * @param triList the triangulation
   * @param triHole triangle which is a hole
   */
  private void removeHole(List<HullTri> triList, HullTri triHole) {
    PriorityQueue<HullTri> queue = new PriorityQueue<HullTri>();
    queue.add(triHole);
    
    while (! queue.isEmpty()) {
      HullTri tri = queue.poll();
      
      if (tri != triHole && isInHull(tri)) 
        break;
      
      if (tri == triHole || isRemovableHole(tri)) {
        //-- the non-null adjacents are now on the border
        HullTri adj0 = (HullTri) tri.getAdjacent(0);
        HullTri adj1 = (HullTri) tri.getAdjacent(1);
        HullTri adj2 = (HullTri) tri.getAdjacent(2);
        
        tri.remove(triList);
        
        //-- add border adjacents to queue
        addBorderTri(adj0, queue);
        addBorderTri(adj1, queue);
        addBorderTri(adj2, queue);
      }
    }
  }
  
  private boolean isRemovableBorder(HullTri tri) {
    /**
     * Tri must have exactly 2 adjacent tris (i.e. a single boundary edge).
     * If it it has only 0 or 1 adjacent then removal would remove a vertex.
     * If it has 3 adjacent then it is not on border.
     */
    if (tri.numAdjacent() != 2) return false;
    /**
     * The tri cannot be removed if it is connecting, because
     * this would create more than one result polygon.
     */
    return ! tri.isConnecting();
  }
  
  private boolean isRemovableHole(HullTri tri) {
    /**
     * Tri must have exactly 2 adjacent tris (i.e. a single boundary edge).
     * If it it has only 0 or 1 adjacent then removal would remove a vertex.
     * If it has 3 adjacent then it is not connected to hole.
     */
    if (tri.numAdjacent() != 2) return false;
    /**
     * Ensure removal does not disconnect hull area.
     * This is a fast check which ensure holes and boundary
     * do not touch at single points.
     * (But it is slightly over-strict, since it prevents
     * any touching holes.)
     */
    return ! tri.hasBoundaryTouch();
  }
  
  private Geometry toGeometry(List<HullTri> triList, GeometryFactory geomFactory) {
    if (! isHolesAllowed) {
      return HullTriangulation.traceBoundaryPolygon(triList, geomFactory);
    }
    //-- in case holes are present use union (slower but handles holes)
    return HullTriangulation.union(triList, geomFactory);
  }
}