做顯微圖像分割的人大概都遇到過同一種窘境灰度直方圖上兩個(gè)峰的間隔還算明顯但閾值一提上去纖維束的走向就斷了一提下來又把相鄰的兩根纖維糊成一片。問題不在于閾值選得不好而在于這類圖像的判別信息根本不長(zhǎng)在灰度上。兩根纖維的灰度可能一模一樣區(qū)別只在方向不同、紋理走向不同。梯度結(jié)構(gòu)張量就是專門對(duì)付這類問題的工具它把鄰域內(nèi)所有梯度向量的外積做加權(quán)平均壓縮成 2×2 的對(duì)稱矩陣再用特征值和特征向量把“這里邊緣有多強(qiáng)”和“邊緣朝哪個(gè)方向”拆成兩個(gè)獨(dú)立的量。有了這兩個(gè)量濾波器就能沿邊緣方向自由擴(kuò)散、垂直邊緣方向幾乎不擴(kuò)散分割算法也可以把“方向突然變了”直接當(dāng)作邊界判據(jù)。這套思路在纖維類材料分析、遙感道路提取、組織切片分割、織物瑕疵定位里都用得上。下面的內(nèi)容適合已經(jīng)寫過基礎(chǔ)圖像處理代碼、但對(duì)張量場(chǎng)還停留在“知道有這個(gè)詞”階段的讀者從公式怎么來的講到參數(shù)怎么調(diào)、坑怎么繞。1. 方向信息比灰度更值得作為分割依據(jù)的原因1.1 灰度閾值在紋理型圖像上的失效路徑先把失效機(jī)制說透。一張纖維、道路、血管或者晶粒這類結(jié)構(gòu)的圖像像素灰度其實(shí)是“結(jié)構(gòu)反射率”和“成像噪聲”疊加的結(jié)果而結(jié)構(gòu)本身的幾何走向幾乎不體現(xiàn)在單個(gè)像素的取值上。相鄰兩根纖維的交界處兩邊的灰度可能都在 120 上下中間那條縫可能反而是 145因?yàn)榻唤缣幋嬖谏⑸浠蛘叻垂?。這時(shí)候你用 Otsu 找一個(gè)全局閾值算法只會(huì)去找直方圖上讓類間方差最大的那個(gè)點(diǎn)它壓根不知道“方向”這回事。更麻煩的是局部對(duì)比度衰減。圖像邊緣在離焦、散射、有限孔徑這些因素作用下會(huì)變成一段漸變帶寬度可能有三到五個(gè)像素。在漸變帶里梯度幅值是很小的但梯度的方向是穩(wěn)定且一致的——垂直于邊緣。這正是結(jié)構(gòu)張量要提取的東西弱梯度不可靠但大量弱梯度方向一致時(shí)統(tǒng)計(jì)上是高度可靠的。我在處理一批取向明顯的材料顯微圖時(shí)做過對(duì)比同一張圖用 Otsu 得到的 Dice 大概在 0.6 出頭主要錯(cuò)誤全部集中在相鄰纖維的粘連處和端點(diǎn)位置而在同一張圖上先跑一遍各向異性擴(kuò)散再分割Dice 能上到 0.85 以上。差別不在于閾值選得多精細(xì)而在于擴(kuò)散過程把“跨邊界的模糊”壓掉了把“沿邊界的連續(xù)性”保住了。這里有個(gè)判斷標(biāo)準(zhǔn)值得記住如果你的分割錯(cuò)誤總是集中在兩條平行結(jié)構(gòu)之間或者總是表現(xiàn)為“某某結(jié)構(gòu)的端點(diǎn)丟失”那基本可以確定灰度不是足夠的判別特征應(yīng)當(dāng)引入方向信息。1.2 結(jié)構(gòu)張量到底壓縮了什么信息可以把它理解成一種“局部梯度協(xié)方差”。在一個(gè)小窗口里你有幾百個(gè)梯度向量每個(gè)都是二維的。直接看這些向量沒法做判斷因?yàn)樗鼈冋?fù)交替、長(zhǎng)短不一。結(jié)構(gòu)張量的做法是把每個(gè)梯度向量和自己做外積得到一個(gè)小矩陣再把這些小矩陣加權(quán)平均起來。外積操作有個(gè)關(guān)鍵性質(zhì)它把方向信息保留下來但把正負(fù)號(hào)抹掉了因?yàn)?(?g)(?g)? 和 g g? 完全一樣。這恰好符合我們的需求邊緣的一側(cè)梯度往左另一側(cè)往右但我們只關(guān)心“這條邊緣的方向”不關(guān)心它朝哪邊。平均之后得到的 2×2 對(duì)稱矩陣有兩個(gè)特征值和兩個(gè)正交特征向量。大特征值對(duì)應(yīng)的特征向量指向梯度最集中的方向也就是垂直于邊緣的方向小特征值對(duì)應(yīng)的特征向量平行于邊緣。兩個(gè)特征值的數(shù)值大小分別刻畫了“垂直于邊緣方向上有多少變化”和“沿邊緣方向上有多少變化”。平坦區(qū)域兩個(gè)都趨近于零角落或者斑點(diǎn)兩個(gè)都比較大這是后續(xù)所有判據(jù)的物理基礎(chǔ)。需要強(qiáng)調(diào)一點(diǎn)結(jié)構(gòu)張量是半正定的特征值天然滿足 λ1 ≥ λ2 ≥ 0不存在負(fù)值。這個(gè)性質(zhì)在數(shù)值實(shí)現(xiàn)上有用可以省掉一些符號(hào)判斷也意味著你不可能因?yàn)樘卣髦档恼?fù)搞混方向。1.3 各向同性與各向異性的直觀差別各向同性擴(kuò)散是熱傳導(dǎo)方程那一套?I/?t ΔI每個(gè)方向擴(kuò)散速率完全一樣。它的效果就是高斯模糊噪聲被抹平了邊緣也被抹平了因?yàn)榭缭竭吘壏较虻臄U(kuò)散和沿著邊緣方向的擴(kuò)散遵循同一套規(guī)則。沾一滴油在水面上往四周均勻攤開這就是各向同性。各向異性擴(kuò)散則是把擴(kuò)散速率寫成一個(gè)張量讓它依賴位置和方向。沿著邊緣方向給一個(gè)大的擴(kuò)散系數(shù)垂直于邊緣方向給一個(gè)小的效果就是“邊緣內(nèi)部越磨越干凈邊緣兩側(cè)互不侵犯”。水滲進(jìn)木頭里沿紋理跑得快、垂直紋理跑得慢這個(gè)類比基本能說清概念。工程上各向異性通常有兩種用法很多人會(huì)混。一種是把它當(dāng)作預(yù)處理用各向異性擴(kuò)散去噪保邊然后交給 Otsu、分水嶺、區(qū)域生長(zhǎng)這些傳統(tǒng)分割器屬于“先把圖修好再分割”。另一種是把張量場(chǎng)的特征相干性、主方向直接作為分割的能量項(xiàng)或者距離度量屬于“讓分割器自己讀懂方向”。前者實(shí)現(xiàn)簡(jiǎn)單、魯棒性好后者上限更高但對(duì)參數(shù)更敏感。后面我會(huì)把兩條路都寫出來。2. 梯度結(jié)構(gòu)張量的數(shù)學(xué)骨架與關(guān)鍵參數(shù)含義2.1 從梯度外積到高斯加權(quán)三步構(gòu)造法構(gòu)造過程拆成三步每一步的意圖都不一樣。第一步是求梯度但不是簡(jiǎn)單差分。直接對(duì)原圖做中心差分會(huì)把噪聲放大得很厲害因?yàn)椴罘炙阕訉?duì)高頻的增益接近 2。標(biāo)準(zhǔn)做法是對(duì)圖像做帶高斯核的偏導(dǎo)等價(jià)于先高斯平滑再求導(dǎo)或者把高斯函數(shù)求導(dǎo)后直接和圖像卷積。這一步引入尺度參數(shù) σ它決定了梯度的噪聲抑制強(qiáng)度也決定了多小的結(jié)構(gòu)會(huì)被濾掉。第二步是做外積。對(duì)每個(gè)像素用 Ix、Iy 三個(gè)組合 Ix2、IxIy、Iy2 構(gòu)成一個(gè)對(duì)稱矩陣的獨(dú)立分量。這一步是純代數(shù)操作不引入?yún)?shù)但它是把“向量”升格成“張量”的關(guān)鍵因?yàn)橥夥e保留了方向分布的各向異性而簡(jiǎn)單的梯度平方和只會(huì)得到一個(gè)標(biāo)量。第三步是加權(quán)平均。用一個(gè)高斯核去平滑這三個(gè)分量核的寬度是尺度參數(shù) ρ。這里有個(gè)容易被忽略的細(xì)節(jié)這個(gè)平滑必須在張量分量上進(jìn)行不能先算特征值再平滑特征值也不能先平滑圖像再算張量。在張量分量上做線性平滑等價(jià)于在張量空間里做加權(quán)平均得到的仍然是合法的張量而在特征值或角度上做平滑會(huì)引入周期性歧義和插值誤差是最常見的實(shí)現(xiàn)錯(cuò)誤之一。用公式把這三步寫出來就是 J G_ρ * (?I_σ ?I_σ?)星號(hào)表示卷積?I_σ 是平滑后的梯度。整個(gè)過程只有兩個(gè)參數(shù)但這兩個(gè)參數(shù)的物理含義完全不同下一節(jié)單獨(dú)說。2.2 特征值分解兩個(gè)數(shù)字各自管什么2×2 對(duì)稱矩陣有解析的特征值公式不需要調(diào)庫這在批量處理時(shí)省很多開銷。設(shè)矩陣分量為 Jxx、Jxy、Jyy先算跡 tr Jxx Jyy 和判別式 d √((Jxx ? Jyy)2 4Jxy2)則 λ1 (tr d)/2λ2 (tr ? d)/2。因?yàn)?J 半正定這個(gè)公式不會(huì)出現(xiàn)負(fù)數(shù)開方數(shù)值上很穩(wěn)。λ1 的物理意義是“垂直于邊緣方向的梯度能量”λ2 是“沿邊緣方向的梯度能量”。對(duì)一條理想單像素寬邊緣λ2 應(yīng)該接近 0λ1 明顯大于 0。對(duì)一塊平坦區(qū)域兩者都接近 0。對(duì)一個(gè)角點(diǎn)或者交叉點(diǎn)兩個(gè)方向都有明顯梯度λ1 和 λ2 都比較大且比值接近 1。這個(gè)判據(jù)在實(shí)際分割里非常有用因?yàn)楹芏嗳蝿?wù)恰恰需要區(qū)分“邊緣”和“角點(diǎn)”比如道路網(wǎng)提取里你要的是連續(xù)的線交叉口要單獨(dú)處理。主特征向量 e1 (cosθ, sinθ)其中 θ 0.5·atan2(2Jxy, Jxx ? Jyy)。這個(gè) θ 是梯度主方向也就是垂直邊緣的方向。邊緣本身的走向是 θ π/2。很多實(shí)現(xiàn)里直接拿 θ 當(dāng)“邊緣角度”用結(jié)果畫出來的方向場(chǎng)整體差了 90 度看起來像所有線條都跟真實(shí)結(jié)構(gòu)垂直這是第 4 節(jié)要細(xì)講的一個(gè)坑。需要提醒的是特征值的大小是相對(duì)的跟圖像的灰度尺度直接掛鉤。你把輸入圖像除以 255 歸一化λ 的量級(jí)就跟著變擴(kuò)散參數(shù) k 必須重新標(biāo)定。這是很多人換數(shù)據(jù)集后效果突然崩掉的真實(shí)原因不是算法不行是量綱變了。2.3 相干性把兩個(gè)特征值壓成一個(gè)可解釋的標(biāo)量單純看 λ1 有幾個(gè)問題它受圖像整體亮度影響也受局部對(duì)比度影響不同區(qū)域的數(shù)值不可比。相干性 C (λ1 ? λ2)/(λ1 λ2) 就是為了解決這個(gè)而設(shè)計(jì)的它的取值范圍是 [0, 1]與亮度尺度無關(guān)。C 接近 1 表示結(jié)構(gòu)高度有序梯度幾乎全部集中在一個(gè)方向上典型情況是清晰筆直的邊緣。C 接近 0 有兩種截然不同的成因必須分開對(duì)待一是各向同性噪聲區(qū)兩個(gè)特征值都小且接近二是角點(diǎn)或者交叉結(jié)構(gòu)兩個(gè)特征值都大且接近。前者應(yīng)該被重度平滑后者必須被保護(hù)。這個(gè)區(qū)別光看 C 是分不出來的還得看 λ1 λ2 的絕對(duì)大小。我通常的做法是把 (λ1 λ2, C) 這兩維做成一張散點(diǎn)圖或者二維直方圖噪聲會(huì)聚在左下角邊緣沿右側(cè)豎條分布角點(diǎn)落在右上角。這張圖能讓你在五分鐘內(nèi)判斷出當(dāng)前參數(shù)下哪些區(qū)域會(huì)被誤處理比盲目試參數(shù)高效得多。數(shù)值上必須加一個(gè)小常數(shù)C (λ1 ? λ2)/(λ1 λ2 ε)。ε 取 λ 量級(jí)的千分之一左右。不加這個(gè) ε平坦區(qū)域會(huì)出現(xiàn) 0/0 產(chǎn)生 NaNNaN 再經(jīng)過擴(kuò)散迭代會(huì)擴(kuò)散到整張圖最后整幅圖變成一片灰這個(gè) bug 極其隱蔽因?yàn)槌绦虿粓?bào)錯(cuò)。2.4 兩個(gè)尺度 σ 與 ρ 的分工與常見誤配σ 和 ρ 分開是有道理的雖然很多簡(jiǎn)化實(shí)現(xiàn)把兩者設(shè)成相等。σ 管的是梯度估計(jì)的質(zhì)量它要足夠大以壓制噪聲又不能大到把窄結(jié)構(gòu)整體模糊掉。ρ 管的是方向統(tǒng)計(jì)的空間范圍它決定了“多大范圍內(nèi)的梯度被認(rèn)為應(yīng)該方向一致”。一般經(jīng)驗(yàn)是 ρ 取 2σ 到 4σ。原因很直接如果 ρ 太小參與平均的獨(dú)立梯度樣本太少方向估計(jì)的方差大方向場(chǎng)會(huì)花如果 ρ 太大跨過相鄰平行結(jié)構(gòu)做平均兩條本來方向不同的纖維會(huì)被平均成一個(gè)中間方向交界處的判別能力直接歸零。當(dāng)你知道目標(biāo)結(jié)構(gòu)的典型寬度 w以像素計(jì)ρ 應(yīng)該控制在 w/2 以內(nèi)。舉個(gè)具體例子一張 1024×1024 的纖維圖像纖維寬度約 12 像素纖維間距約 6 像素。σ 取 1.0 壓制傳感器噪聲ρ 取 3.0 可以在約 19 個(gè)像素的窗口內(nèi)統(tǒng)計(jì)方向此時(shí)窗口半徑約 9 像素還不到纖維寬度不會(huì)跨結(jié)構(gòu)。如果把 ρ 提到 6窗口直徑 24 像素穩(wěn)穩(wěn)跨過間距交界處的方向估計(jì)就廢了。還有一個(gè)容易忽略的點(diǎn)這兩次高斯平滑的邊界處理方式。默認(rèn)的反射填充在圖像邊緣會(huì)制造虛假的對(duì)稱結(jié)構(gòu)導(dǎo)致最外圈幾行的方向估計(jì)完全錯(cuò)誤。如果你的目標(biāo)結(jié)構(gòu)貼著圖像邊緣務(wù)必用裁剪后再處理的方式或者干脆把邊緣區(qū)域的張量置為無效、不參與后續(xù)決策。3. 從張量場(chǎng)到各向異性分割的完整實(shí)現(xiàn)3.1 整體流程與兩條技術(shù)路線先把兩條路線的骨架擺清楚再逐段給代碼。路線一工程上更常用讀圖 → 灰度歸一化 → 計(jì)算梯度結(jié)構(gòu)張量 → 由張量構(gòu)造擴(kuò)散張量 → 各向異性擴(kuò)散迭代 → 在擴(kuò)散結(jié)果上做分水嶺或區(qū)域生長(zhǎng) → 后處理去除小塊、孔洞填充→ 輸出標(biāo)簽圖。特點(diǎn)是每一步都有成熟實(shí)現(xiàn)調(diào)試友好失敗時(shí)能快速定位到是擴(kuò)散沒做好還是分割器沒調(diào)好。路線二上限更高讀圖 → 計(jì)算張量 → 提取相干性圖和方向場(chǎng) → 構(gòu)造測(cè)地距離或者邊界停止函數(shù) → 活動(dòng)輪廓/水平集演化 → 輸出標(biāo)簽圖。特點(diǎn)是不需要先做一次可能損傷細(xì)節(jié)的平滑直接在原始信息上做決策但參數(shù)敏感性高收斂性依賴于初始輪廓。實(shí)際項(xiàng)目里我通常先做路線一因?yàn)榭煽康幕€比炫技重要。等基線跑通、評(píng)估指標(biāo)穩(wěn)定之后再用路線二在難點(diǎn)區(qū)域比如粘連嚴(yán)重的地方做局部增強(qiáng)兩者用同一套評(píng)估腳本對(duì)比避免自我感覺良好。下面按路線一展開路線二在 3.4 節(jié)單獨(dú)給出核心代碼。3.2 張量計(jì)算與方向場(chǎng)可視化的可復(fù)現(xiàn)代碼先給張量計(jì)算用 SciPy 就夠了不需要額外的圖像處理庫。import numpy as np from scipy.ndimage import gaussian_filter def structure_tensor(img, sigma1.0, rho3.0): img img.astype(np.float64) # 帶高斯導(dǎo)數(shù)的偏導(dǎo)等價(jià)于先平滑再求導(dǎo) Ix gaussian_filter(img, sigma, order(0, 1), modereflect) Iy gaussian_filter(img, sigma, order(1, 0), modereflect) # 在張量分量上做第二尺度平滑 Jxx gaussian_filter(Ix * Ix, rho, modereflect) Jxy gaussian_filter(Ix * Iy, rho, modereflect) Jyy gaussian_filter(Iy * Iy, rho, modereflect) return Jxx, Jxy, Jyy def tensor_eig(Jxx, Jxy, Jyy, eps1e-10): tr Jxx Jyy diff Jxx - Jyy disc np.sqrt(diff * diff 4.0 * Jxy * Jxy) l1 0.5 * (tr disc) l2 0.5 * (tr - disc) # 主特征向量對(duì)應(yīng) l1與 x 軸的夾角即梯度主方向 theta 0.5 * np.arctan2(2.0 * Jxy, diff) coherence (l1 - l2) / (l1 l2 eps) return l1, l2, theta, coherence注意order(0, 1)的含義是沿第二個(gè)軸列方向也就是 x求一階導(dǎo)第一個(gè)軸不導(dǎo)order(1, 0)是沿第一個(gè)軸行方向y求導(dǎo)。搞反了不影響 λ 的分布但會(huì)讓方向場(chǎng)整體旋轉(zhuǎn) 90 度屬于典型的靜默錯(cuò)誤。想把方向場(chǎng)畫出來檢查最直觀的方式是抽稀后畫線段疊加到灰度圖上。import matplotlib.pyplot as plt def plot_orientation(img, theta, step16, length8): plt.figure(figsize(8, 8)) plt.imshow(img, cmapgray) h, w img.shape for y in range(step // 2, h, step): for x in range(step // 2, w, step): t theta[y, x] np.pi / 2.0 # 轉(zhuǎn)成邊緣走向 dx np.cos(t) * length / 2 dy np.sin(t) * length / 2 plt.plot([x - dx, x dx], [y - dy, y dy], r-, lw0.8) plt.axis(off) plt.show()抽稀步長(zhǎng) step 是必須的逐像素畫線會(huì)糊成一片紅。步長(zhǎng)取 12 到 20 之間通常剛好能看清趨勢(shì)又不遮擋原圖。這張圖在整個(gè)調(diào)試過程中價(jià)值極高方向場(chǎng)對(duì)了后面基本就順了方向場(chǎng)不對(duì)再怎么調(diào)擴(kuò)散參數(shù)都是白費(fèi)。3.3 各向異性擴(kuò)散的離散實(shí)現(xiàn)與穩(wěn)定性約束擴(kuò)散張量由張量場(chǎng)的特征向量和特征值共同構(gòu)造。核心思路是沿邊緣方向θ π/2給大擴(kuò)散系數(shù) d_edge沿梯度方向θ給小擴(kuò)散系數(shù) d_grad后者隨 λ1 增大而減小。def anisotropic_diffusion(img, l1, theta, k20.0, dt0.15, iters60): I img.astype(np.float64).copy() # 梯度方向抑制系數(shù)λ1 越大抑制越強(qiáng) d_grad 1.0 / (1.0 l1 / (k * k)) d_edge np.ones_like(d_grad) c, s np.cos(theta), np.sin(theta) D11 d_grad * c * c d_edge * s * s D12 (d_grad - d_edge) * c * s D22 d_grad * s * s d_edge * c * c for _ in range(iters): Ix 0.5 * (np.roll(I, -1, axis1) - np.roll(I, 1, axis1)) Iy 0.5 * (np.roll(I, -1, axis0) - np.roll(I, 1, axis0)) fx D11 * Ix D12 * Iy fy D12 * Ix D22 * Iy div (0.5 * (np.roll(fx, -1, axis1) - np.roll(fx, 1, axis1)) 0.5 * (np.roll(fy, -1, axis0) - np.roll(fy, 1, axis0))) I I dt * div return I這段代碼用np.roll意味著周期性邊界對(duì)絕大多數(shù)圖像是不對(duì)的因?yàn)閳D像的左邊緣和右邊緣毫無關(guān)系。生產(chǎn)代碼里應(yīng)該用顯式的切片索引實(shí)現(xiàn) Neumann 邊界或者簡(jiǎn)單地在擴(kuò)散前對(duì)外圍 1 像素做鏡像擴(kuò)展擴(kuò)散完成后再裁掉。用 roll 的版本跑出來的圖左右邊緣會(huì)出現(xiàn)莫名其妙的粘連如果你只關(guān)注中間區(qū)域可能看不出來但一旦做全圖統(tǒng)計(jì)就會(huì)暴露。穩(wěn)定性是必須算清楚的一件事。顯式格式的二維擴(kuò)散方程穩(wěn)定條件近似為 dt ≤ 1/(4·max(D))。本例中 d_edge 1所以 max(D) 最大為 1dt ≤ 0.25。取 0.15 留了約 40% 的余量比較安全。如果把 d_edge 設(shè)成 2 或者更高dt 必須相應(yīng)縮小否則迭代幾輪之后數(shù)值爆炸圖像出現(xiàn)棋盤狀條紋再迭代下去全是 NaN。迭代次數(shù)不是越多越好。經(jīng)驗(yàn)值是 40 到 100 之間具體看圖像尺寸和目標(biāo)結(jié)構(gòu)的尺度。判斷標(biāo)準(zhǔn)不是固定次數(shù)而是看連續(xù)兩次迭代之間的變化量當(dāng) max|I_new ? I_old| 小于灰度范圍千分之一時(shí)就可以停這時(shí)繼續(xù)迭代只會(huì)延長(zhǎng)運(yùn)行時(shí)間并輕微模糊邊界。可以加一個(gè)簡(jiǎn)單的收斂檢查。關(guān)于 k 的取值它是擴(kuò)散的“對(duì)比度閾值”物理意義是“λ1 達(dá)到多大時(shí)認(rèn)為這是一條需要保護(hù)的邊緣”。最靠譜的選法是統(tǒng)計(jì) λ1 在整幅圖上的分布取 80% 到 90% 分位數(shù)作為初值。如果全圖大部分是平坦區(qū)域分位數(shù)會(huì)偏低此時(shí)擴(kuò)散會(huì)很保守如果圖里到處都是紋理分位數(shù)偏高擴(kuò)散會(huì)很激進(jìn)。兩種情況都得靠肉眼掃一遍若干張樣本的擴(kuò)散結(jié)果來最終確認(rèn)。3.4 把張量特征直接喂給分割器擴(kuò)散做完之后分割這一步反而簡(jiǎn)單了因?yàn)檎尺B已經(jīng)被削弱。分水嶺是首選但直接對(duì)灰度做分水嶺容易過分割需要可靠的標(biāo)記?？梢杂镁植繕O大值做標(biāo)記也可以先把擴(kuò)散結(jié)果做一遍形態(tài)學(xué)開運(yùn)算再取標(biāo)記。from scipy import ndimage as ndi from skimage.feature import peak_local_max from skimage.segmentation import watershed def segment_by_watershed(enhanced, min_distance12): # 前景掩膜用 Otsu 在擴(kuò)散結(jié)果上取值 thresh enhanced (enhanced.min() enhanced.max()) / 2.0 dist ndi.distance_transform_edt(thresh) coords peak_local_max(dist, min_distancemin_distance, labelsthresh) markers np.zeros(enhanced.shape, dtypenp.int32) for i, (y, x) in enumerate(coords, start1): markers[y, x] i labels watershed(-dist, markers, maskthresh) return labels路線二的做法是用相干性構(gòu)造邊界代價(jià)。相干性高說明這里是一條可靠邊緣交叉它應(yīng)該付出大代價(jià)。把這個(gè)代價(jià)轉(zhuǎn)成測(cè)地距離再做區(qū)域生長(zhǎng)就不需要先平滑圖像。def geodesic_cost(coherence, alpha6.0): # 相干性越高穿越代價(jià)越大 return np.exp(-alpha * np.clip(coherence, 0, 1))alpha 控制對(duì)邊界的敏感度。alpha 太大代價(jià)圖接近二值弱邊界完全被忽略alpha 太小處處都是軟邊界區(qū)域會(huì)漏出去。經(jīng)驗(yàn)區(qū)間是 4 到 10先用中等值跑一遍看代價(jià)圖的對(duì)比度再?zèng)Q定往哪邊調(diào)。還有一種更徹底的用法把 (λ1, λ2, C, θ) 四個(gè)量做特征標(biāo)準(zhǔn)化后送進(jìn) k-means 或者超像素聚類。這條路在紋理類別多、且不同類別方向分布差異明顯的場(chǎng)景下很好用比如把織物圖像分成經(jīng)紗、緯紗、背景三類。缺點(diǎn)是完全丟棄了空間連通性聚類結(jié)果需要再做一次連通域處理小碎塊非常多。我的習(xí)慣是先用超像素把空間約束加進(jìn)去再在超像素特征上聚類能省掉一半后處理工作。3.5 參數(shù)選擇的計(jì)算過程與一次完整標(biāo)定把參數(shù)標(biāo)定做成一個(gè)固定流程比每次憑感覺調(diào)高效得多。第一步估算噪聲水平。取圖像中一小塊明顯平坦的區(qū)域算標(biāo)準(zhǔn)差 s。這一步?jīng)Q定了 σ 的下限。第二步測(cè)量目標(biāo)結(jié)構(gòu)寬度。在圖像上量出典型結(jié)構(gòu)的像素寬度 w以及相鄰結(jié)構(gòu)之間的間距 g。這一步?jīng)Q定了 ρ 的上限。第三步按約束定參數(shù)。σ 取 max(1.0, 0.5·s) 到 2.0 之間實(shí)測(cè)中 1.0 到 1.5 覆蓋了大多數(shù)情況。ρ 取 2σ 到 3σ同時(shí)強(qiáng)制滿足 ρ ≤ min(w, g)/2。第四步定 k。跑一遍張量計(jì)算取 λ1 的 85% 分位數(shù)。如果這個(gè)值明顯低于噪聲區(qū)的 λ1 水平說明圖像對(duì)比度低需要先做一次對(duì)比度歸一化。第五步定迭代次數(shù)和 dt。dt 固定 0.15迭代從 40 起觀察收斂曲線的拐點(diǎn)。舉個(gè)具體數(shù)字某批圖像量出來 s ≈ 48 位灰度w ≈ 14 像素g ≈ 7 像素。則 σ 取 1.2ρ 取 3.0滿足 3.0 ≤ 7/2 3.5卡得比較緊但可行k 取 λ1 的 85% 分位數(shù)約 180在歸一化前的灰度尺度下dt 0.15迭代 60 次。這套參數(shù)在該批數(shù)據(jù)上穩(wěn)定復(fù)現(xiàn)換批次時(shí)只需要重新測(cè) s 和 w。4. 實(shí)操中的典型坑與排查技巧實(shí)錄4.1 方向場(chǎng)整體偏轉(zhuǎn) 90 度這是出現(xiàn)頻率最高的錯(cuò)誤癥狀是畫出來的方向線全部垂直于真實(shí)結(jié)構(gòu)。原因幾乎總是把梯度主方向和邊緣走向搞混了。θ 0.5·atan2(2Jxy, Jxx ? Jyy) 得到的是對(duì)應(yīng) λ1 的特征向量角度而 λ1 是垂直于邊緣方向的梯度能量所以 θ 是法向邊緣走向必須是 θ π/2。這個(gè)錯(cuò)誤在純擴(kuò)散的代碼里可能不表現(xiàn)出來因?yàn)槿绻阍跇?gòu)造 D 的時(shí)候同時(shí)對(duì)兩個(gè)方向都做了相應(yīng)處理可能剛好抵消。但只要你在任何一處用了 θ 去表示邊緣方向就會(huì)出現(xiàn)問題。防御性做法是在代碼里顯式命名兩個(gè)變量normal_dir和tangent_dir永遠(yuǎn)不要用裸的theta傳遞方向。4.2 尺度不匹配導(dǎo)致交界處方向被平均癥狀是相鄰兩條平行結(jié)構(gòu)之間的邊界模糊分割結(jié)果在這個(gè)位置系統(tǒng)性偏移。根本原因是 ρ 太大方向統(tǒng)計(jì)窗口跨過了間距。診斷方法很直接在你懷疑出問題的交界處把 ρ 從小到大掃一遍看方向場(chǎng)在那個(gè)位置是怎么變化的。如果小 ρ 時(shí)兩條結(jié)構(gòu)的方向清晰可辨大 ρ 時(shí)變成同一個(gè)中間角度就坐實(shí)了這個(gè)問題。修法有兩種。一是把 ρ 降到 min(w, g)/2 以下二是改用自適應(yīng)的 ρ在 λ1 高的區(qū)域用小的 ρ在平坦區(qū)域用大的 ρ。第二種做法實(shí)現(xiàn)上稍麻煩需要兩套不同尺度的張量場(chǎng)并做加權(quán)混合但效果確實(shí)更好尤其是在結(jié)構(gòu)尺度變化大的圖像上。4.3 圖像邊緣的虛假結(jié)構(gòu)反射邊界填充會(huì)讓圖像最外圈出現(xiàn)對(duì)稱的鏡像結(jié)構(gòu)梯度在這里是虛假的。癥狀是圖像四周出現(xiàn)一圈方向紊亂的區(qū)域且這個(gè)紊亂在擴(kuò)散后會(huì)往內(nèi)部傳播幾行。解決方式在實(shí)現(xiàn)層面很簡(jiǎn)單擴(kuò)散和張量計(jì)算都用modenearest或者干脆用零填充然后在所有后續(xù)統(tǒng)計(jì)和可視化里把最外圈大約 2ρ 寬度的一圈裁掉。如果你的目標(biāo)結(jié)構(gòu)本身就貼著圖像邊界那就只能在采集階段解決多留一點(diǎn)余量。4.4 平坦區(qū)的相干性數(shù)值噪聲在幾乎沒有梯度的區(qū)域λ1 和 λ2 都是浮點(diǎn)誤差量級(jí)的小數(shù)它們的差值同樣是誤差量級(jí)。此時(shí) C (λ1 ? λ2)/(λ1 λ2 ε) 的結(jié)果完全由 ε 決定可能呈現(xiàn)出隨機(jī)的、看起來很像結(jié)構(gòu)的高相干值如果用它做邊界代價(jià)平坦區(qū)會(huì)莫名其妙地被分割成碎片。修法是給相干性加一個(gè)有效性掩膜只在 λ1 λ2 超過某個(gè)絕對(duì)閾值的位置認(rèn)為 C 可信其余位置直接置零或者填一個(gè)中性值。這個(gè)閾值取噪聲區(qū) λ1 的 99% 分位數(shù)用一個(gè)平坦區(qū)域采樣就能估出來。4.5 常見問題速查表現(xiàn)象最可能的原因排查動(dòng)作修正方向方向線全部垂直真實(shí)結(jié)構(gòu)梯度方向與邊緣走向混用打印 θ 分布并疊加原圖θ 加 π/2 得到邊緣走向平行結(jié)構(gòu)交界處分割偏移ρ 跨過結(jié)構(gòu)間距掃 ρ 看交界方向場(chǎng)ρ 降到間距一半以下四周方向紊亂并向內(nèi)蔓延邊界填充方式導(dǎo)致的虛假梯度檢查最外圈換 nearest 填充并裁邊圖像出現(xiàn)棋盤狀條紋或 NaNdt 過大導(dǎo)致數(shù)值發(fā)散檢查迭代中的極值dt ≤ 0.25或按 max(D) 重算平坦區(qū)被切碎相干性在零梯度處失真看低梯度區(qū)的 C 直方圖用 λ1λ2 做有效性掩膜換數(shù)據(jù)集后效果崩掉灰度歸一化改變導(dǎo)致 k 失效比較兩批數(shù)據(jù) λ1 分位數(shù)重新標(biāo)定 k擴(kuò)散后整圖均勻變灰分母無 ε 產(chǎn)生 NaN 并擴(kuò)散檢查 NaN 數(shù)量分母加 ε迭代很多次邊界仍模糊d_grad 未隨 λ1 有效下降打印 d_grad 直方圖調(diào)小 k 或改抑制函數(shù)5. 效果評(píng)估、場(chǎng)景適配與性能優(yōu)化5.1 用什么指標(biāo)說明“方向用得對(duì)不對(duì)”分割任務(wù)常用的 Dice 和 IoU 當(dāng)然要報(bào)但它們只反映區(qū)域重疊反映不出方向信息有沒有真正起作用。我習(xí)慣再加兩個(gè)診斷性指標(biāo)。一是邊界法向一致性在參考邊界上采樣取每個(gè)采樣點(diǎn)的切向量和該位置張量場(chǎng)給出的邊緣走向做夾角統(tǒng)計(jì)。如果夾角中位數(shù)在 10 度以內(nèi)說明方向場(chǎng)是準(zhǔn)的分割失敗一定出在別的地方如果夾角分布很散方向場(chǎng)本身就有問題調(diào)分割器沒意義。二是區(qū)域內(nèi)的平均相干性。真實(shí)的結(jié)構(gòu)區(qū)域相干性應(yīng)該顯著高于背景如果某個(gè)分割出來的區(qū)域平均相干性接近背景水平那這個(gè)區(qū)域大概率是誤檢即使它的 Dice 貢獻(xiàn)看起來還行。這兩個(gè)指標(biāo)的好處是能定位問題層次。Dice 差可能是方向場(chǎng)錯(cuò)也可能是分割器參數(shù)錯(cuò)但加上這兩個(gè)診斷量就能立刻區(qū)分。5.2 不同場(chǎng)景下的參數(shù)傾向不同應(yīng)用領(lǐng)域?qū)?shù)的偏好差別很大強(qiáng)行用一套參數(shù)跨領(lǐng)域是不現(xiàn)實(shí)的。場(chǎng)景類型σ 傾向ρ 傾向k 傾向迭代次數(shù)主要關(guān)注點(diǎn)纖維/織物顯微圖1.0 到 1.52σ 到 3σ受間距強(qiáng)約束λ1 的 85% 分位50 到 80交界處不粘連遙感道路提取1.5 到 2.0偏大4σ 也可接受偏小保護(hù)弱邊緣80 到 120弱對(duì)比路段的連續(xù)性醫(yī)學(xué)組織切片1.0 到 2.0中等看細(xì)胞尺度λ1 的 80% 分位40 到 60不破壞細(xì)小結(jié)構(gòu)材料晶粒邊界1.0 到 1.2偏小晶界很窄偏大抑制晶內(nèi)紋理60 到 100晶界銳度金屬表面缺陷1.5 到 2.5偏大中30 到 50缺陷方向與非缺陷區(qū)分這張表不是標(biāo)準(zhǔn)答案是起點(diǎn)。我建議的做法是先用表格里的一行跑通再針對(duì)自己數(shù)據(jù)里最差的那 10% 樣本單獨(dú)調(diào)把調(diào)好的參數(shù)寫進(jìn)配置文件按數(shù)據(jù)批次管理不要硬編碼在腳本里。5.3 性能與工程落地高斯濾波在 SciPy 里是可分離的復(fù)雜度是 O(N·r) 而不是 O(N·r2)但 ρ 大的時(shí)候仍然慢。如果 ρ 超過 8用盒式濾波近似能快好幾倍代價(jià)是方向場(chǎng)的平滑性稍差做方向場(chǎng)可視化時(shí)能看出塊狀。批量處理時(shí)角度和相干性的計(jì)算都是逐像素的算術(shù)運(yùn)算可以用 NumPy 向量化一次完成不要寫循環(huán)。擴(kuò)散迭代是性能瓶頸。60 次迭代、每次若干次數(shù)組操作對(duì) 1024×1024 的圖像大概在秒級(jí)可以接受對(duì) 4K 圖像就明顯吃力了。幾個(gè)提速方向一是把擴(kuò)散放到金字塔的粗層先跑再上采樣作為細(xì)層初始值總迭代次數(shù)能省一半以上二是擴(kuò)散張量 D 在整個(gè)迭代過程中通常是固定的可以預(yù)先算好并緩存不必每輪重算三是如果必須實(shí)時(shí)轉(zhuǎn)向基于 CUDA 的實(shí)現(xiàn)或者把張量計(jì)算和擴(kuò)散都寫成批處理的張量運(yùn)算。還有一個(gè)工程細(xì)節(jié)值得提前設(shè)計(jì)把整個(gè)流程封裝成“參數(shù)對(duì)象 處理函數(shù)”的形式參數(shù)對(duì)象序列化成 JSON 和數(shù)據(jù)版本綁定。這套東西在單次實(shí)驗(yàn)里顯得多余但當(dāng)你需要回頭復(fù)現(xiàn)三個(gè)月前某個(gè)結(jié)果的時(shí)候會(huì)發(fā)現(xiàn)省下的時(shí)間遠(yuǎn)超當(dāng)初寫它的時(shí)間。最后分享一個(gè)我在實(shí)際使用中養(yǎng)成的小習(xí)慣。每次調(diào)參之前先挑三張有代表性的樣本一張結(jié)構(gòu)清晰、一張?jiān)肼曋亍⒁粡堈尺B嚴(yán)重。三張一起跑看方向場(chǎng)疊加圖和相干性直方圖只根據(jù)這三張的反饋調(diào)參數(shù)絕不在全量數(shù)據(jù)上反復(fù)試。全量跑一輪的反饋周期太長(zhǎng)你會(huì)失去對(duì)參數(shù)敏感度的直覺。三張樣本的快速循環(huán)跑上二三十輪對(duì)這套方法的脾性基本就摸透了之后再上新數(shù)據(jù)往往第一組參數(shù)就能到可用的水平。