From d33f08388c8406f1cdaf5568c8c534554a1b6b58 Mon Sep 17 00:00:00 2001
From: Martin Davis <mtnclimb@gmail.com>
Date: Tue, 4 Jan 2022 09:50:04 -0800
Subject: [PATCH] ConcaveHull for Points (#823)

* ConcaveHull for Points, with Length, Length Factor and Area Ratio target criteria

Signed-off-by: Martin Davis <mtnclimb@gmail.com>
---
 .../function/ConstructionFunctions.java       |  55 ++
 .../jts/algorithm/hull/ConcaveHull.java       | 769 ++++++++++++++++++
 .../org/locationtech/jts/geom/Triangle.java   |  25 +-
 .../locationtech/jts/triangulate/tri/Tri.java |  63 +-
 .../triangulate/tri/TriangulationBuilder.java |   4 +-
 .../jts/algorithm/hull/ConcaveHullTest.java   | 159 ++++
 .../locationtech/jts/hull/ConcaveHull.java    |  53 --
 .../jts/hull/ConcaveHullTest.java             |  45 -
 8 files changed, 1066 insertions(+), 107 deletions(-)
 create mode 100644 modules/core/src/main/java/org/locationtech/jts/algorithm/hull/ConcaveHull.java
 create mode 100644 modules/core/src/test/java/org/locationtech/jts/algorithm/hull/ConcaveHullTest.java
 delete mode 100644 modules/lab/src/main/java/org/locationtech/jts/hull/ConcaveHull.java
 delete mode 100644 modules/lab/src/test/java/org/locationtech/jts/hull/ConcaveHullTest.java

diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/ConstructionFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/ConstructionFunctions.java
index 13df15dc05..e3301c8258 100644
--- a/modules/app/src/main/java/org/locationtech/jtstest/function/ConstructionFunctions.java
+++ b/modules/app/src/main/java/org/locationtech/jtstest/function/ConstructionFunctions.java
@@ -16,6 +16,7 @@
 import org.locationtech.jts.algorithm.MinimumDiameter;
 import org.locationtech.jts.algorithm.construct.LargestEmptyCircle;
 import org.locationtech.jts.algorithm.construct.MaximumInscribedCircle;
+import org.locationtech.jts.algorithm.hull.ConcaveHull;
 import org.locationtech.jts.densify.Densifier;
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.Geometry;
@@ -130,4 +131,58 @@ public static Geometry circleByRadiusLine(Geometry radiusLine,
     return radiusLine.getFactory().createPolygon(circlePts);
   }
  
+  public static Geometry concaveHullByLen(Geometry geom, 
+      @Metadata(title="Length")
+      double maxLen) {
+    return ConcaveHull.concaveHullByLength(geom, maxLen);
+  }
+  
+  public static Geometry concaveHullWithHolesByLen(Geometry geom, 
+      @Metadata(title="Length")
+      double maxLen) {
+    return ConcaveHull.concaveHullByLength(geom, maxLen, true);
+  }
+  
+  public static Geometry concaveHullByLenFactor(Geometry geom, 
+      @Metadata(title="Length factor")
+      double maxLen) {
+    return ConcaveHull.concaveHullByLengthFactor(geom, maxLen);
+  }
+  
+  public static Geometry concaveHullWithHolesByLenFactor(Geometry geom, 
+      @Metadata(title="Length factor")
+      double maxLen) {
+    return ConcaveHull.concaveHullByLengthFactor(geom, maxLen, true);
+  }
+  
+  public static Geometry concaveHullByArea(Geometry geom, 
+      @Metadata(title="Area ratio")
+      double minAreaPct) {
+    return ConcaveHull.concaveHullByArea(geom, minAreaPct);
+  }
+  
+  public static double concaveHullLenGuess(Geometry geom) {
+    return ConcaveHull.uniformGridEdgeLength(geom);
+  }
+  
+  /**
+   * A concaveness measure defined in terms of the perimeter length
+   * relative to the convex hull perimeter.
+   * <pre>
+   * C = ( P(geom) - P(CH) ) / P(CH)
+   * </pre>
+   * Concaveness values are >= 0.  
+   * A convex polygon has C = 0. 
+   * A higher concaveness indicates a more concave polygon.
+   * <p>
+   * Originally defined by Park & Oh, 2012.
+   * 
+   * @param geom a polygonal geometry
+   * @return the concaveness measure of the geometry
+   */
+  public static double concaveness(Geometry geom) {
+    double convexLen = geom.convexHull().getLength();
+    return (geom.getLength() - convexLen) / convexLen;
+  }
+  
 }
diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/hull/ConcaveHull.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/hull/ConcaveHull.java
new file mode 100644
index 0000000000..991ac27c7c
--- /dev/null
+++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/hull/ConcaveHull.java
@@ -0,0 +1,769 @@
+/*
+ * 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.ArrayDeque;
+import java.util.ArrayList;
+import java.util.Deque;
+import java.util.List;
+import java.util.PriorityQueue;
+
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateList;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.Triangle;
+import org.locationtech.jts.operation.overlayng.CoverageUnion;
+import org.locationtech.jts.triangulate.DelaunayTriangulationBuilder;
+import org.locationtech.jts.triangulate.quadedge.QuadEdge;
+import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
+import org.locationtech.jts.triangulate.quadedge.TriangleVisitor;
+import org.locationtech.jts.triangulate.tri.Tri;
+import org.locationtech.jts.triangulate.tri.TriangulationBuilder;
+import org.locationtech.jts.util.Assert;
+
+/**
+ * Constructs a concave hull of a set of points.
+ * The hull is constructed by removing the longest outer edges 
+ * of the Delaunay Triangulation of the points
+ * until certain target criteria are reached.
+ * 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 Factor</b> - determine the Maximum Edge Length 
+ * as 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-independent.
+ * <li><b>Maximum Area Ratio</b> - the ratio of the concave hull area to the convex hull area 
+ * will be no larger than this value.
+ * </ul>
+ * Usually only a single criteria is specified, but both may be provided.
+ * The preferred criteria is the <b>Maximum Edge Length Factor</b>, since it is 
+ * scale-independent, and local (so that no assumption needs to be made about the 
+ * total amount of concavity present).
+ * Other length criteria can be used by setting the Maximum Edge Length.
+ * 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 criteria.
+ * <p>
+ * Optionally the concave hull can be allowed to contain holes.
+ * Note that this may result in substantially slower computation,
+ * and it can produce results of low quality.
+ * 
+ * @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 the concave hull of the vertices in a geometry
+   * using the target criteria 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 the concave hull of the vertices in a geometry
+   * using the target criteria 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 the concave hull of the vertices in a geometry
+   * using the target criteria of maximum edge length factor.
+   * 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 lengthFactor the target edge length factor
+   * @return the concave hull
+   */
+  public static Geometry concaveHullByLengthFactor(Geometry geom, double lengthFactor) {
+    return concaveHullByLengthFactor(geom, lengthFactor, false);
+  }
+  
+  /**
+   * Computes the concave hull of the vertices in a geometry
+   * using the target criteria 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 concaveHullByLengthFactor(Geometry geom, double lengthFactor, boolean isHolesAllowed) {
+    ConcaveHull hull = new ConcaveHull(geom);
+    hull.setMaximumEdgeLengthFactor(lengthFactor);
+    hull.setHolesAllowed(isHolesAllowed);
+    return hull.getHull();
+  }
+  
+  /**
+   * Computes the concave hull of the vertices in a geometry
+   * using the target criteria of maximum area ratio.
+   * 
+   * @param geom the input geometry
+   * @param areaRatio the target maximum area ratio
+   * @return the concave hull
+   */
+  public static Geometry concaveHullByArea(Geometry geom, double areaRatio) {
+    ConcaveHull hull = new ConcaveHull(geom);
+    hull.setMaximumAreaRatio(areaRatio);
+    return hull.getHull();
+  }
+  
+  private Geometry inputGeometry;
+  private double maxEdgeLength = 0.0;
+  private double maxEdgeLengthFactor = -1;
+  private double maxAreaRatio = 0.0; 
+  private boolean isHolesAllowed = false;
+  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.maxEdgeLength = edgeLength;
+    maxEdgeLengthFactor = -1;
+  }
+  
+  /**
+   * Sets the target maximum edge length factor for the concave hull.
+   * The edge length factor 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 edgeLengthFactor a length factor value between 0 and 1
+   */
+  public void setMaximumEdgeLengthFactor(double edgeLengthFactor) {
+    if (edgeLengthFactor < 0 || edgeLengthFactor > 1)
+      throw new IllegalArgumentException("Edge length ratio must be in range [0,1]e");
+    this.maxEdgeLengthFactor = edgeLengthFactor;
+  }
+  
+  /**
+   * Sets the target maximum concave hull area as a ratio of the convex hull area.
+   * It is a value in the range 0 to 1. 
+   * <ul>
+   * <li>The value 0.0 produces a concave hull with the smallest area
+   * that is still connected.
+   * <li>The value 1.0 produces the convex hull 
+   * (unless a maximum edge length is also specified).
+   * </ul>
+   * 
+   * @param areaRatio a ratio value between 0 and 1
+   */
+  public void setMaximumAreaRatio(double areaRatio) {
+    if (areaRatio < 0 || areaRatio > 1)
+      throw new IllegalArgumentException("Area ratio must be in range [0,1]");
+    this.maxAreaRatio = areaRatio;
+  }
+  
+  /**
+   * 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 = createDelaunayTriangulation(inputGeometry);
+    if (maxEdgeLengthFactor >= 0) {
+      maxEdgeLength = computeTargetEdgeLength(triList, maxEdgeLengthFactor);
+    }
+    if (triList.isEmpty())
+      return inputGeometry.convexHull();
+    computeHull(triList);
+    Geometry hull = toPolygon(triList, geomFactory);
+    return hull;
+  }
+
+  private static double computeTargetEdgeLength(List<? extends Tri> triList, 
+      double edgeLengthFactor) {
+    if (edgeLengthFactor == 0) return 0;
+    double maxEdgeLen = -1;
+    double minEdgeLen = -1;
+    for (Tri tri : triList) {
+      for (int i = 0; i < 3; i++) {
+        double len = tri.getCoordinate(i).distance(tri.getCoordinate(Tri.next(i)));
+        if (len > maxEdgeLen) 
+          maxEdgeLen = len;
+        if (minEdgeLen < 0 || len < minEdgeLen)
+          minEdgeLen = len;
+      }
+    }
+    //-- ensure all edges are included
+    if (edgeLengthFactor == 1) return 2 * maxEdgeLen;
+    return edgeLengthFactor * (maxEdgeLen - minEdgeLen) + minEdgeLen;
+  }
+
+  private void computeHull(List<HullTri> triList) {
+    //-- used if area is the threshold criteria
+    double areaConvex = Tri.area(triList);
+    double areaConcave = areaConvex;
+    
+    PriorityQueue<HullTri> queue = initQueue(triList);
+    // remove tris in order of decreasing size (edge length)
+    while (! queue.isEmpty()) {
+      if (isBelowAreaThreshold(areaConcave, areaConvex))
+        break;
+
+      HullTri tri = queue.poll();
+      
+      if (isBelowLengthThreshold(tri)) 
+        break;
+      
+      if (isRemovable(tri, triList)) {
+        //-- 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);
+        
+        //-- remove tri
+        tri.remove();
+        triList.remove(tri);
+        areaConcave -= tri.getArea();
+        
+        //-- if holes not allowed, add new border adjacents to queue
+        if (! isHolesAllowed) {
+          addBorderTri(adj0, queue);
+          addBorderTri(adj1, queue);
+          addBorderTri(adj2, queue);
+        }
+      }
+    }
+  }
+
+  private PriorityQueue<HullTri> initQueue(List<HullTri> triList) {
+    PriorityQueue<HullTri> queue = new PriorityQueue<HullTri>();
+    for (HullTri tri : triList) {
+      if (! isHolesAllowed) {
+        //-- add only border triangles which could be eroded
+        // (if tri has only 1 adjacent it can't be removed because that would isolate a vertex)
+        if (tri.numAdjacent() != 2)
+          continue;
+        tri.setSizeToBorder();
+      }
+      queue.add(tri);
+    }
+    return queue;
+  }
+
+  /**
+   * Adds a Tri to the queue.
+   * Only add tris with a single border edge.
+   * The ordering size is the length of the border edge.
+   * 
+   * @param tri the Tri to add
+   * @param queue the priority queue
+   */
+  private void addBorderTri(HullTri tri, PriorityQueue<HullTri> queue) {
+    if (tri == null) return;
+    if (tri.numAdjacent() != 2) return;
+    tri.setSizeToBorder();
+    queue.add(tri);
+  }
+    
+  private boolean isBelowAreaThreshold(double areaConcave, double areaConvex) {
+    return areaConcave / areaConvex <= maxAreaRatio;
+  }
+
+  private boolean isBelowLengthThreshold(HullTri tri) {
+    double len = 0;
+    if (isHolesAllowed) {
+      len = tri.lengthOfLongestEdge();
+    }
+    else {
+      len = tri.lengthOfBorder();
+    }
+    return len < maxEdgeLength;
+  }
+
+  /**
+   * Tests whether a Tri can be removed while preserving 
+   * the connectivity of the hull.
+   * 
+   * @param tri the Tri to test
+   * @param triList 
+   * @return true if the Tri can be removed
+   */
+  private boolean isRemovable(HullTri tri, List<HullTri> triList) {
+    if (isHolesAllowed) {
+      /**
+       * Don't remove if that would separate a single vertex
+       */
+      if (hasVertexSingleAdjacent(tri, triList))
+        return false;
+      return HullTri.isConnected(triList, tri);
+    }
+    
+    //-- compute removable for no holes allowed
+    int numAdj = tri.numAdjacent();
+    /**
+     * Tri must have exactly 2 adjacent tris.
+     * 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 (numAdj != 2) return false;
+    /**
+     * The tri cannot be removed if it is connecting, because
+     * this would create more than one result polygon.
+     */
+    return ! isConnecting(tri);
+  }
+
+  private static boolean hasVertexSingleAdjacent(HullTri tri, List<HullTri> triList) {
+    for (int i = 0; i < 3; i++) {
+      if (degree(tri.getCoordinate(i), triList) <= 1)
+        return true;
+    }
+    return false;
+  }
+
+  /**
+   * The degree of a Tri vertex is the number of tris containing it.
+   * This must be done by searching the entire triangulation, 
+   * since the containing tris may not be adjacent or edge-connected. 
+   * 
+   * @param v a vertex coordinate
+   * @param triList a triangulation
+   * @return the degree of the vertex
+   */
+  private static int degree(Coordinate v, List<HullTri> triList) {
+    int degree = 0;
+    for (HullTri tri : triList) {
+      for (int i = 0; i < 3; i++) {
+        if (v.equals2D(tri.getCoordinate(i)))
+          degree++;
+      }
+    }
+    return degree;
+  }
+
+  /**
+   * Tests if a tri is the only one connecting its 2 adjacents.
+   * Assumes that the tri is on the border of the triangulation
+   * and that the triangulation does not contain holes
+   * 
+   * @param tri the tri to test
+   * @return true if the tri is the only connection
+   */
+  private static boolean isConnecting(Tri tri) {
+    int adj2Index = adjacent2VertexIndex(tri);
+    boolean isInterior = isInteriorVertex(tri, adj2Index);
+    return ! isInterior;
+  }
+
+  /**
+   * A vertex of a triangle is interior if it 
+   * is fully surrounded by triangles.
+   * 
+   * @param tri a tri containing the vertex
+   * @param index the vertex index
+   * @return true if the vertex is interior
+   */
+  private static boolean isInteriorVertex(Tri triStart, int index) {
+    Tri curr = triStart;
+    int currIndex = index;
+    do {
+      Tri adj = curr.getAdjacent(currIndex);
+      if (adj == null) return false;
+      int adjIndex = adj.getIndex(curr);
+      curr = adj;
+      currIndex = Tri.next(adjIndex);
+    }
+    while (curr != triStart);
+    return true;
+  }
+
+  private static int adjacent2VertexIndex(Tri tri) {
+    if (tri.hasAdjacent(0) && tri.hasAdjacent(1)) return 1;
+    if (tri.hasAdjacent(1) && tri.hasAdjacent(2)) return 2;
+    if (tri.hasAdjacent(2) && tri.hasAdjacent(0)) return 0;
+    return -1;
+  }
+  
+  private static class HullTri extends Tri 
+      implements Comparable<HullTri> 
+  {
+    private double size;
+    private boolean isMarked = false;
+    
+    public HullTri(Coordinate p0, Coordinate p1, Coordinate p2) {
+      super(p0, p1, p2);
+      this.size = lengthOfLongestEdge();
+    }
+
+    public double getSize() {
+      return size;
+    }
+    
+    /**
+     * Sets the size to be the length of the border edges.
+     * This is used when constructing hull without holes,
+     * by erosion from the triangulation border.
+     */
+    public void setSizeToBorder() {
+      size = lengthOfBorder();
+    }
+    
+    public boolean isMarked() {
+      return isMarked;
+    }
+    
+    public void setMarked(boolean isMarked) {
+      this.isMarked = isMarked;
+    }
+
+    public boolean isBorder() {
+      return isBorder(0) || isBorder(1) || isBorder(2);
+    }
+    
+    public boolean isBorder(int index) {
+      return ! hasAdjacent(index);
+    }
+    
+    public int borderIndex() {
+      if (isBorder(0)) return 0;
+      if (isBorder(1)) return 1;
+      if (isBorder(2)) return 2;
+      return -1;
+    }
+    
+    /**
+     * Gets the most CCW border edge index.
+     * This assumes there is at least one non-border edge.
+     * 
+     * @return the CCW border edge index
+     */
+    public int borderIndexCCW() {
+      int index = borderIndex();
+      int prevIndex = prev(index);
+      if (isBorder(prevIndex)) {
+        return prevIndex;
+      }
+      return index;
+    }
+    
+    /**
+     * Gets the most CW border edge index.
+     * This assumes there is at least one non-border edge.
+     * 
+     * @return the CW border edge index
+     */
+    public int borderIndexCW() {
+      int index = borderIndex();
+      int nextIndex = next(index);
+      if (isBorder(nextIndex)) {
+        return nextIndex;
+      }
+      return index;
+    }
+    
+    public double lengthOfLongestEdge() {
+      return Triangle.longestSideLength(p0, p1, p2);
+    }
+    
+    private double lengthOfBorder() {
+      double len = 0.0;
+      for (int i = 0; i < 3; i++) {
+        if (! hasAdjacent(i)) {
+          len += getCoordinate(i).distance(getCoordinate(Tri.next(i)));
+        }
+      }
+      return len;
+    }
+    
+    public HullTri nextBorderTri() {
+      HullTri tri = this;
+      //-- start at first non-border edge CW
+      int index = next(borderIndexCW());
+      //-- scan CCW around vertex for next border tri
+      do {
+        HullTri adjTri = (HullTri) tri.getAdjacent(index);
+        if (adjTri == this)
+          throw new IllegalStateException("No outgoing border edge found");
+        index = next(adjTri.getIndex(tri));
+        tri = adjTri;
+      }
+      while (! tri.isBorder(index));
+      return (tri);
+    }
+
+    /**
+     * PriorityQueues sort in ascending order.
+     * To sort with the largest at the head,
+     * smaller sizes must compare as greater than larger sizes.
+     * (i.e. the normal numeric comparison is reversed).
+     * If the sizes are identical (which should be an infrequent case),
+     * the areas are compared, with larger areas sorting before smaller.
+     * (The rationale is that larger areas indicate an area of lower point density,
+     * which is more likely to be in the exterior of the computed shape.)
+     * This improves the determinism of the queue ordering. 
+     */
+    @Override
+    public int compareTo(HullTri o) {
+      /**
+       * If size is identical compare areas to ensure a (more) deterministic ordering.
+       * Larger areas sort before smaller ones.
+       */
+      if (size == o.size) {
+        return -Double.compare(this.getArea(), o.getArea());
+      }
+      return -Double.compare(size, o.size);
+    }
+    
+    public static boolean isConnected(List<HullTri> triList, HullTri exceptTri) {
+      if (triList.size() == 0) return false;
+      clearMarks(triList);
+      HullTri triStart = findTri(triList, exceptTri);
+      if (triStart == null) return false;
+      markConnected(triStart, exceptTri);
+      exceptTri.setMarked(true);
+      return isAllMarked(triList);
+    }
+    
+    public static void clearMarks(List<HullTri> triList) {
+      for (HullTri tri : triList) {
+        tri.setMarked(false);
+      }
+    }
+    
+    public static HullTri findTri(List<HullTri> triList, Tri exceptTri) {
+      for (HullTri tri : triList) {
+        if (tri != exceptTri) return tri;
+      }
+      return null;
+    }
+    
+    public static boolean isAllMarked(List<HullTri> triList) {
+      for (HullTri tri : triList) {
+        if (! tri.isMarked())
+          return false;
+      }
+      return true;
+    }
+    
+    public static void markConnected(HullTri triStart, Tri exceptTri) {
+      Deque<HullTri> queue = new ArrayDeque<HullTri>();
+      queue.add(triStart);
+      while (! queue.isEmpty()) {
+        HullTri tri = queue.pop();
+        tri.setMarked(true);
+        for (int i = 0; i < 3; i++) {
+          HullTri adj = (HullTri) tri.getAdjacent(i);
+          //-- don't connect thru this tri
+          if (adj == exceptTri)
+            continue;
+          if (adj != null && ! adj.isMarked() ) {
+            queue.add(adj);
+          }
+        }
+      }
+    }
+  }
+  
+  private static List<HullTri> createDelaunayTriangulation(Geometry geom) {
+    //TODO: implement a DT on Tris directly?
+    DelaunayTriangulationBuilder dt = new DelaunayTriangulationBuilder();
+    dt.setSites(geom);
+    QuadEdgeSubdivision subdiv = dt.getSubdivision();
+    List<HullTri> triList = toTris(subdiv);
+    return triList;
+  }
+  
+  private static List<HullTri> toTris(QuadEdgeSubdivision subdiv) {
+    HullTriVisitor visitor = new HullTriVisitor();
+    subdiv.visitTriangles(visitor, false);
+    List<HullTri> triList = visitor.getTriangles();
+    TriangulationBuilder.build(triList);
+    return triList;
+  }
+  
+  private static class HullTriVisitor implements TriangleVisitor {
+    private List<HullTri> triList = new ArrayList<HullTri>();
+
+    public HullTriVisitor() {
+    }
+
+    public void visit(QuadEdge[] triEdges) {
+      Coordinate p0 = triEdges[0].orig().getCoordinate();
+      Coordinate p1 = triEdges[1].orig().getCoordinate();
+      Coordinate p2 = triEdges[2].orig().getCoordinate();
+      HullTri tri;
+      if (Triangle.isCCW(p0, p1, p2)) {
+        tri = new HullTri(p0, p2, p1);
+      }
+      else {
+        tri = new HullTri(p0, p1, p2);
+      }
+      triList.add(tri);
+    }
+    
+    public List<HullTri> getTriangles() {
+      return triList;
+    }
+  }
+
+  private Geometry toPolygon(List<HullTri> triList, GeometryFactory geomFactory) {
+    if (! isHolesAllowed) {
+      return extractPolygon(triList, geomFactory);
+    }
+    //-- in case holes are present use union (slower but handles holes)
+    return union(triList, geomFactory);
+  }
+
+  private Geometry extractPolygon(List<HullTri> triList, GeometryFactory geomFactory) {
+    if (triList.size() == 1) {
+      Tri tri = triList.get(0);
+      return tri.toPolygon(geomFactory);
+    }
+    Coordinate[] pts = traceBorder(triList);
+    return geomFactory.createPolygon(pts);
+  }
+
+  private static Geometry union(List<? extends Tri> triList, GeometryFactory geomFactory) {
+    List<Polygon> polys = new ArrayList<Polygon>();
+    for (Tri tri : triList) {
+      Polygon poly = tri.toPolygon(geomFactory);
+      polys.add(poly);
+    }
+    return CoverageUnion.union(geomFactory.buildGeometry(polys));
+  }
+  
+  /**
+   * Extracts the coordinates along the border of a triangulation,
+   * by tracing CW around the border triangles.
+   * Assumption: there are at least 2 tris, they are connected,
+   * and there are no holes.
+   * So each tri has at least one non-border edge, and there is only one border.
+   * 
+   * @param triList the triangulation
+   * @return the border of the triangulation
+   */
+  private static Coordinate[] traceBorder(List<HullTri> triList) {
+    HullTri triStart = findBorderTri(triList);
+    CoordinateList coordList = new CoordinateList();
+    HullTri tri = triStart;
+    do {
+      int borderIndex = tri.borderIndexCCW();
+      //-- add border vertex
+      coordList.add(tri.getCoordinate(borderIndex).copy(), false);
+      int nextIndex = Tri.next(borderIndex);
+      //-- if next edge is also border, add it and move to next
+      if (tri.isBorder(nextIndex)) {
+        coordList.add(tri.getCoordinate(nextIndex).copy(), false);
+        borderIndex = nextIndex;
+      }
+      //-- find next border tri CCW around non-border edge
+      tri = tri.nextBorderTri();
+    } while (tri != triStart);
+    coordList.closeRing();
+    return coordList.toCoordinateArray();
+  }
+  
+  public static HullTri findBorderTri(List<HullTri> triList) {
+    for (HullTri tri : triList) {
+      if (tri.isBorder()) return tri;
+    }
+    Assert.shouldNeverReachHere("No border triangles found");
+    return null;
+  }
+}
diff --git a/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java b/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java
index 8cdaa1d16f..8188a166b7 100644
--- a/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java
+++ b/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java
@@ -290,6 +290,19 @@ public static Coordinate centroid(Coordinate a, Coordinate b, Coordinate c)
     return new Coordinate(x, y);
   }
 
