第 9 章・形態學影像處理 · Pixels → Perception

你將學到

如何把二值影像看成像素座標的集合，並用一個小小的結構元素（structuring element）去探測它
兩個最基本的運算：侵蝕（erosion）與膨脹（dilation），以及它們組合成的開運算（opening）與閉運算（closing）
擊中擊不中轉換（hit-or-miss transform）如何找出角點、孤立點這類精確的像素圖樣
由這些基本運算組成的經典演算法：邊界、填洞、連通元件、凸包、細化、骨架與修剪
形態學重建（morphological reconstruction）：移除物件，卻不讓留下來的物件變形
灰階形態學：形態學梯度、用頂帽轉換校正不均勻照明，以及粒度測定（granulometry）

先看全貌

前幾章的濾波器大多是線性的：每個輸出像素都是輸入像素的加權和。形態學不一樣，它問的是關於形狀的是非題：「這個小形狀放得進這裡嗎？」「它有碰到這裡的東西嗎？」所有答案只靠最小值、最大值、交集和聯集算出來，完全不需要乘法。

為什麼重要：形態學是清理二值遮罩（不論來自第 10 章的閾值化，或來自神經網路）的標準工具，也用來量測物件（計數、量大小、找骨架）以及校正不均勻的背景。它的理論基礎數學形態學（mathematical morphology）源自 Georges Matheron 與 Jean Serra 對隨機集合與多孔材料幾何的研究 [2]、[3]。Haralick、Sternberg 與 Zhuang 的教學論文把它推廣到整個影像分析社群 [4]，而 Soille 的書至今仍是最實用的標準參考 [5]。本章依照 Gonzalez 與 Woods [1] 的主題順序，程式碼使用 scikit-image [12] 與 SciPy。

預備知識

白話說：二值影像就是一串「亮著」的像素位置；結構元素則是一個標了中心點的小模板，我們拿它在這串位置上滑動。

把影像當成集合

設一張二值影像中，前景像素值為 1、背景像素值為 0。我們用下面的集合描述前景：

A = \{ (x, y) \in \mathbb{Z}^2 \mid I(x, y) = 1 \},

其中 $I$ 是影像， $(x, y)$ 是整數像素座標。於是集合運算都有了影像上的意義：補集 $A^c$ 是背景，聯集 $A \cup C$ 是逐像素 OR，交集 $A \cap C$ 是逐像素 AND，差集 $A - C = A \cap C^c$ 則是把 $C$ 的部分從 $A$ 裡拿掉。

結構元素

結構元素（SE） $B$ 是另一個由位移量組成的小集合，並指定一個原點（通常在中心）。把它想成探針。常見的選擇有 $3 \times 3$ 方塊、十字（4 鄰域）、圓盤或線段。結構元素的形狀決定了運算會對哪些特徵起反應：圓盤不偏任何方向，水平線段只對水平結構有反應。scikit-image 把結構元素稱為 footprint。實作上結構元素就是一個小小的布林陣列，True 的位置就是它的成員；圖示中標成「不在乎」的位置，其實就是不屬於這個集合。

反射與平移

接下來每個定義都會用到兩個集合運算。 $B$ 的反射（reflection）是

\hat{B} = \{ w \mid w = -b,\ b \in B \},

也就是把每個位移 $b$ 換成 $-b$ ：結構元素繞原點旋轉 180°。 $B$ 以向量 $z = (z_1, z_2)$ 做的平移（translation）是

(B)_z = \{ c \mid c = b + z,\ b \in B \},

意思是把結構元素的原點放到像素 $z$ 上。圓盤、方塊這類對稱的結構元素滿足 $\hat{B} = B$ ，反射看不出差別；不對稱的結構元素就要小心。

import numpy as np
from skimage import morphology

B = morphology.disk(2)            # 5x5 disk; origin at the centre (2, 2)
B_hat = B[::-1, ::-1]             # reflection: flip both axes
L = np.array([[1, 1, 0],
              [0, 1, 0],
              [0, 0, 0]], bool)   # an asymmetric element
print(L[::-1, ::-1].astype(int))  # its reflection points the other way

侵蝕與膨脹

白話說：侵蝕只在「整個模板都放得進形狀裡」的地方保留像素；膨脹則只要「模板碰到形狀」就把像素點亮。

侵蝕

$A$ 被 $B$ 侵蝕的結果是

A \ominus B = \{ z \mid (B)_z \subseteq A \},

也就是平移後的結構元素完全落在前景內的所有位置 $z$ 。等價寫法是 $A \ominus B = \{ z \mid (B)_z \cap A^c = \emptyset \}$ ：結構元素不能和任何背景像素重疊。侵蝕會讓物件大約縮小一個結構元素半徑，刪掉比結構元素更細的部分，並把被細橋連在一起的物件切開。若 $B$ 包含原點，侵蝕是反擴張（anti-extensive）的： $A \ominus B \subseteq A$ 。

膨脹

$A$ 被 $B$ 膨脹的結果是

A \oplus B = \{ z \mid (\hat{B})_z \cap A \neq \emptyset \},

也就是反射後的結構元素至少碰到一個前景像素的所有位置。另一個等價寫法是把 $B$ 的複本放在每個前景像素上再取聯集： $A \oplus B = \bigcup_{a \in A} (B)_a$ 。膨脹讓物件變大、接起比結構元素窄的縫隙、補上小洞。若 $B$ 包含原點，膨脹是擴張（extensive）的： $A \subseteq A \oplus B$ 。和侵蝕不同，膨脹滿足交換律： $A \oplus B = B \oplus A$ 。

膨脹公式中的反射，讓膨脹成為真正的 Minkowski 和，也和卷積的「翻轉」一致。對稱結構元素則完全沒差。

對偶性

侵蝕與膨脹在補集與反射之下互為對偶（duality）：

(A \ominus B)^c = A^c \oplus \hat{B}, \qquad (A \oplus B)^c = A^c \ominus \hat{B}.

侵蝕前景，等於用反射後的結構元素膨脹背景，再把結果反過來。證明只要一行： $z \notin A \ominus B$ 恰好表示 $(B)_z$ 碰到 $A^c$ ，而這正是 $z \in A^c \oplus \hat{B}$ 的定義（因為 $\hat{\hat{B}} = B$ ）。寫程式時，必須一致地處理影像邊界，對偶性才會精確成立：對侵蝕而言影像外面要當成前景，對補集的膨脹而言外面才會是背景。

