做圖像處理的人遲早會(huì)撞上“細(xì)化骨架”這道坎手里一張二值圖線條粗細(xì)不均、有毛刺、有斷口想量寬度得按列一行行掃想找交點(diǎn)得靠模板硬懟稍微轉(zhuǎn)個(gè)角度規(guī)則全廢。細(xì)化骨架thinning / skeletonization就是把這些“胖線條”壓成單像素寬的中心線讓寬度歸一、拓?fù)浔Ａ?、后續(xù)的測量和識(shí)別有統(tǒng)一的坐標(biāo)基準(zhǔn)。它不改變圖形的大致形狀只剝掉多余的外層像素剩下的那條“脊梁”承載了長度、走向、分支、端點(diǎn)這些真正有用的信息。這篇內(nèi)容我會(huì)從算法內(nèi)核講到工程落地包含可直接跑的代碼、參數(shù)怎么定、坑在哪適合已經(jīng)會(huì)寫幾行 OpenCV、但被細(xì)化結(jié)果里滿屏毛刺和斷裂折磨過的同學(xué)也適合剛接觸形態(tài)學(xué)、想弄清楚“為什么刪這個(gè)像素不刪那個(gè)”的新手。1. 先把“細(xì)化骨架”這件事想明白1.1 骨架到底在替我們解決什么問題原始二值圖里的線條是有寬度的寬度本身攜帶噪聲打印抖一點(diǎn)、掃描糊一點(diǎn)、閾值偏一點(diǎn)同一條筆畫這頭 7 像素那頭 11 像素。你拿這樣的圖去做長度統(tǒng)計(jì)誤差直接跟著寬度走去做模板匹配字體一放大就崩去做直線擬合邊緣兩側(cè)的像素互相拉扯。骨架化的本質(zhì)是降維——把二維的面狀目標(biāo)壓成一維的線狀目標(biāo)寬度信息被丟棄只留中心位置。這里有個(gè)容易被忽略的點(diǎn)骨架不是“隨便找條中線”。數(shù)學(xué)上嚴(yán)格定義的骨架來自中軸變換某個(gè)像素屬于骨架當(dāng)且僅當(dāng)以它為圓心的最大內(nèi)切圓至少與目標(biāo)邊界有兩個(gè)接觸點(diǎn)。這個(gè)定義決定了骨架的三個(gè)天然屬性單像素寬、與原目標(biāo)同倫該連的還連著該有的洞還在、位置居中。工程上我們用的迭代模板法只是這個(gè)定義的近似實(shí)現(xiàn)所以才會(huì)出現(xiàn)毛刺、斜線鋸齒、末端多出一小截這類“不完美”后面所有后處理都是在給這個(gè)近似擦屁股。我見過不少項(xiàng)目在這上面翻車原因往往是期望值錯(cuò)了以為細(xì)化完就是一條干凈漂亮的線。實(shí)際拿到的第一版結(jié)果十個(gè)里有八個(gè)是“主干挺好邊上全是小刺”。接受這個(gè)前提后面的工作才好安排。1.2 三條主流路線選錯(cuò)了后面全是補(bǔ)丁把線條壓成單像素業(yè)界有三條路各自的脾氣完全不同選之前先想清楚下游要什么。路線代表實(shí)現(xiàn)核心思路優(yōu)點(diǎn)短板適用場景迭代模板細(xì)化Zhang-Suen、Guo-Hall、Hilditch反復(fù)掃描按 8 鄰域模板逐層剝皮結(jié)果拓?fù)浞€(wěn)定、實(shí)現(xiàn)簡單、可控迭代次數(shù)隨線寬增長、斜線有鋸齒、毛刺多字符筆畫、走線、裂紋、通用二值圖形態(tài)學(xué)腐蝕骨架開運(yùn)算/腐蝕組合用結(jié)構(gòu)元反復(fù)腐蝕再取殘差速度快、純形態(tài)學(xué)算子斷裂嚴(yán)重、位置偏移、幾乎不可用于測量快速預(yù)覽、粗定位距離變換脊線medial_axis、中軸算每個(gè)前景點(diǎn)到背景的距離取局部極大保留半徑信息、管徑可直接讀有孔洞時(shí)生成環(huán)、結(jié)果對噪聲極敏感血管管徑、顆粒尺寸、需要半徑的場合我的默認(rèn)選擇是第一條。原因很實(shí)在絕大多數(shù)下游任務(wù)要的是“一條能追蹤的路徑”而不是“精確的內(nèi)切圓半徑”。迭代模板法在拓?fù)浔３稚献龅米詈脝蜗袼貙挾纫灿斜ＷC唯一的代價(jià)是慢而慢可以用 ROI 裁剪和索引表優(yōu)化掉。形態(tài)學(xué)腐蝕那條路我只在一種情況下用需要極快的粗篩比如先看看圖里有沒有長條狀目標(biāo)有的話再走精細(xì)流程。它的斷裂問題不是調(diào)參能救的腐蝕到只剩一線的時(shí)候稍微細(xì)一點(diǎn)的位置早就斷了。距離變換那條路適合“測量”類任務(wù)。舉個(gè)具體例子血管分割后要報(bào)管徑medial_axis 返回的骨架每個(gè)點(diǎn)帶著到邊界的距離直徑就是兩倍距離值一步到位。但你得有心理準(zhǔn)備只要目標(biāo)內(nèi)部有個(gè)小孔骨架就會(huì)繞著孔長出一個(gè)閉環(huán)后面又得寫代碼去環(huán)。1.3 我給“好骨架”定的四條驗(yàn)收標(biāo)準(zhǔn)調(diào)參之前先立標(biāo)準(zhǔn)不然永遠(yuǎn)在“看著差不多”和“好像又差一點(diǎn)”之間搖擺。我一般用這四條卡單像素性骨架里不應(yīng)該存在 2×2 的實(shí)心塊任意骨架像素的 8 鄰域內(nèi)前景點(diǎn)數(shù)量應(yīng)盡量不超過 2分支點(diǎn)除外。連通性一致細(xì)化前后的 8 連通分量數(shù)量應(yīng)當(dāng)相同。細(xì)化把一個(gè)連通域拆成兩半說明算法在這一段判錯(cuò)了。分支合理分支點(diǎn)的數(shù)量和位置要能和原圖對得上。原圖是個(gè)十字骨架就該有個(gè)四叉原圖是一條直線骨架就不該冒出分叉。端點(diǎn)收斂端點(diǎn)數(shù)量不能爆炸。一條直線細(xì)化出七八個(gè)端點(diǎn)說明毛刺沒清干凈。這四條里第二條最硬。前三條可以靠肉眼第二條必須用代碼數(shù)因?yàn)槿搜蹖Α皵嗔?3 個(gè)像素”這種事實(shí)在不敏感。2. 算法內(nèi)核模板怎么判、連通性怎么保2.1 八鄰域編碼是所有細(xì)化算法的公共語言所有迭代細(xì)化算法都建立在同一套記號(hào)上。對當(dāng)前像素 P1按順時(shí)針把周圍八個(gè)鄰居編號(hào) P2 到 P9P9 P2 P3 P8 P1 P4 P7 P6 P5基于這套編號(hào)有兩個(gè)量反復(fù)出現(xiàn)。一個(gè)是B(P1)八個(gè)鄰居里前景像素的個(gè)數(shù)它衡量這個(gè)點(diǎn)周圍有多“厚”。另一個(gè)是A(P1)沿 P2→P3→…→P9→P2 這個(gè)環(huán)形序列走一圈出現(xiàn)“0 變 1”的次數(shù)它衡量這個(gè)點(diǎn)是不是連接樞紐。A 的物理意義值得多說兩句。A 等于 1說明繞這個(gè)點(diǎn)轉(zhuǎn)一圈前景只起了一次頭這個(gè)點(diǎn)是個(gè)“過路點(diǎn)”刪掉它不會(huì)把兩側(cè)斷開。A 大于 1說明前景起了好幾回頭這個(gè)點(diǎn)是多條線的交匯處或者拐角刪了就會(huì)撕裂。所以幾乎所有算法的核心判據(jù)都是A 1這一條就是拓?fù)浔３值谋ｋU(xiǎn)絲。至于位權(quán)編碼是把八個(gè)鄰居按1,2,4,8,16,32,64,128打包成一個(gè) 0 到 255 的整數(shù)。這個(gè)技巧在索引表法里是主角后面細(xì)說。2.2 Zhang-Suen 兩遍掃描為什么非要分兩遍Zhang-Suen 是最常被拿來入門的版本它的規(guī)則可以壓縮成一句話在滿足 A1 且 2≤B≤6 的前提下分兩個(gè)子迭代各砍掉一部分邊界點(diǎn)。第一子迭代的附加條件是 P2·P4·P60 且 P4·P6·P80第二子迭代換成 P2·P4·P80 且 P2·P6·P80。這兩組條件看起來像天書其實(shí)是在限制“只能砍東南方向”和“只能砍西北方向”。為什么非要分兩遍因?yàn)閱伪閽呙钑?huì)出問題。假設(shè)有一條寬 2 像素的斜線如果一次性把所有滿足 A1 的邊界點(diǎn)全刪掉整條線可能被削成一條對角線上的孤立點(diǎn)甚至直接消失。分兩遍相當(dāng)于給每次刪除加了方向約束第一遍只削一側(cè)第二遍再削另一側(cè)剩余的部分還能撐住結(jié)構(gòu)。這個(gè)設(shè)計(jì)帶來一個(gè)副作用結(jié)果是順序相關(guān)的。同一張圖正著掃和反著掃結(jié)果會(huì)有細(xì)微差別通常表現(xiàn)為末端的形狀不同。做像素級(jí)對比測試的時(shí)候別拿這個(gè)當(dāng) bug。B 的下界設(shè)成 2 也有講究。B1 說明這個(gè)點(diǎn)本身就是端點(diǎn)刪了就少一截B0 是孤立點(diǎn)。上界 6 是為了保護(hù)粗線內(nèi)部太厚的點(diǎn)周圍鄰居多還沒輪到它。2.3 索引表法把逐點(diǎn)判斷壓成一次查表純 Python 寫 Zhang-Suen一張 2000×2000 的圖能跑到分鐘級(jí)原因是每個(gè)像素每輪都要做十幾次條件判斷而迭代輪數(shù)又跟線寬成正比。優(yōu)化思路很直接八個(gè)鄰居的組合只有 256 種那就預(yù)先算好哪 256 種組合可以刪運(yùn)行時(shí)只做一次位運(yùn)算拼索引、一次數(shù)組查表。import numpy as np def build_zs_lut(): lut np.zeros(256, dtypenp.uint8) for code in range(256): # 按 P2..P9 順序取位位權(quán)遞增 p [(code i) 1 for i in range(8)] B sum(p) if B 2 or B 6: continue seq p [p[0]] A sum(1 for i in range(8) if seq[i] 0 and seq[i 1] 1) if A ! 1: continue p2, p3, p4, p5, p6, p7, p8, p9 p # 第一子迭代條件 if p2 * p4 * p6 0 and p4 * p6 * p8 0: lut[code] | 1 # 第二子迭代條件 if p2 * p4 * p8 0 and p2 * p6 * p8 0: lut[code] | 2 return lut建好表之后每輪循環(huán)只做卷積式的鄰域取值和查表用 numpy 切片就能向量化速度能提一到兩個(gè)數(shù)量級(jí)。我實(shí)測過一張 1600×1200 的工程圖紙純循環(huán)版本跑 40 秒往上索引表加向量化壓到 1 秒出頭差距主要在這一步。注意拼索引的時(shí)候位序必須和建表時(shí)完全一致。P2 到 P9 順時(shí)針位權(quán) 1 到 128寫錯(cuò)一位結(jié)果就是“看起來還在細(xì)化實(shí)際刪錯(cuò)一堆像素”而且很難從輸出圖上直接看出來只能靠連通性測試抓。2.4 Guo-Hall 與 Hilditch什么時(shí)候值得換算法Zhang-Suen 的毛病是斜線鋸齒明顯。一條 45 度的寬線細(xì)化完骨架上會(huì)出現(xiàn)規(guī)律的臺(tái)階。Guo-Hall 用四個(gè)條件替換判據(jù)更細(xì)斜向的支持更好結(jié)果通常更平滑而且子迭代次數(shù)固定為兩次收斂更快。Hilditch 是更老的經(jīng)典方法它引入“并行”和“串行”的混合判定考慮鄰居在本次迭代中是否已被標(biāo)記刪除。這個(gè)特性讓它在某些手工筆畫上表現(xiàn)更好但實(shí)現(xiàn)復(fù)雜度高而且結(jié)果對掃描順序敏感我在工程里用得不多。實(shí)際選擇上我的經(jīng)驗(yàn)是這樣字符識(shí)別、筆畫分析這類對骨架平滑度敏感的優(yōu)先試 Guo-Hall通用二值目標(biāo)做長度和分支統(tǒng)計(jì)Zhang-Suen 夠用圖像里既有粗線又有大量細(xì)弱連接先別急著換算法去檢查預(yù)處理八成是二值化把細(xì)節(jié)吃掉了。換算法是最后的手段因?yàn)閾Q完之后后處理的參數(shù)全要重調(diào)。3. 實(shí)操從一張臟圖到一條干凈骨架3.1 預(yù)處理決定八成成敗細(xì)化是放大器輸入有多臟輸出就有多臟。噪聲點(diǎn)被細(xì)化成短線段邊緣毛刺被細(xì)化成分叉小孔洞被細(xì)化成閉環(huán)。所以順序上所有清理工作都必須放在細(xì)化之前細(xì)化之后再補(bǔ)救成本翻好幾倍。我固定的預(yù)處理鏈?zhǔn)沁@樣灰度化 → 中值濾波 → 背景校正可選→ 二值化 → 形態(tài)學(xué)開運(yùn)算 → 面積過濾 → 細(xì)化。中值濾波核取 3 或 5視噪聲顆粒大小。這里要提醒一句別用高斯高斯會(huì)把邊緣模糊二值化之后再細(xì)化末端會(huì)縮成圓頭長度統(tǒng)計(jì)偏短。中值對椒鹽噪聲的效果好而且保邊。開運(yùn)算的核大小按你預(yù)估的最小線寬來定一般是線寬的三分之一左右。線寬 9 像素核取 3×3 就夠核一大細(xì)線直接消失。開運(yùn)算的目的是去掉孤立的凸起和細(xì)小的粘連讓后續(xù)細(xì)化少長幾根刺。面積過濾這一步很多人省掉我覺得不能省。二值化之后邊緣常有一堆十幾像素的小碎塊細(xì)化后每條碎塊都變成一根小短棍全混進(jìn)結(jié)果里分支點(diǎn)統(tǒng)計(jì)直接廢掉。用cv2.connectedComponentsWithStats按面積篩一遍閾值取預(yù)估最小目標(biāo)面積的十分之一效果立竿見影。3.2 二值化參數(shù)怎么定別迷信一種方法方法關(guān)鍵參數(shù)適用條件坑固定閾值threshold 值光照穩(wěn)定、背景純凈換一張圖就廢別用在批量流程Otsu無直方圖雙峰明顯前景占比極小5%時(shí)閾值偏移嚴(yán)重自適應(yīng)閾值blockSize奇數(shù)、C光照不均、背景漸變blockSize 過小會(huì)引入塊狀噪聲Sauvola 類窗口、k文檔、低對比度計(jì)算量大參數(shù)需按圖調(diào)我最常用的是“自適應(yīng) Otsu 兜底”的組合先算 Otsu看前景占比是否落在 3% 到 40% 這個(gè)區(qū)間落在里面就用 Otsu落外面說明這張圖光照或?qū)Ρ榷扔袉栴}切自適應(yīng)。自適應(yīng)閾值的 blockSize 取圖中典型線寬的 3 到 5 倍并強(qiáng)制為奇數(shù)C 取 5 到 10實(shí)測在多數(shù)工業(yè)圖上都穩(wěn)。這里有個(gè)細(xì)節(jié)值得說二值化后要做一次連通性方向的一致性檢查。如果目標(biāo)的細(xì)線在二值化后出現(xiàn)大量單像素?cái)嗫诩?xì)化后就是一堆斷頭骨架。補(bǔ)救辦法是在細(xì)化前做一次 3×3 的閉運(yùn)算把斷口先接上。寧可接錯(cuò)幾個(gè)像素也別留斷口因?yàn)閿嗫谠诠羌茈A段幾乎無法區(qū)分“真斷”和“假斷”。3.3 細(xì)化實(shí)現(xiàn)三行調(diào)用與手寫版本如果只是要結(jié)果OpenCV 的 contrib 模塊和 scikit-image 都能直接用import cv2 import numpy as np from skimage.morphology import skeletonize, medial_axis # 前提binary 為 uint8前景 255背景 0 binary cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU)[1] # 方案一OpenCV contrib速度快返回 0/255 skel_cv cv2.ximgproc.thinning(binary, thinningTypecv2.ximgproc.THINNING_ZHANGSUEN) # 方案二scikit-image輸入需為 bool返回 bool skel_sk skeletonize(binary 0, methodzhang) # 方案三需要半徑信息時(shí)用 medial_axis skel_ma, dist medial_axis(binary 0, return_distanceTrue)三個(gè)方案的差異值得記一下ximgproc.thinning支持 Guo-Hall速度最快skeletonize的接口干凈返回的 bool 數(shù)組做后續(xù)邏輯運(yùn)算方便medial_axis多返回一個(gè)距離圖做管徑、顆粒尺寸測量時(shí)直接用它省掉一步。手寫版本用前面建的索引表向量化實(shí)現(xiàn)大概是這樣def zhang_suen_fast(binary, max_iter200): img (binary 0).astype(np.uint8) img np.pad(img, 1, modeconstant) lut build_zs_lut() for _ in range(max_iter): changed False for phase in (0, 1): # 兩個(gè)子迭代 p2 img[:-2, 1:-1]; p3 img[:-2, 2:]; p4 img[1:-1, 2:] p5 img[2:, 2:]; p6 img[2:, 1:-1]; p7 img[2:, :-2] p8 img[1:-1, :-2]; p9 img[:-2, :-2] p1 img[1:-1, 1:-1] code (p2.astype(np.uint16) | (p3.astype(np.uint16) 1) | (p4 2) | (p5 3) | (p6 4) | (p7 5) | (p8 6) | (p9 7)) removable (lut[code] (1 phase)) 0 hit removable (p1 1) if hit.any(): p1[hit] 0 changed True if not changed: break return (img[1:-1, 1:-1] * 255).astype(np.uint8)max_iter設(shè) 200 是保險(xiǎn)值。理論上輪數(shù)約等于最粗處的半寬實(shí)際因?yàn)閮蓚€(gè)子迭代一輪各刪一層粗 40 像素的線大概 20 輪內(nèi)就收斂。留足余量是為了防止某張圖有異常粗的塊導(dǎo)致提前退出。提示ximgproc需要裝opencv-contrib-python只裝基礎(chǔ)包會(huì)報(bào)module cv2 has no attribute ximgproc。這是新手最常見的第一個(gè)攔路虎。3.4 后處理剪枝、去毛刺、補(bǔ)斷口細(xì)化出來的第一版通常長這樣主干清楚但每隔幾個(gè)像素就伸出一根兩三個(gè)像素的小刺。這些小刺的來源是原圖邊緣的微小起伏本質(zhì)上是中軸變換對邊界噪聲的正常反應(yīng)不是算法錯(cuò)了。剪枝的標(biāo)準(zhǔn)做法是迭代刪除端點(diǎn)找 8 鄰域只有 1 個(gè)前景點(diǎn)的像素刪掉重復(fù) N 次。N 取多大會(huì)直接決定結(jié)果——N 太小刺還在N 太大會(huì)把真實(shí)的短枝也吃掉。def prune(skel, max_len5): skel (skel 0).astype(np.uint8) for _ in range(max_len): kernel np.array([[1, 1, 1], [1, 10, 1], [1, 1, 1]], dtypenp.uint8) # 鄰域和等于 11 表示中心為前景且只有一個(gè)前景鄰居端點(diǎn) nb cv2.filter2D(skel, -1, kernel) endp ((skel 1) (nb 11)).astype(np.uint8) if endp.sum() 0: break skel[endp 1] 0 return skel那個(gè) 10 的技巧值得解釋一下中心權(quán)重設(shè) 10八個(gè)鄰居各 1中心為 1 時(shí)總和至少 11。鄰域和正好 11說明中心是前景且周圍只有一個(gè)前景點(diǎn)那就是端點(diǎn)。這個(gè)方法比手動(dòng)數(shù)鄰居快得多。N 的取值我一般按線寬來線寬的 0.5 到 1 倍。線寬 10 像素N 取 5 到 8。但別全局一刀切如果圖里既有粗線又有細(xì)線最好按連通域分別處理細(xì)線的刺短用同一個(gè) N 會(huì)把細(xì)線的真實(shí)末端削掉。補(bǔ)斷口是另一件事。骨架斷裂的根源在二值化但細(xì)化之后也能救找所有端點(diǎn)兩兩配對距離小于閾值且方向大致相反的就連一條直線。閾值取預(yù)估最大斷口寬度的 2 倍方向判據(jù)用端點(diǎn)處的局部走向夾角超過 60 度就別連。def find_endpoints_and_branches(skel): s (skel 0).astype(np.uint8) kernel np.ones((3, 3), np.uint8) nb cv2.filter2D(s, -1, kernel) - s # 只數(shù)鄰居 endpoints np.argwhere((s 1) (nb 1)) branches np.argwhere((s 1) (nb 3)) return endpoints, branches端點(diǎn)數(shù)量是個(gè)很好的健康指標(biāo)。我一般會(huì)把它打到日志里同一批圖里某張圖的端點(diǎn)數(shù)突然翻三倍基本可以判定那張圖的二值化出了問題回去查比在骨架階段硬修劃算得多。3.5 骨架轉(zhuǎn)路徑從像素集合到可用的線骨架是像素集合下游要的往往是“一條從 A 到 B 的點(diǎn)序列”或者“線段列表”。轉(zhuǎn)換的關(guān)鍵是處理分支點(diǎn)遇到分支就切段。做法是從端點(diǎn)出發(fā)做深度優(yōu)先遍歷走到分支點(diǎn)就停把走過的路徑記下來然后從分支點(diǎn)繼續(xù)往未訪問的鄰居走。每段路徑存成點(diǎn)列表附帶起止點(diǎn)坐標(biāo)和長度。這套結(jié)構(gòu)拿去做直線擬合、做長度統(tǒng)計(jì)、做拓?fù)鋱D分析都方便。def trace_paths(skel): s (skel 0).astype(np.uint8) kernel np.ones((3, 3), np.uint8) nb cv2.filter2D(s, -1, kernel) - s ends [tuple(p) for p in np.argwhere((s 1) (nb 1))] visited np.zeros_like(s) paths [] for start in ends: if visited[start] 1: continue path [start] visited[start] 1 cur start while True: y, x cur nxt None for dy in (-1, 0, 1): for dx in (-1, 0, 1): if dy 0 and dx 0: continue ny, nx y dy, x dx if 0 ny s.shape[0] and 0 nx s.shape[1] \ and s[ny, nx] 1 and visited[ny, nx] 0: nxt (ny, nx) break if nxt: break if nxt is None: break path.append(nxt) visited[nxt] 1 cur nxt if nb[cur] 3: # 分支點(diǎn)切段 paths.append(path) path [cur] if len(path) 1: paths.append(path) return paths這段代碼沒有處理“分支點(diǎn)重復(fù)入?！钡倪吔缜闆r實(shí)際工程里建議加一個(gè)分支點(diǎn)訪問計(jì)數(shù)允許一個(gè)分支點(diǎn)被多條路徑共享。我第一版就是漏了這個(gè)結(jié)果交叉的十字只追出一條邊就停了排查了半天。4. 性能與工程化大圖、批量、實(shí)時(shí)4.1 復(fù)雜度從哪來加速從哪下手細(xì)化算法的復(fù)雜度是O(N × W)N 是像素總數(shù)W 是最大半寬對應(yīng)的迭代輪數(shù)。這個(gè)式子說明兩件事圖越大越慢線越粗越慢。加速就三個(gè)方向。降像素。如果目標(biāo)只占圖的一小塊先做 ROI 裁剪再細(xì)化。工程圖紙這種大面積留白的場景裁完能省 80% 的時(shí)間。裁剪前記得把 ROI 邊界外擴(kuò)幾個(gè)像素避免目標(biāo)被切斷。降輪數(shù)。線寬如果普遍偏粗先做一次腐蝕把線削細(xì)再細(xì)化能顯著減少迭代輪數(shù)。但要小心腐蝕會(huì)改拓?fù)浼?xì)的地方直接斷。只有當(dāng)線寬普遍大于 6 像素、且細(xì)線不是關(guān)注對象時(shí)才敢這么干。向量化與索引表。前面講過把逐像素判斷換成 numpy 鄰域運(yùn)算加 256 項(xiàng)查表速度差一個(gè)數(shù)量級(jí)。這是收益最高的一步做完之后一般不需要再考慮 C 重寫。我做過一個(gè)粗略對比同一張 2400×1800 的走線圖純 Python 雙循環(huán)約 55 秒numpy 向量化約 3.5 秒向量化加索引表約 1.2 秒OpenCV 的ximgproc.thinning約 0.4 秒?？梢娋幾g好的實(shí)現(xiàn)還是最快自研版本的價(jià)值在于能改條件、能加約束比如只允許特定方向的刪除。4.2 批處理與內(nèi)存別一次全讀進(jìn)來批量跑圖的時(shí)候最常見的翻車是內(nèi)存。幾百張圖一次性讀成數(shù)組加上中間變量幾百兆就上去了。我的做法是單張?zhí)幚硗炅⒖提尫沤Y(jié)果只保留骨架的坐標(biāo)列表和統(tǒng)計(jì)量原始數(shù)組不保留。中間變量的類型也要注意。astype(np.uint8)和astype(np.uint16)在拼接鄰域索引時(shí)必須有區(qū)別索引會(huì)超過 255用 uint8 會(huì)靜默溢出表現(xiàn)是“某些像素莫名被刪”。這個(gè) bug 我踩過排查了兩天才定位到。還有一點(diǎn)OpenCV 的filter2D對 bool 數(shù)組不支持得先轉(zhuǎn) uint8。scikit-image 的輸入輸出都是 bool跟 OpenCV 混用的時(shí)候類型轉(zhuǎn)換別忘了一忘就報(bào)維度不匹配。4.3 對接下游不同任務(wù)對骨架的要求不一樣字符與筆畫分析。傳統(tǒng) OCR 流程里細(xì)化過骨架用來統(tǒng)計(jì)筆畫方向和寬度變換。現(xiàn)在深度學(xué)習(xí)方案大多直接用灰度圖細(xì)化反而會(huì)丟信息尤其在低分辨率下。但如果要做模板匹配、要算筆畫的細(xì)節(jié)點(diǎn)骨架依然是基礎(chǔ)。走線與裂紋檢測。骨架直接給出“長度”這個(gè)指標(biāo)不用再折算寬度。分支點(diǎn)就是分叉位置端點(diǎn)就是斷口位置。裂紋檢測里一條裂紋的骨架總長就是它的長度分支數(shù)就是分叉數(shù)比在原圖上做區(qū)域統(tǒng)計(jì)穩(wěn)得多。血管與管狀目標(biāo)。這類任務(wù)必須用距離變換那條路因?yàn)橐獔?bào)管徑。medial_axis的返回值里帶著每個(gè)骨架點(diǎn)的內(nèi)切半徑直徑取兩倍位置取骨架坐標(biāo)一次搞定。用迭代模板法也能做但需要額外算距離圖去查多一步。顆粒與纖維。纖維長度、彎曲度都是骨架上的量。這里要注意細(xì)化前必須做足去噪纖維邊緣的一點(diǎn)毛毛細(xì)化后就是一根分叉直接影響彎曲度的計(jì)算。5. 踩坑實(shí)錄與排查速查表5.1 典型故障對照表現(xiàn)象最可能的原因排查動(dòng)作處理方式骨架大面積斷裂二值化閾值偏高細(xì)線先沒了回看二值圖量細(xì)線處前景寬度降閾值或改自適應(yīng)細(xì)化前閉運(yùn)算主干上密集短刺邊緣噪聲、未做開運(yùn)算看原圖邊緣是否毛糙加大中值核加開運(yùn)算剪枝 N 調(diào)大交叉處變成糊團(tuán)交叉區(qū)域太粗迭代不收斂看交叉處原始寬度加大 max_iter或?qū)徊?ROI 單獨(dú)處理出現(xiàn)閉環(huán)原圖無孔線寬不一致細(xì)處斷開后繞回檢查是否存在環(huán)形細(xì)結(jié)構(gòu)別急著去環(huán)先看是不是真實(shí)孔洞整圖消失或只剩點(diǎn)輸出了 0/1 被當(dāng)成 0/255 顯示打印像素值范圍細(xì)化后統(tǒng)一乘 255端點(diǎn)數(shù)量異常多噪點(diǎn)碎塊被保留統(tǒng)計(jì)連通域面積分布細(xì)化前加面積過濾結(jié)果和上次不一樣算法本身順序相關(guān)或輸入有變化固定輸入做回歸測試記錄算法版本與參數(shù)別裸跑提速后結(jié)果變差位序拼錯(cuò)或類型溢出與參考實(shí)現(xiàn)對拍逐像素 diff檢查 uint8/uint16 與位權(quán)順序5.2 剪枝閾值到底怎么定這是被問得最多的問題。我的答案永遠(yuǎn)是先量線寬再定閾值。具體做法是在二值圖上隨機(jī)取若干行統(tǒng)計(jì)每段前景的連續(xù)長度取中位數(shù)當(dāng)線寬。剪枝長度取線寬的 0.5 到 1 倍這是經(jīng)驗(yàn)值原理是真實(shí)分支的長度通常大于線寬而毛刺的長度通常小于線寬。但如果圖里線寬差異大一刀切必然誤傷。這時(shí)候按連通域分別統(tǒng)計(jì)線寬各自剪各自的。代碼上就是connectedComponentsWithStats加循環(huán)不復(fù)雜但能把誤剪率降下來。還有一個(gè)更聰明的做法不做長度剪枝改成按分支的“深度”剪。從每個(gè)端點(diǎn)往回走如果走不到一個(gè)長度超過閾值的主干就停說明這是根懸空的短枝整根刪掉如果能走到主干說明是真分支保留。這個(gè)判據(jù)比單純數(shù)端點(diǎn)準(zhǔn)確得多代價(jià)是要多寫幾十行追蹤邏輯。5.3 幾條不那么容易查到的經(jīng)驗(yàn)第一條先在灰度域做背景校正別在二值域硬修。光照不均導(dǎo)致的粗細(xì)不均二值化怎么調(diào)都調(diào)不好用大核形態(tài)學(xué)做一次背景估計(jì)再相減效果比調(diào)十次閾值強(qiáng)。第二條別用骨架去算面積。骨架是線面積是面這倆沒有對應(yīng)關(guān)系。要面積就老老實(shí)實(shí)在二值圖上數(shù)前景像素。第三條細(xì)化前后都要存中間結(jié)果。出問題的時(shí)候能一眼看出是二值化錯(cuò)了還是細(xì)化錯(cuò)了省掉大量重跑時(shí)間。我現(xiàn)在的流程固定輸出四張圖原圖、二值圖、細(xì)化圖、剪枝后圖命名帶批次號(hào)。第四條回歸測試用固定的幾張“刁鉆圖”。我手上常備五張一條 45 度斜線、一個(gè)十字交叉、一條 2 像素細(xì)線、一個(gè)帶小孔的區(qū)域、一張光照不均的實(shí)拍。改了參數(shù)先跑這五張肉眼過了再上批量。這五張圖能覆蓋八成的問題。第五條滑動(dòng)窗口的處理別忘補(bǔ)邊。大圖切塊處理時(shí)塊與塊的接縫處骨架必?cái)唷＝鉀Q辦法是塊之間重疊線寬的 3 倍以上拼回來看重疊區(qū)重合的保留不重合的按距離就近連。這個(gè)細(xì)節(jié)不做整張圖會(huì)呈現(xiàn)規(guī)整的網(wǎng)格狀斷裂非常難看。6. 驗(yàn)收與量化評估別只靠眼睛6.1 單像素寬度怎么量最直接的方法是統(tǒng)計(jì)每個(gè)骨架像素的 8 鄰域前景數(shù)。理想情況下直線段上的點(diǎn)鄰域數(shù)為 2端點(diǎn)數(shù)為 1分支點(diǎn)數(shù)大于等于 3。定義單像素率為鄰域數(shù)小于等于 2 的骨架點(diǎn)占比這個(gè)值在 0.85 以上算合格。低于 0.8 說明存在大量 2×2 殘留塊通常是算法沒收斂或者max_iter設(shè)小了。還有一個(gè)指標(biāo)是 2×2 塊的數(shù)量。用四個(gè)方向的 2×2 窗口卷積一遍統(tǒng)計(jì)全 1 的窗口數(shù)理想值是 0。實(shí)際工程里能壓到個(gè)位數(shù)就很好了因?yàn)槟承┙徊纥c(diǎn)確實(shí)會(huì)殘留。6.2 拓?fù)浔３衷趺礈y連通分量數(shù)量用cv2.connectedComponents(binary, connectivity8)數(shù)細(xì)化前后對比。差值大于 0 說明斷了小于 0 說明連了不該連的。這兩個(gè)方向都要報(bào)警?？锥磾?shù)量用歐拉數(shù)或者簡單地對背景做連通分量統(tǒng)計(jì)。原圖有幾個(gè)洞骨架就應(yīng)該有幾個(gè)環(huán)。這條在裂紋、血管這類任務(wù)里特別重要因?yàn)榄h(huán)對應(yīng)的是真實(shí)的組織間隙或者血管分叉。分支點(diǎn)的數(shù)量不需要絕對一致但變化率應(yīng)該控制在 20% 以內(nèi)。超過這個(gè)數(shù)回去看剪枝是不是剪狠了或者預(yù)處理把細(xì)小結(jié)構(gòu)吃掉了。6.3 一張參數(shù)自查表參數(shù)作用典型取值調(diào)整方向中值核去噪保邊3 或 5噪聲大往大調(diào)細(xì)線多往小調(diào)開運(yùn)算核去凸起線寬/3 的奇數(shù)毛刺多調(diào)大細(xì)線消失調(diào)小面積閾值去碎塊最小目標(biāo)面積/10目標(biāo)多且小就調(diào)小自適應(yīng) blockSize局部閾值窗口線寬×3 到 5奇數(shù)光照漸變快調(diào)小塊狀噪聲調(diào)大max_iter細(xì)化輪數(shù)上限200異常粗塊多就調(diào)大剪枝長度 N去毛刺線寬×0.5 到 1刺殘留調(diào)大真分支被削調(diào)小斷口連接距離補(bǔ)斷最大斷口寬×2連錯(cuò)就調(diào)小分支夾角門限斷口方向校驗(yàn)60 度交叉結(jié)構(gòu)多就調(diào)大這張表我貼在自己的便簽上每次開新項(xiàng)目先按表填一遍初值再根據(jù)輸出微調(diào)。比起從零試參數(shù)能省不少時(shí)間。提醒一句參數(shù)之間是有耦合的。開運(yùn)算核調(diào)大線寬統(tǒng)計(jì)就變小剪枝閾值跟著要調(diào)小自適應(yīng)窗口調(diào)小二值圖噪聲變多面積閾值就得跟著調(diào)大。別孤立地看單個(gè)參數(shù)改一個(gè)至少回歸測三項(xiàng)。我個(gè)人在實(shí)際項(xiàng)目里的體會(huì)是細(xì)化骨架這活兒的難點(diǎn)從來不在算法本身那幾十行代碼網(wǎng)上一搜一堆。真正花時(shí)間的是預(yù)處理鏈的設(shè)計(jì)和后處理閾值的拿捏而這兩塊的答案只能從自己那批圖里長出來。所以別急著抄參數(shù)先把手上圖的線寬分布、噪聲水平、斷口比例這幾個(gè)數(shù)摸清楚剩下的就是按圖填空。我現(xiàn)在的習(xí)慣是先拿十張有代表性的圖跑一遍全流程把每個(gè)環(huán)節(jié)的中間結(jié)果和統(tǒng)計(jì)量都打出來形成一份“基線檔案”后面所有改動(dòng)都對著它比。這份檔案比任何參數(shù)表都值錢。