+  /**
+   * Compute the length of the perimeter of a triangle
+   * 
+   * @param a a vertex of the triangle
+   * @param b a vertex of the triangle
+   * @param c a vertex of the triangle
+   * @return the length of the triangle perimeter
+   */
+  public static double length(Coordinate a, Coordinate b, Coordinate c)
+  {
+    return a.distance(b) + b.distance(c) + c.distance(a);
+  }
+  
   /**
    * Computes the length of the longest side of a triangle
    * 
@@ -430,7 +443,7 @@ public static double area3D(Coordinate a, Coordinate b, Coordinate c)
 
     return area3D;
   }
-
+  
   /**
    * Computes the Z-value (elevation) of an XY point on a three-dimensional
    * plane defined by a triangle whose vertices have Z-values. The defining
@@ -563,6 +576,16 @@ public Coordinate centroid()
     return centroid(p0, p1, p2);
   }
 
+  /**
+   * Computes the length of the perimeter of this triangle.
+   * 
+   * @return the length of the perimeter
+   */
+  public double length()
+  {
+    return length(p0, p1, p2);
+  }
+  
   /**
    * Computes the length of the longest side of this triangle
    * 
diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java
index 337905e815..e2af1e2d12 100755
--- a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java
+++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java
@@ -19,6 +19,7 @@
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.GeometryFactory;
 import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.Triangle;
 import org.locationtech.jts.io.WKTWriter;
 import org.locationtech.jts.util.Assert;
 
@@ -49,6 +50,20 @@ public static Geometry toGeometry(List<Tri> triList, GeometryFactory geomFact) {
     return geomFact.createGeometryCollection(geoms);
   }
 
+  /**
+   * Computes the area of a set of Tris.
+   * 
+   * @param triList a set of Tris
+   * @return the total area of the triangles
+   */
+  public static double area(List<? extends Tri> triList) {
+    double area = 0;
+    for (Tri tri : triList) {
+      area += tri.getArea();
+    }
+    return area;
+  }
+  
   /**
    * Validates a list of Tris.
    * 
@@ -84,17 +99,17 @@ public static Tri create(Coordinate[] pts) {
     return new Tri(pts[0], pts[1], pts[2]);
   }
   
-  private Coordinate p0;
-  private Coordinate p1;
-  private Coordinate p2;
+  protected Coordinate p0;
+  protected Coordinate p1;
+  protected Coordinate p2;
   
   /**
    * triN is the adjacent triangle across the edge pN - pNN.
    * pNN is the next vertex CW from pN.
    */
