3

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.

Stuart Berg
  • 17,026
  • 12
  • 67
  • 99
galactica
  • 1,753
  • 2
  • 26
  • 36
  • 1
    Not my area of expertise, but [scipy.spatial](https://docs.scipy.org/doc/scipy/reference/spatial.html#delaunay-triangulation-convex-hulls-and-voronoi-diagrams) should be probably your first lib to check out (not numpy, sorry). ConvexHull there is based on [Qhull](http://www.qhull.org/) and it seems it's one of the more problematic code-sections in scipy if i remember correctly. Skimming through the code of your skimage-function, it's seem to be also based on scipy.spatial! – sascha Sep 19 '17 at 22:31
  • thanks! yes, scipy.spatial.ConvexHull() seems to handle 3D points correctly but it'd be great if skimage.morphology.convex_hull_image() could handle 3D images. Not sure if other libs are more complete on this. – galactica Sep 20 '17 at 00:14
  • sklearn's function is based on Delaunay instead of ConvexHull. So maybe there is some potential route to go to generalize their algorithm. But again: i don't know much here and this is only my first impression from 10 secs looking into their code (while wondering why a convex_hull function does not use spatial.ConvexHull, but Delaunay). **Remark:** i edited the tags and sacrificed the 3d-tag for computational-geometry to attract [experts](https://stackoverflow.com/tags/computational-geometry/topusers). – sascha Sep 20 '17 at 00:20
  • Also: [very relevant SO question (using Delaunay again)](https://stackoverflow.com/questions/24733185/volume-of-convex-hull-with-qhull-from-scipy), [other link 1](https://math.stackexchange.com/questions/60279/how-to-calculate-volume-of-3d-convex-hull), [other link 2](https://math.stackexchange.com/questions/600193/volume-of-a-convex-hull-in-n-dimensions). **Edit:** Oh, according to the newer answer in firstlink scipy calculates what you want.Check out the available attributes [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html#scipy.spatial.ConvexHull). – sascha Sep 20 '17 at 00:32
  • thanks for the pointers but not quite, i'm looking for a list of points that are within the convex hull, not just the volume or area. – galactica Sep 20 '17 at 01:23
  • Your title reads different, but ok. Let's wait for experts. But this sound pretty easy if one is motivated to do the [math](https://en.wikipedia.org/wiki/Half-space_(geometry)). – sascha Sep 20 '17 at 01:24
  • Not long after this question was asked, scikit-image added support for 3D convex hull images. But I think the answer below is still better. I've edited your question to reflect the current state of affairs. – Stuart Berg Aug 18 '20 at 02:43

1 Answers1

9

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.

Daniel F
  • 13,620
  • 2
  • 29
  • 55
  • Thanks! The idea of flood filling is cool! It's even better seeing this elegant numpy implementation (new in numpy)! Comparing to the slice by slice version, i saw at least a speed-up of 2x. – galactica Sep 20 '17 at 15:35
  • Well, it's not an actual flood fill implementation, but the output is the same. Basically it turns the hull into Delaunay simplices and looks for which simplex each point in the field is in. If it's not in a simplex, it returns -1. Add one and you can catch the points in the hull with `np.nonzero`. – Daniel F Sep 20 '17 at 15:43
  • Thanks for the further explanations! Looks like no other solutions are available at the moment. – galactica Sep 21 '17 at 00:45
  • This is a great answer. Although scikit-image now supports the 3D case out-of-the-box, this answer is still faster. You can make it even faster (slightly) by being stingy with RAM. That would detract from the elegance of this post, but here's one modified version: https://gist.github.com/stuarteberg/8982d8e0419bea308326933860ecce30 – Stuart Berg Aug 18 '20 at 02:42
  • Still much faster than scikit-image in 2023. :) If your initial image is a solid binary object, you can get an additional (slight) speedup if you give this function an image of just the outer countour points, which you can do efficiently with: `flood_fill_hull(scipy.ndimage.laplace(image, mode='constant'))` – Andrew Aug 25 '23 at 08:59