做显微图像分割的人大概都遇到过同一种窘境：灰度直方图上两个峰的间隔还算明显，但阈值一提上去，纤维束的走向就断了，一提下来又把相邻的两根纤维糊成一片。问题不在于阈值选得不好，而在于这类图像的判别信息根本不长在灰度上。两根纤维的灰度可能一模一样，区别只在方向不同、纹理走向不同。梯度结构张量就是专门对付这类问题的工具，它把邻域内所有梯度向量的外积做加权平均，压缩成 2×2 的对称矩阵，再用特征值和特征向量把“这里边缘有多强”和“边缘朝哪个方向”拆成两个独立的量。有了这两个量，滤波器就能沿边缘方向自由扩散、垂直边缘方向几乎不扩散，分割算法也可以把“方向突然变了”直接当作边界判据。这套思路在纤维类材料分析、遥感道路提取、组织切片分割、织物瑕疵定位里都用得上。下面的内容适合已经写过基础图像处理代码、但对张量场还停留在“知道有这个词”阶段的读者，从公式怎么来的讲到参数怎么调、坑怎么绕。

1. 方向信息比灰度更值得作为分割依据的原因

1.1 灰度阈值在纹理型图像上的失效路径

先把失效机制说透。一张纤维、道路、血管或者晶粒这类结构的图像，像素灰度其实是“结构反射率”和“成像噪声”叠加的结果，而结构本身的几何走向几乎不体现在单个像素的取值上。相邻两根纤维的交界处，两边的灰度可能都在 120 上下，中间那条缝可能反而是 145，因为交界处存在散射或者反光。这时候你用 Otsu 找一个全局阈值，算法只会去找直方图上让类间方差最大的那个点，它压根不知道“方向”这回事。

更麻烦的是局部对比度衰减。图像边缘在离焦、散射、有限孔径这些因素作用下会变成一段渐变带，宽度可能有三到五个像素。在渐变带里，梯度幅值是很小的，但梯度的方向是稳定且一致的——垂直于边缘。这正是结构张量要提取的东西：弱梯度不可靠，但大量弱梯度方向一致时，统计上是高度可靠的。

我在处理一批取向明显的材料显微图时做过对比，同一张图用 Otsu 得到的 Dice 大概在 0.6 出头，主要错误全部集中在相邻纤维的粘连处和端点位置；而在同一张图上先跑一遍各向异性扩散再分割，Dice 能上到 0.85 以上。差别不在于阈值选得多精细，而在于扩散过程把“跨边界的模糊”压掉了，把“沿边界的连续性”保住了。

这里有个判断标准值得记住：如果你的分割错误总是集中在两条平行结构之间，或者总是表现为“某某结构的端点丢失”，那基本可以确定灰度不是足够的判别特征，应当引入方向信息。

1.2 结构张量到底压缩了什么信息

可以把它理解成一种“局部梯度协方差”。在一个小窗口里，你有几百个梯度向量，每个都是二维的。直接看这些向量没法做判断，因为它们正负交替、长短不一。结构张量的做法是把每个梯度向量和自己做外积，得到一个小矩阵，再把这些小矩阵加权平均起来。外积操作有个关键性质：它把方向信息保留下来，但把正负号抹掉了，因为 (−g)(−g)ᵀ 和 g gᵀ 完全一样。这恰好符合我们的需求，边缘的一侧梯度往左，另一侧往右，但我们只关心“这条边缘的方向”，不关心它朝哪边。

平均之后得到的 2×2 对称矩阵有两个特征值和两个正交特征向量。大特征值对应的特征向量指向梯度最集中的方向，也就是垂直于边缘的方向；小特征值对应的特征向量平行于边缘。两个特征值的数值大小，分别刻画了“垂直于边缘方向上有多少变化”和“沿边缘方向上有多少变化”。平坦区域两个都趋近于零，角落或者斑点两个都比较大，这是后续所有判据的物理基础。

