how to compute convex hull image/volume in 3d numpy arrays

Question

I wonder if there are any numpy based tools that can:

1. given a binary input numpy image in 3D, find its convex hull;
2. and return a list of indices or similar of the voxels (3D pixels) that are within this 3D convex hull.

One possibility is to use skimage.morphology.convex_hull_image(), but this only supports 2D images, so then i have to call this function slice by slice (in the z-axis), which is slow. [Edit: See note below.]

I definitely prefer to a more efficient way. For example, scipy.spatial.ConvexHull() could take a list of points in N-dimensional space, and return a convex hull object which seems not support finding its convex hull image/volume.

points = np.transpose(np.where(image))
hull = scipy.spatial.ConvexHull(points)
# but now wonder an efficient way of finding the list of 3D voxels that are inside this convex hull

Any ideas how? Please note efficiency is important to my application. Thanks!

---

Update: In the meantime, convex_hull_image() has been extended to support ND images, but it is quite slow for moderately sized data. The accepted answer below is much faster.

Answer 1

You should be able to do this:

def flood_fill_hull(image):    
    points = np.transpose(np.where(image))
    hull = scipy.spatial.ConvexHull(points)
    deln = scipy.spatial.Delaunay(points[hull.vertices]) 
    idx = np.stack(np.indices(image.shape), axis = -1)
    out_idx = np.nonzero(deln.find_simplex(idx) + 1)
    out_img = np.zeros(image.shape)
    out_img[out_idx] = 1
    return out_img, hull

It may not be the fastest, but short of a ready-made function it should work.

Testing:

points = tuple(np.rint(10 * np.random.randn(3,100)).astype(int) + 50)
image = np.zeros((100,)*3)
image[points] = 1

%timeit flood_fill_hull(image)
10 loops, best of 3: 96.8 ms per loop

out, h = flood_fill_hull(image)

plot.imshow(out[50])

Can't upload pictures, but it seems to do the trick.