用一个凹多边形区域对一组点进行三角剖分

Question

设置

假设在一个凸包内有一些节点，假设该区域包含一个或多个凹陷的区域：

蓝点是点，黑线表示该域。假设这些点被保存为长度为n的2D数组points，其中 n 是点对数。

然后我们使用类似于scipy.spatial中的Delaunay方法来进行三角剖分：

作为您所看到的，一个人可能会经历穿过域的三角形的创建。

问题

什么是消除跨越域外的任何三角形的良好算法方法？理想情况下，但不一定需要，简单形状仍然保留域形状（即，在删除三角形的地方没有主要间隙）。

---

由于我的问题似乎一直有相当多的活动，我想跟进一下我当前正在使用的应用程序。

假设您已定义了边界，则可以使用射线投射算法来确定多边形是否在域内。

要做到这一点：

1. Take the centroid of each polygon as C_i = (x_i,y_i).
2. Then, imagine a line L = [C_i,(+inf,y_i)]: that is, a line that spans east past the end of your domain.
3. For each boundary segment s_i in boundary S, check for intersections with L. If yes, add +1 to an internal counter intersection_count; else, add nothing.

4. After the count of all intersections between L and s_i for i=1..N are calculated:

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

如果您的边界没有明确定义，我建议将形状“映射”到一个布尔数组上，并使用邻居跟踪算法来定义它。请注意，这种方法假定了一个坚实的领域，如果领域中有“孔洞”，则需要使用更复杂的算法。

- Daniel R. Livingston

---

这几乎不是一个Python问题。 - Mauricio

尝试在BOOST中使用polygon包中的算法。 - Prune

你熟悉alpha hulls / alpha shapes吗？https://en.wikipedia.org/wiki/Alpha_shape - Rethunk

@Rethunk 没有，但感谢您提供的链接/信息！ - Daniel R. Livingston

我会在任何凸多边形中使用Marching Square算法。它旨在轻松查找边界。 - asdfasdf

Answer 1

我这里有一些Python代码，可以实现你想要的功能。

首先，构建alpha形状（参见我的上一个答案）：

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)

现在，在计算了点集的alpha形状的外边界的边缘之后，以下函数将确定一个点(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

在上图所示的输入中（摘自我之前的答案），调用is_inside(1.5, 0.0, points, edges)将返回True，而is_inside(1.5, 3.0, points, edges)将返回False。 请注意，上面的is_inside函数不处理退化情况。我添加了两个断言来检测这种情况（您可以定义适合您应用程序的任何epsilon值）。在许多应用程序中，这已经足够了，但如果没有遇到这些边缘情况，它们需要单独处理。 例如，请参见此处有关实现几何算法时鲁棒性和精度问题的信息。

- Iddo Hanniel

Answer 2

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

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

这可能还需要插入和修剪凹壳。

- Aki Suihkonen

Answer 3

由于我的问题似乎仍然得到了相当多的关注，我想跟进一下我目前正在使用的应用程序。

假设您已经定义了边界，您可以使用射线投射算法来确定多边形是否在域内。

要做到这一点：

1. Take the centroid of each polygon as C_i = (x_i,y_i).
2. Then, imagine a line L = [C_i,(+inf,y_i)]: that is, a line that spans east past the end of your domain.
3. For each boundary segment s_i in boundary S, check for intersections with L. If yes, add +1 to an internal counter intersection_count; else, add nothing.

4. After the count of all intersections between L and s_i for i=1..N are calculated:

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

如果您的边界没有明确定义，我建议您将其映射到一个布尔数组上，并使用邻居跟踪算法neighbor tracing algorithm来定义它。请注意，这种方法假定领域是实心的，对于具有“孔”的领域，您需要使用更复杂的算法。

- Daniel R. Livingston

Answer 4

你可以尝试使用 Sloan 算法或 CGAL 库等算法实现 约束性 Delaunay 算法。

[1] 一种暴力约束 Delaunay 三角剖分方法？

- Micromega

Answer 5

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

要了解更多信息，请查看以下帖子：https://gis.stackexchange.com/a/352442 请注意，shapely中也实现了三角剖分，即使是MultiPoint对象也可以实现。

我已经使用过它了，性能非常惊人，代码只有五行。

- magraf

Answer 6

使用此算法计算三角形的重心并检查其是否在多边形内。

- abenci

---

不确定那会起作用；考虑多边形上可能包含重心但不包含三角形其余部分的“尖锐部位”。 - meowgoesthedog

最终的DT不包含这样一个三角形，即使中间的镶嵌可以包含它。 - Aki Suihkonen