基于数字图像处理的小目标计数（二）

存在接壤现象的细胞计数

因为前面（一）的中的图片每个波点之间间隔较小，这边的细胞图中存在多个目标接壤问题，所以怎样在存在接壤现象的细胞中准确数出个数是这里较前面不同的地方。

解决思路：

• 进行锐化和平滑滤波操作
• 根据灰度图进行阈值分割
• 形态学处理
• 分水岭算法
• 修正

---

如果使用Python，可以选择的图像处理库有
1、工业级的图像处理库：Opencv
2、PIL
3、Skimage
我使用的是Skimage，因为它和numpy联和的比较紧密，使用较为方便，能暂时满足我基本需求。

进行锐化和平滑滤波操作

我选取的图片是一张细胞图，背景和前景的照片亮度信息差别很大，因此只需要根据灰度值，来进行前景和背景的分离。

# 锐化+平滑滤波操作
img_path1 = '细菌图.jpg'

from PIL import Image, ImageFilter
im = Image.open(img_path1)
# 没找到skimage的锐化api...所以用了PIL的
im.filter(ImageFilter.SHARPEN).save('细菌图_temp.jpg')

sharpen_data = io.imread('细菌图_temp.jpg')
img = color.rgb2gray(sharpen_data)
dst = median_each(img)

plt.figure('filters',figsize=(14,14))

plt.subplot(121)
plt.title('origin image')
orgin_img = io.imread(img_path1)
plt.imshow(orgin_img)

plt.subplot(122)
plt.title('sharpen and Smoothed image')
plt.imshow(dst,plt.cm.gray)

根据灰度图进行阈值分割

# 直方图信息
plt.figure("hist")
arr=dst.flatten()
plt.hist(arr, bins=256, normed=1,facecolor='b',edgecolor='b',hold=1)
plt.show()

# 手动查看鼠标对应H值
%matplotlib qt5
loop_num = 5
plt.imshow(dst,plt.cm.gray)

pos=plt.ginput(loop_num)
for i in range(loop_num):
    x,y = int(pos[i][1]),int(pos[i][0])
    print '第%d个点击的 x,y：' % int(i+1) ,'(', x , y,')'
    print '对应的灰度值为：',dst[x,y],'\n'

# 手动选取阈值
%matplotlib inline
rows,cols=dst.shape
labels=np.zeros([rows,cols])
for i in range(rows):
    for j in range(cols):
        if(dst[i,j] < 135):
            labels[i,j]=1
        else:
            labels[i,j]=0
io.imshow(labels)

形态学处理

先膨胀填充孔洞，然后进行开操作去除孤立的小点

# 膨胀,开操作
labels= sm.dilation(labels,sm.square(3))
labels = sm.opening(labels,sm.disk(5))

分水岭算法

偶然发现了分水岭算法，它是一种基于区域的图像分割算法，很适合把这里的接壤细胞区分开来。在分水岭算法第一步会使用一个关键的api——peak_local_max，他可以找到局部的“山峰”，相当于找到了我们要找的细胞的中心！
分水岭算法也在skimage的形态学包里，具体的使用可以查看官网。

from scipy import ndimage as ndi
# labels = dst
distance = ndi.distance_transform_edt(labels) #距离变换
# min_distance：最小的像素在2×min_distance + 1区分离（即峰峰数至少min_distance分隔）。找到峰值的最大数量，使用min_distance = 1。
# exclude_border：不排除峰值在图像的边界
# indices：False会返回和数组相同大小的布尔数组，为True时，会返回峰值的坐标
local_maxi =feature.peak_local_max(distance, exclude_border = 0,min_distance = 12,indices=False,
                                   footprint=np.ones((10, 10)),labels=labels) #寻找峰值
markers = ndi.label(local_maxi)[0] #初始标记点
label_ =morphology.watershed(-distance, markers, mask=labels) #基于距离变换的分水岭算法

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 12))
axes = axes.ravel()
ax0, ax1, ax2, ax3 = axes

ax0.imshow(labels, cmap=plt.cm.gray)#, interpolation='nearest')
ax0.set_title("Original")
ax1.imshow(-distance, cmap=plt.cm.jet, interpolation='nearest')
ax1.set_title("Distance")
ax2.imshow(sm.dilation(markers,sm.square(10)), cmap=plt.cm.spectral, interpolation='nearest')
ax2.set_title("Markers")
ax3.imshow(label_, cmap=plt.cm.spectral, interpolation='nearest')
ax3.set_title("Segmented")

for ax in axes:
    ax.axis('off')

fig.tight_layout()
plt.show()

修正

可能是参数没有调好，我发现了很多找出的“峰值”靠的很近（8领域内）。我使用欧式距离剔除了选出的重复点。

# 修正
import math
err = []
for i in range(x_y.shape[0]):
    h1,w1 =  x_y[i][0],x_y[i][1]
    if i in err:
        continue
    for j in range(i+1,x_y.shape[0]):
        h2,w2 =  x_y[j][0],x_y[j][1]
        ab = math.sqrt(math.pow(abs(h2-h1), 2) + math.pow(abs(w2-w1), 2))
        if ab <= 10:
#             print 'error:' , x_y[i],' and ', x_y[j],'i,j = ',i,j
            err.append(j)
new_x_y = []
for i in range(len(x_y)):
    if i not in err:
        new_x_y.append(x_y[i])
print '一共有',len(new_x_y),'个圈'

---

本思路在钢筋图像上的效果：

---

总结

其实从这几张测试结果图可以看出这个算法存在的问题，因为核心是通过直方图的阈值分割来进行前景背景的分离。在面对这种比较现实的图像时背景比较复杂，因此对背景的识别很混乱。而且现实图片多变，需要调试的参数和优化的地方还有很多。可行改进方案是进行感兴趣区域提取，在这个区域进行统计。当然这里介绍的整个方法基于数字图像处理，并没有识别的步骤，要完成计数也可以采用机器学习的监督学习。