搞图像处理的人早晚都会撞上距离变换Distance Transform这个词。我第一次觉得它要命是在做分割后处理的时候一张 512×512 的图片前景像素占一半暴力算欧式距离变换屏幕上光标转圈一杯咖啡喝完还没出结果。后来换成 Felzenszwalb 算法同样一张图几十毫秒就出来了效果还完全一致那一刻我意识到这个算法演进背后的数学值得好好讲一讲。本文用最快的方式梳理欧式距离变换EDT的前世今生从最直白的暴力算法开始到一维分解的思路再到 Felzenszwalb 的线性复杂度算法。如果你用过cv2.distanceTransform或scipy.ndimage.distance_transform_edt但不知道里面发生了什么这篇文章适合你如果你需要自己实现一个高效 EDT这篇文章更是可以直接抄作业。1. 先搞清楚 EDT 是什么1.1 二值图与距离图的关系欧式距离变换做的事情一句话就能说清给一张二值图白色前景是目标区域黑色背景是其他部分EDT 输出一张灰度图每个像素点的灰度值等于它到最近前景像素的欧氏距离前景本身距离为 0。这时候很多人会问这不就是查表或者扫一遍的事吗还真不是。难点在于每个背景像素都必须找到最近的“那个”前景像素而前景像素可能非常多像素之间的空间关系又是连续的。反映在输出上距离图自然形成一片连续过渡的灰度场边缘处距离小远离前景的地方距离大最后在形状的中心会出现一条“脊线”——这条脊线恰好是形状的中轴这也是后面要讲到的骨架提取的基础。1.2 为什么只用平方距离不开根号也能算严格来说欧氏距离应该带根号d sqrt(dx² dy²)。但算法里几乎全程只维护 dx² dy²直到最后输出才开一次根号。原因是平方距离保留了所有大小关系a² b² 等价于 a b找最近像素只需要比较大小不需要真实距离。省掉 sqrt 就是省掉大量 CPU 周期尤其在几百万像素的图上差别是秒级的。另外用平方距离还有一个更深的好处能把二维问题拆成两个一维问题。因为 (dx² dy²) 恰好是横向平方距离和纵向平方距离之和这种可分离性在后面是核心突破口。1.3 距离变换不是冷门玩具它很有用我接触 EDT 主要是三块。第一是骨架提取/中轴变换距离图的局部极大值组成形状的中轴线常用于笔迹识别、血管分析第二是语义分割的后处理把网络输出的 mask 转成距离场再做分水岭就能把贴在一起的实例分开第三是机器人路径规划里的代价地图离障碍物越远代价越小导航路径自然偏向安全区域。除此之外字体渲染、点云配准、医学图像测量都会用到它。理解了 EDT很多看似无关的模块其实是同一个数学内核。2. 暴力算法原理直白但慢到怀疑人生2.1 思路和代码一样直白暴力算法的思路别想复杂了遍历每个背景像素遍历每个前景像素算距离取最小。伪代码 10 行正确性一眼能看出来。import numpy as np def edt_bruteforce(binary): h, w binary.shape fg np.argwhere(binary 0) # 所有前景坐标形状 (M, 2) out np.zeros((h, w), dtypenp.float64) if fg.shape[0] 0: return out for y in range(h): for x in range(w): if binary[y, x] 0: continue dy fg[:, 0] - y dx fg[:, 1] - x out[y, x] np.min(dy * dy dx * dx) return np.sqrt(out)注意我全程用平方距离只在最后开根号。这段代码在 64×64 小图上调试很舒服但当你想认真跑一张真实图片时它会教你做人。2.2 复杂度从还算能跑到彻底跑不动设图像像素总数为 N H×W前景像素数为 M。暴力算法的复杂度是 O(N·M)。听起来还好但你代入数字感受一下512×512 的图像N ≈ 26 万前景如果有 10 万相乘是 260 亿次距离计算每次还要做浮点乘法和开方。用 Python 跑大概要几分钟到几十分钟即使用 C 裸写也是秒级往上走。生产环境里 512×512 算一张都嫌慢更别说视频流或者医学影像那种千万像素的大图。暴力算法最大的问题不是图画得大而是前景越密越慢。图像内容越多M 越大耗时越差完全不可扩展。所以后来出现的优化算法核心目标就是让计算量和前景数量解耦最好只和总像素数 N 成线性关系。2.3 暴力版的价值当基准别急着扔掉暴力版。我做算法移植时有一个习惯永远留一份暴力实现当“标尺”。Felzenszwalb 算法再花哨最终结果也必须和暴力版在数值上对上。后面第 5 章我会给对拍代码你先记住暴力版不是摆设它是验证一切优化的试金石。3. 先降维一维分解让问题豁然开朗3.1 关键观察二维欧氏距离可以分两步到这里如果只看暴力算法会觉得这个问题无解。换个角度欧氏距离平方 (x - x)² (y - y)² 是两个独立方向之和。于是可以这样做第一步对每一行做一维变换行内的每个像素记录它到该行最近前景像素的水平距离平方。前景是 0背景是一个正的平方数。这一步相当于把二维问题暂时压成“行内的水平最近距离”。第二步对每一列做一维变换列上每个像素拿到上面那一步的结果再加上纵向的 (y - y)²再取最小。合并之后每个点得到的就是完整的欧氏距离平方。为什么这样拆开是精确的因为全局最近前景像素 (x*, y*) 一定落在某一行的某个位置。先做行变换时该行每个 x 已经标注了到 x*或同行更近前景的水平距离再做列变换时把纵向距离加上就能在最小化过程中取到真正的全局最优同时也不会产生比全局最优更小的结果。这个“分两步但结果不变”的结论是很多精确欧氏距离算法的数学基石。3.2 一维问题的数学表达式于是核心问题变成给定长度为 n 的一维数组 f其中部分是有效值0 或者上一轮的结果要求输出d[i] min_{j 有效} ( f[j] (i - j)² )这是一个标准的一维二次距离变换。现在暴力算法在一维上就已经够简单了遍历 i 和 j取最小复杂度 O(n²)。但别急着满足这个一维形式有一个了不起的几何解释——把每个有效点 j 想象成一条抛物线。3.3 演进岔路倒角、波前和线性精确算法讲 Felzenszwalb 之前先花两分钟看看这条演进路上出现过哪些方案避免你误以为只有暴力到线性这两条路。时期方法复杂度精确性代表1968数字距离变换O(N)非欧氏Rosenfeld倒角模板1980序贯扫描近似 EDTO(N)近似Danielsson2000分段线性精确 EDTO(N)精确Meijster 等2004/2012抛物线包络线性 EDTO(N)精确Felzenszwalb Huttenlocher最早实用的是倒角距离变换Chamfer Distance Transform它不追求欧氏距离而是用局部模板两次扫描近似比如 3×3 模板跑一遍前向、一遍后向就出结果O(N) 且实现极简单但结果是曼哈顿/切比雪夫的混合方向上有明显菱形痕迹后来做精确测量就不够用了。后来又出现基于波前传播的快速行进法Fast Marching Method复杂度 O(N log N)精度好一些但实现复杂且本质上是数值求解不是严格解析解。真正的转折点是 Meijster 等人在 2000 年左右提出线性时间精确 EDT 算法Felzenszwalb 和 Huttenlocher 又把这个思路用非常干净的抛物线包络方法讲清楚了配合两行关键推导代码量砍到 40 行以内。现在 OpenCV 和 scipy 里的精确欧氏距离变换底层基本都是这类线性算法。4. Felzenszwalb 算法把距离变换变成求抛物线下包络4.1 把每个点看成一条抛物线一维问题的每个有效源点 j 其实定义了一条抛物线P_j(x) (x - j)² f[j]而 d[i] 就是所有抛物线在 x i 处的最小值。整个一维距离变换可以理解为在一排无穷延伸的抛物线里找到每个位置最低的那个点。这个最低曲线叫下包络lower envelope。关键在于两条抛物线最多相交一次因为它们都是凸二次函数差值是一次函数。这意味着下包络一定可以被分成若干区间每个区间由一条抛物线占据而且这些抛物线的顺序和位置天然有序——不会出现抛物线 A 统治一段、B 统治一段、然后 A 又回来统治一段的情况。于是问题转化为如何高效地构造这个下包络答案是维护一个栈从左到右扫描。4.2 用栈维护下包络Felzenszwalb 的做法从左到右扫描所有有效位置维护一个栈栈里存放的是当前可能出现在下包络里的抛物线编号。每来一条新抛物线 q就计算它和栈顶抛物线的交点 s。如果 s 比栈中记录的区间边界还靠左说明栈顶抛物线在成为最小点之前就被新抛物线压过去它永远没有机会出场弹出。这个过程一直持续到交点落在正确位置然后把新抛物线压栈。这里面的数学就是解一个一次方程。两条抛物线 P_i(x) 和 P_j(x) 的交点 s 满足(s - i)² f[i] (s - j)² f[j]展开后二次项抵消剩下s ( f[j] j² - f[i] - i² ) / (2j - 2i)这个公式是整个算法的心脏代码里就那么一行但如果你没看懂它后面所有调试都会抓瞎。我自己初学时在这卡了一下午后来用手画了两条抛物线的交点和下包络才真正明白所谓弹出栈顶就是发现旧抛物线的“统治区间”已经被新抛物线完全吞并留着只会碍事。4.3 完整可运行的代码我贴一份我生产中用过的实现改成通用版本二维任意长宽都能直接调用。import numpy as np def dt_1d(f): 一维二次距离变换。 输入 f长度 n 的数组无效点用 np.inf 表示。 输出 dd[i] min_{j 有效} (f[j] (i - j)^2) n f.shape[0] positions np.flatnonzero(np.isfinite(f)) m positions.size if m 0: return np.full(n, np.inf) v np.zeros(m, dtypenp.int64) # 栈中抛物线的下标源位置 z np.zeros(m 1, dtypenp.float64) # 相邻抛物线的交点 k 0 v[0] positions[0] z[0] -np.inf z[1] np.inf for idx in range(1, m): q positions[idx] # 当前栈顶抛物线 v[k] s ((f[q] q * q) - (f[v[k]] v[k] * v[k])) / (2.0 * q - 2.0 * v[k]) while s z[k]: k - 1 s ((f[q] q * q) - (f[v[k]] v[k] * v[k])) / (2.0 * q - 2.0 * v[k]) k 1 v[k] q z[k] s z[k 1] np.inf # 回填结果 d np.empty(n, dtypenp.float64) k 0 for x in range(n): while z[k 1] x: k 1 d[x] (x - v[k]) ** 2 f[v[k]] return d def edt_2d(binary): 二维欧式距离变换精确。 输入 binary二维 0/1 数组1 表示前景。 输出与输入同形状的浮点数组每个元素是到最近前景的欧氏距离。 h, w binary.shape f np.where(binary.astype(bool), 0.0, np.inf).astype(np.float64) # 先逐行一维变换 for y in range(h): f[y, :] dt_1d(f[y, :]) # 再逐列一维变换 for x in range(w): f[:, x] dt_1d(f[:, x]) return np.sqrt(f)这里有三个变量要特别说明。v 数组存的是“下包络当前第 k 段由哪个源点主导”z[k] 存的是 v[k-1] 和 v[k] 的交点也就是第 k 段区间的左边界z[k1] 是右边界。初始化时 z[0] -inf、z[1] inf保证第一段区间覆盖整个数轴。扫描完所有源点后下包络就固定了最后回填时用一个指针 k 跟着 x 走在每个区间内直接取对应 v[k] 算距离这段是 O(n)。4.4 复杂度与二维扩展每个源点最多入栈一次、出栈一次回填每个 x 也最多移动 k 一次所以一维复杂度是 O(n)。二维要做 h w 次一维变换总复杂度 O(h·w) O(N)。这个线性复杂度跟前景数量完全无关是暴力算法无法企及的。我还想强调二维精确的由来行变换能给出水平方向的最小平方距离列变换再把纵向距离叠加。因为每个二维候选点 (x*, y*) 都被唯一地对应到一行和一列两次最小化依次作用后总能取到全局最优。这个过程在数学上是有严格保证的不需要任何近似。换句话说Felzenszwalb 算法不是“算得挺快但差不多”它是和暴力版完全等价的高效实现这也是它能被大量工程库采纳的根本原因。5. 实操经验验证、踩坑与性能优化5.1 如何验证你没有写错我刚写完这个算法的时候心里其实没底毕竟栈和交点的逻辑第一眼不够直觉。验证方法很简单拿随机二值图和暴力版对照。np.random.seed(0) binary (np.random.rand(64, 64) 0.7).astype(np.uint8) ref edt_bruteforce(binary) fast edt_2d(binary) print(np.max(np.abs(ref - fast)))在我机器上这个值一般在 1e-12 量级属于浮点误差。如果出现明显差异优先检查输入是不是被当成 0/255 了以及是否忘记把前景/背景反转。测试图不要只用一张至少覆盖全前景应该全 0、全背景应该全 inf、单像素前景、两个分离前景等边界情况。我在实际对拍时还习惯把测试图放大到 128×128 以上因为太小的图栈结构过于简单很多边界逻辑根本走不到。5.2 常见问题速查表症状可能原因解决方法输出比真实距离小前景/背景反了输入被当成 0/255 导致 255 也算前景判断改成binary 255统一用 0/1 输入结果全是 inf全背景没有任何源点提前判空返回全 0结果有 NaNinf 参与计算导致 inf - inf无效点用 np.inf过滤时只用np.isfinite栈越界z[k1] 没有初始化为 inf每次入栈后设置z[k1] np.inf暴力版和快速版不一致浮点累计误差容忍 1e-8 量级差异用np.allclose对比还有一条非常实用的建议不要在中间步骤开根号。整个流程中始终使用平方距离只有最后输出时才 sqrt。很多人一上来就把距离转成真实值后面再比较数值精度和性能都会下降。5.3 工程上的选型建议如果你只是想拿结果不搞研究我建议直接用成熟库scipy.ndimage.distance_transform_edt内部就是精确欧氏距离变换OpenCV 的cv2.distanceTransform搭配DIST_L2和DIST_MASK_PRECISE也是精确欧氏的。自己实现这个算法的意义在于一是理解原理遇到结果不对能判断是库的用法问题还是算法适配问题二是特殊场景比如要在 GPU 上写自定义核、要处理视频流需要把行/列变换拆到不同线程、需要拿到中间的行变换结果做别的用途——这些只有掌握原理才能灵活改。实测上Python 版 Felzenszwalb 在 512×512 随机图上大概几十毫秒暴力版跑不出来纯 C 的话同样的图可以做到几毫秒。如果行变换并行化还能再压一截。我在实际工程里习惯把行列变换都写到循环里配合 OpenMP 或 TBB 对行做并行瓶颈基本就只剩下内存带宽。最后分享一个个人体会。我最初理解 Felzenszwalb 算法时老想着“为什么不用排序/四叉树/BFS”后来发现都是走弯路。这个算法最妙的地方是把一个离散优化问题翻译成了连续几何问题每个源点变成一条抛物线问题变成求下包络。一旦接受了这个视角代码里每一个数组、每一个 while 循环都变得顺理成章。弄懂这套之后我再遇到类似的最近邻居、下包络、单调最优化问题第一反应都是想想能不能用同样的模型去套。距离变换不是终点它是一把很好用的钥匙。如果你准备动手实现建议按这个顺序来先跑暴力版再跑一维版最后再拼成二维版。每层都留一个能对拍验证的函数踩坑会少很多。