企业官网建设流程全解析

1. EDT解决什么问题：从"最近背景有多远"到距离场

先花三十秒想一个场景：你做图像分割，模型输出了一张二值图，里面是目标区域，外面是背景。接下来你要给每个目标像素算一个值——它离最近的背景像素有多远。这个值不是随便算着玩的，它直接决定了后续的骨架提取、形态学操作、边缘光滑度评估、路径规划中的障碍物距离计算。这个操作就叫欧式距离变换（Euclidean Distance Transform，简称EDT）。

距离变换的思想并不复杂：把二值图当成一个"标记场"，目标像素记为1，背景像素记为0，然后对每个1像素，找到离它最近的0像素，记录下这段欧氏距离。把所有结果拼在一起，就得到一张灰度图，背景是黑（0），离背景越远的前景像素越亮。这张灰度图就是距离场，它把"位置关系"翻译成了"数值关系"，后面所有算法都能在这个数值场上做文章。

我在早期做医学图像里的血管中心线提取时，最依赖的就是这个距离场。血管在分割结果里是一条条粗细不一的连通区域，直接做形态学细化容易产生毛刺，但先对分割图做EDT，中心线附近的像素一定是局部距离极大值，沿着极大值连起来就是血管骨架。同样的思路可以用在任何"找中心"的任务上，比如道路中心线、纸上笔迹的骨架、零件的中轴。可以说，EDT是图像分析里一个不起眼但底座级的基础操作，理解了它，后面看SDF（符号距离场）、看形状匹配、看路径规划里的势场法都会顺很多。

这篇文章的路线很直接：先从最笨的暴力算法讲起，然后看中间那些"便宜但不够精确"的折中方案，最后深入Felzenszwalb在2012年发表的线性时间精确算法。看完你会发现，这个所谓的高端算法，核心其实是一堆抛物线求下包络的几何问题，代码写出来不到四十行。三分钟读完不敢保证，但十分钟内理清全貌是没问题的。

2. 暴力算法：定义最直观，代价也最直观

2.1 暴力版的实现与算账

如果要给EDT写一个最符合人类直觉的实现，那一定是这样：先把所有背景像素的坐标收集到一个列表里，然后遍历每个前景像素，计算它到列表中每个背景像素的欧氏距离，取最小值。伪代码长这样：

import numpy as np from scipy.spatial.distance import cdist def brute_force_edt(binary_img): # binary_img: ndarray, 1为前景, 0为背景 H, W = binary_img.shape ys, xs = np.where(binary_img == 0) # 背景坐标 backgrounds = np.stack([ys, xs], axis=1) result = np.zeros((H, W), dtype=np.float64) for y in range(H): for x in range(W): if binary_img[y, x] == 0: continue dist2 = np.min(np.sum((backgrounds - np.array([y, x])) ** 2, axis=1)) result[y, x] = np.sqrt(dist2) return result

这个实现很容易理解，但你要真的拿一张512×512的图像去跑，会等到怀疑人生。问题不在代码写法上，而在算法本身的算账方式：假设图像里有一半像素是背景，一半是前景，那么512×512图像里约有13万个背景像素和13万个前景像素。每个前景像素都要和所有背景像素算一次距离，总计算量是13万乘以13万，约169亿次欧氏距离计算。这个量级在C++里勉强能忍，在Python里基本属于"不可运行"级别。

复杂度的完整写法是：设图像总像素数为N=H×W，背景比例q，前景比例p=1-q，则暴力算法需要计算p·q·N²次距离。就算p和q都是0.5，也是0.25N²；如果背景和前景各占一半，N=10万像素的图就是25亿次计算。注意这里的N是像素总数，图像尺寸翻一倍，N翻四倍，计算量直接翻十六倍。这就是为什么暴力算法只能用在"演示原理"的层面，工程上没人这么干。

2.2 复杂度到底炸在哪：O(N²)的真相

很多人第一次看到复杂度分析会问：不就是对每个像素找一个最近点吗，为什么是O(N²)而不是O(N)？关键在于"找最近点"这个动作没有免费的索引结构。暴力做法里，每个前景像素的点查询成本是O(背景像素数)，而前景像素又有O(N)个，所以总成本是O(N·背景数)。如果背景数和前景数同量级，就是O(N²)。也就是说，复杂度不是炸在"有多少像素要做变换"，而是炸在"每个像素都要问一圈所有背景像素"。