需要强调一点，结构张量是半正定的，特征值天然满足 λ1 ≥ λ2 ≥ 0，不存在负值。这个性质在数值实现上有用，可以省掉一些符号判断，也意味着你不可能因为特征值的正负搞混方向。

1.3 各向同性与各向异性的直观差别

各向同性扩散是热传导方程那一套，∂I/∂t = ΔI，每个方向扩散速率完全一样。它的效果就是高斯模糊，噪声被抹平了，边缘也被抹平了，因为跨越边缘方向的扩散和沿着边缘方向的扩散遵循同一套规则。沾一滴油在水面上往四周均匀摊开，这就是各向同性。

各向异性扩散则是把扩散速率写成一个张量，让它依赖位置和方向。沿着边缘方向给一个大的扩散系数，垂直于边缘方向给一个小的，效果就是“边缘内部越磨越干净，边缘两侧互不侵犯”。水渗进木头里沿纹理跑得快、垂直纹理跑得慢，这个类比基本能说清概念。

工程上各向异性通常有两种用法，很多人会混。一种是把它当作预处理，用各向异性扩散去噪保边，然后交给 Otsu、分水岭、区域生长这些传统分割器，属于“先把图修好再分割”。另一种是把张量场的特征（相干性、主方向）直接作为分割的能量项或者距离度量，属于“让分割器自己读懂方向”。前者实现简单、鲁棒性好，后者上限更高但对参数更敏感。后面我会把两条路都写出来。

2. 梯度结构张量的数学骨架与关键参数含义

2.1 从梯度外积到高斯加权：三步构造法

构造过程拆成三步，每一步的意图都不一样。

第一步是求梯度，但不是简单差分。直接对原图做中心差分会把噪声放大得很厉害，因为差分算子对高频的增益接近 2。标准做法是对图像做带高斯核的偏导，等价于先高斯平滑再求导，或者把高斯函数求导后直接和图像卷积。这一步引入尺度参数 σ，它决定了梯度的噪声抑制强度，也决定了多小的结构会被滤掉。

第二步是做外积。对每个像素，用 Ix、Iy 三个组合 Ix²、IxIy、Iy² 构成一个对称矩阵的独立分量。这一步是纯代数操作，不引入参数，但它是把“向量”升格成“张量”的关键，因为外积保留了方向分布的各向异性，而简单的梯度平方和只会得到一个标量。

第三步是加权平均。用一个高斯核去平滑这三个分量，核的宽度是尺度参数 ρ。这里有个容易被忽略的细节：这个平滑必须在张量分量上进行，不能先算特征值再平滑特征值，也不能先平滑图像再算张量。在张量分量上做线性平滑，等价于在张量空间里做加权平均，得到的仍然是合法的张量；而在特征值或角度上做平滑，会引入周期性歧义和插值误差，是最常见的实现错误之一。

用公式把这三步写出来就是 J = G_ρ * (∇I_σ ∇I_σᵀ)，星号表示卷积，∇I_σ 是平滑后的梯度。整个过程只有两个参数，但这两个参数的物理含义完全不同，下一节单独说。

2.2 特征值分解：两个数字各自管什么

2×2 对称矩阵有解析的特征值公式，不需要调库，这在批量处理时省很多开销。设矩阵分量为 Jxx、Jxy、Jyy，先算迹 tr = Jxx + Jyy 和判别式 d = √((Jxx − Jyy)² + 4Jxy²)，则 λ1 = (tr + d)/2，λ2 = (tr − d)/2。因为 J 半正定，这个公式不会出现负数开方，数值上很稳。

λ1 的物理意义是“垂直于边缘方向的梯度能量”，λ2 是“沿边缘方向的梯度能量”。对一条理想单像素宽边缘，λ2 应该接近 0，λ1 明显大于 0。对一块平坦区域，两者都接近 0。对一个角点或者交叉点，两个方向都有明显梯度，λ1 和 λ2 都比较大且比值接近 1。这个判据在实际分割里非常有用，因为很多任务恰恰需要区分“边缘”和“角点”，比如道路网提取里你要的是连续的线，交叉口要单独处理。

