The method from the first thread you reference can be adopted to the concave case using the alpha-shape (sometimes called the concave hull) concept, which is what the answer from your second reference suggests.
The alpha-shape is a subset of triangles of the Delaunay triangulation, where each triangle satisfies a circumscribing radius condition.
The following code is modified from my previous answer to compute the set of Delaunay triangles in the alpha-shape. Once the Delaunay triangulation and alpha-shape mask are computed, the fast method you reference can be adopted to the alpha-shape as I'll explain below.
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
The method from your first reference can then be adjusted to check not only if the point is in one of the Delaunay triangles (in which case it is in the convex hull), but also whether it is in one of the alpha-shape triangles.
The following function does this:
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
This method is very fast since it relies on the efficient search implementation of the Delaunay find_simplex() function.
Running it (with alpha=2) on the example data points from your post with the code below gives the results in the following figure, which I believe are not what you wished for...
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)
The reason for the result above is that, since there are large gaps between your input points, the alpha-shape of your data does not follow the polygon from your points. Increasing the alpha parameter won't help either since it will cut concave corners in other places. If you add more dense sample points then this alpha-shape method can be well-suited for your task. If not, then below I propose another solution.
Since your original polygon is not suited for the alpha-shape method, you need an implementation of a function that returns whether point(s) are inside a given polygon. The following function implements such an algorithm based on accumulating inner/outer angles (see here for an explanation).
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
Calling it with the code below results in the following figure, which I believe is what you were looking for.
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)
