Source code for alphashape.alphashape
"""
Tools for working with alpha shapes.
"""
__all__ = ['alphashape']
import itertools
import logging
from shapely.ops import unary_union, polygonize
from shapely.geometry import MultiPoint, MultiLineString
from scipy.spatial import Delaunay
import numpy as np
from typing import Union, Tuple, List
try:
import geopandas
USE_GP = True
except ImportError:
USE_GP = False
[docs]def circumcenter(points: Union[List[Tuple[float]], np.ndarray]) -> np.ndarray:
"""
Calculate the circumcenter of a set of points in barycentric coordinates.
Args:
points: An `N`x`K` array of points which define an (`N`-1) simplex in K
dimensional space. `N` and `K` must satisfy 1 <= `N` <= `K` and
`K` >= 1.
Returns:
The circumcenter of a set of points in barycentric coordinates.
"""
points = np.asarray(points)
num_rows, num_columns = points.shape
A = np.bmat([[2 * np.dot(points, points.T),
np.ones((num_rows, 1))],
[np.ones((1, num_rows)), np.zeros((1, 1))]])
b = np.hstack((np.sum(points * points, axis=1),
np.ones((1))))
return np.linalg.solve(A, b)[:-1]
[docs]def circumradius(points: Union[List[Tuple[float]], np.ndarray]) -> float:
"""
Calculte the circumradius of a given set of points.
Args:
points: An `N`x`K` array of points which define an (`N`-1) simplex in K
dimensional space. `N` and `K` must satisfy 1 <= `N` <= `K` and
`K` >= 1.
Returns:
The circumradius of a given set of points.
"""
points = np.asarray(points)
return np.linalg.norm(points[0, :] - np.dot(circumcenter(points), points))
[docs]def alphasimplices(points: Union[List[Tuple[float]], np.ndarray]) -> \
Union[List[Tuple[float]], np.ndarray]:
"""
Returns an iterator of simplices and their circumradii of the given set of
points.
Args:
points: An `N`x`M` array of points.
Yields:
A simplex, and its circumradius as a tuple.
"""
coords = np.asarray(points)
tri = Delaunay(coords)
for simplex in tri.simplices:
simplex_points = coords[simplex]
try:
yield simplex, circumradius(simplex_points)
except np.linalg.LinAlgError:
logging.warn('Singular matrix. Likely caused by all points '
'lying in an N-1 space.')
[docs]def alphashape(points: Union[List[Tuple[float]], np.ndarray],
alpha: Union[None, float] = None):
"""
Compute the alpha shape (concave hull) of a set of points. If the number
of points in the input is three or less, the convex hull is returned to the
user. For two points, the convex hull collapses to a `LineString`; for one
point, a `Point`.
Args:
points (list or ``shapely.geometry.MultiPoint`` or \
``geopandas.GeoDataFrame``): an iterable container of points
alpha (float): alpha value
Returns:
``shapely.geometry.Polygon`` or ``shapely.geometry.LineString`` or
``shapely.geometry.Point`` or ``geopandas.GeoDataFrame``: \
the resulting geometry
"""
# If given a geodataframe, extract the geometry
if USE_GP and isinstance(points, geopandas.GeoDataFrame):
crs = points.crs
points = points['geometry']
else:
crs = None
# If given a triangle for input, or an alpha value of zero or less,
# return the convex hull.
if len(points) < 4 or (alpha is not None and not callable(
alpha) and alpha <= 0):
if not isinstance(points, MultiPoint):
points = MultiPoint(list(points))
result = points.convex_hull
if crs:
gdf = geopandas.GeoDataFrame(geopandas.GeoSeries(result)).rename(
columns={0: 'geometry'}).set_geometry('geometry')
gdf.crs = crs
return gdf
else:
return result
# Determine alpha parameter if one is not given
if alpha is None:
try:
from optimizealpha import optimizealpha
except ImportError:
from .optimizealpha import optimizealpha
alpha = optimizealpha(points)
# Convert the points to a numpy array
if USE_GP and isinstance(points, geopandas.geoseries.GeoSeries):
coords = np.array([point.coords[0] for point in points])
else:
coords = np.array(points)
# Create a set to hold unique edges of simplices that pass the radius
# filtering
edges = set()
# Create a set to hold unique edges of perimeter simplices.
# Whenever a simplex is found that passes the radius filter, its edges
# will be inspected to see if they already exist in the `edges` set. If an
# edge does not already exist there, it will be added to both the `edges`
# set and the `permimeter_edges` set. If it does already exist there, it
# will be removed from the `perimeter_edges` set if found there. This is
# taking advantage of the property of perimeter edges that each edge can
# only exist once.
perimeter_edges = set()
for point_indices, circumradius in alphasimplices(coords):
if callable(alpha):
resolved_alpha = alpha(point_indices, circumradius)
else:
resolved_alpha = alpha
# Radius filter
if circumradius < 1.0 / resolved_alpha:
for edge in itertools.combinations(
point_indices, r=coords.shape[-1]):
if all([e not in edges for e in itertools.combinations(
edge, r=len(edge))]):
edges.add(edge)
perimeter_edges.add(edge)
else:
perimeter_edges -= set(itertools.combinations(
edge, r=len(edge)))
if coords.shape[-1] > 3:
return perimeter_edges
elif coords.shape[-1] == 3:
import trimesh
result = trimesh.Trimesh(vertices=coords, faces=list(perimeter_edges))
trimesh.repair.fix_normals(result)
return result
# Create the resulting polygon from the edge points
m = MultiLineString([coords[np.array(edge)] for edge in perimeter_edges])
triangles = list(polygonize(m))
result = unary_union(triangles)
# Convert to pandas geodataframe object if that is what was an input
if crs:
gdf = geopandas.GeoDataFrame(geopandas.GeoSeries(result)).rename(
columns={0: 'geometry'}).set_geometry('geometry')
gdf.crs = crs
return gdf
else:
return result