离散点列上计算曲率，本质上是对位置序列做二阶差分。假设轮廓坐标 x(n)、y(n) 带有零均值、标准差为 σ 的随机定位噪声，那么中心差分公式会把噪声项按 1/h² 的比例放大。采样越密，h 越小，噪声放大越严重。因此直接拿原始坐标求曲率，往往得到剧烈振荡的毛刺曲线，根本无法用于角点检测、缺陷识别或运动规划。解决思路是在曲率估计之前对轮廓进行合理平滑，或对曲率序列本身做滤波后处理。
需要注意的是，平滑不是越强越好。窗口过大或低通截止频率过低，会抹掉真实的小曲率变化，让尖锐拐角变成圆弧。因此需要根据噪声幅值和目标特征尺度来选方法。下面先从数学上说明噪声为什么会被曲率计算强烈放大，再讨论可用的平滑与滤波策略。
一、曲率噪声的来源与数学特性
对于参数曲线 r(t)=(x(t),y(t))，曲率通常定义为 k=|r′×r″|/|r′|³。对于离散点列，工程上常用差分代替导数。以 x 坐标为例，一阶中心差分为 dx(i)≈(x(i+1)−x(i−1))/(2h)，二阶中心差分为 ddx(i)≈(x(i+1)−2x(i)+x(i−1))/h²。如果每个采样点存在独立随机误差 ε(i)，其标准差为 σ，那么二阶差分后的误差项为 (ε(i+1)−2ε(i)+ε(i−1))/h²。该误差的方差是 (1²+2²+1²)σ²/h⁴ = 6σ²/h⁴，标准差约为 2.45σ/h²。
这意味着当采样步长 h 较小时，坐标噪声会被急剧放大。例如图像轮廓提取中，边缘点定位误差通常在 0.3 到 1 个像素之间，如果 h 取 2 像素，二阶导数噪声标准差可能达到 0.18 到 0.61 像素⁻¹。很多真实轮廓的曲率值本身就小于 0.1，因此直接计算出的曲率信噪比非常低。除了定位误差，图像离散化、边缘检测算子的亚像素估计误差、传感器抖动以及轮廓追踪过程中的舍入误差也会引入高频扰动。
更麻烦的是，曲率计算中的分母 |r′|³ 对速度波动也很敏感。如果轮廓点分布不均匀，某些区域点距明显偏小或偏大，即使坐标没有噪声，曲率值也会出现上下跳动。因此稳定的曲率估计通常需要同时处理坐标噪声和参数化不均匀问题。平滑处理与滤波可以抑制坐标序列中的高频分量，但不能完全替代重采样或弧长参数化。理解这一点后，就可以更有针对性地选择方法。
二、基于局部拟合的平滑处理
局部拟合类方法假设轮廓点在邻域内近似满足低阶多项式或局部平滑模型。通过在窗口内拟合模型，再取中心位置的平滑值或导数值，可以有效降低随机噪声的影响。这类方法实现简单，参数直观，是曲率计算前处理中最常用的手段。
移动平均和高斯平滑是最基础的两类。移动平均给窗口内所有点相同权重，计算量低，但边界效应明显，且对窗口内所有频率分量做同样的抑制，容易造成轮廓收缩。高斯平滑使用高斯核作为权重，离中心越近权重越大，能更好地保持轮廓中心位置，同时保留更多低频形状信息。对于闭合轮廓，高斯卷积可以直接通过循环卷积实现，避免端点不连续问题。高斯核的标准差 σ 控制平滑强度，σ 越大，保留的低频形状越少，角点也会变得更圆。
Savitzky-Golay 滤波在局部窗口内做多项式最小二乘拟合，直接用拟合多项式在中心点处的一阶、二阶导数作为平滑后的导数估计。与简单平滑后差分不同，它同时完成了平滑和求导，能在一定程度上避免“先平滑再差分”造成两次误差累积。窗口长度通常取奇数，多项式阶数一般取 2 到 4。阶数过低会过度平滑，阶数过高则容易保留噪声。对于曲率估计，推荐用窗口长度 7 到 15、三次多项式，具体取决于轮廓点密度。
import numpy as np
from scipy.ndimage import gaussian_filter1d
from scipy.signal import savgol_filter
def curvature(x, y):
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
dx = np.gradient(x)
dy = np.gradient(y)
ddx = np.gradient(dx)
ddy = np.gradient(dy)
speed = np.sqrt(dx ** 2 + dy ** 2)
speed[speed < 1e-9] = 1e-9
return (dx * ddy - dy * ddx) / speed ** 3
# 示例：对轮廓先做 Savitzky-Golay 平滑，再计算曲率
x_sg = savgol_filter(x, window_length=15, polyorder=3)
y_sg = savgol_filter(y, window_length=15, polyorder=3)
k_sg = curvature(x_sg, y_sg)
# 示例：使用高斯平滑后再计算曲率
x_g = gaussian_filter1d(x, sigma=2.0, mode='wrap')
y_g = gaussian_filter1d(y, sigma=2.0, mode='wrap')
k_g = curvature(x_g, y_g)
样条平滑是另一类重要方法。利用三次样条对轮廓点进行拟合，可以直接得到连续且二阶可导的参数表达式，然后用解析公式计算曲率。样条拟合通过平滑参数 s 控制逼近程度：s 越大，拟合越平滑但可能偏离原始点；s 越小，插值越精确但噪声会重新出现。样条方法适合需要连续曲率输出的场景，例如 CNC 轨迹规划或机器人运动控制。缺点是当轮廓点非常密集或噪声较大时，样条拟合可能产生过冲现象，需要额外做节点选择或最小二乘样条逼近。
三、滤波算法在曲率计算中的应用
如果不想改变原始轮廓坐标，也可以先计算出曲率序列，再对曲率序列做滤波。不过更推荐的做法是对坐标序列先滤波再计算曲率，因为曲率是二阶量，直接对曲率滤波需要更大的窗口，容易把真实曲率峰值一起平滑掉。对坐标低频化后再求曲率，通常能得到更稳定的结果。
低通滤波在频域直接截断高频分量，对闭合轮廓特别有效。将 x(n)、y(n) 分别做一维 FFT，保留低频部分，将高频系数置零，再逆变换回空间域。这样能去除大部分高频噪声，同时保持轮廓的整体闭合性。但直接截断会在频域产生矩形窗效应，逆变换后可能出现振铃。可以通过加汉宁窗或使用高斯频域衰减减弱振铃。对于非闭合轮廓，FFT 方法需要先做延拓或镜像填充，否则端点处会出现明显不连续。
def lowpass_contour(x, y, keep_ratio=0.1):
n = len(x)
fx = np.fft.rfft(x)
fy = np.fft.rfft(y)
k = int(len(fx) * keep_ratio)
if k < 1:
k = 1
fx[k:] = 0
fy[k:] = 0
x_lp = np.fft.irfft(fx, n)
y_lp = np.fft.irfft(fy, n)
return x_lp, y_lp
中值滤波对孤立脉冲噪声和离群点非常有效。如果轮廓中存在少量错误跟踪的点，中值滤波可以快速剔除这些点的影响，而不会像均值滤波那样把异常值扩散到邻域。但中值滤波会移动真实角点的位置，尤其是窗口较大时，角点可能被削平或偏移。因此中值滤波更适合作为前处理步骤去除粗大误差，不建议直接用于精细曲率估计。
双边滤波在空间距离权重之外引入值域差异权重，可以在平滑的同时保持边缘或角点。对二维图像中的轮廓坐标，也可以把 x、y 看作两个通道做双边滤波。它的缺点是需要同时调节空间窗口和值域带宽两个参数，计算成本也比线性滤波高。如果轮廓曲线没有明显的角点需要保留，双边滤波的优势并不明显。小波去噪则适合非平稳噪声，通过阈值化小波系数去除局部高频成分，能得到较自然的平滑效果，但参数选择相对复杂。
四、参数选择、验证与工程建议
参数选择需要与轮廓的最小特征尺度匹配。如果希望保留的角点半径约为 R，那么平滑窗口对应的空间尺度应明显小于 R。对于图像轮廓，窗口长度通常选择 5 到 15 个像素，高斯 σ 取 1 到 3 个像素，FFT 低频保留比例取 0.05 到 0.2。若轮廓点总数较多，可以先做等弧长重采样，使点距均匀，然后再固定窗口长度，这样参数在不同缩放尺度下更一致。闭合轮廓建议统一采用循环边界条件，避免端点处产生人为曲率突变。
验证平滑效果时，建议构造已知几何形状的合成轮廓。例如生成一个标准圆或椭圆，加入不同强度的随机噪声，分别计算平滑前后曲率序列与理论曲率之间的均方根误差。也可以标记出真实角点位置，统计平滑后角点偏移量和漏检率。实测中常见的情况是：高斯平滑计算快但角点变圆；Savitzky-Golay 对角点和峰值的保持更好；中值滤波能有效去除脉冲但会让角点偏移；FFT 低通在闭合周期轮廓上效果稳定，但对非闭合轮廓端点处理较麻烦。
工程上比较稳妥的流程是：先做中值滤波去除孤立粗差，再用 Savitzky-Golay 或高斯平滑对坐标序列进行低通处理，然后计算曲率。如果对速度不敏感，可以直接计算曲率后再做一次轻量高斯平滑以去除残余毛刺。对曲率峰值点应用非极大值抑制可以避免同一角点周围产生多个候选项。最后，边界区域和采样稀疏段应裁剪或单独处理，防止这些区域污染整体判断。
稳定提取曲率特征没有一套参数可以适应所有场景。噪声幅值、轮廓复杂程度、角点锐利程度以及后续任务对曲率误差的容忍度都会影响最优方法。关键是根据实际数据做小样本实验，观察平滑前后曲率曲线是否保留了真实拐点，同时抑制了高频抖动。只有在噪声模型和特征尺度之间找到平衡，才能让曲率估计真正服务于角点检测、形状匹配和运动规划等任务。