-  private Tri tri0;
-  private Tri tri1;
-  private Tri tri2;
+  protected Tri tri0;
+  protected Tri tri1;
+  protected Tri tri2;
 
   /**
    * Creates a triangle with the given vertices.
@@ -246,6 +261,24 @@ private void replace(Tri triOld, Tri triNew) {
     }
   }
   
+  /**
+   * Removes this triangle from a triangulation.
+   * All adjacent references and the references to this
+   * Tri in the adjacent Tris are set to <code>null</code.
+   */
+  public void remove() {
+    remove(0);
+    remove(1);
+    remove(2);
+  }
+
+  private void remove(int index) {
+    Tri adj = getAdjacent(index);
+    if (adj == null) return;
+    adj.setTri(adj.getIndex(this), null);
+    setTri(index, null);
+  }
+  
   /**
    * Gets the triangles adjacent to the quadrilateral
    * formed by this triangle and an adjacent one.
@@ -528,6 +561,24 @@ public Coordinate midpoint(int edgeIndex) {
     return new Coordinate(midX, midY);
   }
   
+  /**
+   * Gets the area of the triangle.
+   * 
+   * @return the area of the triangle
+   */
+  public double getArea() {
+    return Triangle.area(p0, p1, p2);
+  }
+  
+  /**
+   * Gets the length of the perimeter of the triangle.
+   * 
+   * @return the length of the perimeter
+   */
+  public double getLength() {
+    return Triangle.length(p0, p1, p2);
+  }
+  
   /**
    * Creates a {@link Polygon} representing this triangle.
    * 
diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java
index a4828fc9d4..5a9fd381fb 100755
--- a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java
+++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java
@@ -30,7 +30,7 @@ public class TriangulationBuilder {
    * 
    * @param triList the list of Tris
    */
