查找位于点云凸壳内的点。

Question

在这个帖子中，提出了一种方法用于屏蔽位于凸包内的点，例如：

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

以下是翻译的结果：

使用这些遮罩点集，看起来像下面这样。但我正在寻找一种使用凹壳（类似于蓝色点所示）的方法。

我找到了这个帖子，其中提供了一些凹边框的功能，但我还不确定它是否适用于我的情况。有人有建议吗？

- azerila

---

你的问题没有明确定义，看起来你正在寻找欧拉路径，其中有(n-1)!个，最大面积？ - Bob

@QuangHoang，凸包是有的，但他正在寻求某种凹壳（尽管下面的图像似乎是凸包，所以我并不确定）。 - CJR

啊，我明白了。我相信所给的点表示一个（凹）多边形，并按照这个顺序排列。那就是 OP 所询问的内容。 - Quang Hoang

是的，蓝色点确实代表一个凹多边形。@QuangHoang - azerila

@Warren，我原本想说“凸包”，但问题表述得更加通用可能对其他人也有帮助。不过你所说的正是我要找的。 - azerila

Answer 1

您所提到的第一个线程中的方法可以使用alpha-shape（有时称为凹壳）概念来采用凹形情况，这也是您所提到的第二个参考答案建议的方法。

Alpha-shape是Delaunay三角剖分的三角形子集，其中每个三角形都满足外接半径条件。以下代码是从我的以前的回答修改而来，用于计算alpha-shape中Delaunay三角形的集合。一旦计算出Delaunay三角形和alpha-shape掩码，就可以采用您所引用的快速方法来应用于alpha-shape，我将在下面解释。

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://dev59.com/GlUL5IYBdhLWcg3wUWit#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-shape三角形中。下面的函数可以实现这个功能：

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-shape不会遵循多边形从您的点。增加alpha参数也无济于事，因为它会在其他地方切割凹角。如果您添加更密集的采样点，则此alpha-shape方法非常适合您的任务。如果不是，则下面我提出另一种解决方案。

由于原始多边形不适用于alpha-shape方法，因此您需要实现一个函数，该函数返回一个给定多边形内的点是否存在。以下函数基于累积内/外角度实现这样的算法（请参见这里进行解释）。

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)

- Iddo Hanniel