read_images calculating mean expression of bounding box, not cell?

[bswhite]

Hello,

If I'm not mistaken aren't the following lines in read_images computing the mean expression of the bounding box surrounding a masked region (e.g., a cell):

for v in measurements['image_intensity']:
            cell_exp = getattr(v, intensity_measure).__call__(axis=(0, 1))

I'm no regionprops expert, but its document says that 'image_intensity' is the image inside the bounding box. Don't we instead want to calculate the mean expression of the region within the bounding box?

Thank you,
Brian

[bswhite]

I believe to calculate the mean intensity only within the region you want to first subset the intensity using that region mask:

for v, m in zip(measurements['image_intensity'], measurements['image']):
    cell_exp = getattr(v[m], intensity_measure).__call__(axis=(0))

Calculating the mean of 'image_intensity' gives the mean of the bounding box. Calculating the mean of 'image_intensity' after it has been subsetted by the region mask 'image' gives the mean of the region within the bounding box.

Here is a toy example with two cells and two channels, showing the current spatialtis implementation (mean of 'image_intensity'), the above version (mean of 'imagine_intensity' after subsetting by mask), and that the above version is consistent with having regionprops calculate the mean intensities:

from skimage.measure import regionprops_table

# a is the mask with two cells
a = np.zeros((10, 10)).astype(int)
a[2:7, 2:7] = 1
a[5:7,5:7]=2

# b is a 2-channel image
nchannels = 2
b = np.zeros((a.shape[0], a.shape[1], nchannels)).astype(np.float32)
b[:,:,0] = a.copy()
b[:,:,1] = a.copy()
b[2,2,0] = 100
b[6:7,6:7,0]=49

# remove the top left corner of cell 1
a[2,2] = 0
print('Mask:')
print(a)
print('\nIntensity for channel 0:')
print(b[:,:,0])

print('\nIntensity for channel 1:')
print(b[:,:,1])

print('\nNote that both channels have intensity values at (2,2).')
print('These are inside the bounding box, but outside the mask, of cell 1.')
print('i.e., cell 1 mean expression should be 1 for both channels')

properties = ['label', 'intensity_mean', 'intensity_min', 'intensity_max', 'image', 'image_intensity', 'slice']
measurements = regionprops_table(a, intensity_image=b,  properties=properties)

# Calculate mean intensities as spatialtis does
intensity_measure = "mean"
print("\nspatialtis-style mean intensities")
for idx, v in enumerate(measurements['image_intensity']):
    print('spatialtis-style mean intensities for cell ' + str(idx) + ' across both channels')
    cell_exp = getattr(v, intensity_measure).__call__(axis=(0, 1))
    print(cell_exp)

print("\nspatialtis-style mean intensities after subsetting by mask")
# This mimics the logic of regionprops mean_intensity here:
# https://github.com/scikit-image/scikit-image/blob/f1892afb76e60b6fe41cf0c73be97f918581aa33/skimage/measure/_regionprops.py#L571
# i.e., it subsets the image intensity by the mask ('image')
for idx, v, m in zip(range(0,len(measurements['image_intensity'])), measurements['image_intensity'], measurements['image']):
    print('spatialtis-style mean intensities after subsetting by mask for cell ' + str(idx) + ' across both channels')
    cell_exp = getattr(v[m], intensity_measure).__call__(axis=(0))
    print(cell_exp)

# Report mean intensities from regionprops
print("")
print('regionprops intensity_mean-0 (i.e., channel 0) across both cells')
print(measurements['intensity_mean-0'])

print('regionprops intensity_mean-1 (i.e., channel 1) across both cells')
print(measurements['intensity_mean-1'])

# Note that there is a difference between subsetting by a mask (which only includes the region within the bounding box)
# and multiplying by the mask (which zeroes out areas outside the mask within the bounding box).
# The means will be different in the two cases
print("\nimage_intensity * image (mask)")
for idx, v, m in zip(range(0,len(measurements['image_intensity'])), measurements['image_intensity'], measurements['image']):
    print("image_intensity * image (mask) for cell " + str(idx) + ' and channel 0')
    print(v[:,:,0] * m)
    print("mean of image_intensity * image (mask) for cell " + str(idx) + ' and channel 0')
    print(np.mean(v[:,:,0] * m, axis=(0,1)))
    print('')

print("\nimage_intensity subsetted by image (mask)")
for idx, v, m in zip(range(0,len(measurements['image_intensity'])), measurements['image_intensity'], measurements['image']):
    print("image_intensity subsetted by (mask) for cell " + str(idx) + ' and channel 0')
    print(v[m,0])
    print("mean of image_intensity subsetted by image (mask) for cell " + str(idx) + ' and channel 0')
    print(np.mean(v[m,0]))
    print('')

Here is the output:

Mask:
[[0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0]
[0 0 0 1 1 1 1 0 0 0]
[0 0 1 1 1 1 1 0 0 0]
[0 0 1 1 1 1 1 0 0 0]
[0 0 1 1 1 2 2 0 0 0]
[0 0 1 1 1 2 2 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]]

Intensity for channel 0:
[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 100. 1. 1. 1. 1. 0. 0. 0.]
[ 0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
[ 0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
[ 0. 0. 1. 1. 1. 2. 2. 0. 0. 0.]
[ 0. 0. 1. 1. 1. 2. 49. 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.]]

Intensity for channel 1:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
[0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
[0. 0. 1. 1. 1. 1. 1. 0. 0. 0.]
[0. 0. 1. 1. 1. 2. 2. 0. 0. 0.]
[0. 0. 1. 1. 1. 2. 2. 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.]]
Note that both channels have intensity values at (2,2).
These are inside the bounding box, but outside the mask, of cell 1.
i.e., cell 1 mean expression should be 1 for both channels

spatialtis-style mean intensities
spatialtis-style mean intensities for cell 0 across both channels
[0.8 0.8]
spatialtis-style mean intensities for cell 1 across both channels
[13.75 2. ]

spatialtis-style mean intensities after subsetting by mask
spatialtis-style mean intensities after subsetting by mask for cell 0 across both channels
[1. 1.]
spatialtis-style mean intensities after subsetting by mask for cell 1 across both channels
[13.75 2. ]

regionprops intensity_mean-0 (i.e., channel 0) across both cells
[ 1. 13.75]
regionprops intensity_mean-1 (i.e., channel 1) across both cells
[1. 2.]

image_intensity * image (mask)
image_intensity * image (mask) for cell 0 and channel 0
[[0. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]
[1. 1. 1. 0. 0.]
[1. 1. 1. 0. 0.]]
mean of image_intensity * image (mask) for cell 0 and channel 0
0.8

image_intensity * image (mask) for cell 1 and channel 0
[[ 2. 2.]
[ 2. 49.]]
mean of image_intensity * image (mask) for cell 1 and channel 0
13.75

image_intensity subsetted by image (mask)
image_intensity subsetted by (mask) for cell 0 and channel 0
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
mean of image_intensity subsetted by image (mask) for cell 0 and channel 0
1.0

image_intensity subsetted by (mask) for cell 1 and channel 0
[[ 2. 2.]
[ 2. 49.]]
mean of image_intensity subsetted by image (mask) for cell 1 and channel 0
13.75

[Mr-Milk]

Sorry, I'm not longer maintained this package, this issue arises from the API changes in skimage.

I will mark this project as Archived. If you want to use it, you are encouraged to fork it to your own repository and start from there.