首先,构建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)

edges 的alpha形状的外边界 points (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 参见,例如, here

至少从提供的图像中,我们可以推导出一些启发式方法,即修剪出凹壳上具有所有顶点的一些三角形。在没有证明的情况下,当顶点按与定义凹壳相同的顺序排序时,要修剪的三角形具有负区域。

这可能需要插入凹面外壳,并将其修剪掉。

由于我的问题似乎继续得到相当数量的活动,我想继续使用我目前使用的应用程序。

假设定义了边界,可以使用 ray casting algorithm

为此:

1. 将每个多边形的质心作为 C_i = (x_i,y_i)
2. 然后,想象一条线 L = [C_i,(+inf,y_i)]
3. 对于每个边界段 s_i 在边界内 S ,检查与 L intersection_count ; 否则,不添加任何内容。

4. 在计算之间的所有交点之后 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

受限

[1] A Brute-Force Constrained Delaunay Triangulation?

一种简单但优雅的方法是在三角形上循环,并检查它们是否在我们的范围内 domain 或者不是。这个 shapely 包裹可以帮你。

有关更多信息,请查看以下帖子: https://gis.stackexchange.com/a/352442

我使用了它,性能非常惊人,代码只有五行。

this algorithm .