Crop area from 3d array

Question

I have a 3D numpy array with shape of (400, 512, 512), and I have bounding box, for example:

[376, 250, 206]
[380, 256, 211]

Now, I want to crop the bounding box area from the array and visualize. Normally in 2D array/image we can crop using the min max value like:

img[y:y + h, x:x + w]

But I don't think it's like that in 3D array, as the z axis is also involved. Can anyone suggest any example how to crop an area from a 3D numpy array using min max bounding box coordinates? In the end I will fit a cylinder in extracted bounding box area.

Edit: I want all the indices and their intensity values inside the bounding box of the main big 3D array. I can do like:

big_3d_array[z1:z2, x1:x2, y1:y2]
where,    
z1 = 376, z2 = 380, x1 = 250, x2 = 256, y1=206, y2= 211

Now problem with this approach is that I get the intensity values only, can not get the indices. Also I am not sure whether I am getting all the coordinates inside the bounding box using this approach. Because in 3d space bounding box should have 8 corner points. Here I only have the min max of the bounding box.

Thanks @MadPhysicist. When I read at the first. It seemed fine. Your answer is correct. But i was looking something like that: [link](https://stackoverflow.com/questions/42352622/finding-points-within-a-bounding-box-with-numpy) Get all the points from a 3d array inside the bbox, using that bbox's min max coordinates. May be your answer does the trick. I might have misunderstood. — Tokai, Oct 06 '20 at 12:16
You really need to make the question clear then. Is the bbox in array indices or in the x-y-z coordinates of the data? Please provide a meaningful example of what you want — Mad Physicist, Oct 06 '20 at 13:06
Why would the link not work for you in 3D? What you're doing in the question is not following the suggestion of the link at all. Either that solution works, or mine does, but it's hard to tell unless you show what you want and what you did. — Mad Physicist, Oct 06 '20 at 13:08
@MadPhysicist, I am beginner. So i might make mistake :/ So, I have two coordinates of a bounding box `[[z1,x1,y1] [z2,x2,y2]]` I have a big 3d array of `z*x*y`. Now I want to get all the points and their indices that is inside the bounding box. As you can understand that bounding box has 8 corner points. And I have the two with min max. I am trying to extract all the indices that is inside this bounding box and also their respective intensity values. I hope that makes things a bit clear. — Tokai, Oct 07 '20 at 07:23
Please edit your question with an MCVE. Show how you decide when a point is in the bounding box. There are a couple of different ways to interpret. The other answer shows you one way, I show another — Mad Physicist, Oct 07 '20 at 11:20
Three pairs of coordinates defines eight corners. Since corners are taken as combinations of min/max for each coordinate, there are 2^n corners for n pairs of bounds — Mad Physicist, Oct 07 '20 at 11:23
What do you mean you don't have the indices? The indices are `np.arange(x1, x2)`, etc... — Mad Physicist, Oct 07 '20 at 11:24

Mad Physicist · Answer 1 · 2020-10-01T02:49:05.810

Assuming that the bounding box array consists of three columns for each dimension, and the rows are start and end coordinates, in that order:

c = np.array([[376, 250, 206]
              [380, 256, 211]])

Let's take a look at a basic index in 3D:

a = np.random.randint(10, size=(400, 512, 512))

You write an index like this:

a[376:380, 250:256, 206:211]

This index represents a function call to a.__getitem__ (actually type(a).__getitem__(a, ...). The argument is a tuple of slice objects. The tuple is specified by commas in the brackets, and slices are created when you use the : character:

a.__getitem__((slice(376, 380), slice(250, 256), slice(206, 211)))

Using this information, you can index the array in a couple of different ways:

a[c[0, 0]:c[1, 0], c[0, 1]:c[1, 1], c[0, 2]:c[1, 2]]
a[tuple(slice(*x) for x in c.T)]

The second way is more scalable in general.

Numpy offers a second way to index arrays, besides slicing, called advanced or fancy indexing. In this method, you provide arrays of coordinates that broadcast together into the shape of the output that you want. For example:

ix = np.arange(*c[:, 0]).reshape(-1, 1, 1)
iy = np.arange(*c[:, 1]).reshape(-1, 1)
iz = np.arange(*c[:, 2])
a[ix, iy, iz]

There is a helper object called np.ogrid which lets you create the index arrays much more seamlessly:

np.ogrid[tuple(slice(*x) for x in c.T)]

In all cases, if your bounding box is inclusive on the upper bound, you will need to increment it before using in an index:

c[1, :] += 1

Crop area from 3d array

1 Answers1