Computing Concave hull/alpha shape for a point cloud

I have a point cloud and I want to find it’s outline(also includes holes). Eg:

A delaunay triangulation of such point cloud might look like
image

I have tried to use existing tools in the ecosystem but most of them don’t seem to give expected results.

GitHub - lstagner/ConcaveHull.jl: Julia package for calculating 2D concave/convex hulls - stack overflow error, looks unmaintained
AlphaShapes.jl Documentation · AlphaShapes.jl - output across different parameters remain similar to convex hull.

Any help would be really appreciated.

1 Like

If I remember correctly the Voronoi/Delaunay tessellation has information that could be used to identify this boundary. I don’t know if libraries provide this information though.

I’m familiar only with the 2D case. There you have at least two options, both of which involve some coding. Option 1 is to write the algorithm youself, which is not that hard and below is a python example, option 2 is to use LibGEOS.jl to calculate the concave hull. Currently the concave hull function is not wrapped, so you have to write the wrapper. But there are plenty of examples in the source code on how to do it.

Python example, assembled from some Stackoverflow some years ago, unfortunately not well documented, so I have to guess like you. But I remember using it, so it works in principle. It uses shapely for the triangulation I think.

def alpha_shape(points, alpha, buffer=0):

    if len(points) < 20:
        return geometry.MultiPoint(points).convex_hull, None

    def add_edge(edges, edge_points, coords, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(coords[ [i, j] ])

    coords = np.array(points)

    try:
        tri = Delaunay(coords)
    except:
        return None, None

    edges = set()
    edge_points = []
    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
        a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
        b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
        c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
        # Semiperimeter of triangle
        s = (a + b + c)/2.0
        try:
            # Area of triangle by Heron's formula
            area = math.sqrt(s*(s-a)*(s-b)*(s-c))
        except ValueError:
            area = 0
        # print(ia, ib, ic, area)
        circum_r = a*b*c/(4.0*area+10**-5)
        # print(circum_r)
        # Here's the radius filter.
        if circum_r < alpha:
            add_edge(edges, edge_points, coords, ia, ib)
            add_edge(edges, edge_points, coords, ib, ic)
            add_edge(edges, edge_points, coords, ic, ia)
    try:
        m = geometry.MultiLineString(edge_points)
        triangles = list(polygonize(m))
        concave_hull = cascaded_union(triangles)
        concave_hull = concave_hull.buffer(buffer)

    except:
        return None, None

    # Lets check, if the resulting polygon contains at least 90% of the points.
    # If not, we return the convex hull.

    points_total = len(points)
    points_inside = 0

    for p in shapely.geometry.MultiPoint(points):
        points_inside += concave_hull.contains(p)

    if points_inside/points_total<0.9:
        return geometry.MultiPoint(points).convex_hull, None
    elif not concave_hull.is_empty:
        return concave_hull, edge_points
    else:
        return None, None

To get you started with the wrapper, here is the documentation of the C API, which you need to wrap. Here and here is how it is done for the convex hull, so you can copy that more or less 1:1.

Thanks alot.
I tried both methods.
The python example script tends to work well, for a circumradius that filters the outline well, but also filters some internal triangles introducing unwanted holes in the polygon later. So there’s not much we can do on that part.

I tried to warp the concave hull function.
But LibGEOS seems to be quite tricky for first time use. I have created a coordseq(MultiPoint is what I need?)
I am running into segfaults if I try to pass that directly to the concave/convex hull method.
What do you think I might be missing here?

Can you post the code you wrote?

I tried something like.

a = LibGEOS.createCoordSeq(Vector{Float64}[[1, 2, 3], [4, 5, 6]])

hull = LibGEOS.convexhull(a)
# or
hull = LibGEOS.concavehull(a, 0.2, false)

GMT.jl wraps and exports it

help?> convexhull
search: convexhull

  convexhull(geom; gdataset=false)

  Parameters
  ––––––––––––

    •  geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

    •  gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

  A new geometry object is created and returned containing the convex hull of the geometry on which the method is invoked.

  Returns
  –––––––––

  A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

Thanks. I am actually looking for concave hull. GMT.jl doesn’t seem to have implemented that.

I just showed it has (well, only showed the doc string but…)

EDIT: Sorry, I didn’t read it right. You said concave. But if it’s wrapped by GDAL it can be done in GMT.jl

Bleeding edge(s) news. Needs GMT.jl master You will need a GDAL built with geos3.11 that is still marked as beta1 but since last December.

Then one can:
(Note: we are about to release GMT (the lib, not the wrapper) in a couple of days. The Windows version will come with a GDAL latest with geos3.11)

using GMT

mat = rand(50,2);
D = concavehull(mat, 0.5);
imshow(D, plot=(data=mat, marker=:star, ms=0.2, mc=:red))

2 Likes

@joa-quim this sounds like gdal builds on geos, is that correct?

If you want to go further, I can maybe guide you a bit. The code you posted only shows how you call the functions, which seemingly fails. Can you post the wrapper code you wrote? Then I can walk you to a working version and you can make a PR in the end :wink: But the GMT approach sounds also nice!

Yes, it wraps geos, PROJ (that wraps part of Gepgraphiclib). That’s my favorite way off accessing those libraries and many of them’s functions are wrapped in GMT.jl

Hmm yeah.
I pushed the changes here GitHub - JuliaGeo/LibGEOS.jl at sv-concavehull

It fails for the following example too

julia> input_ = LibGEOS._readgeom(
               "MULTIPOINT (130 240, 130 240, 130 240, 570 240, 570 240, 570 240, 650 240)",
           )
Ptr{Nothing} @0x000055a9ef17c710

julia> output = LibGEOS.concavehull(input_)
ERROR: could not load symbol "GEOSConcaveHull_r":
/home/user/.julia/artifacts/a795a1524a0f75ee944e3a050eeed19a4ce01458/lib/libgeos_c.so: undefined symbol: GEOSConcaveHull_r
Stacktrace:
 [1] GEOSConcaveHull_r
   @ ~/.julia/dev/LibGEOS/src/libgeos_api.jl:596 [inlined]
 [2] concavehull (repeats 2 times)
   @ ~/.julia/dev/LibGEOS/src/geos_functions.jl:501 [inlined]
 [3] top-level scope
   @ REPL[84]:1

LibGEOS.jl comes with GEOS 3.10, the latest release. GEOS 3.11 is in beta and includes the concave hull, see geos/NEWS.md at main · libgeos/geos · GitHub. I don’t think it’s a good idea to start shipping LibGEOS.jl with beta releases so we’ll have to wait for the release first.

Ahh I completely missed that point! Of course one can’t wrap sth that is not there. Let’s wait for the release then.