年前在狗東新買了一個茶濾，店家號稱有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:
continue
pixel_data[x, y] = WHITE
pixels = [(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 = tmp
connected.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 Image
BLACK = 0
WHITE = 1
def 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] = WHITE
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
def 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] = WHITE
area.append((x1, y1))
return area
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:
continue
pixel_data[x, y] = WHITE
pixels = [(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 = tmp
connected.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。
