从 Delaunay 三角剖分计算 alpha 形状的边界多边形

Posted 2023-02-19

【中文标题】从 Delaunay 三角剖分计算 alpha 形状的边界多边形

【英文标题】：Calculate bounding polygon of alpha shape from the Delaunay triangulation

【发布时间】：2014-05-29 04:26:15

【问题描述】：

给定平面上的一组点，对于给定的正数 alpha，alpha 形状的概念是通过找到 Delaunay 三角剖分并删除至少一条边长度超过 alpha 的任何三角形来定义的。这是一个使用 d3 的示例：

http://bl.ocks.org/gka/1552725

问题是，当有数千个点时，简单地绘制所有内部三角形对于交互式可视化来说太慢了，所以我只想找到边界多边形。这不是那么简单，因为从该示例中您可以看到，有时可能有两个这样的多边形。

为简化起见，假设已经执行了一些聚类，以确保每个三角剖分都有一个唯一的边界多边形。找到这个边界多边形的最佳方法是什么？特别是，边缘的顺序必须一致，并且必须支持“孔”的可能性（想想圆环或甜甜圈形状——这可以在 GeoJSON 中表达）。

【参考方案1】：

我知道这是一个延迟的答案，但由于各种原因，此处发布的方法对我不起作用。

提到的 Alphashape 包通常很好，但它的缺点是它使用 Shapely 的cascade_union 和三角测量方法为您提供最终的凹壳，这对于大型数据集和低 alpha 值（高精度）来说非常慢。在这种情况下，Iddo Hanniel 发布的代码（以及 Harold 的修订版）将在内部保留大量边缘，消除它们的唯一方法是使用上述 cascade_union 和 Shapely 的三角测量。

我通常使用复杂的表格，下面的代码运行良好且速度很快（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

【讨论】：

您的解决方案是我在不使用 shapely 的情况下获得边缘最快的解决方案（我在尝试运行它时不断出错），itertools 组合对我来说很慢，所以我最终用 numpy 切片和时间减少了近 50%

【参考方案2】：

这里有一些 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 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

如果您使用以下测试代码运行它，您将获得 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()

【讨论】：

感谢您的代码。如果我们的 alpha 形状返回不连贯的区域，我们如何将边缘分解为单独的形状/多边形？

在另一个答案中，我添加了缝合边缘的代码。看看这个答案的结尾***.com/a/50714300/9702190，我相信它会做你想要的。

【参考方案3】：

现在有一个python包alphashape非常好用，可以通过pip或conda安装。

主函数的输入与@Iddo Hanniel 给出的答案相似，调整第二个位置参数将为您提供所需的轮廓。 或者，您可以将第二个位置参数留空，该函数将为您优化该参数，为您提供最佳凹壳。注意，如果让函数优化值，计算时间会大大增加。

【参考方案4】：

稍微修改 Hanniel 对 3d 点案例（四面体）的回答。

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

【讨论】：

你能用这个绘制一个形状吗？这个算法给我带来了麻烦。

【参考方案5】：

在@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

【参考方案6】：

Alpha 形状定义为边缘不超过 Alpha 的 delaunay 三角剖分。首先删除所有内部三角形，然后删除所有超过 alpha 的边。

【讨论】：

这是错误的。看我回答中的图片。有许多边界边比内部边长。

您的回答似乎表明您可以开始反复删除最长的 Delaunay 边，直到剩下一堆多边形。那是行不通的。 “问号”形状在其边界周围有很多比大多数边缘都长的边缘，但那些不应被删除。此外，形状内部的边缘比大多数边缘都短，应该被删除。 - 也许你的回答是想说一些不同的东西？您可以添加更多解释。

【参考方案7】：

事实证明，TopoJSON 有一个合并算法来执行这个任务：https://github.com/mbostock/topojson/wiki/API-Reference#merge

甚至还有一个例子展示了它的实际效果：http://bl.ocks.org/mbostock/9927735

就我而言，生成 TopoJSON 数据对我来说很容易，这个库函数完美地为我完成了任务。

【参考方案8】：

创建一个图，其中节点对应于 Delaunay 三角形，其中两个三角形之间存在图边当且仅当它们共享两个顶点。

计算图的连通分量。

对于每个连通分量，找出所有相邻节点少于三个的节点（即度数为 0、1 或 2 的节点）。这些对应于边界三角形。我们将不与另一个三角形共享的边界三角形的边定义为该边界三角形的边界边。

例如，我在您的示例“问号”Delaunay 三角剖分中突出显示了这些边界三角形：

根据定义，每个边界三角形最多与另外两个边界三角形相邻。边界三角形的边界边形成循环。您可以简单地遍历这些循环来确定边界的多边形形状。如果您在实现中牢记它们，这也适用于带孔的多边形。

【讨论】：

删除所有共享边并重新连接剩余边以形成边界多边形不是更容易吗？

@juniper- 这与我所描述的有何不同？请记住，我描述的方法允许算法保留边界的拓扑结构（例如，气泡字母 O 有两个边界，一个在另一个内部）。

困难的部分是以非凸的方式执行第 1 步。