import numpy as np
from scipy import ndimage as ndi
from skimage import data, morphology

A = ~data.horse()                 # True = horse pixels
B = morphology.disk(3)

eroded  = ndi.binary_erosion(A, structure=B)
dilated = ndi.binary_dilation(A, structure=B)
print(A.sum(), eroded.sum(), dilated.sum())     # erosion shrinks, dilation grows

# Duality check. border_value=1 treats outside pixels as foreground for the
# erosion, which matches background-outside for the dilation of ~A.
lhs = ~ndi.binary_erosion(A, structure=B, border_value=1)
rhs = ndi.binary_dilation(~A, structure=B[::-1, ::-1])
print("duality holds:", np.array_equal(lhs, rhs))   # True

開運算與閉運算

白話說：開運算是「先侵蝕、再膨脹回來」，能去掉小斑點和細長部分，大形狀則幾乎不變。閉運算是「先膨脹、再侵蝕回來」，能補上小洞和窄裂縫。

定義

$A$ 以 $B$ 做開運算，就是用同一個結構元素先侵蝕再膨脹：

A \circ B = (A \ominus B) \oplus B,

閉運算則是先膨脹再侵蝕：

A \bullet B = (A \oplus B) \ominus B.

單獨做侵蝕會丟掉細節，也會讓所有東西縮小；第二步則把倖存部分的大小補回來。這兩者也互為對偶： $(A \bullet B)^c = A^c \circ \hat{B}$ 、 $(A \circ B)^c = A^c \bullet \hat{B}$ 。對前景做閉運算，就是對背景做開運算。

幾何解釋：滾球

開運算有一個很直觀的圖像：它等於所有能放進 $A$ 的 $B$ 平移複本的聯集，

A \circ B = \bigcup \{ (B)_z \mid (B)_z \subseteq A \}.

想像一顆球（結構元素）在前景裡面滾動，而且必須整顆待在裡面。球滾得到的像素都保留；比球還尖的角、細頸和比球小的斑點都會消失。閉運算則是同樣的圖像換到外面：讓球在背景中滾動，凡是球到不了的背景像素（窄縫、小洞、又深又細的凹口）都變成前景。

性質

開運算與閉運算在嚴格意義上都是形態學濾波器（morphological filter）：它們是遞增的，而且是冪等的 [5]。

開運算是反擴張的： $A \circ B \subseteq A$ ，只會移除像素。
閉運算是擴張的： $A \subseteq A \bullet B$ ，只會加入像素。
遞增性： 若 $C \subseteq D$ ，則 $C \circ B \subseteq D \circ B$ ，閉運算亦同。輸入變大，輸出絕不會變小。
冪等性： $(A \circ B) \circ B = A \circ B$ ， $(A \bullet B) \bullet B = A \bullet B$ 。同樣的開運算做兩次不會再改變任何東西，因為第一次留下來的部分本來就都是球放得進去的地方。

冪等性正是開、閉運算和侵蝕、膨脹的分水嶺：重複侵蝕會一直縮下去，重複開運算則會停下來。圖 9.1 把四種運算套用在一張帶雜訊的剪影上。開運算去掉背景中的白點；閉運算去掉馬身上的黑點；先開再閉則兩種都去掉，這是非常常用的形態學雜訊濾波器。

六個面板：帶雜訊的馬剪影，以及以半徑 2 圓盤做的侵蝕、膨脹、開運算、閉運算、先開後閉。 — 圖 9.1 — 對翻轉了 4% 像素的馬剪影，以半徑 2 的圓盤做二值侵蝕、膨脹、開運算與閉運算。先開後閉可以同時去除兩種斑點，並保留輪廓。

import numpy as np
from skimage import data, morphology

A = ~data.horse()
B = morphology.disk(4)
op = morphology.binary_opening(A, B)
cl = morphology.binary_closing(A, B)
print("opening is anti-extensive:", not (op & ~A).any())
print("closing is extensive:     ", not (A & ~cl).any())
print("opening is idempotent:    ", np.array_equal(morphology.binary_opening(op, B), op))

擊中擊不中轉換

白話說：擊中擊不中轉換（HMT）就是二值圖樣的模板比對。你指定哪些像素一定要是前景、哪些一定要是背景，它就只在兩個條件同時成立的位置亮起來。

使用兩個不相交的結構元素： $B_1$ 是必須落在前景的像素， $B_2$ 是必須落在背景的像素。 $A$ 的 HMT 為

A \circledast B_{1,2} = \{ z \mid (B_1)_z \subseteq A \ \text{and}\ (B_2)_z \subseteq A^c \} = (A \ominus B_1) \cap (A^c \ominus B_2).

其中 $B_{1,2}$ 代表這一對 $(B_1, B_2)$ ； $A \ominus B_1$ 找出前景圖樣吻合的位置， $A^c \ominus B_2$ 找出背景圖樣吻合的位置。既不在 $B_1$ 也不在 $B_2$ 的像素是不在乎的位置。習慣上會把這一對畫成一個小網格：1 代表必須是前景、0 代表必須是背景、空白代表不在乎。

幾個實用的模板：

孤立點： $B_1$ 是中心像素， $B_2$ 是它的 8 個鄰居。
左上角： $B_1$ 是中心加上右邊與下面的鄰居， $B_2$ 是上面與左邊的鄰居。
單像素線段的端點： 中心與恰好一個鄰居是前景，其餘都是背景。每個方向需要一個模板（共 8 個旋轉）。

圖 9.2 套用了前兩種模板。注意數位圓盤的階梯狀輪廓在像素層級上確實有好幾個「左上角」：HMT 是逐字比對，不是幾何意義上的角。

三個面板：含矩形、L 形、圓盤與兩個單像素的二值場景；以紅圈標出偵測到的左上角；以橘圈標出偵測到的孤立像素。 — 圖 9.2 — 把擊中擊不中轉換當成精確的圖樣偵測器。中：左上角（包含數位圓盤上的小階梯）。右：孤立的前景像素。

