给定平面中的一组点，对于给定的正数α，α形状的概念通过找到Delaunay三角剖分并且删除至少一个边长度超过α的任何三角形来定义。这是使用d3的示例：
http://bl.ocks.org/gka/1552725
问题是，当有数千个点时，简单地绘制所有内部三角形对于交互式可视化来说太慢了，所以我想找到边界多边形。这不是那么简单，因为从这个例子可以看出，有时可能会有两个这样的多边形。
作为简化，假设已经执行了一些聚类，因此保证每个三角测量的唯一边界多边形。找到这个边界多边形的最佳方法是什么？特别是，边缘必须一致地排序，它必须支持“洞”的可能性（想想圆环或圆环形状 - 这在GeoJSON中是可表达的）。
举个例子，我在你的例子“问号”Delaunay三角剖分中突出显示了这些边界三角形：
根据定义，每个边界三角形最多与其他两个边界三角形相邻。边界三角形的边界边界形成循环。您可以简单地遍历这些循环以确定边界的多边形形状。如果您在实现中记住它们，这也适用于带孔的多边形。
这是一些Python代码，可以满足您的需求。我修改了here的alpha形状（凹壳）计算，使其不插入内边（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's 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.vertices:
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
如果您使用以下测试代码运行它，您将获得this figure：
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()
事实证明，TopoJSON有一个合并算法，它只执行这个任务：https://github.com/mbostock/topojson/wiki/API-Reference#merge
甚至有一个例子显示它在行动：http://bl.ocks.org/mbostock/9927735
在我的例子中，我很容易生成TopoJSON数据，这个库函数完美地完成了我的任务。
在@Timothy的答案的基础上，我使用以下算法来计算Delaunay三角剖分的边界环。
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形状被定义为没有边缘超过alpha的delaunay三角剖分。首先删除所有内部三角形，然后删除所有超出alpha的边。