-  public static void build(List<Tri> triList) {
+  public static void build(List<? extends Tri> triList) {
     new TriangulationBuilder(triList);
   }
   
@@ -41,7 +41,7 @@ public static void build(List<Tri> triList) {
    * 
    * @param triList the list of Tris
    */
-  private TriangulationBuilder(List<Tri> triList) {
+  private TriangulationBuilder(List<? extends Tri> triList) {
     triMap = new HashMap<TriEdge, Tri>();
     for (Tri tri : triList) {
       add(tri);
diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/hull/ConcaveHullTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/hull/ConcaveHullTest.java
new file mode 100644
index 0000000000..6d8f05ba03
--- /dev/null
+++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/hull/ConcaveHullTest.java
@@ -0,0 +1,159 @@
+/*
+ * 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 org.locationtech.jts.geom.Geometry;
+
+import junit.textui.TestRunner;
+import test.jts.GeometryTestCase;
+
+public class ConcaveHullTest extends GeometryTestCase {
+  
+  public static void main(String args[]) {
+    TestRunner.run(ConcaveHullTest.class);
+  }
+
+  public ConcaveHullTest(String name) { super(name); }
+  
+  public void testLengthEmpty() {
+    checkHullByLength("MULTIPOINT EMPTY", 
+       70, "POLYGON EMPTY" );
+  }
+
+  public void testLengthPoint() {
+    checkHullByLength("MULTIPOINT ((10 10), (10 10))", 
+       70, "POINT (10 10)" );
+  }
+
+  public void testLengthCollinear() {
+    checkHullByLength("LINESTRING (10 10, 20 20, 30 30))", 
+       70, "LINESTRING (10 10, 30 30)" );
+  }
+
+  public void testLengthTriangle() {
+    checkHullByLength("MULTIPOINT ((10 10), (90 10), (30 70))", 
+       70, "POLYGON ((10 10, 30 70, 90 10, 10 10))" );
+  }
+
+  public void testLengthChevron() {
+    checkHullByLength("MULTIPOINT ((10 10), (90 10), (30 70), (70 70), (50 60))", 
+       70, "POLYGON ((30 70, 70 70, 90 10, 50 60, 10 10, 30 70))" );
+  }
+
+  public void testLengthZero() {
+    checkHullByLength("MULTIPOINT ((10 10), (90 10), (70 70), (50 60), (50 90), (40 70), (30 30))", 
+       0, "POLYGON ((10 10, 40 70, 50 90, 70 70, 90 10, 50 60, 30 30, 10 10))" );
+  }
+
+  public void testLengthConvex() {
+    checkHullByLength("MULTIPOINT ((10 10), (90 10), (70 70), (50 60), (50 90), (40 70), (30 30))", 
+       100, "POLYGON ((10 10, 40 70, 50 90, 70 70, 90 10, 10 10))" );
+  }
+
+  public void testLengthCShape() {
+    checkHullByLength("MULTIPOINT ((70 80), (80 90), (90 70), (50 80), (30 70), (20 40), (30 20), (50 10), (90 20), (40 50), (40 30), (41 67))", 
+       50, "POLYGON ((30 70, 50 80, 80 90, 90 70, 70 80, 40 50, 40 30, 90 20, 50 10, 30 20, 20 40, 30 70))" );
+  }
+
+  public void testLengthSShape() {
+    checkHullByLength("MULTIPOINT ((0 81), (65 86), (70 71), (80 59), (92 49), (107 44), (122 41), (137 40), (152 41), (167 42), (182 47), (195 55), (203 68), (201 83), (188 92), (173 97), (158 100), (143 103), (128 106), (113 109), (98 112), (83 115), (68 120), (53 125), (40 133), (28 143), (18 155), (13 170), (12 185), (16 200), (26 213), (38 223), (51 231), (66 236), (81 240), (96 243), (111 245), (126 245), (141 245), (156 245), (171 244), (186 241), (201 238), (216 233), (229 225), (242 216), (252 204), (259 190), (262 175), (194 171), (189 186), (178 197), (164 203), (149 205), (134 206), (119 205), (104 203), (89 198), (77 188), (80 173), (93 165), (108 160), (123 157), (138 154), (153 151), (168 149), (183 146), (198 142), (213 138), (227 132), (241 126), (253 116), (263 104), (269 90), (271 75), (270 60), (264 46), (254 34), (243 23), (229 16), (215 10), (200 6), (185 3), (170 1), (155 0), (139 0), (123 0), (108 1), (93 3), (78 5), (63 10), (49 16), (35 23), (23 33), (13 45), (6 59), (16 82), (32 83), (48 84), (245 174), (228 173), (211 172), (131 128), (63 148), (222 207), (127 230), (154 131), (240 82), (72 220), (210 32), (90 22), (206 208), (57 202), (195 117), (55 166), (246 55), (201 101), (224 73), (211 192), (42 176), (152 228), (172 113), (24 61), (76 33), (92 216), (46 69), (118 138), (169 23), (213 118), (221 56), (44 192), (118 22), (224 40), (56 57), (192 32), (179 220), (34 44), (145 18), (239 194), (40 155), (92 136), (231 106), (40 207), (108 228), (256 81), (28 185), (54 33), (74 205), (172 132), (221 93), (249 96), (69 47), (78 146), (155 115), (202 223))", 
+       20, "POLYGON ((16 200, 26 213, 38 223, 51 231, 66 236, 81 240, 96 243, 111 245, 126 245, 141 245, 156 245, 171 244, 186 241, 201 238, 216 233, 229 225, 242 216, 252 204, 259 190, 262 175, 245 174, 228 173, 211 172, 194 171, 189 186, 178 197, 164 203, 149 205, 134 206, 119 205, 104 203, 89 198, 77 188, 80 173, 93 165, 108 160, 123 157, 138 154, 153 151, 168 149, 183 146, 198 142, 213 138, 227 132, 241 126, 253 116, 263 104, 269 90, 271 75, 270 60, 264 46, 254 34, 243 23, 229 16, 215 10, 200 6, 185 3, 170 1, 155 0, 139 0, 123 0, 108 1, 93 3, 78 5, 63 10, 49 16, 35 23, 23 33, 13 45, 6 59, 0 81, 16 82, 32 83, 48 84, 65 86, 70 71, 80 59, 92 49, 107 44, 122 41, 137 40, 152 41, 167 42, 182 47, 195 55, 203 68, 201 83, 188 92, 173 97, 158 100, 143 103, 128 106, 113 109, 98 112, 83 115, 68 120, 53 125, 40 133, 28 143, 18 155, 13 170, 12 185, 16 200))" );
+  }
+
+  //------------------------------------------------
+
+  public void testLengthFactorZero() {
+    checkHullByLengthFactor("MULTIPOINT ((10 90), (10 10), (90 10), (90 90), (40 40), (60 30), (30 70), (40 60), (60 50), (60 72), (47 66), (90 60))", 
+       0, "POLYGON ((30 70, 10 90, 60 72, 90 90, 90 60, 90 10, 60 30, 10 10, 40 40, 60 50, 47 66, 40 60, 30 70))" );
+  }
+
+  public void testLengthFactorP5() {
+    checkHullByLengthFactor("MULTIPOINT ((10 90), (10 10), (90 10), (90 90), (40 40), (60 30), (30 70), (40 60), (60 50), (60 72), (47 66), (90 60))", 
+       0.5, "POLYGON ((30 70, 10 90, 60 72, 90 90, 90 60, 90 10, 60 30, 10 10, 40 40, 30 70))" );
+  }
+
+  public void testLengthFactorOne() {
+    checkHullByLengthFactor("MULTIPOINT ((10 90), (10 10), (90 10), (90 90), (40 40), (60 30), (30 70), (40 60), (60 50), (60 72), (47 66), (90 60))", 
+       1, "POLYGON ((10 10, 10 90, 90 90, 90 60, 90 10, 10 10))" );
+  }
+
+  public void testLengthFactorXYZChevronP5() {
+    checkHullByLengthFactorXYZ("MULTIPOINT Z ((10 10 1), (90 10 2), (30 70 3), (70 70 4), (50 60 5))", 
+       0.5, "POLYGON Z ((30 70 3, 70 70 4, 90 10 2, 50 60 5, 10 10 1, 30 70 3))" );
+  }
+  
+  //------------------------------------------------
+  
+  public void testLengthHolesCircle() {
+    checkHullWithHolesByLength("MULTIPOINT ((90 20), (80 10), (45 5), (10 20), (20 10), (21 30), (40 20), (11 60), (20 70), (20 90), (40 80), (70 80), (80 60), (90 70), (80 90), (56 95), (95 45), (80 40), (70 20), (15 45), (5 40), (40 96), (60 15))", 
+       40, "POLYGON ((20 90, 40 96, 56 95, 80 90, 90 70, 95 45, 90 20, 80 10, 45 5, 20 10, 10 20, 5 40, 11 60, 20 90), (20 70, 15 45, 40 20, 70 20, 80 40, 80 60, 70 80, 40 80, 20 70))" );
+  }
+
+  public void testLengthHolesCircle0() {
+    checkHullWithHolesByLength("MULTIPOINT ((90 20), (80 10), (45 5), (10 20), (20 10), (21 30), (40 20), (11 60), (20 70), (20 90), (40 80), (70 80), (80 60), (90 70), (80 90), (56 95), (95 45), (80 40), (70 20), (15 45), (5 40), (40 96), (60 15))", 
+       0, "POLYGON ((20 90, 40 96, 56 95, 80 90, 90 70, 95 45, 90 20, 80 10, 60 15, 45 5, 20 10, 10 20, 5 40, 11 60, 15 45, 21 30, 40 20, 70 20, 80 40, 80 60, 70 80, 40 80, 20 70, 20 90))" );
+  }
+
+  //------------------------------------------------
+  
+  public void testAreaSimple() {
+    checkHullByArea("MULTIPOINT ((10 10), (90 10), (30 70), (70 70), (50 60))", 
+       .5, "POLYGON ((30 70, 70 70, 90 10, 50 60, 10 10, 30 70))" );
+  }
+
+  public void testAreaZero() {
+    checkHullByArea("MULTIPOINT ((10 10), (90 10), (70 70), (50 60), (50 90), (40 70), (30 30))", 
+       0, "POLYGON ((10 10, 40 70, 50 90, 70 70, 90 10, 50 60, 30 30, 10 10))" );
+  }
+
+  public void testAreaConvex() {
+    checkHullByArea("MULTIPOINT ((10 10), (90 10), (70 70), (50 60), (50 90), (40 70), (30 30))", 
+       1, "POLYGON ((10 10, 40 70, 50 90, 70 70, 90 10, 10 10))" );
+  }
+
+  //==========================================================================
+  
+  private void checkHullByLengthFactor(String wkt, double threshold, String wktExpected) {
+    Geometry geom = read(wkt);
+    Geometry actual = ConcaveHull.concaveHullByLengthFactor(geom, threshold);
+    Geometry expected = read(wktExpected);
+    checkEqual(expected, actual);
+  }
+  
+  private void checkHullByLengthFactorXYZ(String wkt, double threshold, String wktExpected) {
+    Geometry geom = read(wkt);
+    Geometry actual = ConcaveHull.concaveHullByLengthFactor(geom, threshold);
+    Geometry expected = read(wktExpected);
+    checkEqualXYZ(expected, actual);
+  }
+  
+  private void checkHullByLength(String wkt, double threshold, String wktExpected) {
+    Geometry geom = read(wkt);
+    Geometry actual = ConcaveHull.concaveHullByLength(geom, threshold);
+    Geometry expected = read(wktExpected);
+    checkEqual(expected, actual);
+  }
+  
+  private void checkHullWithHolesByLength(String wkt, double threshold, String wktExpected) {
+    Geometry geom = read(wkt);
+    Geometry actual = ConcaveHull.concaveHullByLength(geom, threshold, true);
+    Geometry expected = read(wktExpected);
+    checkEqual(expected, actual);
+  }
+  
+  private void checkHullByArea(String wkt, double threshold, String wktExpected) {
+    Geometry geom = read(wkt);
+    Geometry actual = ConcaveHull.concaveHullByArea(geom, threshold);
+    Geometry expected = read(wktExpected);
+    checkEqual(expected, actual);
+  }
+}
diff --git a/modules/lab/src/main/java/org/locationtech/jts/hull/ConcaveHull.java b/modules/lab/src/main/java/org/locationtech/jts/hull/ConcaveHull.java
deleted file mode 100644
index da239803b6..0000000000
--- a/modules/lab/src/main/java/org/locationtech/jts/hull/ConcaveHull.java
+++ /dev/null
@@ -1,53 +0,0 @@
-/*
- * Copyright (c) 2016 Vivid Solutions.
- *
- * 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.hull;
-
-import java.util.List;
-
-import org.locationtech.jts.geom.Geometry;
-import org.locationtech.jts.triangulate.DelaunayTriangulationBuilder;
-import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
-import org.locationtech.jts.triangulate.quadedge.QuadEdgeTriangle;
-
-public class ConcaveHull {
-  
-  private Geometry geom;
-  private double tolerance;
-
-  public ConcaveHull(Geometry geom, double tolerance) {
-    this.geom = geom;
-    this.tolerance = tolerance;
-  }
-  
-  public Geometry getResult() {
-    QuadEdgeSubdivision subdiv = buildDelaunay();
-    List tris = extractTriangles(subdiv);
-    Geometry hull = computeHull(tris);
-    return hull;
-  }
-
-  private List extractTriangles(QuadEdgeSubdivision subdiv) {
-    List qeTris = QuadEdgeTriangle.createOn(subdiv);
-    return qeTris;
-  }
-
-  private Geometry computeHull(List tris) {
-    return null;
-    
-  }
-
-  private QuadEdgeSubdivision buildDelaunay() {
-    DelaunayTriangulationBuilder builder = new DelaunayTriangulationBuilder();
-    builder.setSites(geom);
-    return builder.getSubdivision();
-  }
-}
diff --git a/modules/lab/src/test/java/org/locationtech/jts/hull/ConcaveHullTest.java b/modules/lab/src/test/java/org/locationtech/jts/hull/ConcaveHullTest.java
deleted file mode 100644
index 3c8adfbb0e..0000000000
--- a/modules/lab/src/test/java/org/locationtech/jts/hull/ConcaveHullTest.java
+++ /dev/null
@@ -1,45 +0,0 @@
-/*
- * Copyright (c) 2016 Vivid Solutions.
- *
- * 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.hull;
-
-import junit.textui.TestRunner;
-
-import org.locationtech.jts.geom.Geometry;
-
-import test.jts.GeometryTestCase;
-
-public class ConcaveHullTest extends GeometryTestCase {
-
-  public static void main(String args[]) {
-    TestRunner.run(ConcaveHullTest.class);
-  }
-  
-  public ConcaveHullTest(String name) {
-    super(name);
-  }
-
-  public void testSimple() {
-    checkHull(
-        "POLYGON ((100 200, 200 180, 300 200, 200 190, 100 200))",
-        150,
-        "POLYGON ((100 200, 200 180, 300 200, 200 190, 100 200))"
-        );
-  }
-  
-  private void checkHull(String inputWKT, double tolerance, String expectedWKT) {
-    Geometry input = read(inputWKT);
-    Geometry expected = read(expectedWKT);
-    ConcaveHull hull = new ConcaveHull(input, tolerance);
-    Geometry actual = hull.getResult();
-    //checkEqual(expected, actual);
-  }
-}