pythonmatplotlibpolygonshapelyboundary

Find a concave hull of an array of integer points on a grid

I have arrays of thousands of points on an integer grid. I'm looking for a fast way of finding a concave hull for the points. Points are adjacent if they are 1 unit away in any cardinal direction. I am ambivalent about whether the boundary should move diagonally (i.e. cutting corners as shown below from [391,4036] -> [392,4037]), and instead prioritise speed of calculation. There are no interior holes. I'm working in python.

an example case using a smaller set of points

My initial thought was to loop through each point, and look up whether it's cardinal neighbours are also in the set of points, and if one of them is not then mark the point as being on the bounds of the shape. I would then need some other algorithm for ordering those points to get the (clockwise or anti-clockwise) boundary.

This would not scale well with number of points, as for each point I need to check it's four cardinal directions against every other point in the set for membership.

Python code for finding boundary points:

boundary_pixels = [
    (row_idx, col_idx)
    for (row_idx, col_idx) in full_pixels
    if not (
        ((row_idx+1, col_idx) in full_pixels) & 
        ((row_idx-1, col_idx) in full_pixels) &
        ((row_idx, col_idx+1) in full_pixels) & 
        ((row_idx, col_idx-1) in full_pixels)
    )
]

I know that finding concave hulls is a difficult problem, but is there a solution for when the points are evenly spaced on a grid?

Solution

  • There is a young repo (concav_hull, that relies on a JavaScript algorithm).

    You might find it interesting. If so, here is a proposition that you can adjust if needed :

    import matplotlib.pyplot as plt
from concave_hull import concave_hull_indexes

idxes = concave_hull_indexes(points[:, :2], length_threshold=-0)

plt.plot(points[:, 0], points[:, 1], "rx")

for f, t in zip(idxes[:-1], idxes[1:]):
    seg = points[[f, t]]
    plt.plot(seg[:, 0], seg[:, 1], "k-.")

    enter image description here

    Used input :

    import numpy as np

op = [
    (389, 4034), (389, 4035), (389, 4036),
    (390, 4035), (390, 4036),
    (391, 4034), (391, 4035), (391, 4036),
    (392, 4032), (392, 4033), (392, 4034), (392, 4035), (392, 4036), (392, 4037),
    (393, 4032), (393, 4033), (393, 4034), (393, 4035), (393, 4036), (393, 4037),
    (394, 4034), (394, 4035), (394, 4036), (394, 4037)
]

points = np.array(op)