import numpy as np
from scipy import ndimage as ndi

A = np.zeros((7, 9), bool)
A[1:5, 1:6] = True                # a rectangle
A[5, 7] = True                    # an isolated pixel
hit  = np.array([[0, 0, 0], [0, 1, 1], [0, 1, 0]], bool)   # must be foreground
miss = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]], bool)   # must be background
corners = ndi.binary_hit_or_miss(A, structure1=hit, structure2=miss)
print(np.argwhere(corners))       # [[1 1]]: the upper-left corner only

HMT 是下面細化、粗化與凸包演算法的基本積木。

一些基本的形態學演算法

白話說：有了侵蝕、膨脹和 HMT，就能組出許多實用工具。大部分都是一個短迴圈：反覆做同一件事，直到沒有任何變化為止。

邊界擷取

一個集合的內邊界，就是侵蝕會拿掉的部分：

\beta(A) = A - (A \ominus B),

其中 $B$ 取 $3 \times 3$ 方塊（得到 8 連通的邊界）或十字（得到較粗、4 連通的邊界）。 $B$ 越大，邊界越粗。圖 9.3 第一格就是馬的 $\beta(A)$ 。

填洞

洞（hole）是被前景完全包圍的背景區域。只要在洞裡給一個種子像素，就能用膨脹讓種子長大，但每一步都把它限制在背景內：

X_k = (X_{k-1} \oplus B) \cap A^c, \qquad k = 1, 2, 3, \ldots

其中 $X_0$ 只含種子， $B$ 是 4 連通的十字，和 $A^c$ 取交集讓生長在洞壁停下。當 $X_k = X_{k-1}$ 時停止， $X_k \cup A$ 就是填好的形狀。這種「先膨脹、再裁切」的模式稱為條件膨脹（conditional dilation），後面的測地膨脹會再看到它。 $B$ 用十字（而不是方塊）很重要：前景牆是 8 連通時，4 連通的背景填充就不會從斜角漏出去。

import numpy as np
from scipy import ndimage as ndi

def fill_hole(A, seed):
    """Grow X from a seed inside one hole: X_k = dilate(X_{k-1}) AND not A."""
    B = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], bool)   # 4-connected cross
    X = np.zeros_like(A); X[seed] = True
    while True:
        X_next = ndi.binary_dilation(X, structure=B) & ~A
        if np.array_equal(X_next, X):
            return A | X
        X = X_next

yy, xx = np.mgrid[-12:13, -12:13]
ring = (np.hypot(yy, xx) <= 10) & (np.hypot(yy, xx) > 6)   # boolean ring
print(ring.sum(), fill_hole(ring, (12, 12)).sum())        # 204 -> 317

擷取連通元件

把角色對調：要擷取 $A$ 中包含某個種子像素的連通元件，就讓種子膨脹，但限制在前景內：

X_k = (X_{k-1} \oplus B) \cap A.

$B$ 取 $3 \times 3$ 方塊時得到種子的 8 連通元件；取十字時得到 4 連通元件。對每個尚未標記的種子重複這個步驟，就能標記所有元件。實際的函式庫使用更快的一次或兩次掃描標記演算法，但結果相同：scipy.ndimage.label(A, structure=np.ones((3, 3))) 會傳回標記影像與元件數。

凸包

若一個集合中任兩點之間的線段都留在集合內，它就是凸的。凸包（convex hull） $C(A)$ 是包含 $A$ 的最小凸集合。純形態學的近似做法使用四個 HMT 模板 $B^i$ （ $i = 1, \ldots, 4$ ），各自偵測「某一側（左、上、右、下）有前景」的背景像素。從 $X_0^i = A$ 開始迭代

X_k^i = (X_{k-1}^i \circledast B^i) \cup X_{k-1}^i

直到收斂得到 $D^i$ ，再取 $C(A) = \bigcup_{i=1}^4 D^i$ 。每個模板會在一個方向上持續補像素，直到該方向不再有凹陷。結果是相對於這四個方向的凸包，可能比真正的歐氏凸包更大，因此實作上常把生長限制在外接矩形內。實務上，skimage.morphology.convex_hull_image 直接計算前景像素座標的多邊形凸包再點陣化，速度更快，也更接近歐氏凸包。

細化

細化（thinning）有選擇地移除邊界像素，但不把形狀拆開，直到只剩一個像素寬的線。用模板 $B$ 做一次細化是

A \otimes B = A - (A \circledast B),

也就是刪掉 HMT 比對到的像素。完整的細化會輪流套用一串模板 $\{B^1, B^2, \ldots, B^n\}$ （通常是「邊緣像素」圖樣的 8 個旋轉），

A \otimes \{B\} = (( \cdots ((A \otimes B^1) \otimes B^2) \cdots ) \otimes B^n),

並重複整個循環，直到沒有像素改變。這些模板的設計保證：刪掉被比對到的像素，不會讓形狀斷開，也不會刪掉線段的端點。Zhang 與 Suen 的平行細化演算法是這個想法廣為使用的高效版本 [6]；skimage.morphology.skeletonize 對 2-D 影像預設就用它，skimage.morphology.thin 則實作了另一種相關的平行細化。

粗化

粗化（thickening）是細化的對偶，它加入被模板比對到的背景像素：

A \odot B = A \cup (A \circledast B).

實務上通常不另外設計粗化模板，而是對背景做細化，再取補集。把背景細化到收斂，會留下一張位於物件之間正中央的細線網，所以粗化是讓物件長大到彼此相接、卻永遠不會合併的經典方法。

骨架

骨架（skeleton） $S(A)$ 是一組沿著形狀中央延伸的細線。一種正式定義使用最大圓盤：若某點是一個放得進 $A$ 、而且不被其他放得進的圓盤包含的圓盤中心，它就在骨架上。這種圓盤至少在兩處碰到邊界。骨架有純形態學的公式：

S(A) = \bigcup_{k=0}^{K} S_k(A), \qquad S_k(A) = (A \ominus kB) - \big( (A \ominus kB) \circ B \big),