主特征向量 e1 = (cosθ, sinθ)，其中 θ = 0.5·atan2(2Jxy, Jxx − Jyy)。这个 θ 是梯度主方向，也就是垂直边缘的方向。边缘本身的走向是 θ + π/2。很多实现里直接拿 θ 当“边缘角度”用，结果画出来的方向场整体差了 90 度，看起来像所有线条都跟真实结构垂直，这是第 4 节要细讲的一个坑。

需要提醒的是，特征值的大小是相对的，跟图像的灰度尺度直接挂钩。你把输入图像除以 255 归一化，λ 的量级就跟着变，扩散参数 k 必须重新标定。这是很多人换数据集后效果突然崩掉的真实原因，不是算法不行，是量纲变了。

2.3 相干性：把两个特征值压成一个可解释的标量

单纯看 λ1 有几个问题：它受图像整体亮度影响，也受局部对比度影响，不同区域的数值不可比。相干性 C = (λ1 − λ2)/(λ1 + λ2) 就是为了解决这个而设计的，它的取值范围是 [0, 1]，与亮度尺度无关。

C 接近 1 表示结构高度有序，梯度几乎全部集中在一个方向上，典型情况是清晰笔直的边缘。C 接近 0 有两种截然不同的成因，必须分开对待：一是各向同性噪声区，两个特征值都小且接近；二是角点或者交叉结构，两个特征值都大且接近。前者应该被重度平滑，后者必须被保护。这个区别光看 C 是分不出来的，还得看 λ1 + λ2 的绝对大小。

我通常的做法是把 (λ1 + λ2, C) 这两维做成一张散点图或者二维直方图，噪声会聚在左下角，边缘沿右侧竖条分布，角点落在右上角。这张图能让你在五分钟内判断出当前参数下哪些区域会被误处理，比盲目试参数高效得多。

数值上必须加一个小常数：C = (λ1 − λ2)/(λ1 + λ2 + ε)。ε 取 λ 量级的千分之一左右。不加这个 ε，平坦区域会出现 0/0 产生 NaN，NaN 再经过扩散迭代会扩散到整张图，最后整幅图变成一片灰，这个 bug 极其隐蔽，因为程序不报错。

2.4 两个尺度 σ 与 ρ 的分工与常见误配

σ 和 ρ 分开是有道理的，虽然很多简化实现把两者设成相等。σ 管的是梯度估计的质量，它要足够大以压制噪声，又不能大到把窄结构整体模糊掉。ρ 管的是方向统计的空间范围，它决定了“多大范围内的梯度被认为应该方向一致”。

一般经验是 ρ 取 2σ 到 4σ。原因很直接：如果 ρ 太小，参与平均的独立梯度样本太少，方向估计的方差大，方向场会花；如果 ρ 太大，跨过相邻平行结构做平均，两条本来方向不同的纤维会被平均成一个中间方向，交界处的判别能力直接归零。

当你知道目标结构的典型宽度 w（以像素计），ρ 应该控制在 w/2 以内。举个具体例子：一张 1024×1024 的纤维图像，纤维宽度约 12 像素，纤维间距约 6 像素。σ 取 1.0 压制传感器噪声，ρ 取 3.0 可以在约 19 个像素的窗口内统计方向，此时窗口半径约 9 像素，还不到纤维宽度，不会跨结构。如果把 ρ 提到 6，窗口直径 24 像素，稳稳跨过间距，交界处的方向估计就废了。

还有一个容易忽略的点：这两次高斯平滑的边界处理方式。默认的反射填充在图像边缘会制造虚假的对称结构，导致最外圈几行的方向估计完全错误。如果你的目标结构贴着图像边缘，务必用裁剪后再处理的方式，或者干脆把边缘区域的张量置为无效、不参与后续决策。

