用80行代码数1000个孔——《Python也可以》之四

年前在狗东新买了一个茶滤，店家号称有921个孔，我觉得好像有点“夸张”，总想证实一下。
肉眼肯定是数不过来了，所以我就想着要不写个程序来数吧，就拍了张照片放着。年前太忙了，然后又是过年，昨天终于有空了，便写了个80行的小脚本数了一下。

预处理

因为是白瓷，所以洞孔最好全黑，于是拍照的时候我就在下面放了个黑色的无纺布袋子，拍出来的效果还可以。

用微信传到电脑上，然后简单地用macOS自带的图像app处理一下，主要是去掉边缘无关元素和调整明暗。

预处理完成以后如下图：

概述

总体来说，这是一个很简单的任务，大概的流程是把图像二值化（黑白），只留下黑色的孔洞和白色的边缘，如下图。
可以看到有一些干扰因素，比如正下方的阴影和右上方的孔洞反光，需要做一下开运算去除掉干扰掉。这时候我们可以发现我们只要统计黑色区域（连通域）的数量就是孔洞的数量了，所以整体的代码框架如下：

# 打开
im = Image.open("WechatIMG1039.png")
w, h = im.size
# 二值化
im = im.convert('1')
pixel_data = im.load()
# 做几次开闭运算提升对比度
for i in range(3):# 用腐蚀算法去除杂点corrode(pixel_data, w, h)# 用膨胀算法填上空白expand(pixel_data, w, h)
# 统计连通域
print(f'holes count: {count_connected_area(pixel_data, w, h)}')

算法

腐蚀

腐蚀算法可以去掉孤点，是去噪的一种方法。在这里简单的把周边8邻点里黑色点不到2个的点定义为孤点，然后设置为白色去掉。

def corrode(pixel_data, w, h):for x in range(w):for y in range(h):if count_black8(pixel_data, x, y, w, h) < 2:# 设置为白色pixel_data[x, y] = WHITE

膨胀

膨胀算法是腐蚀算法的反操作，在这里简单的把周边8邻点里黑色点超过4个的点定义为设置为黑色，这样可以填上孔洞和扩展黑色区域。

def expand(pixel_data, w, h):for x in range(w):for y in range(h):if count_black8(pixel_data, x, y, w, h) > 4:# 设置为黑色pixel_data[x, y] = BLACK

统计连通域

所谓连通域就是一个黑色区块，我们从左上角开始遍历，发现一个黑色就扩展它的8邻点，然后再一步一步扩展开来，最后统计它的点阵数量。

def count_connected_area(pixel_data, w, h):connected = []for x in range(w):for y in range(h):if pixel_data[x, y] != BLACK:continuepixel_data[x, y] = WHITEpixels = [(x, y)]area = probe(pixel_data, x, y, w, h)while area:pixels.extend(area)tmp = []for px, py in area:tmp.extend(probe(pixel_data, px, py, w, h))area = tmpconnected.append(len(pixels))return reduce(lambda s, x: s + 1 if 30 < x < 500 else s, connected, 0)

我在这里使用的是Seed-Filling方法，比较慢，但胜在简单明了，反正我只是写个小脚本玩的，不在乎性能，就无所谓了。还有一种Two-Pass方法，比较快，但复杂些；也可以直接用OpenCV的算法，比自己手写可能要快个一二十倍。

全部代码

小脚本很简单，全部代码80行。

from functools import reduce
from PIL import ImageBLACK = 0
WHITE = 1def extend8(x, y, w, h):return ((x1, y1) for x1, y1 in ((x-1, y-1), (x, y-1), (x+1, y-1),(x-1, y),             (x+1, y),(x-1, y+1), (x, y+1), (x+1, y+1),) if x1 >= 0 and y1 >= 0 and x1 < w and y1 < h)def count_black8(pixel_data, x, y, w, h):return reduce(lambda s, x: s + 1 if pixel_data[x] == BLACK else s, extend8(x, y, w, h), 0)def corrode(pixel_data, w, h):for x in range(w):for y in range(h):if count_black8(pixel_data, x, y, w, h) < 2:# 设置为白色pixel_data[x, y] = WHITEdef expand(pixel_data, w, h):for x in range(w):for y in range(h):if count_black8(pixel_data, x, y, w, h) > 4:# 设置为黑色pixel_data[x, y] = BLACKdef probe(pixel_data, x, y, w, h):area = []for x1, y1 in extend8(x, y, w, h):if pixel_data[x1, y1] == BLACK:pixel_data[x1, y1] = WHITEarea.append((x1, y1))return areadef count_connected_area(pixel_data, w, h):connected = []for x in range(w):for y in range(h):if pixel_data[x, y] != BLACK:continuepixel_data[x, y] = WHITEpixels = [(x, y)]area = probe(pixel_data, x, y, w, h)while area:pixels.extend(area)tmp = []for px, py in area:tmp.extend(probe(pixel_data, px, py, w, h))area = tmpconnected.append(len(pixels))return reduce(lambda s, x: s + 1 if 30 < x < 500 else s, connected, 0)# 打开
im = Image.open("WechatIMG1039.png")
w, h = im.size
# 二值化
im = im.convert('1')
pixel_data = im.load()
# 做几次开运算提升对比度
for i in range(3):# 用腐蚀算法去除杂点corrode(pixel_data, w, h)# 用膨胀算法填上空白expand(pixel_data, w, h)
# 统计连通域
print(f'holes count: {count_connected_area(pixel_data, w, h)}')
# 目视查看是不是图像全都白了。
im.show()

最后计算出来的结果是922个孔，居然比店家说的还多1个，看来是个诚信店家，可以买它！

其它

这种方法主要的理论基础是数学形态学，是一门很通用的学科，不仅可以用来数孔洞，还可以用来数血小板、细菌等等，应用广泛。
另，发现我十几年前写的《用Python做图像处理》还是最简明好用的PIL教程，我用它重新学习了PIL。