Skimage Watershed And Particles Size Detection

I have the following image. I was able to use watershed to detect all the particles using the code below. However, now I need to calculate the size of each particles in the figure

Solution 1:

I am sharing an approach with watershed and regionprops

from skimage import io
import numpy as np
import matplotlib.pyplot as plt
from skimage.feature import peak_local_max
from skimage.measure import regionprops
from skimage.morphology import watershed
from scipy.ndimage.morphology import binary_erosion, binary_dilation, distance_transform_edt
from scipy.ndimage import label

import pandas as pd


img = io.imread('obvvX.jpg')

a = gaussian(img, sigma=5)
a = np.sum(a, axis=2)
a_thr = a < 1
plt.imshow(a)

# clean up specks
a_thr = binary_erosion(a_thr, iterations = 5)
a_thr = binary_dilation(a_thr, iterations = 5)

# do distance transform as prepartion for watershed
distances = distance_transform_edt(a_thr)

# find watershed seeds
seeds = peak_local_max(distances, indices =False, min_distance=20, footprint=np.ones((3,3)))
seeds = label(seeds)[0]

# watershed
ws = watershed(a, seeds, mask=a_thr)

plt.imshow(ws, cmap='tab20c')

So, the scale bar is also recognised as objects. We can now use regionprops to get the areas:

# compute region properties
props = regionprops(ws)

# exclude the bar on the bottom left:
props = [p for p in props if p['centroid'][0]<950 and p['centroid'][1]>400]

# get the sizes for each of the remaining objects and store in dataframe
entries = []
for p in props:
    entry = [p['label'], p['area'], p['perimeter'], *p['centroid']]
    entries.append(entry)


df = pd.DataFrame(entries, columns= ['label', 'area', 'perimeter', 'y', 'x'])

The dataframe has some entries that are too small to be actual objects. These can be deleted by setting a lower size threshold:

df=df[df['area']>40]labelareaperimeteryx0143275.01219317.0486111182.2361111249079.25483448.781633679.4387762358086.08326198.012069851.2603453460189.740115116.3826961047.9434284572998.911688126.149520972.5541845659588.669048226.092437663.6739506766494.325902263.8087351018.5602417813643.313708323.875000756.86764789382107.012193332.437173764.95811511126936.041631359.4202901028.507246121338670.426407475.4145081498.5466321415576117.876154503.248264481.036458181914660.656854524.890411484.308219192041589.597980532.655422492.6674702021580114.118795533.4086211383.151724222469596.568542581.5856121038.273381232528871.976659605.1145831522.27083324267732.485281611.6103901529.7792212628666124.704581634.734234676.509009272920552.769553696.9219511083.165854283055584.426407719.8126131220.690090293160588.669048745.538843743.3041323133637119.497475762.742543931.612245323449179.254834784.340122410.175153333570097.154329793.7357141179.764286343671296.911688846.039326987.450843353752889.740115932.549242984.071970

Solution 2:

Here is one way to do it using blobs in Python/OpenCV.

• Read the image
• Convert to grayscale
• Gaussian smooth the image to reduce noise
• Apply adaptive thresholding
• Use Simple Blob Detector with appropriate limits on characteristics to get key points and their size and locations

Input:

import numpy as np
import cv2
import math

# read image
img = cv2.imread("particles.jpg")

# convert to grayscale
gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

# apply Gaussian Blur
smoothed = cv2.GaussianBlur(gray, (0,0), sigmaX=9, sigmaY=9, borderType = cv2.BORDER_DEFAULT)

# do adaptive threshold on gray image
thresh = cv2.adaptiveThreshold(smoothed, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 65, 10)

cv2.imshow("Threshold", thresh)
cv2.waitKey(0)
cv2.destroyAllWindows()

# Set up the SimpleBlobdetector with default parameters.
params = cv2.SimpleBlobDetector_Params()

# Change thresholds
params.minThreshold = 0
params.maxThreshold = 256# Filter by Area.
params.filterByArea = True
params.minArea = 30
params.maxArea = 10000# Filter by Color (black=0)
params.filterByColor = True
params.blobColor = 0# Filter by Circularity
params.filterByCircularity = True
params.minCircularity = 0.5
params.maxCircularity = 1# Filter by Convexity
params.filterByConvexity = True
params.minConvexity = 0.5
params.maxConvexity = 1# Filter by InertiaRatio
params.filterByInertia = True
params.minInertiaRatio = 0
params.maxInertiaRatio = 1# Distance Between Blobs
params.minDistBetweenBlobs = 0# Do detecting
detector = cv2.SimpleBlobDetector_create(params)

# Get keypoints
keypoints = detector.detect(thresh)

print(len(keypoints))
print('')

# Get keypoint locations and radiusfor keypoint in keypoints:
   x = int(keypoint.pt[0])
   y = int(keypoint.pt[1])
   s = keypoint.size
   r = int(math.floor(s/2))
   print (x,y,r)
   #cv2.circle(img, (x, y), r, (0, 0, 255), 2)# Draw blobs
