geos::algorithm::hull::ConcaveHull::concaveHullByLength segfault

[orca-hydromancer]

In libgeos v3.11.0 the following code results in segfault.

auto g = geos::io::WKTReader().read("MULTIPOINT( 1139294.6389832512941211 8201313.5346954688429832, 1139360.8549531854223460 8201271.1898052766919136, 1139497.5995843114797026 8201199.9955425458028913, 1139567.7837303513661027 8201163.3485335074365139, 1139635.3942210066597909 8201119.9025274068117142 )");
geos::algorithm::hull::ConcaveHull::concaveHullByLength( g.get(), 240, false );

Callstack:

geos::triangulate::tri::Tri::getIndex(const geos::triangulate::tri::Tri * const this, const geos::triangulate::tri::Tri * tri) (/externals/geos/src/triangulate/tri/Tri.cpp:312)
geos::algorithm::hull::HullTriangulation::nextBorderTri(geos::algorithm::hull::HullTri * triStart) (/externals/geos/src/algorithm/hull/HullTriangulation.cpp:161)
geos::algorithm::hull::HullTriangulation::traceBoundary(geos::triangulate::tri::TriList<geos::algorithm::hull::HullTri> & triList) (/externals/geos/src/algorithm/hull/HullTriangulation.cpp:129)
geos::algorithm::hull::HullTriangulation::traceBoundaryPolygon(geos::triangulate::tri::TriList<geos::algorithm::hull::HullTri> & triList, const geos::geom::GeometryFactory * factory) (/externals/geos/src/algorithm/hull/HullTriangulation.cpp:105)
geos::algorithm::hull::ConcaveHull::toGeometry(geos::algorithm::hull::ConcaveHull * const this, geos::triangulate::tri::TriList<geos::algorithm::hull::HullTri> & triList, const geos::geom::GeometryFactory * factory) (/externals/geos/src/algorithm/hull/ConcaveHull.cpp:375)
geos::algorithm::hull::ConcaveHull::getHull(geos::algorithm::hull::ConcaveHull * const this) (/externals/geos/src/algorithm/hull/ConcaveHull.cpp:152)
geos::algorithm::hull::ConcaveHull::concaveHullByLength(const geos::geom::Geometry * geom, double maxLength, bool holesAllowed) (/externals/geos/src/algorithm/hull/ConcaveHull.cpp:80)

This is what the multipoint looks like

[orca-hydromancer]

geosop -a "MULTIPOINT( 1139294.6389832512941211 8201313.5346954688429832, 1139360.8549531854223460 8201271.1898052766919136, 1139497.5995843114797026 8201199.9955425458028913, 1139567.7837303513661027 8201163.3485335074365139, 1139635.3942210066597909 8201119.9025274068117142 )" concaveHull -f wkt

[dr-jts]

Fails in JTS too. It looks like the underlying Delaunay Triangulation is incorrect (it contains only 2 triangles and is not convex).

[dr-jts]

Here's a simpler reproducer:

MULTIPOINT ((0 194), (66 151), (203 80), (273 43), (340 0))

It looks like the triangular "frame" that is created around the points to contain the evolving triangulation is too small, interfering with correct construction of the DT. Increasing the size of the frame fixes the problem. The open question is: how much of an increase is needed to handle all cases?

The criteria for determining a safe frame size is that none of the frame vertices lie in a circumcircle of a triangle in the resulting DT. This means that the frame size is influenced by the nature of the triangulated points. In this case the triangulated points are almost collinear, which means the DT triangles are quite flat, and thus have large circumcircles. This is the source of the problem.

Most (all?) Incremental DT algorithms gloss over this point. So this is a yet another example of how published algorithms often trivialize details which are key for producing a fully robust algorithm.