给定一组平面点，对于给定的正数alpha，定义了alpha-shape的概念，通过找到Delaunay三角剖分并删除至少有一个边长超过alpha的三角形来实现。以下是使用d3的示例：
http://bl.ocks.org/gka/1552725
问题在于，当存在成千上万的点时，简单地绘制所有内部三角形对于交互式可视化来说太慢了，因此我想只找到边界多边形。这并不简单，因为从那个示例中可以看出，有时可能会有两个这样的多边形。
作为简化，假设已经执行了一些聚类，以保证每个三角剖分都有唯一的边界多边形。找到这个边界多边形的最佳方法是什么？特别是，必须一致地对边进行排序，并且它必须支持“孔”的可能性（考虑环面或甜甜圈形状-这可以用GeoJSON表示）。
这里有一些Python代码可以满足您的需求。我修改了来自这里的alpha-shape（凸包）计算，以便它不插入内部边缘（only_outer参数）。我还使它自包含，因此它不依赖外部库。
from scipy.spatial import Delaunay
import numpy as np
def alpha_shape(points, alpha, only_outer=True):
"""
Compute the alpha shape (concave hull) of a set of points.
:param points: np.array of shape (n,2) points.
:param alpha: alpha value.
:param only_outer: boolean value to specify if we keep only the outer border
or also inner edges.
:return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
the indices in the points array.
"""
assert points.shape[0] > 3, "Need at least four points"
def add_edge(edges, 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
assert (j, i) in edges, "Can't go twice over same directed edge right?"
if only_outer:
# if both neighboring triangles are in shape, it is not a boundary edge
edges.remove((j, i))
return
edges.add((i, j))
tri = Delaunay(points)
edges = set()
# Loop over triangles:
# ia, ib, ic = indices of corner points of the triangle
for ia, ib, ic in tri.simplices:
pa = points[ia]
pb = points[ib]
pc = points[ic]
# Computing radius of triangle circumcircle
# www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
s = (a + b + c) / 2.0
area = np.sqrt(s * (s - a) * (s - b) * (s - c))
circum_r = a * b * c / (4.0 * area)
if circum_r < alpha:
add_edge(edges, ia, ib)
add_edge(edges, ib, ic)
add_edge(edges, ic, ia)
return edges
from matplotlib.pyplot import *
# Constructing the input point data
np.random.seed(0)
x = 3.0 * np.random.rand(2000)
y = 2.0 * np.random.rand(2000) - 1.0
inside = ((x ** 2 + y ** 2 > 1.0) & ((x - 3) ** 2 + y ** 2 > 1.0) & ((x - 1.5) ** 2 + y ** 2 > 0.09))
points = np.vstack([x[inside], y[inside]]).T
# Computing the alpha shape
edges = alpha_shape(points, alpha=0.25, only_outer=True)
# Plotting the output
figure()
axis('equal')
plot(points[:, 0], points[:, 1], '.')
for i, j in edges:
plot(points[[i, j], 0], points[[i, j], 1])
show()
例如，在您的示例“问号”Delaunay三角剖分中，我已突出显示这些边界三角形：
根据定义，每个边界三角形最多与另外两个边界三角形相邻。边界三角形的边界边形成循环。您可以简单地遍历这些循环以确定边界的多边形形状。如果考虑到多边形内部的空洞，则此实现也适用于具有孔的多边形。
我知道这是迟来的回答，但这里发布的方法出于各种原因对我都不起作用。
提到的 Alphashape 包通常很好，但它的缺点是使用 Shapely 的 cascade_union 和三角剖分方法来生成最终的凹多边形壳体，对于大数据集和低 alpha 值（高精度）非常慢。在这种情况下，Iddo Hanniel（和 Harold 的修订版）发布的代码将保留大量内部边缘，唯一溶解它们的方法是使用前面提到的 Shapely 的cascade_union和三角剖分。
我通常处理复杂形状，下面的代码运行良好且速度快（100K 2D 点只需 2 秒）：
import numpy as np
from shapely.geometry import MultiLineString
from shapely.ops import unary_union, polygonize
from scipy.spatial import Delaunay
from collections import Counter
import itertools
def concave_hull(coords, alpha): # coords is a 2D numpy array
# i removed the Qbb option from the scipy defaults.
# it is much faster and equally precise without it.
# unless your coords are integers.
# see http://www.qhull.org/html/qh-optq.htm
tri = Delaunay(coords, qhull_options="Qc Qz Q12").vertices
ia, ib, ic = (
tri[:, 0],
tri[:, 1],
tri[:, 2],
) # indices of each of the triangles' points
pa, pb, pc = (
coords[ia],
coords[ib],
coords[ic],
) # coordinates of each of the triangles' points
a = np.sqrt((pa[:, 0] - pb[:, 0]) ** 2 + (pa[:, 1] - pb[:, 1]) ** 2)
b = np.sqrt((pb[:, 0] - pc[:, 0]) ** 2 + (pb[:, 1] - pc[:, 1]) ** 2)
c = np.sqrt((pc[:, 0] - pa[:, 0]) ** 2 + (pc[:, 1] - pa[:, 1]) ** 2)
s = (a + b + c) * 0.5 # Semi-perimeter of triangle
area = np.sqrt(
s * (s - a) * (s - b) * (s - c)
) # Area of triangle by Heron's formula
filter = (
a * b * c / (4.0 * area) < 1.0 / alpha
) # Radius Filter based on alpha value
# Filter the edges
edges = tri[filter]
# now a main difference with the aforementioned approaches is that we dont
# use a Set() because this eliminates duplicate edges. in the list below
# both (i, j) and (j, i) pairs are counted. The reasoning is that boundary
# edges appear only once while interior edges twice
edges = [
tuple(sorted(combo)) for e in edges for combo in itertools.combinations(e, 2)
]
count = Counter(edges) # count occurrences of each edge
# keep only edges that appear one time (concave hull edges)
edges = [e for e, c in count.items() if c == 1]
# these are the coordinates of the edges that comprise the concave hull
edges = [(coords[e[0]], coords[e[1]]) for e in edges]
# use this only if you need to return your hull points in "order" (i think
# its CCW)
ml = MultiLineString(edges)
poly = polygonize(ml)
hull = unary_union(list(poly))
hull_vertices = hull.exterior.coords.xy
return edges, hull_vertices
pip或conda安装。与@Iddo Hanniel所提供的答案类似，主要函数具有类似的输入，调整第二个位置参数即可获得所需轮廓。 或者，您可以将第二个位置参数留空，该函数将为您优化该参数以给出最佳凹壳。请注意，如果让函数优化该值，则计算时间会大大增加。对于三维点情况（四面体），Hanniel的答案 进行了轻微修改。
from scipy.spatial import Delaunay
import numpy as np
import math
def alpha_shape(points, alpha, only_outer=True):
"""
Compute the alpha shape (concave hull) of a set of points.
:param points: np.array of shape (n, 3) points.
:param alpha: alpha value.
:param only_outer: boolean value to specify if we keep only the outer border
or also inner edges.
:return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
the indices in the points array.
"""
assert points.shape[0] > 3, "Need at least four points"
def add_edge(edges, i, j):
"""
Add a line between the i-th and j-th points,
if not in the set already
"""
if (i, j) in edges or (j, i) in edges:
# already added
if only_outer:
# if both neighboring triangles are in shape, it's not a boundary edge
if (j, i) in edges:
edges.remove((j, i))
return
edges.add((i, j))
tri = Delaunay(points)
edges = set()
# Loop over triangles:
# ia, ib, ic, id = indices of corner points of the tetrahedron
print(tri.vertices.shape)
for ia, ib, ic, id in tri.vertices:
pa = points[ia]
pb = points[ib]
pc = points[ic]
pd = points[id]
# Computing radius of tetrahedron Circumsphere
# http://mathworld.wolfram.com/Circumsphere.html
pa2 = np.dot(pa, pa)
pb2 = np.dot(pb, pb)
pc2 = np.dot(pc, pc)
pd2 = np.dot(pd, pd)
a = np.linalg.det(np.array([np.append(pa, 1), np.append(pb, 1), np.append(pc, 1), np.append(pd, 1)]))
Dx = np.linalg.det(np.array([np.array([pa2, pa[1], pa[2], 1]),
np.array([pb2, pb[1], pb[2], 1]),
np.array([pc2, pc[1], pc[2], 1]),
np.array([pd2, pd[1], pd[2], 1])]))
Dy = - np.linalg.det(np.array([np.array([pa2, pa[0], pa[2], 1]),
np.array([pb2, pb[0], pb[2], 1]),
np.array([pc2, pc[0], pc[2], 1]),
np.array([pd2, pd[0], pd[2], 1])]))
Dz = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], 1]),
np.array([pb2, pb[0], pb[1], 1]),
np.array([pc2, pc[0], pc[1], 1]),
np.array([pd2, pd[0], pd[1], 1])]))
c = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], pa[2]]),
np.array([pb2, pb[0], pb[1], pb[2]]),
np.array([pc2, pc[0], pc[1], pc[2]]),
np.array([pd2, pd[0], pd[1], pd[2]])]))
circum_r = math.sqrt(math.pow(Dx, 2) + math.pow(Dy, 2) + math.pow(Dz, 2) - 4 * a * c) / (2 * abs(a))
if circum_r < alpha:
add_edge(edges, ia, ib)
add_edge(edges, ib, ic)
add_edge(edges, ic, id)
add_edge(edges, id, ia)
add_edge(edges, ia, ic)
add_edge(edges, ib, id)
return edges
from matplotlib.tri import Triangulation
import numpy as np
def get_boundary_rings(x, y, elements):
mpl_tri = Triangulation(x, y, elements)
idxs = np.vstack(list(np.where(mpl_tri.neighbors == -1))).T
unique_edges = list()
for i, j in idxs:
unique_edges.append((mpl_tri.triangles[i, j],
mpl_tri.triangles[i, (j+1) % 3]))
unique_edges = np.asarray(unique_edges)
ring_collection = list()
initial_idx = 0
for i in range(1, len(unique_edges)-1):
if unique_edges[i-1, 1] != unique_edges[i, 0]:
try:
idx = np.where(
unique_edges[i-1, 1] == unique_edges[i:, 0])[0][0]
unique_edges[[i, idx+i]] = unique_edges[[idx+i, i]]
except IndexError:
ring_collection.append(unique_edges[initial_idx:i, :])
initial_idx = i
continue
# if there is just one ring, the exception is never reached,
# so populate ring_collection before returning.
if len(ring_collection) == 0:
ring_collection.append(np.asarray(unique_edges))
return ring_collection
Alpha shapes被定义为没有边超过alpha的Delaunay三角剖分。首先删除所有内部三角形，然后删除所有超过alpha的边。