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)