Find points that lie in a concave hull of a point cloud

Question

in this thread a method is suggested for masking out point that lie in a convex hull for example:

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

with which the set of masked points look like the following. However I am looking for a way that does so with a concave hull (similar to what the blue points suggest)

I found this thread that suggest some functionality for a concave border but am not sure yet if it is applicable in my case. Does anyone has a suggestion?

Answer 1

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)