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