Skimage分水岭和粒径检测 [英] 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.

但是，现在我需要计算图中每个粒子的大小，如果使用标签"图像，由于某些原因，我将无法使用功能cv2.findContours.

However, now I need to calculate the size of each particles in the figure and if I use the "labels" image, for some reasons I am not capable of using the function cv2.findContours.

有人愿意分享一些想法吗?如果您提出一些代码，请加入说明，因为我是初学者. :)

Anyone willing to share some ideas? If you propose some code, please include explanation because I am a beginner. :)

非常感谢！

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

#-------------------------------------------------------------------------------------------#
# 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) # Figure two has subplots 2 raw, 2 columns, and this is plot 1
plt.gca().set_title('Gaussian_Blur')
plt.xticks([]), plt.yticks([]) # To hide axes
plt.imshow(Gaussian_Blur, cmap='gray')

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

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


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

# Now we want to separate the two objects in image
# Generate the markers as local maxima of the distance to the background
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, footprint=np.ones((3, 3)), labels=img_mop)
# indices: if True, the output will be an array representing peak coordinates. If False, the output will be a boolean
# array shaped as image.shape with peaks present at True elements.
# If footprint == 1 represents the local region within which to search for peaks at every point in image.
# labels: if provided, each unique region labels == value represents a unique region to search for peaks. Zero is
# reserved for background.
# returns an array of boolean with True on max points
print('local_maxi lenght: ', local_maxi.shape)
print('local_maxi: ', local_maxi[0])
markers = ndi.label(local_maxi)[0]
print('markers lenght: ', markers.shape)
print('markers: ', markers[0])
labels = watershed(-distance, markers, mask=img_mop)
print('labels lenght: ', labels.shape)
print('labels: ', labels[0])


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) # Figure two has subplots 2 raw, 2 columns, and this is plot 1
plt.gca().set_title('local_maxi')
plt.xticks([]), plt.yticks([]) # To hide axes
plt.imshow(local_maxi, cmap='gray')

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

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

plt.show()

#-------------------------------------------------------------------------------------------#
# DATA ANALYSIS ---- WORK IN PROGRESS

cnts, _ = cv2.findContours(labels, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
img = cv2.drawContours(img, cnts, -1, (0, 255, 255), 2) # To print all contours
cv2.imshow('Contours',  cv2.resize(img, dsize=(0, 0), fx=0.3, fy=0.3))
print('\nCnts length: ', len(cnts), '\n') # 11 objects (10 nanoparticles + scale barr)





# Divide the cnts array into scalebar and nanoparticles
# Get bounding rectangles for the scale and the particles from detailed contour determine on line 32.
# cv2.boundingRect() outputs: x, y of starting point (top left corner), and width and height of rectangle.
# Find contours. For more info see: https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_imgproc/py_contours/py_contour_features/py_contour_features.html
# cv2.contourArea() outputs the area of each detailed contour, does not work on rectangle generated by cv2.boundingRect.
thr_size = 5000
for cnt in cnts:
    if cv2.contourArea(cnt) > thr_size:
        scale = [cv2.boundingRect(cnt)] # returns x, y, w, h

img = cv2.rectangle(img, (scale[0][0], scale[0][1]), (scale[0][0] + scale[0][2], scale[0][1] + scale[0][3]), (255, 255, 0), 2)
print('Scale is: ', scale) #only one box (object) = scalebar
print("scale[0][1] is scalebar's width of {} pixels".format(scale[0][2]), '\n')


# 8. MINIMUM ENCLOSING CIRCLE
i = 1
for cnt in cnts:
    if cv2.contourArea(cnt) < thr_size:
        # Find min enclosing circle and get xy of centre
        (x, y), radius = cv2.minEnclosingCircle(cnt)
        center = (int(x), int(y))

        # Get radius average method
        #rx, ry, w, h = cv2.boundingRect(cnt)
        #radius = int((((w+h)/2))*1.5)
        img = cv2.circle(img, center, radius, (255, 0, 255), 3)

        cv2.putText(img, str(i), (int(x), int(y)-20), cv2.FONT_HERSHEY_COMPLEX, 1, (0, 255, 0), 2)
        print('Particle ' + str(i) + ' | Horizontal diameter: ' + '{:.2f}'.format((radius/ scale[0][2] * 200)*2) + ' nm')
        i=i+1
cv2.imshow('img',  cv2.resize(img, dsize=(0, 0), fx=0.3, fy=0.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.

我仍然有一些疑问: 1)分水岭分割发现的区域超出预期.例如，如果您仔细检查纳米粒子的这些二元簇之一，它会发现3-4个不同的区域，而不是2个.这些区域通常很小，我使用大小阈值将它们除去，如翘曲所示.但是，是否可以对分水岭进行微调以某种方式合并这些区域并获得更准确的结果?

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)我更喜欢使用cv2.imshow()来显示图像.但是由于某些原因，我无法使用该命令绘制分水岭结果(变量名称:标签)，这就是为什么在代码的第一部分中使用matplotlib的原因.有人对此有解释和修正吗?

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 = 6000
for 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'] > 15 and 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()

这篇关于Skimage分水岭和粒径检测的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！