ssm ssm - 1 year ago 106
Python Question

numpy extract 3d cubes from 3d array

I would like to extract the 3D cubes (3x3x3) from a boolean 3D array of (180x180x197). This is similar to scipy.ndimage.measurements.label but need to fixed size of (3x3x3).
Is there a fast way of doing this than using for loops.

Answer Source

In your special case, I suggest to use ndimage.minimum_filter Let's say your array is called ``a''. The following:

centers = ndimage.minimum_filter(a, 3, mode="constant")

will only contain ones where your array contained such box of True's Then you can use scipy.ndimage.measurements.label with the default structure to classify the boxes and maybe identify connected boxs. To locate them you can use ndimage.measurements.find_objects


The way above, you will correctly get the centers of all cubes in the array. To be clear, I think that is the answer to your initial question. In the comments, it turned out, that indeed only non-overlapping cubes are needed. Therefore, one needs to analyze the output of minimum_filter, where I can imagine many approaches.

One can use the following to obtain only one cube per cluster:

s = ndimage.generate_binary_structure(3,3)
labels, num = ndimage.measurements.label(centers, s)
locations = ndimage.measurements.find_objects(labels)
locations = map(lambda slices: [slices[i].start for i in xrange(3)], locations)

Now the problem occurs that cubes are lost which are not overlapping but share a face. Actually one can imagine quite complicated structures of non-overlapping, facesharing cubes. Certainly, there are several solutions (sets of non-overlapping cubes) that can be found for this problem. So it is a completely new task to choose a set of cubes from the found centers and I think you will have to find one which is ideal for you.

One way could be to iterate through all solutions and set each found cube to False:

get_starting_point = numpy.vectorize(lambda sl: sl.start) #to be applied on slices
s = ndimage.generate_binary_structure(3,3)
result = []

while True:
    labels, num = ndimage.measurements.label(centers, s)
    if not num:
    locations = ndimage.measurements.find_objects(labels)
    sp = get_starting_point(locations)
    for p in sp:
        centers[p[0]-1:p[0]+2, p[1]-1:p[1]+2, p[2]-1:p[2]+2] = False

numiter = len(results)
results = numpy.vstack(results)

I guess only very few iterations will be necessary.

I hope this is what you were looking for