How to find the diameter of objects using image processing in Python?

Question

Given an image with some irregular objects in it, I want to find their individual diameter.

Thanks to this answer, I know how to identify the objects. However, is it possible to measure the maximum diameter of the objects shown in the image?

I have looked into the scipy-ndimage documentation and haven't found a dedicated function.

Code for object identification:

import numpy as np
from scipy import ndimage
from matplotlib import pyplot as plt

# generate some lowpass-filtered noise as a test image
gen = np.random.RandomState(0)
img = gen.poisson(2, size=(512, 512))
img = ndimage.gaussian_filter(img.astype(np.double), (30, 30))
img -= img.min()
img /= img.max()

# use a boolean condition to find where pixel values are > 0.75
blobs = img > 0.75

# label connected regions that satisfy this condition
labels, nlabels = ndimage.label(blobs)

# find their centres of mass. in this case I'm weighting by the pixel values in
# `img`, but you could also pass the boolean values in `blobs` to compute the
# unweighted centroids.
r, c = np.vstack(ndimage.center_of_mass(img, labels, np.arange(nlabels) + 1)).T

# find their distances from the top-left corner
d = np.sqrt(r*r + c*c)

# plot
fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10, 5))
ax[0].imshow(img)
ax[1].hold(True)
ax[1].imshow(np.ma.masked_array(labels, ~blobs), cmap=plt.cm.rainbow)
for ri, ci, di in zip(r, c, d):
    ax[1].annotate('', xy=(0, 0), xytext=(ci, ri),
                   arrowprops={'arrowstyle':'<-', 'shrinkA':0})
    ax[1].annotate('d=%.1f' % di, xy=(ci, ri),  xytext=(0, -5),
                   textcoords='offset points', ha='center', va='top',
                   fontsize='x-large')
for aa in ax.flat:
    aa.set_axis_off()
fig.tight_layout()
plt.show()

Image:

Answer 1

You could use skimage.measure.regionprops to determine the bounding box of all the regions in your image. For roughly circular blobs the diameter of the minimum enclosing circle can be approximated by the largest side of the bounding box. To do so you just need to add the following snippet at the end of your script:

from skimage.measure import regionprops

properties = regionprops(labels)
print 'Label \tLargest side'
for p in properties:
    min_row, min_col, max_row, max_col = p.bbox
    print '%5d %14.3f' % (p.label, max(max_row - min_row, max_col - min_col))

fig = plt.figure()
ax = fig.add_subplot(111)    
ax.imshow(np.ma.masked_array(labels, ~blobs), cmap=plt.cm.gist_rainbow) 
ax.set_title('Labeled objects')
plt.xticks([])
plt.yticks([])
for ri, ci, li in zip(r, c, range(1, nlabels+1)):
    ax.annotate(li, xy=(ci, ri), fontsize=24)
plt.show()

And this is the output you get:

Label   Largest side
    1        106.000
    2         75.000
    3         79.000
    4         56.000
    5        161.000
    6         35.000
    7         47.000  

Answer 2

I would propose using a distance transform. So once you've got your labeled image you do:

dt = ndimage.distance_transform_edt(blobs)
slices = ndimage.find_objects(input=labels)
radii = [np.amax(dt[s]) for s in slices]

This gives the largest inscribed circle (or sphere in 3D). The find_objects function is quite handy. It returns a list of Python slice objects, which you can use to index into the image at the specific locations containing the blobs. These slices can of course be used to index into the distance transform image. Thus the largest value of the distance transform inside the slice is the radius you're looking for.

There is one potential gothcha of the above code: the slice is a square (or cubic) section so might contain small pieces of other blobs if they are close together. You can get around this with a bit more complicated logic as follows:

radii = [np.amax(dt[slices[i]]*(labels[slices[i]] == (i+1))) for i in range(nlabels)]

The above version of the list comprehension masks the distance transform with the blob that is supposed to be indexed by the slice, thereby removing any unwanted interference from neighboring blobs.