extract 40 mm cube from 3d volume data

Question

I have created a 3d numpy array of data (shape: 133 x 512 x 512) from dicom format. Given a center point of a nodule location, how would i extract a 3d volume of size 40mm x 40mm x 40mm. I am not sure how would the pixel to mm conversion happens? I have attached the data here. The location of nodule is (317, 363, 89) which is (x,y,z) where z denotes the number of slice. So for this example, the nodule is on slice 89. Here is the sample data. It is in nrrd format. In the dicom header information, the slice thickness was given 2.5 and pixel spacing as ['0.703125', '0.703125'].

Answer 1

It looks like the NumPy array coordinates are reversed to (z, y, x) aka (plane, row, column), see this scikit-image documentation for more details on this. Additionally, let's take the voxel spacing to be equal to [2.5, 0.703125, 0.703125], though I couldn't find this information with pynrrd. But, visualising it, I these measurements looked right. The units of this are mm/pixel.

Now, you want to convert 40mm, or rather half that, 20mm, to pixels. You do this by dividing by the mm/pixel, so you get a box size of:

>>> 20 / np.array([2.5, 0.703125, 0.703125])
array([ 8.        , 28.44444444, 28.44444444])

That means that the top of the box should be at:

>>> np.array([89, 363, 317]) - [8, 28, 28]
array([ 81, 335, 289])

and the bottom at:

>>> np.array([89, 363, 317]) + [8, 29, 29]
array([ 97, 392, 346])

(I've made it 29 on the other side to get a centre pixel and also because 28.4 * 2 approximates to 57, not 56.)

So, to get your box, you would do:

>>> box40mm = data[81:97, 335:392, 289:346]

See the NumPy indexing documentation for more details about NumPy indexing.

To put it all together in a function, which assumes all the inputs are NumPy arrays:

def get_box(data, centre_pixel, box_size, pixel_spacing):
    box_size_pixels = box_size / pixel_spacing
    box_size_left = np.round(box_size_pixels / 2).astype(int)
    box_size_right = np.round(box_size_pixels - box_size_left).astype(int)
    starts = centre_pixel - box_size_left
    ends = centre_pixel + box_size_right
    slices = tuple(
        slice(start, end)
        for start, end in zip(starts, ends)
    )
    return data[slices]

Note that this code does not account for hitting the boundaries of the volume.