"""
Computation of alpha shape algorithm in 2-D based on original implementation
by Tim Kittel (@timkittel) available at:
https://github.com/timkittel/alpha-shapes
Author(s):
Dani Arribas-Bel daniel.arribas.bel@gmail.com
Levi John Wolf levi.john.wolf@gmail.com
"""
import numpy as np
import scipy.spatial as spat
from packaging.version import Version
from scipy import sparse
from ..common import HAS_JIT, jit, requires
if not HAS_JIT:
from warnings import warn
NUMBA_WARN = (
"Numba not imported, so alpha shape construction may be slower than expected."
)
try:
import shapely
assert Version(shapely.__version__) >= Version("2")
HAS_SHAPELY = True
except (ModuleNotFoundError, AssertionError):
HAS_SHAPELY = False
EPS = np.finfo(float).eps
__all__ = ["alpha_shape", "alpha_shape_auto"]
@jit(nopython=True)
def nb_dist(x, y):
"""numba implementation of distance between points `x` and `y`
Parameters
----------
x : ndarray
Coordinates of point `x`
y : ndarray
Coordinates of point `y`
Returns
-------
dist : float
Distance between `x` and `y`
Examples
--------
>>> x = np.array([0, 0])
>>> y = np.array([1, 1])
>>> dist = nb_dist(x, y)
>>> dist
1.4142135623730951
"""
sum_ = 0
for x_i, y_i in zip(x, y): # noqa: B905
sum_ += (x_i - y_i) ** 2
dist = np.sqrt(sum_)
return dist
@jit(nopython=True)
def r_circumcircle_triangle_single(a, b, c):
"""Computation of the circumcircle of a single triangle
Parameters
----------
a : ndarray
(2,) Array with coordinates of vertex `a` of the triangle
b : ndarray
(2,) Array with coordinates of vertex `b` of the triangle
c : ndarray
(2,) Array with coordinates of vertex `c` of the triangle
Returns
-------
r : float
Circumcircle of the triangle
Notes
-----
Source for equations:
> https://www.mathopenref.com/trianglecircumcircle.html
[Last accessed July 11th. 2018]
Examples
--------
>>> a = np.array([0, 0])
>>> b = np.array([0.5, 0])
>>> c = np.array([0.25, 0.25])
>>> r = r_circumcircle_triangle_single(a, b, c)
>>> r
0.2500000000000001
"""
ab = nb_dist(a, b)
bc = nb_dist(b, c)
ca = nb_dist(c, a)
num = ab * bc * ca
den = np.sqrt((ab + bc + ca) * (bc + ca - ab) * (ca + ab - bc) * (ab + bc - ca))
if den == 0:
return np.array([ab, bc, ca]).max() / 2.0
else:
return num / den
@jit(nopython=True)
def r_circumcircle_triangle(a_s, b_s, c_s):
"""Computation of circumcircles for a series of triangles
Parameters
----------
a_s : ndarray
(N, 2) array with coordinates of vertices `a` of the triangles
b_s : ndarray
(N, 2) array with coordinates of vertices `b` of the triangles
c_s : ndarray
(N, 2) array with coordinates of vertices `c` of the triangles
Returns
-------
radii : ndarray
(N,) array with circumcircles for every triangle
Examples
--------
>>> a_s = np.array([[0, 0], [2, 1], [3, 2]])
>>> b_s = np.array([[1, 0], [5, 1], [2, 4]])
>>> c_s = np.array([[0, 7], [1, 3], [4, 2]])
>>> rs = r_circumcircle_triangle(a_s, b_s, c_s)
>>> rs
array([3.53553391, 2.5 , 1.58113883])
"""
len_a = len(a_s)
r2 = np.zeros((len_a,))
for i in range(len_a):
r2[i] = r_circumcircle_triangle_single(a_s[i], b_s[i], c_s[i])
return r2
@jit(nopython=True)
def get_faces(triangle):
"""Extract faces from a single triangle
Parameters
----------
triangles : ndarray
(3,) array with the vertex indices for a triangle
Returns
-------
faces : ndarray
(3, 2) array with a row for each face containing the indices of the two
points that make up the face
Examples
--------
>>> triangle = np.array([3, 1, 4], dtype=np.int32)
>>> faces = get_faces(triangle)
>>> faces
array([[3., 1.],
[1., 4.],
[4., 3.]])
"""
faces = np.zeros((3, 2))
for i, (i0, i1) in enumerate([(0, 1), (1, 2), (2, 0)]):
faces[i] = triangle[i0], triangle[i1]
return faces
@jit(nopython=True)
def build_faces(faces, triangles_is, num_triangles, num_faces_single):
"""Build facing triangles
Parameters
----------
faces : ndarray
(num_triangles * num_faces_single, 2) array of zeroes in int form
triangles_is : ndarray
(D, 3) array, where D is the number of Delaunay triangles, with the
vertex indices for each triangle
num_triangles : int
Number of triangles
num_faces_single : int
Number of faces a triangle has (i.e. 3)
Returns
-------
faces : ndarray
Two dimensional array with a row for every facing segment containing
the indices of the coordinate points
Examples
--------
>>> import scipy.spatial as spat
>>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
>>> triangulation = spat.Delaunay(pts)
>>> triangulation.simplices
array([[3, 1, 4],
[1, 2, 4],
[2, 1, 0]], dtype=int32)
>>> num_faces_single = 3
>>> num_triangles = triangulation.simplices.shape[0]
>>> num_faces = num_triangles * num_faces_single
>>> faces = np.zeros((num_faces, 2), dtype=np.int_)
>>> mask = np.ones((num_faces,), dtype=np.bool_)
>>> faces = build_faces(
... faces, triangulation.simplices, num_triangles, num_faces_single
... )
>>> faces
array([[3, 1],
[1, 4],
[4, 3],
[1, 2],
[2, 4],
[4, 1],
[2, 1],
[1, 0],
[0, 2]])
"""
for i in range(num_triangles):
from_i = num_faces_single * i
to_i = num_faces_single * (i + 1)
faces[from_i:to_i] = get_faces(triangles_is[i])
return faces
@jit(nopython=True)
def nb_mask_faces(mask, faces):
"""Run over each row in `faces`, if the face in the following row is the
same, then mark both as False on `mask`
Parameters
----------
mask : ndarray
One-dimensional boolean array set to True with as many observations as
rows in `faces`
faces : ndarray
Sorted sequence of faces for all triangles (ie. triangles split by each
segment)
Returns
-------
masked : ndarray
Sequence of outward-facing faces
Examples
--------
>>> import numpy as np
>>> faces = np.array(
... [
... [0, 1], [0, 2], [1, 2], [1, 2], [1, 3], [1, 4], [1, 4], [2, 4], [3, 4]
... ]
... )
>>> mask = np.ones((faces.shape[0], ), dtype=np.bool_)
>>> masked = nb_mask_faces(mask, faces)
>>> masked
array([[0, 1],
[0, 2],
[1, 3],
[2, 4],
[3, 4]])
"""
for k in range(faces.shape[0] - 1):
if mask[k] and np.all(faces[k] == faces[k + 1]):
mask[k] = False
mask[k + 1] = False
return faces[mask]
def get_single_faces(triangles_is):
"""Extract outward facing edges from collection of triangles
Parameters
----------
triangles_is : ndarray
(D, 3) array, where D is the number of Delaunay triangles, with the
vertex indices for each triangle
Returns
-------
single_faces : ndarray
Examples
--------
>>> import scipy.spatial as spat
>>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
>>> alpha = 0.33
>>> triangulation = spat.Delaunay(pts)
>>> triangulation.simplices
array([[3, 1, 4],
[1, 2, 4],
[2, 1, 0]], dtype=int32)
>>> get_single_faces(triangulation.simplices)
array([[0, 1],
[0, 2],
[1, 3],
[2, 4],
[3, 4]])
"""
num_faces_single = 3
num_triangles = triangles_is.shape[0]
num_faces = num_triangles * num_faces_single
faces = np.zeros((num_faces, 2), dtype=np.int_)
mask = np.ones((num_faces,), dtype=np.bool_)
faces = build_faces(faces, triangles_is, num_triangles, num_faces_single)
orderlist = [f"x{i}" for i in range(faces.shape[1])]
dtype_list = [(el, faces.dtype.str) for el in orderlist]
# Arranging each face so smallest vertex is first
faces.sort(axis=1)
# Arranging faces in ascending way
faces.view(dtype_list).sort(axis=0)
# Masking
single_faces = nb_mask_faces(mask, faces)
return single_faces
@requires("geopandas", "shapely")
def _alpha_geoms(alpha, triangles, radii, xys):
"""Generate alpha-shape polygon(s) from `alpha` value, vertices of
`triangles`, the `radii` for all points, and the points themselves
Parameters
----------
alpha : float
Alpha value to delineate the alpha-shape
triangles : ndarray
(D, 3) array, where D is the number of Delaunay triangles, with the
vertex indices for each triangle
radii : ndarray
(N,) array with circumcircles for every triangle
xys : ndarray
(N, 2) array with one point per row and coordinates structured as X and Y
Returns
-------
geoms : GeoSeries
Polygon(s) resulting from the alpha shape algorithm, in a GeoSeries.
The output is a GeoSeries even if only a single polygon is returned. There is
no CRS included in the returned GeoSeries.
Examples
--------
>>> import scipy.spatial as spat
>>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
>>> alpha = 0.33
>>> triangulation = spat.Delaunay(pts)
>>> triangles = pts[triangulation.simplices]
>>> triangles
array([[[6, 7],
[3, 5],
[9, 3]],
<BLANKLINE>
[[3, 5],
[4, 1],
[9, 3]],
<BLANKLINE>
[[4, 1],
[3, 5],
[0, 1]]])
>>> a_pts = triangles[:, 0, :]
>>> b_pts = triangles[:, 1, :]
>>> c_pts = triangles[:, 2, :]
>>> radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
>>> geoms = alpha_geoms(alpha, triangulation.simplices, radii, pts)
>>> geoms
0 POLYGON ((0.00000 1.00000, 3.00000 5.00000, 4....
dtype: geometry
"""
from geopandas import GeoSeries
from shapely.geometry import LineString
from shapely.ops import polygonize
triangles_reduced = triangles[radii < 1 / alpha]
outer_triangulation = get_single_faces(triangles_reduced)
face_pts = xys[outer_triangulation]
geoms = GeoSeries(list(polygonize(list(map(LineString, face_pts)))))
return geoms
[docs]
@requires("geopandas", "shapely")
def alpha_shape(xys, alpha):
"""Alpha-shape delineation (Edelsbrunner, Kirkpatrick & Seidel, 1983)
from a collection of points
Parameters
----------
xys : ndarray
(N, 2) array with one point per row and coordinates structured as X and Y
alpha : float
Alpha value to delineate the alpha-shape
Returns
-------
shapes : GeoSeries
Polygon(s) resulting from the alpha shape algorithm. The GeoSeries
object remains so even if only a single polygon is returned. There is
no CRS included in the object. Note that the returned shape(s) may
have holes, as per the definition of the shape in Edselbrunner et al. (1983)
Examples
--------
>>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
>>> alpha = 0.1
>>> poly = alpha_shape(pts, alpha)
>>> poly
0 POLYGON ((0.00000 1.00000, 3.00000 5.00000, 6....
dtype: geometry
>>> poly.centroid
0 POINT (4.69048 3.45238)
dtype: geometry
References
----------
Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
a set of points in the plane. IEEE Transactions on information theory,
29(4), 551-559.
"""
if not HAS_JIT:
warn(NUMBA_WARN, stacklevel=2)
if xys.shape[0] < 4:
from shapely import geometry as geom
from shapely import ops
return ops.unary_union([geom.Point(xy) for xy in xys]).convex_hull.buffer(0)
triangulation = spat.Delaunay(xys)
triangles = xys[triangulation.simplices]
a_pts = triangles[:, 0, :]
b_pts = triangles[:, 1, :]
c_pts = triangles[:, 2, :]
radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
del triangles, a_pts, b_pts, c_pts
geoms = _alpha_geoms(alpha, triangulation.simplices, radii, xys)
geoms = _filter_holes(geoms, xys)
return geoms
def _valid_hull(geoms, points):
"""Sanity check within ``alpha_shape_auto()`` to verify the generated alpha
shape actually contains the original set of points (xys).
Parameters
----------
geoms : GeoSeries
See alpha_geoms()
points : list
xys parameter cast as shapely.geometry.Point objects
Returns
-------
flag : bool
Valid hull for alpha shape [True] or not [False]
"""
# if there is not exactly one polygon
if geoms.shape[0] != 1:
return False
# if any (xys) points do not intersect the polygon
if HAS_SHAPELY:
return shapely.intersects(geoms[0], points).all()
else:
return all(point.intersects(geoms[0]) for point in points)
[docs]
@requires("geopandas", "shapely")
def alpha_shape_auto(
xys, step=1, verbose=False, return_radius=False, return_circles=False
):
"""Computation of alpha-shape delineation with automated selection of alpha.
This method uses the algorithm proposed by Edelsbrunner, Kirkpatrick &
Seidel (1983) to return the tightest polygon that contains all points in
`xys`. The algorithm ranks every point based on its radius and iterates
over each point, checking whether the maximum alpha that would keep the
point and all the other ones in the set with smaller radii results in a
single polygon. If that is the case, it moves to the next point;
otherwise, it retains the previous alpha value and returns the polygon
as `shapely` geometry. Note that this geometry may have holes.
Parameters
----------
xys : ndarray
Nx2 array with one point per row and coordinates structured as X and Y
step : int
[Optional. Default=1] Number of points in `xys` to jump ahead after
checking whether the largest possible alpha that includes the point and
all the other ones with smaller radii
verbose : Boolean
[Optional. Default=False] If True, it prints alpha values
being tried at every step.
Returns
-------
poly : shapely.Polygon
Tightest alpha-shape polygon containing all points in `xys`
Examples
--------
>>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
>>> poly = alpha_shape_auto(pts)
>>> poly.bounds
(0.0, 1.0, 9.0, 7.0)
>>> poly.centroid.x, poly.centroid.y
(4.690476190476191, 3.4523809523809526)
References
----------
Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
a set of points in the plane. IEEE Transactions on information theory,
29(4), 551-559.
"""
if not HAS_JIT:
warn(NUMBA_WARN, stacklevel=2)
from shapely import geometry as geom
if return_circles:
return_radius = True
if xys.shape[0] < 4:
if xys.shape[0] == 3:
multipoint = geom.MultiPoint(xys)
alpha_shape = multipoint.convex_hull
else:
alpha_shape = geom.Polygon([])
if xys.shape[0] == 1:
if return_radius:
if return_circles:
out = [alpha_shape, 0, alpha_shape]
return alpha_shape, 0
return alpha_shape
elif xys.shape[0] == 2:
if return_radius:
r = spat.distance.euclidean(xys[0], xys[1]) / 2
if return_circles:
circle = _construct_centers(xys[0], xys[1], r)
return [alpha_shape, r, circle]
return [alpha_shape, r]
return alpha_shape
elif return_radius: # this handles xys.shape[0] == 3
radius = r_circumcircle_triangle_single(xys[0], xys[1], xys[2])
if return_circles:
circles = construct_bounding_circles(alpha_shape, radius)
return [alpha_shape, radius, circles]
return [alpha_shape, radius]
return alpha_shape
triangulation = spat.Delaunay(xys)
triangles = xys[triangulation.simplices]
a_pts = triangles[:, 0, :]
b_pts = triangles[:, 1, :]
c_pts = triangles[:, 2, :]
radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
radii[np.isnan(radii)] = 0 # "Line" triangles to be kept for sure
del triangles, a_pts, b_pts, c_pts
radii_sorted_i = radii.argsort()
triangles = triangulation.simplices[radii_sorted_i][::-1]
radii = radii[radii_sorted_i][::-1]
geoms_prev = _alpha_geoms((1 / radii.max()) - EPS, triangles, radii, xys)
points = shapely.points(xys) if HAS_SHAPELY else [geom.Point(pnt) for pnt in xys]
if verbose:
print("Step set to %i" % step)
for i in range(0, len(radii), step):
radi = radii[i]
alpha = (1 / radi) - EPS
if verbose:
print(f"{(i + 1) / radii.shape[0]:.2f}% | Trying a = {alpha:f}")
geoms = _alpha_geoms(alpha, triangles, radii, xys)
if _valid_hull(geoms, points):
geoms_prev = geoms
radi_prev = radi
else:
break
if verbose:
print(geoms_prev.shape)
if return_radius:
out = [geoms_prev[0], radi_prev]
if return_circles:
out.append(construct_bounding_circles(out[0], radi_prev))
return out
# Return a shapely polygon
return geoms_prev[0]
def construct_bounding_circles(alpha_shape, radius):
"""Construct the bounding circles for an alpha shape, given the radius
computed from the `alpha_shape_auto` method.
Parameters
----------
alpha_shape : shapely.Polygon
An alpha-hull with the input radius.
radius : float
The radius of the input alpha_shape.
Returns
-------
center : numpy.ndarray of shape (n,2)
The centers of the circles defining the alpha_shape.
"""
coordinates = list(alpha_shape.boundary.coords)
n_coordinates = len(coordinates)
centers = []
for i in range(n_coordinates - 1):
a, b = coordinates[i], coordinates[i + 1]
centers.append(_construct_centers(a, b, radius))
return centers
@jit(nopython=True)
def _construct_centers(a, b, radius):
midpoint_x = (a[0] + b[0]) * 0.5
midpoint_y = (a[1] + b[1]) * 0.5
d = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5
if b[0] - a[0] == 0:
m = np.inf
axis_rotation = np.pi / 2
else:
m = (b[1] - a[1]) / (b[0] - a[0])
axis_rotation = np.arctan(m)
# altitude is perpendicular bisector of AB
interior_angle = np.arccos(0.5 * d / radius)
chord = np.sin(interior_angle) * radius
dx = chord * np.sin(axis_rotation)
dy = chord * np.cos(axis_rotation)
up_x = midpoint_x - dx
up_y = midpoint_y + dy
down_x = midpoint_x + dx
down_y = midpoint_y - dy
# sign gives us direction of point, since
# shapely shapes are clockwise-defined
sign = np.sign((b[0] - a[0]) * (up_y - a[1]) - (b[1] - a[1]) * (up_x - a[0]))
if sign == 1:
return up_x, up_y
else:
return down_x, down_y
def _filter_holes(geoms, points): # noqa: ARG001
"""
Filter hole polygons using a computational geometry solution
"""
import geopandas
if (geoms.interiors.apply(len) > 0).any():
from shapely.geometry import Polygon
# Extract the "shell", or outer ring of the polygon.
shells = geoms.exterior.apply(Polygon)
# Compute which original geometries are within each shell, self-inclusive
if Version(geopandas.__version__) >= Version("0.13"):
inside, outside = shells.sindex.query(geoms, predicate="within")
else:
inside, outside = shells.sindex.query_bulk(geoms, predicate="within")
# Now, create the sparse matrix relating the inner geom (rows)
# to the outer shell (cols) and take the sum.
# A z-order of 1 means the polygon is only inside if its own exterior.
# This means it's not a hole.
# A z-order of 2 means the polygon is inside of exactly one other exterior.
# Because the hull generation method is restricted to be planar, this
# means the polygon is a hole.
# In general, an even z-order means that the polygon is always exactly
# matched to one exterior, plus some number of intermediate
# exterior-hole pairs. Therefore, the polygon is a hole.
# In general, an odd z-order means that there is an uneven number of exteriors.
# This means the polygon is not a hole.
zorder = sparse.csc_matrix((np.ones_like(inside), (inside, outside))).sum(
axis=1
)
zorder = np.asarray(zorder).flatten()
# Keep only the odd z-orders
to_include = (zorder % 2).astype(bool)
geoms = geoms[to_include]
return geoms
if __name__ == "__main__":
import time
import geopandas as gpd
import matplotlib.pyplot as plt
plt.close("all")
xys = np.random.random((1000, 2))
t0 = time.time()
geoms = alpha_shape_auto(xys, 1)
t1 = time.time()
print("%.2f Seconds to run algorithm" % (t1 - t0))
f, ax = plt.subplots(1)
gpd.GeoDataFrame({"geometry": [geoms]}).plot(ax=ax, color="orange", alpha=0.5)
ax.scatter(xys[:, 0], xys[:, 1], s=0.1)
plt.show()