跳轉到

8.3 分水嶺與區域成長分割#

概念說明#

二值化與連通標記處理的是「已經分好黑白」的影像，但真實照片需要直接從灰階/彩色分割出區域。本節介紹兩個經典的自動分割演算法。

分水嶺 (watershed) 把灰階影像想像成地形：亮度 = 海拔。從每個局部極小值（盆地）注水，水位上升時不同盆地的水會在某條稜線相遇——這些分水嶺線就是區域邊界。實作上有兩種等價路線：沉浸法 (immersion, Vincent–Soille 1991) 按海拔層逐級淹水；優先佇列法用 (海拔, 先進先出) 的 priority queue 從種子向外擴張。結果是把影像切成互不重疊的集水區，邊界恰好落在梯度最大的稜線上。

過度分割 (over-segmentation) 是分水嶺的著名問題：影像裡有幾十個微小極小值，就會淹出幾十個區域。標準解法：標記控制 (marker-controlled) 分水嶺——不用所有極小值當水源，而是人工/自動指定少量「確定是物件」的種子（前景標記）與「確定是背景」的標記，水只從這些標記淹起。種子品質決定分割品質，這使分水嶺從「自動」變成「半自動」，實務上更可靠。得到種子的常用方法：對距離轉換做局部極大值（物件的中心）、閾值後的侵蝕核心等。

區域成長 (region growing) 走另一條路：從種子點出發，反覆檢查種子區域的鄰居，若滿足同質性準則（如與區域平均的灰階差 < 門檻）就併入，直到無人可加。它是貪婪演算法：簡單、直覺，但結果依賴 (1) 種子位置（不同種子長出不同區域）、(2) 同質性門檻（太鬆會溢出邊界、太緊會碎片化）、(3) 成長順序。改進方向：用區域平均而非種子值做比較、結合梯度資訊禁止跨越強邊緣成長。

兩者的關係：分水嶺其實是「同時從所有種子成長」的區域成長——差別在成長順序由海拔（梯度）決定，且邊界處精確停住。把區域成長的「一個種子」推廣到「多個種子按高度競爭」，就得到分水嶺。

演算法步驟#

標記控制分水嶺（優先佇列法）： 1. 計算梯度大小圖 G（第 4 章的方法）——地形是 G，盆地 = 平坦區內部。 2. 準備標記圖 M（標記處為 1,2,…，其餘 0）。 3. 所有標記像素推入優先佇列（key = G 值）。 4. 彈出 G 最小的像素，檢查其鄰居：未標記者賦予相同標記並入隊。 5. 直到佇列空。每個像素歸於「最先淹到它的盆地」。

區域成長： 1. 選種子點（如使用者點擊或亮度極值）。 2. 種子入堆疊；記錄區域平均 μ。 3. 彈出像素，檢查 4/8 鄰居：|I(鄰) − μ| < T 則併入、更新 μ、入堆疊。 4. 直到無新像素。

Python 程式碼範例#

import numpy as np
import heapq

# 地形：兩個亮圓（低梯度=盆地）被高梯度的邊界隔開
yy, xx = np.mgrid[0:40, 0:40]
img = np.full((40, 40), 100.0)
img[(xx-12)**2+(yy-12)**2 < 64] = 180
img[(xx-29)**2+(yy-29)**2 < 49] = 220
gy, gx = np.gradient(img)
G = np.hypot(gx, gy)
print("地形梯度範圍: %.1f ~ %.1f" % (G.min(), G.max()))

def watershed(G, markers):
    H, W = G.shape
    lab = markers.copy()
    pq = [(G[y, x], x, y, markers[y, x])
          for y, x in zip(*np.where(markers > 0))]
    heapq.heapify(pq)
    while pq:
        g, x, y, m = heapq.heappop(pq)
        for dy, dx in ((-1,0),(1,0),(0,-1),(0,1)):
            ny, nx = y+dy, x+dx
            if 0 <= ny < H and 0 <= nx < W and lab[ny, nx] == 0:
                lab[ny, nx] = m
                heapq.heappush(pq, (G[ny, nx], nx, ny, m))
    return lab

markers = np.zeros((40, 40), np.int32)
markers[12, 12] = 1; markers[29, 29] = 2      # 兩個種子
lab = watershed(G, markers)
print("盆地 1 面積:", (lab == 1).sum(), " 盆地 2 面積:", (lab == 2).sum())
print("邊界處（分水嶺）標記像素數:",
      int(((lab[20, :] == 1) | (lab[20, :] == 2)).sum()), "/ 40")

# 區域成長
def region_grow(I, seed, T):
    lab = np.zeros(I.shape, np.int32)
    stack = [seed]; lab[seed] = 1
    members = [I[seed]]; 
    while stack:
        y, x = stack.pop()
        mu = np.mean(members)
        for dy, dx in ((-1,0),(1,0),(0,-1),(0,1)):
            ny, nx = y+dy, x+dx
            if 0 <= ny < I.shape[0] and 0 <= nx < I.shape[1] \
               and lab[ny, nx] == 0 and abs(I[ny, nx] - mu) < T:
                lab[ny, nx] = 1; stack.append((ny, nx))
                members.append(I[ny, nx])
    return lab

T_test = region_grow(img, (12, 12), 30)
print("區域成長（T=30）從圓1種子成長的面積:", T_test.sum(), "（理論約 200）")
T_bad = region_grow(img, (12, 12), 200)
print("區域成長（T=200）面積:", T_bad.sum(), "（門檻太鬆，溢出到整張圖）")

執行後：分水嶺把兩個圓精確分開，邊界落在梯度高的環上；區域成長在 T=30 時收在第一個圓內（面積接近圓的 π·8² ≈ 201），T=200 時則吞併整張影像——生動展示門檻敏感性，也說明為何實務常把「與種子值的差」改成「與區域平均的差」並配合梯度牆。

複雜度分析#

分水嶺（優先佇列）：每像素至多入隊一次，堆操作 O(log N)，總計 O(N log N)。沉浸法（按灰階桶排序）可達 O(N + K)，K 為灰階數。空間 O(N)。區域成長：每像素至多被併入一次、檢查常數個鄰居，O(N)，但注意「每次併入都重算區域平均」若用樸素方式是 O(N²)——需維護 running mean（總和與計數的增量更新）才能線性。實務提醒：分水嶺對極小值數量極敏感（先做 h-minima 變換抑制淺盆地）；區域成長對種子與門檻敏感（多種子 + 條件式成長更穩）。兩者都屬於「區域型分割」，相對於基於邊緣（第 4 章）與基於聚類（k-means、第 10 章的語義分割前身）的另兩大路線。

練習#

  1. 不給種子、把「每個 G 的局部極小值」都當標記跑分水嶺，數一下區域數，驗證過度分割；再加入 h-minima（忽略深度 < h 的盆地）觀察區域數下降。
  2. 對 8.2 節的連通標記結果做「距離轉換 + 局部極大值」自動生成種子，跑標記控制分水嶺，分開兩個互相接觸的圓形物件。
  3. 把區域成長的同質性準則從「與區域平均比」改成「與種子值比」，對漸層影像測試，解釋為何後者會漏掉漸變區。