用凹域对一组点进行三角剖分

Posted 2023-02-14

【中文标题】用凹域对一组点进行三角剖分

【英文标题】：Triangulate a set of points with a concave domain

【发布时间】：2017-12-26 20:55:56

【问题描述】：

设置

给定凸包内的一组节点，假设域包含一个或多个凹区域：

其中蓝点是点，黑线表示域。假设点被保存为长度为n 的二维数组points，其中n 是点对的数量。

然后让我们使用 scipy.spatial 中的 Delaunay 方法对这些点进行三角测量：

如您所见，您可能会体验到创建穿过域的三角形。

问题

有什么好的算法可以去除任何跨越域的三角形？理想情况下但不一定，单纯形边缘仍保留域形状（即，三角形被删除的地方没有大间隙）。

---

由于我的问题似乎继续获得大量活动，我想跟进我目前正在使用的应用程序。

假设您定义了边界，您可以使用ray casting algorithm 来确定多边形是否在域内。

为此：

以每个多边形的质心为C_i = (x_i,y_i)。 然后，想象一条线 L = [C_i,(+inf,y_i)]：即向东跨越您的域末端的线。 对于边界S 中的每个边界段s_i，检查与L 的交点。如果是，将 +1 添加到内部计数器intersection_count；否则，什么都不加。

计算L和s_i for i=1..N之间的所有交集的计数后：

if intersection_count % 2 == 0:
    return True # triangle outside convex hull
else:
    return False # triangle inside convex hull

如果您的边界没有明确定义，我发现将形状“映射”到布尔数组并使用neighbor tracing algorithm 来定义它很有帮助。请注意，这种方法假设一个实体域，您需要对其中有“洞”的域使用更复杂的算法。

【问题讨论】：

这不是一个 python 问题

在 BOOST 中尝试 polygon 包中的算法。

您熟悉 alpha hulls / alpha 形状吗？ en.wikipedia.org/wiki/Alpha_shape

@Rethunk 一点也不，但谢谢你的链接/信息！

我会对任何凸多边形使用行进平方算法。它可以轻松找到边界

【参考方案1】：

这里有一些 Python 代码可以满足您的需求。

首先，构建 alpha 形状（参见 my previous answer）：

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

要计算 alpha 形状的外边界的边缘，请使用以下示例调用：

edges = alpha_shape(points, alpha=alpha_value, only_outer=True)

现在，在计算points的alpha-shape的外边界edges之后，下面的函数将判断一个点(x,y)是否在外边界内。

def is_inside(x, y, points, edges, eps=1.0e-10):
    intersection_counter = 0
    for i, j in edges:
        assert abs((points[i,1]-y)*(points[j,1]-y)) > eps, 'Need to handle these end cases separately'
        y_in_edge_domain = ((points[i,1]-y)*(points[j,1]-y) < 0)
        if y_in_edge_domain:
            upper_ind, lower_ind = (i,j) if (points[i,1]-y) > 0 else (j,i)
            upper_x = points[upper_ind, 0] 
            upper_y = points[upper_ind, 1]
            lower_x = points[lower_ind, 0] 
            lower_y = points[lower_ind, 1]

            # is_left_turn predicate is evaluated with: sign(cross_product(upper-lower, p-lower))
            cross_prod = (upper_x - lower_x)*(y-lower_y) - (upper_y - lower_y)*(x-lower_x)
            assert abs(cross_prod) > eps, 'Need to handle these end cases separately'
            point_is_left_of_segment = (cross_prod > 0.0)
            if point_is_left_of_segment:
                intersection_counter = intersection_counter + 1
    return (intersection_counter % 2) != 0

在上图中显示的输入中（取自我之前的answer），调用is_inside(1.5, 0.0, points, edges) 将返回True，而is_inside(1.5, 3.0, points, edges) 将返回False。

请注意，上面的is_inside 函数不处理退化情况。我添加了两个断言来检测这种情况（您可以定义任何适合您的应用程序的 epsilon 值）。在许多应用程序中，这已经足够了，但如果没有，并且您遇到这些最终情况，则需要单独处理它们。 例如，请参阅here，了解实施几何算法时的稳健性和精度问题。

【参考方案2】：

其中一种经典 DT 算法首先生成一个边界三角形，然后添加所有按 x 排序的新三角形，然后修剪掉所有在超三角形中有顶点的三角形。

至少从提供的图像中，我们可以得出启发式方法，也可以剪除一些在凹壳上具有所有顶点的三角形。如果没有证明，当它们的顶点按照定义凹壳的顺序排序时，要修剪的三角形的面积为负。

这可能还需要插入凹形外壳，并被修剪掉。

【参考方案3】：

由于我的问题似乎继续获得大量活动，我想跟进我目前正在使用的应用程序。

假设您定义了边界，您可以使用ray casting algorithm 来确定多边形是否在域内。

为此：

以每个多边形的质心为C_i = (x_i,y_i)。 然后，想象一条线 L = [C_i,(+inf,y_i)]：即向东跨越您的域末端的线。 对于边界S 中的每个边界段s_i，检查与L 的交点。如果是，将 +1 添加到内部计数器intersection_count；否则，什么都不加。

计算L和s_i for i=1..N之间的所有交集的计数后：

if intersection_count % 2 == 0:
    return True # triangle outside convex hull
else:
    return False # triangle inside convex hull

如果您的边界没有明确定义，我发现将形状“映射”到布尔数组并使用neighbor tracing algorithm 来定义它很有帮助。请注意，这种方法假设一个实体域，您需要对其中有“洞”的域使用更复杂的算法。

【参考方案4】：

您可以尝试使用 约束 delaunay 算法，例如使用 sloan 算法或 cgal 库。

[1]A Brute-Force Constrained Delaunay Triangulation?

【参考方案5】：

一个简单而优雅的方法是遍历三角形并检查它们是否在我们的domain 内。 shapely 包可以为您解决问题。

有关更多信息，请查看以下帖子：https://gis.stackexchange.com/a/352442 请注意，即使对于 MultiPoin 对象，也实现了 shapely 三角测量。

我用过，性能惊人，代码只有五行。

【参考方案6】：

使用this algorithm 计算三角形质心并检查它是否在多边形内。

【讨论】：

不确定是否可行；考虑多边形上的一个“尖位”，它可能恰好包含质心，但不包含三角形的其余部分

最终的 DT 不包含这样的三角形，即使中间镶嵌可以。