# Computing Concave hull/alpha shape for a point cloud

**URL:** <https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746>\
**Category:** Visualization\
**Tags:** gis, geometry, archgdal, libgeos, convex-hull\
**Created:** [June 14, 2022, 1:17pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746 "2022-06-14T13:17:01Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![arsh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arsh/32/9073_2.png) [@arsh](https://discourse.julialang.org/u/arsh)\
**Post date:** [June 14, 2022, 1:17pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/1 "2022-06-14T13:17:01Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/a/6/a6d216718beeaf3f10d6a95688f2f16a82bc4da6.png)

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](https://github.com/lstagner/ConcaveHull.jl) - stack overflow error, looks unmaintained  
[AlphaShapes.jl Documentation · AlphaShapes.jl](https://harveydevereux.github.io/AlphaShapes.jl/dev/#Index) - output across different parameters remain similar to convex hull.

Any help would be really appreciated.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 14, 2022, 1:25pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/2 "2022-06-14T13:25:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![maxfreu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxfreu/32/17468_2.png) [@maxfreu](https://discourse.julialang.org/u/maxfreu)\
**Post date:** [June 14, 2022, 1:29pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/3 "2022-06-14T13:29:13Z")

</div>

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.

```python
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](https://libgeos.org/doxygen/geos__c_8h.html#aaceadb7351dda799a583edca329931ae) is the documentation of the C API, which you need to wrap. [Here](https://github.com/JuliaGeo/LibGEOS.jl/blob/430c17f5c8f04d1d28387d45f2fb6f4e95e6bcc9/src/geos_functions.jl#L482-L488) and [here](https://github.com/JuliaGeo/LibGEOS.jl/blob/841339d3ee04330cbc7673a29a544ea42f2e877b/src/libgeos_api.jl#L587-L592) is how it is done for the convex hull, so you can copy that more or less 1:1.

---

<div class="post-metadata">

**Author:** ![arsh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arsh/32/9073_2.png) [@arsh](https://discourse.julialang.org/u/arsh)\
**Post date:** [June 14, 2022, 3:56pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/4 "2022-06-14T15:56:54Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![maxfreu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxfreu/32/17468_2.png) [@maxfreu](https://discourse.julialang.org/u/maxfreu)\
**Post date:** [June 14, 2022, 4:04pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/5 "2022-06-14T16:04:20Z")

</div>

Can you post the code you wrote?

---

<div class="post-metadata">

**Author:** ![arsh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arsh/32/9073_2.png) [@arsh](https://discourse.julialang.org/u/arsh)\
**Post date:** [June 14, 2022, 4:07pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/6 "2022-06-14T16:07:03Z")

</div>

I tried something like.

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

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

```

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 14, 2022, 5:10pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/7 "2022-06-14T17:10:42Z")

</div>

> [@maxfreu](#):
>
> 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.

GMT.jl wraps and exports it

```julia
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

```

---

<div class="post-metadata">

**Author:** ![arsh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arsh/32/9073_2.png) [@arsh](https://discourse.julialang.org/u/arsh)\
**Post date:** [June 14, 2022, 5:41pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/8 "2022-06-14T17:41:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 14, 2022, 5:44pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/9 "2022-06-14T17:44:21Z")

</div>

> [@arsh](#):
>
> 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

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 14, 2022, 9:39pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/10 "2022-06-14T21:39:06Z")

</div>

Bleeding edge(s) news. Needs GMT.jl master You will need a GDAL built with geos3.11 that is still marked as [beta1](https://libgeos.org/usage/download/) 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)

```julia
using GMT

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

```

 ![GMTjl_tmp](https://global.discourse-cdn.com/julialang/original/3X/6/9/69836ed8d383d0ed9e8fae09aa8c61c925ffef36.png)

---

<div class="post-metadata">

**Author:** ![maxfreu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxfreu/32/17468_2.png) [@maxfreu](https://discourse.julialang.org/u/maxfreu)\
**Post date:** [June 15, 2022, 8:35pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/11 "2022-06-15T20:35:53Z")

</div>

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

> [@arsh](#):
>
> I tried to warp the concave hull function.

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 😉 But the GMT approach sounds also nice!

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [June 15, 2022, 9:29pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/12 "2022-06-15T21:29:57Z")

</div>

> [@maxfreu](#):
>
> @joa-quim this sounds like gdal builds on geos, is that correct?

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

---

<div class="post-metadata">

**Author:** ![arsh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arsh/32/9073_2.png) [@arsh](https://discourse.julialang.org/u/arsh)\
**Post date:** [June 16, 2022, 7:16am UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/13 "2022-06-16T07:16:27Z")

</div>

Hmm yeah.  
I pushed the changes here [GitHub - JuliaGeo/LibGEOS.jl at sv-concavehull](https://github.com/JuliaGeo/LibGEOS.jl/tree/sv-concavehull)

It fails for the following example too

```julia
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

```

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [June 16, 2022, 8:52am UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/14 "2022-06-16T08:52:16Z")

</div>

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](https://github.com/libgeos/geos/blob/main/NEWS.md). 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.

---

<div class="post-metadata">

**Author:** ![maxfreu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxfreu/32/17468_2.png) [@maxfreu](https://discourse.julialang.org/u/maxfreu)\
**Post date:** [June 16, 2022, 2:56pm UTC](https://discourse.julialang.org/t/computing-concave-hull-alpha-shape-for-a-point-cloud/82746/15 "2022-06-16T14:56:14Z")

</div>

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