CODEWITHSUNDEEP
juliacomputational-geometryconcave-hullpycall
RM-
RM-

Reputation: 1008

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?

Upvotes: 0

Views: 935

Answers (1)

hckr
hckr

Reputation: 5583

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)

Upvotes: 2

Related Questions