3. 从张量场到各向异性分割的完整实现

3.1 整体流程与两条技术路线

先把两条路线的骨架摆清楚，再逐段给代码。

路线一（工程上更常用）：读图 → 灰度归一化 → 计算梯度结构张量 → 由张量构造扩散张量 → 各向异性扩散迭代 → 在扩散结果上做分水岭或区域生长 → 后处理（去除小块、孔洞填充）→ 输出标签图。特点是每一步都有成熟实现，调试友好，失败时能快速定位到是扩散没做好还是分割器没调好。

路线二（上限更高）：读图 → 计算张量 → 提取相干性图和方向场 → 构造测地距离或者边界停止函数 → 活动轮廓/水平集演化 → 输出标签图。特点是不需要先做一次可能损伤细节的平滑，直接在原始信息上做决策，但参数敏感性高，收敛性依赖于初始轮廓。

实际项目里我通常先做路线一，因为可靠的基线比炫技重要。等基线跑通、评估指标稳定之后，再用路线二在难点区域（比如粘连严重的地方）做局部增强，两者用同一套评估脚本对比，避免自我感觉良好。

下面按路线一展开，路线二在 3.4 节单独给出核心代码。

3.2 张量计算与方向场可视化的可复现代码

先给张量计算，用 SciPy 就够了，不需要额外的图像处理库。

import numpy as np from scipy.ndimage import gaussian_filter def structure_tensor(img, sigma=1.0, rho=3.0): img = img.astype(np.float64) # 带高斯导数的偏导，等价于先平滑再求导 Ix = gaussian_filter(img, sigma, order=(0, 1), mode='reflect') Iy = gaussian_filter(img, sigma, order=(1, 0), mode='reflect') # 在张量分量上做第二尺度平滑 Jxx = gaussian_filter(Ix * Ix, rho, mode='reflect') Jxy = gaussian_filter(Ix * Iy, rho, mode='reflect') Jyy = gaussian_filter(Iy * Iy, rho, mode='reflect') return Jxx, Jxy, Jyy def tensor_eig(Jxx, Jxy, Jyy, eps=1e-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) # 主特征向量（对应 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）求一阶导，第一个轴不导；order=(1, 0)是沿第一个轴（行方向，y）求导。搞反了不影响 λ 的分布，但会让方向场整体旋转 90 度，属于典型的静默错误。

想把方向场画出来检查，最直观的方式是抽稀后画线段叠加到灰度图上。

import matplotlib.pyplot as plt def plot_orientation(img, theta, step=16, length=8): plt.figure(figsize=(8, 8)) plt.imshow(img, cmap='gray') 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 # 转成边缘走向 dx = np.cos(t) * length / 2 dy = np.sin(t) * length / 2 plt.plot([x - dx, x + dx], [y - dy, y + dy], 'r-', lw=0.8) plt.axis('off') plt.show()

抽稀步长 step 是必须的，逐像素画线会糊成一片红。步长取 12 到 20 之间通常刚好能看清趋势又不遮挡原图。这张图在整个调试过程中价值极高，方向场对了，后面基本就顺了；方向场不对，再怎么调扩散参数都是白费。

3.3 各向异性扩散的离散实现与稳定性约束

扩散张量由张量场的特征向量和特征值共同构造。核心思路是：沿边缘方向（θ + π/2）给大扩散系数 d_edge，沿梯度方向（θ）给小扩散系数 d_grad，后者随 λ1 增大而减小。

def anisotropic_diffusion(img, l1, theta, k=20.0, dt=0.15, iters=60): I = img.astype(np.float64).copy() # 梯度方向抑制系数：λ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, axis=1) - np.roll(I, 1, axis=1)) Iy = 0.5 * (np.roll(I, -1, axis=0) - np.roll(I, 1, axis=0)) fx = D11 * Ix + D12 * Iy fy = D12 * Ix + D22 * Iy div = (0.5 * (np.roll(fx, -1, axis=1) - np.roll(fx, 1, axis=1)) + 0.5 * (np.roll(fy, -1, axis=0) - np.roll(fy, 1, axis=0))) I = I + dt * div return I

