Mask (ROI) as NumPy array

I have a 3D nii image file. I have drawn a 3D label (ROI) in ITK snap. Now I want to analyze this ROI beyond the limited options in the labelIntensityStatistics. For example, I want to calculate the 10th, 20th, … percentiles. One idea is to convert this ROI into numpy array with the real values from the original image but I can’t find a solution.

I tried:
import SimpleITK as sitk
import numpy as np

ADC600 = sitk.ReadImage( ‘/maps/ADC600.nii’)
ROI = sitk.ReadImage( ‘/ROI2.nii’)
shape = sitk.LabelShapeStatisticsImageFilter()
bounds = shape.GetBoundingBox(1)

new_im = sitk.RegionOfInterest(ADC600,bounds)
new_im = sitk.GetArrayFromImage(newim)
mean = np.mean(new_im)

But the value is totally off

Any help?

The Regions in SimpleITK are defined here:

The RegionOfInterest filter is expecting a region, so try GetRegion. Also consider using the LabelStatisticsImageFilter over the LabelShapeStatisticsImageFilter for this purpose.

Also inspect the “bounds”, and the resulting new_im from your cropping operations to give hints on wha the problems are.

I made the changes but still the same result.

Is there a way to get the individual pixel values from an ROI? then I can do the first order statistics.

Using ITK regions and iterators seem like an overkill for getting intensity statistics. Getting n-th percentile of grayscale voxel values corresponding to a specific label value in the mask volume is a one-liner in Python (assuming both are stored as a numpy array of the same size, corresponding to the same image geometry).

For example, for grayscale image image you can get 10th percentile voxel values of a structure labelled as 18 in the corresponding mask volume by calling:

np.percentile(image[mask==18], 10)

This appears to be a cross posting, with the OP accepting the answer on GitHub.

The proposed solution, which uses a numpy mask, or numpy.extract have the additional memory and computational time on the order of the number of pixels.

One future possibility is to add a method to sitk::LabelIntensityStatisticsImageFilter::GetIntensities which returns an array of the values for a given label. For cases where there are a large number of small labels this would be significantly more computationally efficient.

1 Like

Thanks all.

Yes I have found numpy.extract option to be useful. I didn’t get time to update the forum here.

For reference, for a 512^3 16-bit volume with a fairly large mask, numpy provides result in a fraction of a second (around 200msec). Interestingly, np.extract is about 50% slower than direct array indexing image[condition].

So, using numpy for image masking is indeed much slower than ITK but still good enough for many use cases.

1 Like

Yes certainly.

To clarify my comment, if you need to do this operations for many different labels for example over a set of organs or a set of cells it is computationally beneficial to use the ITK infrastructure which works with multiple labels.

1 Like