Convex hull of 3D image with sparsely dispersed particles


image

These are 2D maximum intensity projections of a 3D image I’m analysing. Since these are max intensity projections, the boundaries of ROI are not as clearly defined as they seem here.

Analysis goals

ROI consists of groups of cells and some noise. I was successfully able to segment individual cells, but to clear the acutal bordering cells, I need to know the boundaries of this ROI consisting of noise & cells.

Challenges

  • I have the above thresholded 3D image of cells. For achieving my goal, I was hoping to calculate 3D convex hull of thresholded image, but doing it directly on 3D image gave a MemoryError due to the size of images I’m using (103x1024x1024).
from skimage.morphology import convex_hull_image

ch = convex_hull_image(thresholded)

MemoryError: Unable to allocate 393. GiB for an array with shape (488, 108003328) and data type float64

  • I tried applying convex hull layer-by-layer, but it gave undersirable results as it couldn’t deal with upward or downward convex ROI. In the following image, I’ve overlapped the convex hull found by processing layer-by-layer.

  • Also, if there’s a better way to achieve my analysis goals, please suggest.

Do you have a full traceback of where the error occurs?

1 Like

Full traceback:

---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-47-ffaefe88f8f4> in <module>
      1 from skimage.morphology import convex_hull_image
----> 2 ch = convex_hull_image(thresholded)
        global ch = array([[[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       ...,

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]],

       [[0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]])
        global convex_hull_image = <function convex_hull_image at 0x0000015D0084A310>
        global thresholded = array([[[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       ...,

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]],

       [[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ...,
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]]])

~\AppData\Local\Programs\Python\Python38\lib\site-packages\skimage\morphology\convex_hull.py in convex_hull_image(image=array([[[False, False, False, ..., False, False,...False, False, False, ..., False, False, False]]]), offset_coordinates=True, tolerance=1e-10)
     87                                 (ndim, -1))
     88         # A point is in the hull if it satisfies all of the hull's inequalities
---> 89         coords_in_hull = np.all(hull.equations[:, :ndim].dot(gridcoords) +
        coords_in_hull = undefined
        global np.all = <function all at 0x0000015D3A74E790>
        hull.equations = array([[ 9.42899317e-01, -2.11533311e-01, -2.57282990e-01,
         8.40162988e+00],
       [ 9.72648329e-01, -2.22843957e-01, -6.55423403e-02,
        -2.95137158e+01],
       [ 9.31985033e-01, -2.97532252e-01, -2.07071140e-01,
         2.17303375e+01],
       ...,
       [ 6.09207699e-01, -6.09207699e-01,  5.07673083e-01,
        -4.79243390e+02],
       [ 5.80697000e-01, -5.80697000e-01,  5.70597922e-01,
        -5.35263772e+02],
       [ 5.80697000e-01, -5.80697000e-01,  5.70597922e-01,
        -5.35263772e+02]])
        ndim.dot = undefined
        gridcoords = array([[   0,    0,    0, ...,  102,  102,  102],
       [   0,    0,    0, ..., 1023, 1023, 1023],
       [   0,    1,    2, ..., 1021, 1022, 1023]])
        ndim = 3
        tolerance = 1e-10
        global axis = undefined
     90                                 hull.equations[:, ndim:] < tolerance, axis=0)
     91         mask = np.reshape(coords_in_hull, image.shape)

MemoryError: Unable to allocate 393. GiB for an array with shape (488, 108003328) and data type float64
1 Like

Hi,

so I can reproduce this with a minimal example:

import numpy as np
from skimage.morphology import convex_hull_image

tmp = np.random.randn(100, 1024, 1024) > 0
ch = convex_hull_image(tmp)

results in

     88         # A point is in the hull if it satisfies all of the hull's inequalities
---> 89         coords_in_hull = np.all(hull.equations[:, :ndim].dot(gridcoords) +
     90                                 hull.equations[:, ndim:] < tolerance, axis=0)
     91         mask = np.reshape(coords_in_hull, image.shape)

MemoryError: Unable to allocate 39.1 GiB for an array with shape (50, 104857600) and data type float64

So I am not sure about the hull object that is returned there and it’s equation methods, but as the error message implies this seems to be quite memory intensive.
I would recommend openingn an issue on https://github.com/scikit-image/scikit-image/ and linking to this thread here.

1 Like

Are you implying that there’s some fault in implementation of skimage.morphology.convex_hull_image. If the implementation could be improved to be less memory intensive, convex hull of such image could be retrieved?

Yup, it’s basically got to build the hull (in 3D, a bunch of planes), and for every plane, for every point’s coordinates, check whether it lives on the inside or outside of the hull according to that one plane. We create the entire array for that structure in one go.

We could check for a memory-error here and if so, split the volume into a few chunks before trying again. Like @VolkerH says, creating an issue with his minimal example in the scikit-image repo would be much appreciated! :blush:

As a side note, you will might see a significant speedup in hull computation if you first do surface = thresholded - binary_erosion(thresholded). I also think that for this application you want to set offset_coordinates=False, which will give you another speedup. See the documentation for convex_hull_image for a description.

1 Like

I think so! Note that the MemoryError happens after the hull itself has been computed. So the hull itself has been found, and the error happens in the creation of the image. We can definitely improve there.

1 Like

I found the right convex hull for my data, but to meet my analysis goals I have to remove all the objects that (I’ve thresholded & labeled) touch the convex hull. Just like the skimage.segmentation.clear_border function, but for me the border is a non-cuboidal convex hull.

And since I want to remove those labels/objects that touch my convex hull, not those which overlap borders of convex hull, I cannot use skimage.segmentation.clear_border function’s mask argument for my goals.

You can use the code for clear_border for inspiration:

You can find the border of your convex hull by doing a binary_erosion and subtracting it from the hull, then use the same logic as used in that function with the rectangular borders. In fact, I just saw that the mask= keyword argument for that function can be used for this! So you just need to compute your mask as above (though its negated, mask=~(hull_image - erosion(hull_image)), because of the definition of mask:

    mask : ndarray of bool, same shape as `image`, optional.
        Image data mask. Objects in labels image overlapping with
        False pixels of mask will be removed. If defined, the
        argument buffer_size will be ignored.
1 Like

Can I just use binary_erosion(hull_image) as mask in clear_borders? Instead of subtracting it from hull image.

You’re absolutely right! :man_facepalming: My solution was putting false pixels exactly on the border, yours puts them on the border and everything outside — which of course doesn’t matter since there are no objects there. Well spotted.

1 Like