其中 $A \ominus kB$ 表示以 $B$ 連續侵蝕 $A$ 共 $k$ 次， $K$ 是侵蝕變成空集合之前的最後一個 $k$ 。每個 $S_k(A)$ 收集侵蝕深度 $k$ 時的「山脊」點：做了 $k$ 次侵蝕後它們還在，但在那個深度做開運算就撐不住。這個分解是可逆的：

A = \bigcup_{k=0}^{K} \big( S_k(A) \oplus kB \big),

所以把每個 $S_k$ 連同索引 $k$ 存起來，就是一種無失真的形狀編碼。缺點是這種骨架不保證連通。細化可以得到連通、一個像素寬的骨架；中軸（medial axis，由距離轉換算出）則在每個骨架點上附帶最大圓盤的半徑，如圖 9.3 最右格。

四個面板：馬的邊界；淺藍色凸包包住深色的馬；紅色的細化骨架疊在馬上；依到邊界距離上色的中軸。 — 圖 9.3 — 由侵蝕與 HMT 組成的演算法：內邊界、凸包、細化骨架，以及以顏色表示局部半寬的中軸。

修剪

骨架對邊界雜訊很敏感：輪廓上一個小凸起就會長出一根短短的雜枝（spur），圖 9.3 中馬尾附近就看得到。修剪（pruning）會移除長度短於 $L$ 的雜枝：

把骨架的端點（恰有一個 8 鄰居的像素）移除 $L$ 次。所有長度不超過 $L$ 的雜枝都會消失，但真正的分支末端也會少掉 $L$ 個像素。
找出縮短後骨架的端點。
從這些端點做 $L$ 次條件膨脹（限制在原骨架內），把主要分支長回原本長度。
取縮短後的骨架與長回來的末端之聯集。

import numpy as np
from scipy import ndimage as ndi
from skimage import data, morphology

def endpoints(S):
    n = ndi.convolve(S.astype(int), np.ones((3, 3), int), mode="constant") - S
    return S & (n == 1)            # foreground pixels with exactly one neighbour

def prune_tips(S, L):
    """Step 1 of pruning: strip L layers of end points."""
    S = S.copy()
    for _ in range(L):
        S &= ~endpoints(S)
    return S

skel = morphology.skeletonize(~data.horse())
print(skel.sum(), prune_tips(skel, 10).sum())

形態學重建

白話說：重建讓你可以說「把含有這顆種子的物件整個留下來」，而且拿回來的是一模一樣的物件，不是被磨圓的複本。它用兩張影像：標記（marker）說明從哪裡開始，遮罩（mask）說明最多能長到哪裡。

測地膨脹與測地侵蝕

設 $F$ 為標記、 $G$ 為遮罩，且 $F \subseteq G$ 。 $F$ 相對於 $G$ 的大小為 1 的測地膨脹（geodesic dilation）是

D_G^{(1)}(F) = (F \oplus B) \cap G,

其中 $B$ 是小的結構元素（通常是 $3 \times 3$ 方塊）。標記長大一步，但不能離開遮罩。大小為 $n$ 就是重複 $n$ 次： $D_G^{(n)}(F) = D_G^{(1)}\big(D_G^{(n-1)}(F)\big)$ 。「測地」指的是沿著留在 $G$ 內的路徑量距離，就像在建築物裡走路，而不是穿牆。

對偶的大小為 1 的測地侵蝕是

E_G^{(1)}(F) = (F \ominus B) \cup G,

此時 $F \supseteq G$ ：標記會縮小，但不會縮到比遮罩還小。

以膨脹重建與以侵蝕重建

把測地膨脹一直迭代到不再變化，就得到遮罩 $G$ 由標記 $F$ 做的以膨脹重建（reconstruction by dilation）：

R_G^D(F) = D_G^{(k)}(F) \quad \text{with } k \text{ such that } D_G^{(k)}(F) = D_G^{(k+1)}(F).

在二值影像中， $R_G^D(F)$ 恰好是 $G$ 裡所有至少含一個標記像素的連通元件之聯集。同理，以侵蝕重建 $R_G^E(F)$ 是把測地侵蝕迭代到穩定。Vincent 的論文給出了標準的高效演算法（以佇列為基礎，每個像素只被拜訪有限次，而不是一次次掃過整張影像），以及許多應用 [7]；skimage.morphology.reconstruction 就是採用這個方法。

以重建做開運算

一般的開運算會移除小物件，但也會把大物件的角磨圓、把細的部分擦掉。以重建做開運算（opening by reconstruction）避開了這個問題：

O_R^{(n)}(F) = R_F^D\big(F \ominus nB\big),

其中 $F$ 是輸入影像， $F \ominus nB$ （以 $B$ 侵蝕 $n$ 次）是標記，輸入影像本身就是遮罩。侵蝕決定哪些物件能留下（至少放得進一個結構元素的那些）；重建則把每個倖存者恢復成原本的形狀。比較圖 9.4 的「一般開運算」與「以重建做開運算」：方塊上的細手臂和尖角都完整回來了。以重建做閉運算是其對偶，由膨脹與以侵蝕重建組成。

自動填洞

重建也能在不給種子的情況下一次填滿所有洞。洞就是從影像邊框走不到的背景區域。因此建立一個標記：在邊框上取影像的補值，其他地方為零，

F(x, y) = \begin{cases} 1 - I(x, y) & \text{if } (x, y) \text{ is on the image border,} \\ 0 & \text{otherwise,} \end{cases}

由它重建背景，再取補集：

H = \big[ R_{I^c}^D(F) \big]^c .

$R_{I^c}^D(F)$ 是從邊框走得到的背景；其餘的部分 $H$ 就是前景加上它的洞。

清除邊界物件

碰到影像邊界的物件通常不完整，量測時應該排除。用影像在邊框上的值當標記，

F(x, y) = \begin{cases} I(x, y) & \text{if } (x, y) \text{ is on the border,} \\ 0 & \text{otherwise,} \end{cases} \qquad X = I - R_I^D(F),

重建會找回所有碰到邊框的物件，從影像中減掉它們，就只剩內部的物件。

六個面板：合成的二值場景；侵蝕後的標記；一般開運算把角磨圓並移除細手臂；以重建做開運算完整恢復大形狀；以紅色標示被填的洞；移除碰到邊界的物件。 — 圖 9.4 — 形態學重建。以半徑 6 圓盤侵蝕得到的標記挑出大物件，重建把它們（包括細手臂）完整恢復。下排：自動填洞與清除邊界物件。

