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