python - 查找位于点云凹包中的点

Question

在this thread建议使用一种方法来屏蔽位于凸包中的点，例如:

x = np.array([0,1,2,3,4,4, 4, 6, 6, 5, 5, 1])
y = np.array([0,1,2,3,4,3, 3.5, 3, 2, 0, 3, 0])

xx = np.linspace(np.min(x)-1, np.max(x)+1, 40)
yy = np.linspace(np.min(y)-1, np.max(y)+1, 40)
xx, yy = np.meshgrid(xx, yy)

plt.scatter(x, y, s=50)
plt.scatter(xx, yy, s=10)

def in_hull(p, hull):
    from scipy.spatial import Delaunay
    if not isinstance(hull, Delaunay):
        hull = Delaunay(hull)
hull1 = np.stack((x,y)).T
p1 = np.stack((xx.ravel(),yy.ravel())).T
cond = in_hull(p1, hull1)
p2 = p1[cond,:]
plt.scatter(x, y)
plt.scatter(p2[:,0],p2[:,1], s=10)
    return hull.find_simplex(p)>=0

掩码点集如下所示。但是，我正在寻找一种使用凹壳的方法(类似于蓝点所暗示的)

我找到了 this thread这暗示了凹形边框的一些功能，但我不确定它是否适用于我的情况。有人有什么建议吗？

Answer 1

您引用的第一个线程中的方法可以采用 alpha 形状(有时称为凹包)概念的凹壳，这是您第二个引用的答案所暗示的。

alpha 形状是 Delaunay 三角剖分的三角形子集，其中每个三角形都满足外接半径条件。以下代码修改为from my previous answer计算 alpha 形状中的 Delaunay 三角形集。一旦计算出 Delaunay 三角剖分和 alpha 形状掩码，您引用的快速方法就可以应用于 alpha 形状，我将在下面解释。

def circ_radius(p0,p1,p2):
    """
    Vectorized computation of triangle circumscribing radii.
    See for example https://www.cuemath.com/jee/circumcircle-formulae-trigonometry/
    """
    a = p1-p0
    b = p2-p0

    norm_a = np.linalg.norm(a, axis=1)
    norm_b = np.linalg.norm(b, axis=1)
    norm_a_b = np.linalg.norm(a-b, axis=1)
    cross_a_b = np.cross(a,b)  # 2 * area of triangles
    return (norm_a*norm_b*norm_a_b) / np.abs(2.0*cross_a_b)


def alpha_shape_delaunay_mask(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set of points and return the Delaunay triangulation and a boolean
    mask for any triangle in the triangulation whether it belongs to the alpha shape.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :return: Delaunay triangulation dt and boolean array is_in_shape, so that dt.simplices[is_in_alpha] contains
    only the triangles that belong to the alpha shape.
    """
    # Modified and vectorized from:
    # https://stackoverflow.com/questions/50549128/boundary-enclosing-a-given-set-of-points/50714300#50714300

    assert points.shape[0] > 3, "Need at least four points"
    dt = Delaunay(points)

    p0 = points[dt.simplices[:,0],:]
    p1 = points[dt.simplices[:,1],:]
    p2 = points[dt.simplices[:,2],:]

    rads = circ_radius(p0, p1, p2)

    is_in_shape = (rads < alpha)

    return dt, is_in_shape

然后可以调整您第一个引用中的方法，不仅检查该点是否位于 Delaunay 三角形之一(在这种情况下它位于凸包中)，而且还检查它是否位于 alpha-形状三角形。以下函数执行此操作:

def in_alpha_shape(p, dt, is_in_alpha):
    simplex_ids = dt.find_simplex(p)
    res = np.full(p.shape[0], False)
    res[simplex_ids >= 0] = is_in_alpha[simplex_ids[simplex_ids >= 0]]  # simplex should be in dt _and_ in alpha
    return res

此方法非常快，因为它依赖于 Delaunay find_simplex() 函数的高效搜索实现。

使用下面的代码在您帖子中的示例数据点上运行它(使用 alpha=2)给出下图中的结果，我相信这不是您想要的...

points = np.vstack([x, y]).T

alpha = 2.
dt, is_in_alpha = alpha_shape_delaunay_mask(points, alpha)

p1 = np.stack((xx.ravel(),yy.ravel())).T
cond = in_alpha_shape(p1, dt, is_in_alpha)
p2 = p1[cond,:]

plt.figure()
plt.scatter(x, y)
plt.scatter(p2[:,0],p2[:,1], s=10)

上述结果的原因是，由于输入点之间存在较大间隙，因此数据的 alpha 形状不遵循点的多边形。增加 alpha 参数也无济于事，因为它会在其他地方切割凹角。如果您添加更密集的样本点，那么这种 alpha 形状方法可能非常适合您的任务。如果没有，那么下面我提出另一种解决方案。

由于您的原始多边形不适合 alpha-shape 方法，您需要实现一个函数来返回点是否在给定的多边形内。以下函数基于累积内/外角度实现了这样的算法(有关解释，请参阅 here)。

def points_in_polygon(pts, polygon):
    """
    Returns if the points are inside the given polygon,

    Implemented with angle accumulation.
    see: 
https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm

    :param np.ndarray pts: 2d points
    :param np.ndarray polygon: 2d polygon
    :return: Returns if the points are inside the given polygon, array[i] == True means pts[i] is inside the polygon.
    """
    polygon = np.vstack((polygon, polygon[0, :]))  # close the polygon (if already closed shouldn't hurt)
    sum_angles = np.zeros([len(pts), ])
    for i in range(len(polygon) - 1):
        v1 = polygon[i, :] - pts
        norm_v1 = np.linalg.norm(v1, axis=1, keepdims=True)
        norm_v1[norm_v1 == 0.0] = 1.0  # prevent divide-by-zero nans
        v1 = v1 / norm_v1
        v2 = polygon[i + 1, :] - pts
        norm_v2 = np.linalg.norm(v2, axis=1, keepdims=True)
        norm_v2[norm_v2 == 0.0] = 1.0  # prevent divide-by-zero nans
        v2 = v2 / norm_v2
        dot_prods = np.sum(v1 * v2, axis=1)
        cross_prods = np.cross(v1, v2)
        angs = np.arccos(np.clip(dot_prods, -1, 1))
        angs = np.sign(cross_prods) * angs
        sum_angles += angs

    sum_degrees = np.rad2deg(sum_angles)
    # In most cases abs(sum_degrees) should be close to 360 (inside) or to 0 (outside).
    # However, in end cases, points that are on the polygon can be less than 360, so I allow a generous margin..
    return abs(sum_degrees) > 90.0

使用下面的代码调用它会产生下图，我相信这正是您要找的。

points = np.vstack([x, y]).T
p1 = np.vstack([xx.ravel(), yy.ravel()]).T
cond = points_in_polygon(p1, points)
p2 = p1[cond,:]

plt.figure()
plt.scatter(x, y)
plt.plot(x, y)
plt.scatter(p2[:,0],p2[:,1], s=10)