import numpy as np
from skimage import morphology, segmentation

def open_by_reconstruction(A, footprint):
    marker = morphology.binary_erosion(A, footprint)
    return morphology.reconstruction(marker.astype(np.uint8), A.astype(np.uint8),
                                     method="dilation").astype(bool)

def fill_holes_auto(A):
    """Background reachable from the image frame is background; the rest is hole."""
    Ac = ~A
    marker = np.zeros_like(Ac)
    marker[0, :], marker[-1, :] = Ac[0, :], Ac[-1, :]
    marker[:, 0], marker[:, -1] = Ac[:, 0], Ac[:, -1]
    reached = morphology.reconstruction(marker.astype(np.uint8), Ac.astype(np.uint8),
                                        method="dilation").astype(bool)
    return ~reached

# Library shortcuts: scipy.ndimage.binary_fill_holes(A), segmentation.clear_border(A)

二值運算總整理

下表整理本章的二值運算。 $A$ 為輸入集合、 $B$ 為結構元素、 $F$ 為標記、 $G$ 為遮罩，並假設 $B$ 包含原點。

運算	公式	作用
平移	$(B)_z = \{ b + z \}$	把結構元素原點移到 $z$
反射	$\hat{B} = \{ -b \}$	結構元素旋轉 180°
補集	$A^c$	前景與背景互換
侵蝕	$A \ominus B$	保留 $B$ 放得進去的位置；縮小、去除細部
膨脹	$A \oplus B$	保留 $\hat{B}$ 碰得到 $A$ 的位置；變大、接起縫隙
開運算	$(A \ominus B) \oplus B$	去除斑點、細頸、尖角；冪等
閉運算	$(A \oplus B) \ominus B$	補上小洞與窄縫；冪等
擊中擊不中	$(A \ominus B_1) \cap (A^c \ominus B_2)$	找出精確的前景／背景圖樣
邊界	$A - (A \ominus B)$	一個像素寬的內輪廓
填洞	$X_k = (X_{k-1} \oplus B) \cap A^c$	填滿含種子的那個洞
連通元件	$X_k = (X_{k-1} \oplus B) \cap A$	擷取含種子的元件
凸包	反覆 HMT 的聯集	填平凹陷
細化	$A - (A \circledast B)$ ，套用一串模板	剝成一個像素寬的骨架，保持連通
粗化	$A \cup (A \circledast B)$	細化的對偶
骨架	$\bigcup_k (A \ominus kB) - (A \ominus kB) \circ B$	中央線；配合索引 $k$ 可還原
修剪	移除端點，再條件式長回	去除短雜枝
測地膨脹	$(F \oplus B) \cap G$	在遮罩內長大一步
測地侵蝕	$(F \ominus B) \cup G$	在遮罩之上縮小一步
以膨脹重建	穩定時的 $D_G^{(k)}(F)$	$G$ 中被 $F$ 碰到的元件
以重建做開運算	$R_F^D(F \ominus nB)$	移除小物件，其他物件保持原樣
自動填洞	$[R_{I^c}^D(F)]^c$	不需種子就填滿所有洞
清除邊界物件	$I - R_I^D(F)$	移除碰到邊界的物件

灰階形態學

白話說：在灰階影像裡，「放得進去」變成「局部最小值」，「碰得到」變成「局部最大值」。把影像想成一片地形，亮度就是高度。侵蝕把每個像素降到模板底下最低的地面；膨脹則把它升到最高處。

灰階侵蝕與膨脹

使用平坦結構元素 $b$ （一組高度都是零的位移）時，影像 $f$ 的灰階侵蝕與膨脹是

[f \ominus b](x, y) = \min_{(s, t) \in b} f(x + s,\ y + t), \qquad [f \oplus b](x, y) = \max_{(s, t) \in b} f(x - s,\ y - t),

其中 $(s, t)$ 走遍結構元素的所有位移。膨脹中的負號一樣是反射。用在二值影像上，這兩個式子會完全退化成集合的定義。侵蝕讓影像變暗、亮的特徵縮小、暗的特徵變寬；膨脹則相反。

非平坦結構元素 $b_N$ 另外帶有高度 $b_N(s, t)$ ：

[f \ominus b_N](x, y) = \min_{(s, t) \in b_N} \big\{ f(x + s,\ y + t) - b_N(s, t) \big\}, \qquad [f \oplus b_N](x, y) = \max_{(s, t) \in b_N} \big\{ f(x - s,\ y - t) + b_N(s, t) \big\}.

實務上很少用非平坦結構元素：結果會受影像的亮度尺度影響，而且侵蝕可能低於影像的最小值。對偶性用灰階補數 $f^c = -f$ （對 $L$ 階影像則為 $L - 1 - f$ ）依然成立： $(f \ominus b)^c = f^c \oplus \hat{b}$ 。

開運算與閉運算

定義不變： $f \circ b = (f \ominus b) \oplus b$ 、 $f \bullet b = (f \oplus b) \ominus b$ 。幾何圖像是一顆球從亮度曲面下方往上推：開運算就是球在不超過 $f$ 的前提下能到達的最高曲面。比球窄的亮峰會被削平，比球寬的部分則不受影響。閉運算是球從上方往下壓：比球窄的暗谷會被填平。各種性質（反擴張／擴張、遞增、冪等）只要把 $\subseteq$ 換成 $\le$ 就照樣成立。

形態學平滑

開運算壓掉小的亮細節，閉運算壓掉小的暗細節，所以用同一個結構元素先開後閉，可以同時去除這兩種小結構，同時保留較大區域的銳利邊緣。這就是圖 9.1 雜訊濾波器的灰階版本。相關的交替序列濾波（alternating sequential filtering）則用逐漸變大的結構元素輪流做開、閉運算，避免一次用大結構元素造成的瑕疵。

形態學梯度

膨脹減去侵蝕，可以衡量局部對比：

g = (f \oplus b) - (f \ominus b).

在均勻區域內最大值與最小值幾乎相等，所以 $g \approx 0$ ；跨過邊緣時兩者相差約一個邊緣高度。和 Sobel 導數不同， $g$ 一定非負，而且（對稱 $b$ 時）不偏任何方向。只取一邊則得到內梯度 $f - (f \ominus b)$ 或外梯度 $(f \oplus b) - f$ ，分別把邊緣響應放在較亮區域的內側或外側。

cameraman 影像的六個面板：輸入、灰階侵蝕、膨脹、先開後閉的平滑、形態學梯度與內梯度。 — 圖 9.5 — 以半徑 3 的平坦圓盤做灰階形態學：侵蝕（局部最小）、膨脹（局部最大）、先開後閉的平滑、形態學梯度，以及內梯度。

from skimage import data, morphology, util

f = util.img_as_float(data.camera())
b = morphology.disk(3)                       # flat structuring element
ero  = morphology.erosion(f, b)              # local minimum
dil  = morphology.dilation(f, b)             # local maximum
grad = dil - ero                             # morphological gradient (>= 0)
smooth = morphology.closing(morphology.opening(f, b), b)

頂帽與底帽轉換

用比亮物件更大的結構元素做開運算，會把物件削掉，但保留緩慢變化的背景。因此把影像減去開運算結果，就只剩下物件。這就是頂帽轉換（top-hat，又稱 white top-hat）：

T_{\text{hat}}(f) = f - (f \circ b).

它的對偶是底帽轉換（bottom-hat，又稱 black top-hat），用來擷取亮背景上的暗物件：

B_{\text{hat}}(f) = (f \bullet b) - f.

兩者都非負。最主要的用途是在閾值化之前校正不均勻照明。在圖 9.6 中，一張從單側打光的合成顆粒影像讓全域 Otsu 閾值（第 10 章）完全失效：亮側的背景比暗側的顆粒還亮。用半徑 8（大於任何顆粒）的圓盤做頂帽轉換後，背景變平，單一閾值就幾乎完美地分出顆粒。唯一的參數是結構元素大小：它必須大於物件，又要小於背景變化的尺度。

六個面板：背景往右變亮的顆粒影像；估計背景的開運算；背景變平的頂帽結果；對輸入做 Otsu 閾值（右側失敗）；對頂帽結果做 Otsu 閾值（乾淨）；以及輸入、開運算、頂帽在同一列的剖面。 — 圖 9.6 — 用頂帽轉換校正不均勻照明。以大於顆粒的圓盤做開運算可估計背景，相減後背景變平。IoU 是和已知的顆粒遮罩比較的結果。

拖動滑桿比較輸入與頂帽結果：

輸入頂帽（圓盤 r=8）

from skimage import data, filters, morphology, util

coins = util.img_as_float(data.coins())
th = morphology.white_tophat(coins, morphology.disk(30))   # bright objects
bh = morphology.black_tophat(coins, morphology.disk(30))   # dark objects
mask = th > filters.threshold_otsu(th)

粒度測定

粒度測定（granulometry）不需要分割，就能量出影像中顆粒的大小分布。用越來越大的結構元素 $r$ 對影像做開運算，每次開運算都會移除比 $b_r$ 小的亮顆粒。記錄剩下的總亮度：

V(r) = \sum_{x, y} [f \circ b_r](x, y), \qquad r = 0, 1, 2, \ldots,

其中 $b_r$ 是半徑 $r$ 的圓盤， $b_0$ 是單一像素（所以 $V(0)$ 就是 $f$ 的總和）。 $V(r)$ 絕不會增加，因為越大的圓盤能放進的位置越少。它的負離散導數

PS(r) = V(r - 1) - V(r)

是一種大小直方圖，稱為圖樣頻譜（pattern spectrum）[8]：在 $r$ 處的值很大，表示有很多影像「質量」位於半徑約為 $r$ 的特徵中。圖 9.7 的場景有兩種顆粒大小，頻譜也對應出現兩個峰。以大小為索引的一族開運算、以及「篩子」的詮釋，可以追溯到 Matheron [2]；Maragos 則把圖樣頻譜發展成多尺度形狀描述子 [8]。

左：暗背景上有大小兩種亮圓盤。中：開運算後總和對半徑的曲線，分兩階段下降。右：在半徑 5 與 12 有峰值的長條圖。 — 圖 9.7 — 半徑 5 與 12 兩種圓盤的粒度測定。粒度曲線下降兩次，圖樣頻譜在兩個半徑處出現峰值。

import numpy as np
from skimage import morphology

def granulometry(f, radii):
    total = f.sum()
    vol = np.array([morphology.opening(f, morphology.disk(r)).sum() / total
                    for r in radii])
    spectrum = -np.diff(vol)                 # loss from radius r-1 to r
    return vol, spectrum

紋理分割

形態學也能把「紋理元素大小不同」的區域分開，例如一片小顆粒旁邊接著一片大顆粒。做法沿用粒度測定的大小邏輯：

用比小顆粒之間縫隙更大的結構元素做閉運算（暗背景上的亮顆粒，閉運算會填平這些暗縫）。小顆粒區變成一塊幾乎均勻的亮區，而大顆粒之間的縫比結構元素寬，仍然分開。
用比大顆粒更大、但比合併後的小顆粒區更小的結構元素對結果做開運算，擦掉孤立的大顆粒，只留下合併後的小顆粒區。
對結果取形態學梯度（或做閾值）。兩種紋理的交界就成為一條明顯的邊緣。

確切的結構元素大小取決於資料中的顆粒與間距大小，可以先用粒度測定量出來。

灰階重建

所有重建的定義都能搬過來，只要把 $\cap$ 換成逐點最小值 $\wedge$ 、 $\cup$ 換成最大值 $\vee$ 。標記影像 $f$ 在遮罩影像 $g$ （ $f \le g$ ）之下的測地膨脹為

D_g^{(1)}(f) = (f \oplus b) \wedge g,

以膨脹重建就是把它迭代到穩定。灰階重建支撐了好幾種實用的濾波器 [7]：

以重建做開運算 $R_f^D(f \ominus nb)$ ：移除比結構元素小的亮細節，其他部分則保持原本形狀。以重建做頂帽 $f - R_f^D(f \ominus nb)$ 比一般頂帽乾淨，因為估計出的背景不會有被磨圓的「肩膀」。
區域極大值與 h-dome 轉換：以 $f - h$ 為標記重建 $f$ 。重建會把每個高度不超過 $h$ 的峰削平；從 $f$ 減去它，就只剩下「圓頂」（峰），每個高度至多 $h$ 。對圓頂做閾值，不論局部背景多亮都能偵測出亮斑。

from skimage import data, morphology, util
import numpy as np

f = util.img_as_float(data.camera())
h = 0.1
background = morphology.reconstruction(np.clip(f - h, 0, None), f, method="dilation")
domes = f - background                     # bright peaks, each at most h tall

現代觀點

形態學雖然歷史悠久，但遠未結束。有五條研究脈絡值得認識。

理論基礎與教學文獻。 Matheron 的隨機集合理論 [2] 與 Serra 的 Image Analysis and Mathematical Morphology [3] 建立了整套代數：侵蝕與膨脹是 Minkowski 運算、開運算是「篩子」，而格論（lattice theory）的觀點後來把一切從集合推廣到灰階函數。Haralick、Sternberg 與 Zhuang 1987 年的教學論文 [4] 至今仍是工程師入門二值與灰階形態學最清楚的文獻之一：它定義運算子、證明基本性質，並連結到實用演算法。Soille 的 Morphological Image Analysis [5] 是實務者的參考書，涵蓋測地轉換、重建、分水嶺、粒度測定與高效實作；函式庫文件寫得太簡略時，就該翻它。

重建與高效演算法。 Vincent 1993 年的論文 [7] 讓重建變得實用。重點在於：天真的「迭代到穩定」迴圈可以換成以佇列為基礎、成本大致隨影像大小線性成長的演算法，於是重建成為填洞、區域極值與 h-dome 的日常工具。

連通運算子與元件樹。 若一個濾波器只能移除或合併平坦區（數值相同的連通區域），而絕不產生新的輪廓，它就是連通（connected）的。以重建做開運算是最簡單的例子。Salembier 與 Wilkinson 的回顧論文 [10] 介紹了整個家族：面積開運算、屬性濾波器、levelings，以及以樹為基礎的濾波。Breen 與 Jones 提出屬性開運算與屬性細化（attribute openings and thinnings）[9]，依據屬性（面積、伸長度、慣性矩）而非「結構元素放不放得進」來決定保留或移除一個連通元件。高效實作把影像表示成最大樹（max-tree）或最小樹（min-tree），也就是所有閾值集合之連通元件構成的樹，再以修剪樹的方式濾波。Carlinet 與 Géraud 在統一框架下比較了多種最大樹建構演算法，並給出挑選演算法的決策樹 [11]。scikit-image 提供了 area_opening、area_closing、diameter_opening 與 max_tree。

大小分布。 Maragos 的圖樣頻譜 [8] 把粒度測定形式化為多尺度形狀描述子，至今仍用於紋理與顆粒大小分析；以屬性為基礎的粒度測定 [9] 則把它推廣到大小以外的形狀屬性。

形態學遇上深度學習。 深度學習從兩個方向改變了形態學。

網路內部的形態學。 最大池化（max pooling）本來就是平坦灰階膨脹再加上降取樣，所以兩者的連結很自然。有好幾個團隊設計了結構元素可以用梯度下降學習的網路層。Mondal 等人研究由可學習的膨脹與侵蝕神經元組成的網路，證明兩個「膨脹與侵蝕後接線性組合」的區塊就能逼近任意連續函數 [13]。Nogueira 等人提出 DeepMorphNet，以可學習的形態學運算取代卷積，並比較它學到的特徵與一般卷積網路的差異 [14]。實務上的難處是 min 與 max 的梯度很稀疏。Kirszenberg 等人 [16] 與 Hermary 等人 [15] 用 max 與 min 的可微近似（銳利程度由可學習參數控制的廣義平均與 soft-max 函數）建構平滑的形態學層，改進了較早的 p-convolution 層，並探討這些層何時真的能學回侵蝕與膨脹。這些仍是研究工具：形態學層並不是主流視覺骨幹網路的標準元件。

網路周邊的形態學。 在分割任務中，形態學一直被大量用作預測遮罩的後處理：用開運算去掉偽陽性斑點、用閉運算或填洞修補破碎的物件、remove_small_objects、清除邊界物件，以及用骨架量測管狀結構。它也出現在損失函數裡：clDice [17] 以反覆的 min 與 max 池化（軟侵蝕與軟膨脹）算出可微的「軟骨架」，再用預測骨架與真實骨架的重疊程度當作損失，鼓勵血管、道路與神經元的分割在拓撲上正確。

沒有改變的是：要清理、量測、校正二值或灰階影像，本章的經典運算子仍然是最快、最可預測的工具，而且它們附帶學習式濾波器無法提供的保證（冪等性、連通性、精確保持形狀）。

重點整理

形態學處理的是形狀，不是頻率。它只用最小值、最大值、聯集與交集，並由一個結構元素控制；結構元素的大小與形狀就是你提出的問題。
侵蝕保留結構元素放得進去的位置；膨脹保留反射後的結構元素碰得到的位置。兩者互為對偶：侵蝕前景等於膨脹背景。
開運算（先侵蝕後膨脹）去除小的亮部；閉運算（先膨脹後侵蝕）補上小的暗縫。兩者都是冪等濾波器，先開後閉是很穩健的雜訊清除法。
擊中擊不中轉換比對精確的前景／背景圖樣，是細化、粗化與凸包的基礎。
重建（標記加遮罩，長到穩定為止）能完整保留整個物件：以重建做開運算、自動填洞與清除邊界物件都由它而來。
在灰階影像中，侵蝕與膨脹就是局部最小值與最大值。梯度、頂帽、底帽分別用於邊緣、背景校正與小物件擷取；粒度測定不需分割就能得到大小分布。
今天形態學仍廣泛用於神經網路遮罩的後處理、clDice 這類損失中的可微積木，也是可學習形態學層的研究主題。

練習

一張二值影像中有一個寬 40 像素、高 10 像素的矩形。分別描述它被 (a) $5 \times 5$ 方塊、(b) 長 15 像素的水平線段、(c) 長 15 像素的垂直線段侵蝕與開運算的結果。哪個結構元素會把矩形完全移除？為什麼？

