I am trying to analyze greyscale TIFF stacks, in which a given frame will look like this. I filter it (using Gaussian blur), and then binarize it (using Otsu's method for threshold).
MATLAB code, which works great:
image_conncomp = bwconncomp(image_binary); # entire stack is held in image_binary
for i=1:image_conncomp.NumObjects
object_size = length(image_conncomp.PixelIdxList{i});
end
Each white spot in the example image is picked up, and its volume (in pixels) is pretty accurately given by object_size.
Python code:
from skimage import measure
labels = measure.label(image_binary, background=1) # same image_binary as above
propsa = measure.regionprops(labels)
for label in propsa:
object_size = len(label.coords)
The Python code seems to work decently... except that most detected objects will have object_size of 1 - 200, and then a couple will have a size of several thousand pixels.
What are these functions doing differently? I would be happy to try another approach in Python to get calculate object sizes, but I struggled to find another one. It'd be great to have a Python version of this code, if I could find a good substitute for Matlab's bwconncomp function.
2 Answers
Something like this?
from skimage.io import imread, imshow
from skimage.filters import gaussian, threshold_otsu
from skimage import measure
import matplotlib.pyplot as plt
original = imread(')
blurred = gaussian(original, sigma=.8)
binary = blurred > threshold_otsu(blurred)
labels = measure.label(binary)
plots = {'Original': original, 'Blurred': blurred,
'Binary': binary, 'Labels': labels}
fig, ax = plt.subplots(1, len(plots))
for n, (title, img) in enumerate(plots.items()):
cmap = plt.cm.gnuplot if n == len(plots) - 1 else plt.cm.gray
ax[n].imshow(img, cmap=cmap)
ax[n].axis('off')
ax[n].set_title(title)
plt.show(fig)
props = measure.regionprops(labels)
for prop in props:
print('Label: {} >> Object size: {}'.format(prop.label, prop.area))
Output:
Label: 1 >> Object size: 37
Label: 2 >> Object size: 66
Label: 3 >> Object size: 1
We could do the same by first applying scipy.ndimage's morphological closing on the thresholded binary image followed by the label() function to merge the connected regions in the binary image, as shown below (size of the regions are a bit different though and will depend upon the size of the morphological kernel):
from scipy.ndimage import label
from scipy.ndimage.morphology import binary_closing
from skimage.filters import threshold_otsu
import matplotlib.pylab as plt
original = plt.imread(')
thres = threshold_otsu(original)
binary = original > thres
binary_closed = binary_closing(binary, structure=np.ones((2,2)))
labeled_image, num_features = label(binary_closed)
feature_areas = np.bincount(labeled_image.ravel())[1:]
print(feature_areas) # if we use 3x3 SE we shall get two regions with areas 24, 46
# [24 42 1]
plt.figure(figsize=(8,7))
plt.gray()
plt.subplot(221), plt.imshow(original), plt.axis('off'), plt.title('original')
plt.subplot(222), plt.imshow(binary), plt.axis('off'), plt.title('binray')
plt.subplot(223), plt.imshow(binary_closed), plt.axis('off'), plt.title('binray closed')
plt.subplot(224), plt.imshow(labeled_image, cmap='inferno'), plt.axis('off'), plt.title('labelled')
plt.show()
to obtain the following output:
