您引用的第一个线程中的方法可以使用 alpha 形状（有时称为凹壳）概念应用于凹面案例，这就是您第二个参考的答案所暗示的。
alpha 形状是 Delaunay 三角剖分的三角形子集，其中每个三角形都满足外接半径条件。以下代码是根据我之前的答案修改的，以计算 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
这种方法非常快，因为它依赖于 Delaunayfind_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 方法，因此您需要一个函数的实现，该函数返回点是否在给定的多边形内。以下函数实现了基于累积内/外角的算法（请参阅此处以获取解释）。
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)