提示

侵蝕保留的是「整個結構元素放得進去」時原點的位置。 $5 \times 5$ 方塊可放的位置構成 $36 \times 6$ 的區塊，而開運算會恢復整個矩形，因為矩形本身就是放得進去的方塊之聯集。水平線段也放得進去，所以開運算結果同樣是整個矩形。垂直線段比矩形還高，永遠放不進去：侵蝕與開運算都是空集合。

證明集合的膨脹滿足交換律與結合律： $A \oplus B = B \oplus A$ 、 $(A \oplus B) \oplus C = A \oplus (B \oplus C)$ 。用第二個等式解釋為什麼用 $3 \times 3$ 方塊膨脹兩次，等於用 $5 \times 5$ 方塊膨脹一次。侵蝕也有同樣的性質嗎？

提示

把膨脹寫成 Minkowski 和 $\{ a + b \mid a \in A, b \in B \}$ ，兩個等式都直接來自向量加法的交換律與結合律。 $3 \times 3$ 方塊加上自己就是 $5 \times 5$ 方塊。侵蝕則用 $(A \ominus B) \ominus C = A \ominus (B \oplus C)$ ，這可由對偶性推得。

設計一組擊中擊不中模板 $(B_1, B_2)$ ，偵測一條 8 連通、一個像素寬、往右延伸到盡頭（線從左邊過來）的線段端點。要抓到所有端點，總共需要幾個旋轉後的模板？用 scipy.ndimage.binary_hit_or_miss 在 skimage.morphology.skeletonize(~skimage.data.horse()) 上驗證你的答案。

提示

中心與左鄰居放在 $B_1$ ，其他七個鄰居放在 $B_2$ 。線也可能從斜方向過來，所以需要 8 個模板（4 個軸向加 4 個斜向）。把你數到的總數和修剪程式碼中 endpoints() 這類「數鄰居」函式的結果比較。

你的分割網路輸出細胞遮罩。許多遮罩有小洞、少數單像素的偽陽性，還有一些細胞被影像邊界切到。用 scikit-image 寫一個五行的後處理流程，並說明步驟順序的理由。

提示

一個合理的順序是：先 remove_small_objects（或小的開運算），避免斑點被填起來或被合併；接著 binary_fill_holes；輪廓參差時再做一次小的閉運算；然後 clear_border；最後用 label 計數。討論一下如果先做閉運算再去斑點會出什麼問題（相鄰的斑點可能合併成一個「細胞」）。

對灰階影像 $f$ 與平坦結構元素 $b$ ，證明頂帽 $T(f) = f - (f \circ b)$ 非負，而且是冪等的： $T(T(f)) = T(f)$ 。

提示

非負性來自 $f \circ b \le f$ 。冪等性則要證明 $T(f) \circ b = 0$ 。任取 $b$ 的一個平移 $W$ ，令 $y$ 為 $W$ 中 $f$ 最小的點。開運算在 $W$ 上處處至少是 $\min_W f$ ，又不超過 $f$ ，所以在 $y$ 處恰等於 $f(y)$ ，因此 $T(f)(y) = 0$ 。於是 $b$ 的每個平移都含有 $T(f)$ 的零點， $T(f)$ 被 $b$ 侵蝕後為零，開運算也是零。

參考文獻

R. C. Gonzalez and R. E. Woods, Digital Image Processing, 4th ed., Pearson, 2018, Ch. 9. publisher page
G. Matheron, Random Sets and Integral Geometry, Wiley, 1975. library record
J. Serra, Image Analysis and Mathematical Morphology, Academic Press, 1982. library record
R. M. Haralick, S. R. Sternberg, and X. Zhuang, “Image Analysis Using Mathematical Morphology,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. PAMI-9, pp. 532–550, 1987. doi
P. Soille, Morphological Image Analysis, Springer, 2004. doi
T. Y. Zhang and C. Y. Suen, “A Fast Parallel Algorithm for Thinning Digital Patterns,” Communications of the ACM, vol. 27, pp. 236–239, 1984. doi
L. Vincent, “Morphological Grayscale Reconstruction in Image Analysis: Applications and Efficient Algorithms,” IEEE Trans. Image Processing, vol. 2, pp. 176–201, 1993. doi
P. Maragos, “Pattern Spectrum and Multiscale Shape Representation,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 11, no. 7, pp. 701–716, 1989. doi
E. J. Breen and R. Jones, “Attribute Openings, Thinnings, and Granulometries,” Computer Vision and Image Understanding, vol. 64, no. 3, pp. 377–389, 1996. doi
P. Salembier and M. H. F. Wilkinson, “Connected Operators: A Review of Region-Based Morphological Image Processing Techniques,” IEEE Signal Processing Magazine, vol. 26, no. 6, pp. 136–157, 2009. doi
E. Carlinet and T. Géraud, “A Comparative Review of Component Tree Computation Algorithms,” IEEE Trans. Image Processing, vol. 23, pp. 3885–3895, 2014. doi
S. van der Walt et al., “scikit-image: Image Processing in Python,” PeerJ, 2014. arXiv
R. Mondal, S. Santra, S. S. Mukherjee, and B. Chanda, “Morphological Network: How Far Can We Go with Morphological Neurons?,” BMVC, 2022. arXiv
K. Nogueira, J. Chanussot, M. Dalla Mura, and J. A. dos Santos, “An Introduction to Deep Morphological Networks,” arXiv:1906.01751, 2019. arXiv
R. Hermary, G. Tochon, É. Puybareau, A. Kirszenberg, and J. Angulo, “Learning Grayscale Mathematical Morphology with Smooth Morphological Layers,” Journal of Mathematical Imaging and Vision, vol. 64, pp. 736–753, 2022. doi
A. Kirszenberg, G. Tochon, É. Puybareau, and J. Angulo, “Going Beyond p-convolutions to Learn Grayscale Morphological Operators,” arXiv:2102.10038, 2021. arXiv
S. Shit et al., “clDice — A Novel Topology-Preserving Loss Function for Tubular Structure Segmentation,” CVPR, 2021. arXiv