那能不能用空间索引把点查询加速到O(log N)甚至O(1)？比如KD-Tree、四叉树、网格哈希？确实可以，但这里有个隐藏问题：EDT要求的不是"少数几个查询点"的最近邻，而是"图像上每一个前景像素"的最近邻。全图所有像素都是查询点，构建索引的代价和查询总代价加起来，往往还是O(N log N)或O(N√N)级别，并且常数不小。更麻烦的是，距离变换的输出是一个密集的距离场，不是稀疏的点对关系，任何基于"点查询"的思路都绕不开"每个像素都要有个结果"这个事实，所以更聪明的方向不是优化查询，而是换一种整体计算的思路——这就是后面所有优化算法的出发点。

在进入中间方案之前，我还想澄清一个常见误解：有些人以为只要把暴力代码里的np.min换成argmin，或者用cdist批量算，速度就上去了。实验确实会快一点，因为cdist底层调了BLAS，但复杂度没变，图像一大照样撑不住。我之前在实习项目里就犯过这个错，用cdist处理一张2000×2000的图，内存直接爆掉，因为cdist要一次性生成一个前景×背景的距离矩阵，2000×2000图如果前后景各占一半，矩阵大小就是200万×200万，光存float64就要32TB，想都不用想。所以暴力算法只能作为正确性参照，不能当产品实现。

3. 中间地带的优化尝试：便宜但不够精确

3.1 倒角算法：一前一后两趟扫描

既然暴力算法用不起，早期研究者想到的第一个优化思路是"局部传播"。一个直观观察是，一个像素到背景的距离，和它邻居到背景的距离之间有关系：如果我知道右边邻居的距离是d，那当前像素的距离大约就是d+1（欧氏距离水平方向）。沿着这个思路，可以用两次扫描（一次从左上到右下，一次从右下到左上）来传播距离信息，每次只查一个小邻域里几个预设方向的邻居。这个算法叫倒角距离变换（Chamfer Distance Transform）。

具体做法是设计一个3×3或5×5的掩码，掩码里每个位置有一个权重。比如经典3×3掩码：

• 水平/垂直邻居权重：1
• 对角邻居权重：√2（约1.414）

第一遍正向扫描时，对每个像素，用掩码左上半区的邻居候选值+对应权重来更新自己；第二遍反向扫描时，用右下半区的邻居候选值+对应权重继续更新。这样每趟扫描只需常数次操作，总复杂度O(N)。听着很美，但代价是：这是近似距离，不是精确的欧氏距离。3×3掩码的最大误差约8%，5×5掩码能压到3%左右，但算法复杂度和内存占用也上去了。

为什么会有误差？因为距离传播本质上是用一个有限方向集合上的最短路径去逼近欧氏距离。想象从前景像素走到背景像素，你只允许走水平、垂直、对角这三种方向，那走出来的路径长度是曼哈顿距离和棋盘距离的混合，而不是直线距离。8%的误差在骨架提取这种对距离精度敏感的任务里很容易造成中心线偏移，所以倒角算法后来基本被精确算法取代，只在一些嵌入式场景里还能见到。

3.2 BFS扩散：距离成了网格脚印

另一种直觉思路是把背景像素当作"源点"，让它们像水波一样往外扩散，每扩散一圈，距离加1。这个思路实现起来就是BFS（广度优先搜索），从所有背景像素出发做多源BFS，访问到每个前景像素时记录步数。复杂度O(N)，实现也简单，但问题同样出在距离定义上：BFS的步数对应的是网格上的最短路径跳数。如果你用四邻域，得到的是曼哈顿距离；用八邻域，得到的是棋盘距离的变体（对角步算一步）。

这两种距离都不是欧氏距离。四邻域BFS在45度方向上的误差高达约41%，八邻域BFS在对角方向误差约8%到15%不等。有人可能会说，既然8邻域误差不大，是不是够用了？在很多早期游戏寻路里确实够用，但在需要精确几何语义的任务里不行。比如你要根据距离场生成模型的等距面，BFS出来的距离场会产生明显的方块状伪影，曲面根本不平滑。SDF（有向距离场）如果这么算，渲染出来的描边效果会直接穿帮。

BFS真正的价值在于它建立了"距离场可以靠传播计算"的心智模型，也启发了一类基于队列的改进算法——通过在队列中加入误差校正项或者用多趟扫描，可以逼近欧氏距离。但这些改进要么引入复杂的优先队列操作（退化成Dijkstra，复杂度O(N log N)），要么仍然存在方向性误差，始终没有一个"既简单又精确"的答案。