这段代码用np.roll意味着周期性边界，对绝大多数图像是不对的，因为图像的左边缘和右边缘毫无关系。生产代码里应该用显式的切片索引实现 Neumann 边界，或者简单地在扩散前对外围 1 像素做镜像扩展，扩散完成后再裁掉。用 roll 的版本跑出来的图，左右边缘会出现莫名其妙的粘连，如果你只关注中间区域可能看不出来，但一旦做全图统计就会暴露。

稳定性是必须算清楚的一件事。显式格式的二维扩散方程，稳定条件近似为 dt ≤ 1/(4·max(D))。本例中 d_edge = 1，所以 max(D) 最大为 1，dt ≤ 0.25。取 0.15 留了约 40% 的余量，比较安全。如果把 d_edge 设成 2 或者更高，dt 必须相应缩小，否则迭代几轮之后数值爆炸，图像出现棋盘状条纹，再迭代下去全是 NaN。

迭代次数不是越多越好。经验值是 40 到 100 之间，具体看图像尺寸和目标结构的尺度。判断标准不是固定次数，而是看连续两次迭代之间的变化量：当 max|I_new − I_old| 小于灰度范围千分之一时就可以停，这时继续迭代只会延长运行时间并轻微模糊边界。可以加一个简单的收敛检查。

关于 k 的取值，它是扩散的“对比度阈值”，物理意义是“λ1 达到多大时认为这是一条需要保护的边缘”。最靠谱的选法是统计 λ1 在整幅图上的分布，取 80% 到 90% 分位数作为初值。如果全图大部分是平坦区域，分位数会偏低，此时扩散会很保守；如果图里到处都是纹理，分位数偏高，扩散会很激进。两种情况都得靠肉眼扫一遍若干张样本的扩散结果来最终确认。

3.4 把张量特征直接喂给分割器

扩散做完之后，分割这一步反而简单了，因为粘连已经被削弱。分水岭是首选，但直接对灰度做分水岭容易过分割，需要可靠的标记。可以用局部极大值做标记，也可以先把扩散结果做一遍形态学开运算再取标记。

from scipy import ndimage as ndi from skimage.feature import peak_local_max from skimage.segmentation import watershed def segment_by_watershed(enhanced, min_distance=12): # 前景掩膜：用 Otsu 在扩散结果上取值 thresh = enhanced > (enhanced.min() + enhanced.max()) / 2.0 dist = ndi.distance_transform_edt(thresh) coords = peak_local_max(dist, min_distance=min_distance, labels=thresh) markers = np.zeros(enhanced.shape, dtype=np.int32) for i, (y, x) in enumerate(coords, start=1): markers[y, x] = i labels = watershed(-dist, markers, mask=thresh) return labels

路线二的做法是用相干性构造边界代价。相干性高说明这里是一条可靠边缘，交叉它应该付出大代价。把这个代价转成测地距离，再做区域生长，就不需要先平滑图像。

def geodesic_cost(coherence, alpha=6.0): # 相干性越高，穿越代价越大 return np.exp(-alpha * np.clip(coherence, 0, 1))

alpha 控制对边界的敏感度。alpha 太大，代价图接近二值，弱边界完全被忽略；alpha 太小，处处都是软边界，区域会漏出去。经验区间是 4 到 10，先用中等值跑一遍看代价图的对比度，再决定往哪边调。