blobs = cv2.drawKeypoints(thresh, keypoints, np.array([]), (0,0,255), cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow("Keypoints", blobs)
cv2.waitKey(0)
cv2.destroyAllWindows()

# Save result
cv2.imwrite("particle_blobs.jpg", blobs)

Results:

25 points:1143 965199969422213192891589205859217987845151180 794154117841593276214743745141221 71913677635151523 606141039 58114211539151383 53314486516211498 47413763330131019 2641466422614973126151048 1161485299146794914

Output Image:

See this example for discussion of arguments

A second approach might be to get the contours in place of the blobs. Then get the bounding boxes of the contours and from that compute the radii and centers.

A third approach might be to use connected components with stats. Again it would get the bounding boxes and areas and centroids from which you could compute the radius and draw circles.

Solution 3:

By following the example of warped, I was able to pretty much solve the problem. You can find the new code below. I though that this might be useful to others.

I still have some questions though: 1) Watershed segmentation finds more areas than expected. For example, if you closely check one of those binary clusters of nanoparticles, it finds 3-4 different areas instead of just 2. These areas are usually small and I got rid of them using a size threshold, as warped suggested. However, is it possible to fine tune the watershed to somehow merge those areas and get a more accurate result?

2) I prefer using cv2.imshow() to show the images. However for some reasons I cannot plot the watershed result (variable name: labels) with that command, that's why I used matplotlib in the first part of the code. Does anyone have an explanation and a fix for this?

import numpy as np
import cv2
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage.morphology import watershed
from skimage.feature import peak_local_max
from skimage.measure import regionprops

#----------------------------------------------------------------------------------------------------------------------## IMAGE PRETREATMENT

img = cv2.imread('Test images/TEM of nanoparticles/NP good 0010.tif')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
Gaussian_Blur = cv2.GaussianBlur(gray,(21, 21), cv2.BORDER_DEFAULT)

# Use fixed threshold to mask black areas
_, thresh = cv2.threshold(Gaussian_Blur, 90, 255, cv2.THRESH_BINARY_INV) # _ = 30# Morphological closing to close holes inside particles; opening to get rid of noise
img_mop1 = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (7, 7)))
img_mop = cv2.morphologyEx(img_mop1, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15)))
tiled_h = np.hstack((img_mop1, img_mop)) # stack images side-by-side

plt.figure('Pretreatment')
plt.subplot(2, 2, 1) # Figure two has subplots 2 raw, 2 columns, and this is plot 1
plt.gca().set_title('Gray')
plt.xticks([]), plt.yticks([]) # To hide axes
plt.imshow(gray, cmap='gray')

plt.subplot(2, 2, 2)
plt.gca().set_title('Gaussian_Blur')
plt.xticks([]), plt.yticks([])
plt.imshow(Gaussian_Blur, cmap='gray')

plt.subplot(2, 2, 3)
plt.gca().set_title('Thresh')
plt.xticks([]), plt.yticks([])
plt.imshow(thresh, cmap='gray')

plt.subplot(2, 2, 4)
plt.gca().set_title('img_mop')
plt.xticks([]), plt.yticks([])
plt.imshow(img_mop, cmap='gray')


#----------------------------------------------------------------------------------------------------------------------## WTERSHED WITH SKIMAGE

distance = ndi.distance_transform_edt(img_mop) # Calculates distance of pixels from background#Find peaks in an image as coordinate list or boolean mask.
local_maxi = peak_local_max(distance, indices=False, min_distance=50, footprint=np.ones((3, 3)), labels=img_mop)
markers = ndi.label(local_maxi)[0]
labels = watershed(-distance, markers, mask=img_mop)

plt.figure('Processing')
plt.subplot(2, 2, 1) # Figure two has subplots 2 raw, 2 columns, and this is plot 1
plt.gca().set_title('Distance trans')
plt.xticks([]), plt.yticks([]) # To hide axes
plt.imshow(distance, cmap='gray')

plt.subplot(2, 2, 2)
plt.gca().set_title('local_maxi')
plt.xticks([]), plt.yticks([])
plt.imshow(local_maxi, cmap='gray')

plt.subplot(2, 2, 3)
plt.gca().set_title('markers')
plt.xticks([]), plt.yticks([])
plt.imshow(markers, cmap='gray')

plt.figure('Watershed')
plt.gca().set_title('Watershed')
plt.xticks([]), plt.yticks([]) # To hide axes
plt.imshow(labels)

plt.show()

#----------------------------------------------------------------------------------------------------------------------## DATA ANALYSIS# Regionprops: Measure properties of labeled image regions. It can give A LOT of properties, see info in:# https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.regionprops
props = regionprops(labels)

# Determine scale bar (largest object) and set the scale.
thr_size = 6000for p in props:
    if p['area'] > thr_size:
        box = p['bbox']
        scale = box[3]-box[1]


# Remove smaller detected areas, and give area and diameter for each of the remaining particles.for p in props:
    if p['equivalent_diameter'] > 15and p['equivalent_diameter'] < 40:
        entry = [p['label'], p['area'], p['equivalent_diameter'], *p['centroid']]
        n = entry[0]
        y = entry[3]
        x = entry[4]-60# so that number shows on the left of particle
        cv2.putText(img, str(n), (int(x), int(y)), cv2.FONT_HERSHEY_COMPLEX, 1, (0, 255, 0), 2)
        print('Particle {} | Area (nm^2): {}; Equivalent diameter (nm): {}'.format(str(n),
                                            str(int(((entry[1]*40000)/(scale**2)))), str(int((entry[2])*200/scale))))

cv2.imshow('img', img)
cv2.waitKey(0)
cv2.destroyAllWindows()