3.3 分块/局部窗口：精度与速度的拉锯

还有一个常见的工程折中是分块暴力：把图像切成小块，每个前景像素只需要在小块内部以及小块边界附近的背景像素里找最近点。如果块大小是w×w，复杂度大约是O(N·w²)。这意味着你可以用一个很小的窗口（比如7×7）跑得飞快，但结果只在窗口范围内局部最优——如果一个前景像素离背景的最近距离超过窗口半径，算出来的就是截断值，不是真实距离。

有人会想，那把w设大一点不就行了？比如32×32窗口，查一次要1024次距离计算，比暴力好很多，但仍然达不到O(N)。而且窗口大了之后，块和块之间的接缝处理、重叠区计算都变麻烦，实际加速比远低于理论值。我在做三维体数据的距离变换时试过这条路，内存占用和代码复杂度都上来了，精度还不能保证，最后果断放弃。所以分块方案更适合"只需要局部小范围距离"的场景，比如碰撞检测里只关心物体表面附近几毫米的距离场，而不是全图精确距离。

这些中间方案告诉我们一件事：精度和速度的拉锯，缺一个数学层面的突破。近似算法速度快但误差不可控，暴力算法精确但复杂度不可接受。我们需要的是"精确且线性"的算法，而Felzenszwalb给出的正是这个。

4. Felzenszwalb的数学内核：下包络与抛物线栈

4.1 把距离平方展开成抛物线

Felzenszwalb和Huttenlocher在论文《Distance Transforms of Sampled Functions》里讨论的问题比二值EDT更一般：给定一个采样函数f(i)，对于每个位置x，要求：

d(x) = min_i [ (x - i)² + f(i) ]

这个式子的含义是：x处的新值等于"某个源位置i的本身值f(i)，加上从i到x的欧氏距离平方"。当f(i)只在背景处为0、在前景处为无穷大时，这个式子就退化成了精确的欧氏距离变换。所以这个一般形式才是算法的真正内核。

为什么写成(x-i)²？因为欧氏距离平方展开后天然是二次的。如果把每个i对应的函数画出来，g_i(x) = (x-i)² + f(i)，这是一条开口向上的抛物线，顶点在x=i处，最低点的值是f(i)。我们需要在每一个x处，从所有这样的抛物线里挑出高度最低的那条。也就是说，d(x) = min(g_1(x), g_2(x), ..., g_n(x))，即这一族抛物线的下包络。

这里有个几何直觉：一堆开口方向和宽度都相同的抛物线，只是顶点位置和最低点高度不同。它们的下包络就像给这个抛物线家族"蒙了一层布"，这层布在每个x处贴在最低的那条上。因为抛物线是光滑凸函数，这层下包络也是一条分段抛物线的组合，而且每条抛物线在下包络上出现的区间是连续的。这个"区间连续"的性质，是后面能用栈实现的关键。

4.2 栈如何在线性时间内维护下包络

现在问题是：怎么高效地构造这组抛物线的下包络？最朴素的想法是枚举每个x，把n条抛物线都算一遍取最小值，复杂度O(n²)。但Felzenszwalb观察到，由于所有抛物线形状相同，只是平移和上下移动，它们的交点结构有非常强的规律性。

想象一条新抛物线q要加入已有的下包络。它不会从中间插进去，只会影响下包络最右侧的那一段。为什么？因为抛物线开口都朝上，且二次项系数相同，两条抛物线的交点的x坐标是唯一的；三条及以上抛物线的下包络，从左到右是一段接一段的，每条抛物线最多出现在一个连续区间里（这个性质叫"下包络的区间单调性"）。于是可以用一个栈来维护当前下包络的"分段信息"：

• 栈里保存的是每条抛物线对应的源位置下标，以及它在x轴上的起始区间左端点。
• 新抛物线q入栈时，和栈顶抛物线算交点横坐标s。如果s落在栈顶抛物线区间的左端点左侧，说明栈顶抛物线在下包络上完全没有露脸机会，直接弹出；否则q的起始区间就从s开始，入栈。
• 由于每条抛物线最多入栈一次、出栈一次，整个过程是O(n)的。

用生活化的类比：想象一队在操场上排队做操的人，每个人都举着一根竖杆，杆顶有一条抛物线形的柔性线。新来一个人要把自己的抛物线线放到最右边，如果发现它和队尾那位的交点在队尾那位已经消失的位置更左边，说明队尾那位从头到尾被新来的完全压住，可以让他直接下场；否则新来的人就排到队尾，从交点处开始接管。整个过程，每个人最多被叫出来一次，所以效率很高。

这个栈维护的是"抛物线源下标序列v"和"每条抛物线的区间起点数组z"。查询阶段，给一个x，从栈里找到区间覆盖x的那条抛物线，计算(x-v[k])²+f(v[k])即为d(x)。由于x是从左往右递增扫描的，区间指针也只会往右走，查询同样是O(n)总耗时。

4.3 一个具体的数值例子

空讲公式容易晕，我拿一个极小的例子走一遍。假设有一列4个像素，经过行方向变换后，某列的值是g = [0, 1, 4, 0]，意思是第0行和第3行是背景（距离为0），第1行和第2行是前景，它们各自的水平最近背景距离平方分别是1和4。现在要在列方向求出最终的二维距离平方。

四条抛物线分别是：

• i=0：g_0(x) = x² + 0
• i=1：g_1(x) = (x-1)² + 1
• i=2：g_2(x) = (x-2)² + 4
• i=3：g_3(x) = (x-3)² + 0

逐步入栈：

1. 抛物线0入栈，区间从负无穷开始。
2. 抛物线1入栈，和抛物线0的交点：解x² = (x-1)²+1，得x=1。交点在抛物线0区间起点右侧，所以0的区间是(-∞, 1]，1的区间从1开始。
3. 抛物线2入栈，和抛物线1的交点：解(x-1)²+1 = (x-2)²+4，得2x = 2? 展开：(x²-2x+1)+1 = (x²-4x+4)+4 → -2x+2 = -4x+8 → 2x=6 → x=3。交点在3，而抛物线1的区间起点是1，所以1的区间变成[1, 3]，抛物线2区间从3开始。
4. 抛物线3入栈，和抛物线2的交点：(x-2)²+4 = (x-3)²+0 → -4x+8 = -6x+9 → 2x=1 → x=0.5。交点0.5小于抛物线2的区间起点3，说明抛物线2在下包络里完全没有露脸机会，弹出抛物线2。继续和现在的栈顶抛物线1求交点：(x-1)²+1 = (x-3)²+0 → -2x+2 = -6x+9 → 4x=7 → x=1.75。这个交点在抛物线1的区间起点1右侧，所以抛物线1的区间更新为[1, 1.75]，抛物线3的区间从1.75开始。

最终下包络分段：x∈(-∞,1]用抛物线0，x∈[1,1.75]用抛物线1，x∈[1.75,+∞)用抛物线3。代入x=0、1、2、3：

• x=0：抛物线0，0
• x=1：抛物线0或1，都是1
• x=2：抛物线3，(2-3)²+0=1
• x=3：抛物线3，0

所以这一列的最终距离平方是[0, 1, 1, 0]。第1行和第2行的像素距离背景的最近距离都是1，这个结果从几何上看也对：第1行向上一步就是背景行，第2行向下一步是背景行，水平方向反而更远。整个过程，栈操作次数是线性的，没有任何逐像素枚举所有源的浪费。

5. 完整实现：两步法把你自己的EDT跑起来

5.1 第一步：行方向的一维变换

二维EDT的分解思路很经典：先对每一行做一维距离变换，得到每个像素到"同一行里最近背景像素"的水平距离平方；然后对每一列，把上一步得到的值当成新的源函数，再做一次一维距离变换，合成二维距离。第一步相对简单，因为二值图在行方向上要么是背景（值0），要么是前景（需要初始化为一个很大的数），可以直接用一趟从左到右、一趟从右到左的扫描完成。

import numpy as np INF = 1e20 def edt_1d_binary(row): """一维二值距离变换：row里0是背景，返回每个位置到最近背景的平方距离""" n = len(row) d = np.full(n, INF, dtype=np.float64) # 从左到右 dist = INF for x in range(n): if row[x] == 0: dist = 0 elif dist < INF: dist += 1 if dist < INF: d[x] = dist * dist # 从右到左 dist = INF for x in range(n - 1, -1, -1): if row[x] == 0: dist = 0 elif dist < INF: dist += 1 if dist < INF and dist * dist < d[x]: d[x] = dist * dist return d

