from shapely.ops import cascaded_union, polygonize
from scipy.spatial import Delaunay
import shapely.geometry as geometry
import numpy as np
import math
def alpha_shape(points, alpha):
"""
Compute the alpha shape (concave hull) of a set
of points.
@param points: Iterable container of points.
@param alpha: alpha value to influence the
gooeyness of the border. Smaller numbers
don't fall inward as much as larger numbers.
Too large, and you lose everything!
"""
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull
coords = np.array([point.coords[0] for point in points])
tri = Delaunay(coords)
triangles = coords[tri.vertices]
a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5
b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5
c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5
s = ( a + b + c ) / 2.0
areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
circums = a * b * c / (4.0 * areas)
filtered = triangles[circums < (1.0 / alpha)]
edge1 = filtered[:,(0,1)]
edge2 = filtered[:,(1,2)]
edge3 = filtered[:,(2,0)]
edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