还有一种更彻底的用法：把 (λ1, λ2, C, θ) 四个量做特征标准化后送进 k-means 或者超像素聚类。这条路在纹理类别多、且不同类别方向分布差异明显的场景下很好用，比如把织物图像分成经纱、纬纱、背景三类。缺点是完全丢弃了空间连通性，聚类结果需要再做一次连通域处理，小碎块非常多。我的习惯是先用超像素把空间约束加进去，再在超像素特征上聚类，能省掉一半后处理工作。

3.5 参数选择的计算过程与一次完整标定

把参数标定做成一个固定流程，比每次凭感觉调高效得多。

第一步，估算噪声水平。取图像中一小块明显平坦的区域，算标准差 s。这一步决定了 σ 的下限。

第二步，测量目标结构宽度。在图像上量出典型结构的像素宽度 w，以及相邻结构之间的间距 g。这一步决定了 ρ 的上限。

第三步，按约束定参数。σ 取 max(1.0, 0.5·s) 到 2.0 之间，实测中 1.0 到 1.5 覆盖了大多数情况。ρ 取 2σ 到 3σ，同时强制满足 ρ ≤ min(w, g)/2。

第四步，定 k。跑一遍张量计算，取 λ1 的 85% 分位数。如果这个值明显低于噪声区的 λ1 水平，说明图像对比度低，需要先做一次对比度归一化。

第五步，定迭代次数和 dt。dt 固定 0.15，迭代从 40 起，观察收敛曲线的拐点。

举个具体数字：某批图像量出来 s ≈ 4（8 位灰度），w ≈ 14 像素，g ≈ 7 像素。则 σ 取 1.2，ρ 取 3.0（满足 3.0 ≤ 7/2 = 3.5，卡得比较紧但可行），k 取 λ1 的 85% 分位数约 180（在归一化前的灰度尺度下），dt 0.15，迭代 60 次。这套参数在该批数据上稳定复现，换批次时只需要重新测 s 和 w。

4. 实操中的典型坑与排查技巧实录

4.1 方向场整体偏转 90 度

这是出现频率最高的错误，症状是画出来的方向线全部垂直于真实结构。原因几乎总是把梯度主方向和边缘走向搞混了。θ = 0.5·atan2(2Jxy, Jxx − Jyy) 得到的是对应 λ1 的特征向量角度，而 λ1 是垂直于边缘方向的梯度能量，所以 θ 是法向，边缘走向必须是 θ + π/2。

这个错误在纯扩散的代码里可能不表现出来，因为如果你在构造 D 的时候同时对两个方向都做了相应处理，可能刚好抵消。但只要你在任何一处用了 θ 去表示边缘方向，就会出现问题。防御性做法是在代码里显式命名两个变量：normal_dir和tangent_dir，永远不要用裸的theta传递方向。

4.2 尺度不匹配导致交界处方向被平均

症状是相邻两条平行结构之间的边界模糊，分割结果在这个位置系统性偏移。根本原因是 ρ 太大，方向统计窗口跨过了间距。

诊断方法很直接：在你怀疑出问题的交界处，把 ρ 从小到大扫一遍，看方向场在那个位置是怎么变化的。如果小 ρ 时两条结构的方向清晰可辨，大 ρ 时变成同一个中间角度，就坐实了这个问题。

修法有两种。一是把 ρ 降到 min(w, g)/2 以下；二是改用自适应的 ρ，在 λ1 高的区域用小的 ρ，在平坦区域用大的 ρ。第二种做法实现上稍麻烦，需要两套不同尺度的张量场并做加权混合，但效果确实更好，尤其是在结构尺度变化大的图像上。

4.3 图像边缘的虚假结构

反射边界填充会让图像最外圈出现对称的镜像结构，梯度在这里是虚假的。症状是图像四周出现一圈方向紊乱的区域，且这个紊乱在扩散后会往内部传播几行。

解决方式在实现层面很简单：扩散和张量计算都用mode='nearest'或者干脆用零填充，然后在所有后续统计和可视化里，把最外圈大约 2ρ 宽度的一圈裁掉。如果你的目标结构本身就贴着图像边界，那就只能在采集阶段解决，多留一点余量。