这段代码里有个容易忽略的细节：对于整行都是前景的情况，dist会一直是INF，距离保持无穷大。实际二值图很少出现整行全前景的情况，但代码要能正确处理，不能变成INF+1导致溢出。这也是为什么我用一个INF=1e20而不是float('inf')，因为后续列方向的抛物线交点计算里，如果出现无穷参与运算，结果会变成NaN，非常难排查。用一个足够大的有限数，可以保证数值稳定。

5.2 第二步：列方向的下包络组合

第二步是Felzenszwalb算法的核心。对每一列，我们把上一步得到的行方向距离平方数组（长度H）当作f，执行一次"任意采样函数的一维距离变换"，也就是第四节里那个抛物线栈算法。这里的f已经不再是0/INF二值，而是各种非负数值，所以不能再用简单的两趟扫描，必须用抛物线栈。

def dt_1d_general(f): """对任意一维数组f做距离变换: d[x] = min_i[(x-i)^2 + f[i]]""" n = len(f) v = [0] * n # 下包络中抛物线对应的源位置下标 z = [0.0] * (n + 1) # 每条抛物线区间的左端点 k = 0 v[0] = 0 z[0] = -INF z[1] = INF def intersect(i, j): # 抛物线i和j的交点横坐标，(x-i)^2+f[i] = (x-j)^2+f[j] # 展开后: 2(i-j)x = f[j]+j^2 - f[i]-i^2 return ((f[i] + i * i) - (f[j] + j * j)) / (2.0 * (i - j)) for q in range(1, n): while True: s = intersect(v[k], q) if s <= z[k]: # 新抛物线把栈顶抛物线完全压住，栈顶出局 k -= 1 else: break k += 1 v[k] = q z[k] = s z[k + 1] = INF # 查询阶段 d = np.empty(n, dtype=np.float64) k = 0 for x in range(n): while z[k + 1] < x: k += 1 diff = x - v[k] d[x] = diff * diff + f[v[k]] return d def felzenszwalb_edt(binary_img): H, W = binary_img.shape g = np.zeros((H, W), dtype=np.float64) # 第一步：行方向二值EDT for r in range(H): g[r, :] = edt_1d_binary(binary_img[r, :].astype(np.float64)) # 第二步：列方向一般EDT d_sq = np.zeros((H, W), dtype=np.float64) for c in range(W): d_sq[:, c] = dt_1d_general(g[:, c].copy()) return np.sqrt(d_sq)

这段代码里最绕的是intersect函数。从公式展开看，(x-i)²+f[i] = x² - 2ix + i² + f[i]，两条抛物线相减后x²消掉，剩线性方程，解出来的交点横坐标就是上面那个((f[i]+i²)-(f[j]+j²))/(2(i-j))。注意分母的符号：因为v[k]总是小于q（源位置是递增入栈的），所以i-j是负数。如果你自己推导时发现符号对不上，多半是展开时漏了交叉项。

查询阶段的while z[k+1] < x也是一个容易踩坑的点。这里严格取小于号，意味着当x恰好等于某条抛物线的区间左端点时，我们仍然使用前一条。这在数值上没有任何影响，因为交点处两条抛物线值相等，选哪条结果一样。但如果你改成<=，在极端情况下可能导致k多走一位，访问越界，所以要小心。

5.3 验证：和scipy结果对齐

算法写完了，不验证等于白写。SciPy里有现成的scipy.ndimage.distance_transform_edt，它返回的是真实欧氏距离（非平方），且背景约定、边界处理都和我们的实现基本一致。我建议你写一段对照代码：

from scipy.ndimage import distance_transform_edt # 生成一张随机二值图，保证至少有一个背景像素 np.random.seed(42) img = np.random.rand(64, 64) > 0.5 img[0, 0] = 0 # 强制一个背景点 my_dt = felzenszwalb_edt(img.astype(np.float64)) scipy_dt = distance_transform_edt(img) print("最大绝对误差:", np.max(np.abs(my_dt - scipy_dt)))

在我的本机测试里，这个最大绝对误差通常在1e-10量级，属于纯浮点误差，可以认为两个结果一致。如果这个值很大，优先检查两个地方：一是背景标记约定，scipy默认distance_transform_edt把非零值当目标，零值当背景，和我们的约定一致，但如果你的输入把前景标成0、背景标成255，那结果会完全反掉；二是INF是否取得足够大，如果距离场里有像素的真实距离超过了你能表示的范围，结果会异常。对于一般图像，INF=1e20在这个算法里是安全的，因为任何中间值都不可能接近这个量级。

