Get alpha-shape (concave hull) of a set of points using Julia

I want to compute the alpha-shape (or even only the concave hull) of a set of points using Julia. In other questions they have solved this problem in python by using Delaunay tesselation Boundary enclosing a given set of points .

This package in Julia can get a Delaunay tesselation https://github.com/JuliaGeometry/VoronoiDelaunay.jl (though I am not sure if it is updated for julia v0.7). I am wondering if there is an implementation already for julia v0.7 that can get eh alpha-shape, or even just the concave hull of a set of points.

Alternatively, is there a way to efficiently call python (scipy.spatial.Delaunay) to do the job?

VoronoiDelaunay.jl works with Julia 1.0 and 1.1. It should also work with Julia 0.7.

VoronoiDelaunay.jl has some numerical restrictions, i.e. (1.0+eps(), 2.0-eps()), on coordinates so you may need to re-scale your data points.

To create a DelaunayTesselation with your own point type, make sure your type is a subtype of AbstractPoint2D, that is <: AbstractPoint2D, and defines getx, and gety methods.

The following example code, I believe, finds what you call the Concave Hull of a set of points using DelaunayTesselation and plots the result. It basically uses the same algorithm in this answer. You may easily edit the code to get the alpha shape.

I did not wrap some code snippets into a function. If you need high performance, please do so. I used === while checking for equality of points which actually checks if two points are the same object (i.e. address in memory). If you somehow end up in a code which breaks this part, you can extend == and use it instead of ===.

using Random, VoronoiDelaunay, Plots

import Base.==

struct MyEdge{T<:AbstractPoint2D}
    _a::T
    _b::T
end

==(e1::MyEdge{T}, e2::MyEdge{T}) where {T<:AbstractPoint2D} = ((e1._a === e2._a) && (e1._b === e2._b)) || ((e1._b === e2._a) && (e2._b === e1._a))

###==(p1::T, p2::T) where {T<:AbstractPoint2D} = (getx(p1) == getx(p2)) && (gety(p1) == gety(p2))

### Create a Delaunay tesselation from random points
tess = DelaunayTessellation2D(46)

for _ in 1:23
    push!(tess, Point2D(rand()+1, rand()+1))
end

edges = MyEdge[]

function add_edge!(edges, edge)
    i = findfirst(e -> e == edge, edges)

    if isnothing(i) # not found
        push!(edges, edge)
    else # found so not outer, remove this edge
        deleteat!(edges, i) 
    end
end

for trig in tess
    a, b, c = geta(trig), getb(trig), getc(trig)
    add_edge!(edges, MyEdge(b, c))
    add_edge!(edges, MyEdge(a, b))
    add_edge!(edges, MyEdge(a, c))
end

### PLOT
x, y = Float64[], Float64[] # outer edges
for edge in edges
    push!(x, getx(edge._a))
    push!(x, getx(edge._b))
    push!(x, NaN)
    push!(y, gety(edge._a))
    push!(y, gety(edge._b))
    push!(y, NaN)
end

xall, yall = getplotxy(delaunayedges(tess)) # all the edges

plot(xall, yall, color=:blue, fmt=:svg, size=(400,400))
plot!(x, y, color=:red, linewidth=3, opacity=0.5)