4.4 平坦区的相干性数值噪声

在几乎没有梯度的区域，λ1 和 λ2 都是浮点误差量级的小数，它们的差值同样是误差量级。此时 C = (λ1 − λ2)/(λ1 + λ2 + ε) 的结果完全由 ε 决定，可能呈现出随机的、看起来很像结构的高相干值，如果用它做边界代价，平坦区会莫名其妙地被分割成碎片。

修法是给相干性加一个有效性掩膜：只在 λ1 + λ2 超过某个绝对阈值的位置认为 C 可信，其余位置直接置零或者填一个中性值。这个阈值取噪声区 λ1 的 99% 分位数，用一个平坦区域采样就能估出来。

4.5 常见问题速查表

5. 效果评估、场景适配与性能优化

5.1 用什么指标说明“方向用得对不对”

分割任务常用的 Dice 和 IoU 当然要报，但它们只反映区域重叠，反映不出方向信息有没有真正起作用。我习惯再加两个诊断性指标。

一是边界法向一致性：在参考边界上采样，取每个采样点的切向量，和该位置张量场给出的边缘走向做夹角统计。如果夹角中位数在 10 度以内，说明方向场是准的，分割失败一定出在别的地方；如果夹角分布很散，方向场本身就有问题，调分割器没意义。

二是区域内的平均相干性。真实的结构区域相干性应该显著高于背景，如果某个分割出来的区域平均相干性接近背景水平，那这个区域大概率是误检，即使它的 Dice 贡献看起来还行。

这两个指标的好处是能定位问题层次。Dice 差，可能是方向场错，也可能是分割器参数错；但加上这两个诊断量，就能立刻区分。

5.2 不同场景下的参数倾向

不同应用领域对参数的偏好差别很大，强行用一套参数跨领域是不现实的。

这张表不是标准答案，是起点。我建议的做法是先用表格里的一行跑通，再针对自己数据里最差的那 10% 样本单独调，把调好的参数写进配置文件按数据批次管理，不要硬编码在脚本里。

5.3 性能与工程落地

高斯滤波在 SciPy 里是可分离的，复杂度是 O(N·r) 而不是 O(N·r²)，但 ρ 大的时候仍然慢。如果 ρ 超过 8，用盒式滤波近似能快好几倍，代价是方向场的平滑性稍差，做方向场可视化时能看出块状。批量处理时，角度和相干性的计算都是逐像素的算术运算，可以用 NumPy 向量化一次完成，不要写循环。

扩散迭代是性能瓶颈。60 次迭代、每次若干次数组操作，对 1024×1024 的图像大概在秒级，可以接受；对 4K 图像就明显吃力了。几个提速方向：一是把扩散放到金字塔的粗层先跑，再上采样作为细层初始值，总迭代次数能省一半以上；二是扩散张量 D 在整个迭代过程中通常是固定的，可以预先算好并缓存，不必每轮重算；三是如果必须实时，转向基于 CUDA 的实现或者把张量计算和扩散都写成批处理的张量运算。

还有一个工程细节值得提前设计：把整个流程封装成“参数对象 + 处理函数”的形式，参数对象序列化成 JSON 和数据版本绑定。这套东西在单次实验里显得多余，但当你需要回头复现三个月前某个结果的时候，会发现省下的时间远超当初写它的时间。

最后分享一个我在实际使用中养成的小习惯。每次调参之前，先挑三张有代表性的样本：一张结构清晰、一张噪声重、一张粘连严重。三张一起跑，看方向场叠加图和相干性直方图，只根据这三张的反馈调参数，绝不在全量数据上反复试。全量跑一轮的反馈周期太长，你会失去对参数敏感度的直觉。三张样本的快速循环跑上二三十轮，对这套方法的脾性基本就摸透了，之后再上新数据，往往第一组参数就能到可用的水平。