验证完之后，我强烈建议你再用一个半圆形的二值图看一眼可视化结果。拿matplotlib画一下距离场的等值线，你会发现等值线是完美的同心圆弧，没有任何方向性偏差——这一步能直观感受到精确EDT和BFS、倒角算法的差距，等值线如果出现方块状锯齿，说明算法里还有近似的成分。

6. 实测量级与选型建议

6.1 不同算法的时间量级对比

我在一张1000×1000、前后景各半的合成二值图上做了不同算法的量级测试，环境是Python 3.11、NumPy 1.24。这里给出的不是精确benchmark，而是让你对数量级有个概念：

这里有个很有意思的点：纯Python实现的Felzenszwalb算法比scipy的C实现慢了约15倍，但它仍然是线性复杂度。如果你自己写一个C扩展，或者用Numba的@njit把那个dt_1d_general函数编译一下，耗时能压到和scipy一个量级。我在项目里就是用Numba包装了这个算法，在三维体数据上跑，2^256的网格大概1.5秒出结果，比调scipy的逐层循环还快一些。

6.2 什么时候用暴力，什么时候用Felzenszwalb

很多人看完复杂度分析会觉得，暴力算法完全没有存在的必要。但实际情况是，暴力算法在两种场景下依然是合理选择：一是图像尺寸极小，比如几十乘几十的图标，暴力算也就几毫秒，没必要引入复杂实现；二是你要处理的是点集而非栅格图像，比如一堆散点找最近邻，这时候暴力 + 空间索引（KD-Tree）才是标准思路，EDT算法反而是为规则栅格量身定制的。

更常见的选择是在"scipy现成函数"和"自己实现Felzenszwalb"之间做权衡。如果你的目标只是拿到距离场结果，直接调scipy就是最优解，又快又稳。但如果你要在距离变换的基础上做定制，比如输出距离平方、同时返回最近背景点的坐标、处理带权重的背景、做三维EDT，那理解Felzenszwalb的实现思路就特别重要了。它并不难，核心就是那几十行代码，掌握了它，你可以轻松地把算法从二维扩展到三维，甚至推广到任意分辨率。

关于三维扩展我再多说一句：二维EDT是"行方向二值+列方向一般"，三维EDT就是"z轴方向一般+每个z平面内的二维EDT"，本质上就是不断套娃调用dt_1d_general。理解了抛物线栈这个公共底座，三维、四维的距离变换都不再是新知识，只是把同一段逻辑调用更多次而已。

6.3 再往前走一步：从二值图到任意采样函数

Felzenszwalb这篇论文真正的价值，远远不止二值图的欧氏距离变换。它处理的对象是"任意采样函数f"，也就是说，f不一定是0/INF二值，可以是任何数组。这带来了一类新的应用场景：把任意标量场转换成一个"带几何意义"的距离场。

举几个例子。在图像分割里，你可以把模型输出的概率图当作f，计算每个像素到"高概率目标中心"的距离，这样得到的距离场天然带有置信度加权的几何信息。在路径规划里，可以把障碍物的代价图当作f，计算出来的距离场不是到障碍物表面的欧氏距离，而是考虑了代价衰减的"广义距离"，路径搜索会更平滑。在三维重建里，用TSDF（截断符号距离场）表示曲面时，对每个体素格点做这种一般化距离变换，就能快速更新整个场的数值。

我自己最常用的一个场景是生成SDF：给定一个二值或灰度的形状图，先算出精确EDT，然后根据像素在形状内外赋予正负号，就得到符号距离场。这个SDF可以直接用于字体描边、阴影生成、形状融合、网格平滑等等。整套流程的耗时瓶颈不在算法本身，而在后续的符号判断和归一化操作，EDT这一步永远是毫秒级。

所以我的建议是：把dt_1d_general这个函数当作一个"算法积木"存进你的工具库。它就像是排序里的快排、图论里的Dijkstra，是一个基础算法零件，遇到"采样函数到距离场"这类问题时，直接拿出来用就行。理解了它的抛物线栈原理，你甚至不需要背代码，现场推一遍交点公式就能写出来——这比背一堆现成库的手册要有用得多。