MATLAB离散曲线曲率计算：差分、样条与平滑实践

简介一套基于MATLAB的曲率计算工具包面向需要处理二维与三维曲线几何特征的工程师与科研人员覆盖曲率、曲率半径求解与可视化可应用于图像处理、几何建模、机器人路径规划等方向。压缩包共6个文件以4个m脚本为主另有1个PDF文档与1个txt许可文件整体仅85KB轻量且便于携带。m脚本中包含曲率核心计算函数、圆心求解辅助函数以及2D/3D演示脚本可直接运行查看效果PDF文档系统讲解三维曲率的理论基础与MATLAB实现细节适合配合源码边读边练。该资源已有1962人学习下载在相关领域具备一定的参考价值。通过示例代码可快速生成曲线并绘制曲率半径向量帮助使用者将抽象的曲率公式转化为直观图形提升对曲线局部弯曲程度的理解与实际应用能力。 做路径规划、图像轮廓分析或者三维点云处理的时候绕不开一个问题怎么量化一条曲线有多“弯”。很多初学者上来就翻公式背下来曲率等于半径的倒数但真拿到一堆离散点数据还是不知道怎么下手。我之前在MATLAB里折腾过一段时间的曲率计算一开始也是踩了不少坑——算出来的结果要么噪声大得没法看要么符号方向完全搞反。这篇文章就把我总结出来的MATLAB曲率计算思路、具体步骤和代码示例整理出来希望能帮你少走一些弯路。1. 曲率这个概念多数人理解得太浅了1.1 先厘清曲率在数学上的定义曲率描述的是曲线在某一点上的弯曲程度。数学上对于一条参数化曲线 \(\mathbf{r}(t) (x(t), y(t))\)曲率 \(\kappa\) 的定义是切线方向角相对于弧长的变化率\[ \kappa \left| \frac{d\theta}{ds} \right| \]其中 \(\theta\) 是切线的倾角\(s\) 是弧长。展开到二维平面直角坐标系下曲率的显式公式是\[ \kappa \frac{|xy - yx|}{(x^2 y^2)^{3/2}} \]这里 \(x, y\) 是一阶导数\(x, y\) 是二阶导数。如果是显式函数曲线 \(y f(x)\)公式可以简化为\[ \kappa \frac{|y|}{(1 y^2)^{3/2}} \]这个公式大家应该都不陌生但要注意它的前提曲线必须是光滑的、导数存在且连续。你在MATLAB里随手画的正弦波 \(y \sin(x)\) 当然满足这个条件但现实工程中的曲线比如道路中心线、零件轮廓线、轨迹数据往往只是一堆离散的测量点。这就引出了第一个核心矛盾连续微分理论和离散数值计算之间存在一个“翻译”过程而这个过程恰恰是误差的主要来源。1.2 从连续导数到离散点工程里的曲率计算没那么理想化当我们手里只有离散点 \((x_i, y_i), i 1,2,\dots,n\) 的时候不能直接套上面那个连续公式了得先想办法估计一阶导数和二阶导数。最朴素的想法是用差分\[ x \approx \frac{x_{i1} - x_i}{\Delta t}, \quad x \approx \frac{x_{i1} - 2x_i x_{i-1}}{\Delta t^2} \]思路没错但问题在于差分会放大噪声。如果原始数据是带测量误差的比如GPS轨迹、传感器采集的运动曲线直接差分算出来的导数会剧烈震荡曲率结果基本不可用。所以在离散曲率计算里真正决定成败的不是那个公式而是你用什么方式处理噪声、怎么选差分窗口、怎么处理边界。这些细节我会在后面展开讲。2. MATLAB里实现曲率计算的三条主线2.1 主线一有显式函数时用符号计算或者匿名函数求导直接算如果你手里现成的就是一条解析曲线比如圆、抛物线、螺旋线那事情简单得很直接用符号数学工具箱求导再代入公式即可。syms x y(x) y sin(x); d1 diff(y, x); % 一阶导 d2 diff(y, x, 2); % 二阶导 kappa_sym abs(d2) / (1 d1^2)^(3/2); % 代入具体数值 kappa_val double(subs(kappa_sym, x, 1)); fprintf(x1处的曲率: %.4f\n, kappa_val);如果你不想用符号计算直接用匿名函数和数值微分也可以。一条圆形参数曲线 \(x(t) R\cos t, y(t) R\sin t\)它的曲率恒等于 \(1/R\)这是检验代码对不对的绝佳用例。我后来写算法第一步永远先用圆做基准测试——如果圆都算不对其他曲线就更不用说了。2.2 主线二离散点用差分法计算关键是微分格式的选择这是实际工程里最常用的方式。对于等间距采样点常见的差分格式有前向差分、后向差分和中心差分三种。三种方法精度不同对噪声的敏感度也不同。% 假设 t 等间距步长 h h t(2) - t(1); % 一阶导中心差分 xp zeros(size(x)); xp(2:end-1) (x(3:end) - x(1:end-2)) / (2*h); xp(1) (-3*x(1) 4*x(2) - x(3)) / (2*h); xp(end) (3*x(end) - 4*x(end-1) x(end-2)) / (2*h); % 二阶导二阶中心差分 xpp zeros(size(x)); xpp(2:end-1) (x(3:end) - 2*x(2:end-1) x(1:end-2)) / h^2; xpp(1) xpp(2); xpp(end) xpp(end-1); % 计算曲率 kappa abs(xp .* ypp - yp .* xpp) ./ (xp.^2 yp.^2).^(3/2);中心差分比前向/后向差分精度高一阶O(h²)对比O(h)所以内部点我优先用中心差分。边界点因为缺少左右邻居只能用单向差分近似误差会大一些但一般控制在两端就可以了。这里要提醒一句上面代码里的xpp(1) xpp(2); xpp(end) xpp(end-1);是一种很粗暴的边界处理只是让计算不报错精度别指望。如果边界曲率对结果很重要建议单独用二次多项式拟合来估算边界点的一二阶导。2.3 主线三先拟合再计算用多项式或者样条曲线平滑数据差分的天生缺陷是对噪声敏感。如果数据点本身质量不高比如手工采集的坐标点我一般不走纯差分路线而是先拟合再做微分或者直接用样条工具箱自带的fnder函数对拟合曲线的微分求值。MATLAB的Curve Fitting Toolbox里有一个非常好用的命令spline做插值然后用fnder求导。% 用三次样条插值平滑曲线 pp csape(t, [x; y], periodic); pp_x fnder(pp(1)); pp_y fnder(pp(2)); % 在细密网格上求导 tt linspace(t(1), t(end), 500); dx fnval(pp_x, tt); dy fnval(pp_y, tt); ddx fnval(fnder(pp_x), tt); ddy fnval(fnder(pp_y), tt); % 曲率 kappa abs(dx .* ddy - dy .* ddx) ./ (dx.^2 dy.^2).^(3/2);样条方法的优势是连续性好、可解析求导得到的一二阶导不再是逐点的离散近似而是整条连续曲线上的解析值。但样条也有自己的问题如果数据噪声太大插值样条会“过拟合”噪声反而比差分更糟糕。这种情况下需要改用平滑样条csaps或者把数据质量先整一整再算。拟合参数的选择我会在第四部分专门讲。3. 离散点的曲率计算真正决定结果好坏的重头戏3.1 差分窗口的取舍窗口越大越平滑但会丢失细节很多人以为差分窗口越大越好其实不是。窗口大参与平均的相邻点变多局部噪声对导数的干扰确实下降了但代价是曲率峰值被“抹平”。比如你有一条在局部有急剧转弯的曲线窗口过大会把那个尖锐的转弯识别成一个缓弯。反过来窗口取最小就是相邻点直接差分对噪声的放大最严重。我实践下来一个比较稳的套路是先用原始数据直接算一遍曲率观察曲率曲线的噪声幅度再逐步拉大差分窗口观察峰值位置是否发生明显偏移。如果拉大窗口后曲率曲线整体形状变化不大说明原数据噪声不严重窗口可以取小一点如果形状明显变了就得配合平滑手段而不是一味加窗口。MATLAB里可以用diff或者gradient配合卷积平滑实现这个思路。gradient函数本身用的是二阶精度中心差分内部点精度不错边界点只做单边差分所以边界还是需要自己处理。3.2 一个隐藏的大坑坐标尺度不一致会让曲率计算失真这个坑我踩过一次印象极深。有一组点数据x方向的范围是0到500米y方向的范围是0到2毫米——因为那条曲线本身就是一条近似水平的细长弧线。直接用上面的代码计算算出来的曲率数值大得离谱而且没有任何物理意义。原因很简单曲率公式对坐标尺度非常敏感。如果x和y的尺度差了好几个量级\(x^2 y^2\) 这一项会被坐标尺度大的那个方向主导小尺度方向的信息被完全淹没。解决办法有两种。第一种是对坐标做标准化/归一化处理把x和y同时缩放到差不多的量级算完曲率再按比例换算回来。第二种是用弧长参数重新参数化曲线。第二种方法物理意义最清晰因为在弧长参数化下\(x^2 y^2 1\)曲率公式会简化成 \(\kappa |xy - yx|\)不仅计算更简单数值稳定性也大幅提升。3.3 实用组合拳平滑预处理 自适应差分在实际处理测量数据时我的标准流程是先做离群点剔除比如用isoutlier然后做轻度平滑smoothdata或者 Savitzky-Golay 滤波再做差分或样条求导最后计算曲率。% 流程示例 x x_orig; y y_orig; % 1. 剔除离群点 idx_out isoutlier(y, movmedian, 10); x(idx_out) []; y(idx_out) []; % 2. Savitzky-Golay 平滑 x_s sgolayfilt(x, 3, 15); y_s sgolayfilt(y, 3, 15); % 3. 后续差分或样条 % ...这里要多说一句sgolayfilt。它本质是局部多项式拟合加卷积滤波对信号形状的保持能力比普通滑动平均好很多特别适合曲率计算这种需要保留导数信息的场景。普通滑动平均虽然平滑但会把峰值削掉直接导致曲率尖峰被磨损。用Savitzky-Golay的时候多项式阶数取3就够窗口长度建议是多项式阶数的5倍以上但不宜超过数据长度的十分之一。4. 曲率可视化光算出来不够还要会看4.1 把曲率画成沿弧长的曲线比在曲线上标注数字直观得多曲率计算的最终目的往往是要定位“弯曲最剧烈”的位置。直接把曲率值标注在点号位置上不太直观我习惯先把点序列转化成弧长坐标然后画 \(\kappa\) 随弧长变化的曲线这样很容易找到峰值的弧长位置再反查是曲线上哪个位置。弧长累计计算用cumtrapz即可% 计算相邻点距离 seg_len sqrt(diff(x).^2 diff(y).^2); s [0; cumtrapz(t, ones(size(t)))]; % 如果 t 等间距这里直接是累计长度 % 上面这行有点绕简洁做法 s cumtrapz(linspace(0, 1, length(x)), sqrt(gradient(x).^2 gradient(y).^2));算了直接给一个最直观的版本dx gradient(x); dy gradient(y); ds hypot(dx, dy); s cumsum(ds); s s / s(end); % 归一化到 [0, 1] figure; plot(s, kappa, linewidth, 1.5); xlabel(归一化弧长); ylabel(曲率 \kappa); title(曲率沿弧长的变化); grid on;用归一化弧长当横轴的好处是不管曲线整体多长你都能在同一坐标系下比较不同曲线的弯曲分布模式。4.2 在曲线上画出曲率圆直观验证计算对不对数值算完一定要做可视化验证。我最喜欢的一个验证方法是任选一个点用算出的曲率半径画一个内切圆看看这个圆和原始曲线在那一小段是否贴合。如果贴合得很自然说明曲率计算基本靠谱如果圆和曲线交叉跑偏那就要回去查导数计算或者平滑参数了。function draw_curvature_circle(x, y, kappa, idx) % 在 idx 处画出曲率圆 r 1 / kappa(idx); xc x(idx) - r * ???; % 需要法线方向 end曲率圆的中心在曲线法线方向上距离为曲率半径。具体实现时要算该点处的法线方向切线方向旋转90度并注意法线方向的正负。这个验证方法虽然简单但非常能抓错误。我调试代码时至少随机抽5个点画曲率圆检查只靠数值比对很难发现局部错误。4.3 一个直观案例圆形与椭圆验证代码是否正确我的首选测试对象永远是圆。一个半径为5的圆理论上每个点的曲率都是0.2。如果你的代码算出来圆心附近是0.2两端就偏到0.15之类那基本可以断定是边界处理或者差分格式的问题。椭圆 \(x a\cos t, y b\sin t\) 则是一个更好的测试用例因为它的曲率随位置显著变化。长轴端点是曲率最大点短轴端点是曲率最小点。具体值可以通过解析公式验证。用这个测试能同时检查幅值精度和位置精度。5. 几个容易踩的坑和我的习惯做法5.1 参数化方式不对得到的结果没法看有一条二维曲线你可以对x做参数 t 1:length(x)也可以直接用相邻点连接的多段线参数化。如果用等间隔索引当参数但原始点并不是等弧长采样的话计算出的导数会偏向密集区曲率结果也会出现伪尖峰。我的经验是优先用弧长参数化也就是把累计弧长当参数t。做了弧长参数化之后导数计算对采样不均匀的抗性会强很多。一个快速判断方法是画 \(x(t)\) 和 \(y(t)\) 对 \(t\) 的曲线如果两段曲线的斜率变化很剧烈说明参数化可能有问题。5.2 闭合曲线比开口曲线多了两个麻烦处理闭合轮廓比如细胞边缘、零件内孔时首尾相接的地方最容易出错。因为差分计算要使用邻点在闭合曲线里“第一个点”的邻居是“最后一个点”很多人初始条件没配对导致接头处曲率跳变。解决方案分两步第一步是确保闭合即把首点复制一份追加到末尾或者用周期样条代码里的periodic选项第二步是算完曲率后把首尾两个点的曲率取平均作为接头处的近似值。5.3 曲率符号什么时候该关心它什么时候可以无视数学公式里曲率常常取绝对值但工程场景有些需要带符号的曲率。比如路径规划中通常定义左转弯为正、右转弯为负这个时候曲率的正负代表转向方向。带符号曲率的计算方式是% 带符号曲率正表示逆时针弯曲左转 kappa_signed (xp .* ypp - yp .* xpp) ./ (xp.^2 yp.^2).^(3/2);如果你只是关心“弯得多厉害”就用绝对值如果要判断“往哪边弯”就必须保留符号。前提是曲线方向定义要统一——是从左往右算的还是从右往左算的会直接影响符号的正负。我习惯在所有代码开头注释清楚“曲线按起点到终点的方向计算”。5.4 数据清洗比算法选择更重要最后说一个我认为最值得分享的体会曲率计算里算法本身占三成数据前处理占七成。拿到一组带离群点的原始数据你花再多精力调差分数值格式和窗口大小都没用离群点造成的伪曲率尖峰会严重干扰你的判断。所以我的标准动作永远是先画散点图肉眼看数据质量再做离群点剔除和平滑最后才轮到曲率计算。这个过程听起来原始但真的能省掉后面大量的排查时间。6. 代码模板一套可以直接拿去用的曲率计算函数这里给出一套我常用的完整函数输入二维曲线点集输出每个点的曲率并附带上简单的内部验证逻辑。特点是内置了平滑、弧长参数化和边界处理你可以直接复制拿去跑自己的数据。function [kappa, s] compute_curvature(x, y, varargin) % COMPUTE_CURVATURE 计算二维离散曲线的曲率 % 输入: % x, y - 曲线点坐标 (n×1 或 1×n) % Smooth - 是否启用 Savitzky-Golay 平滑, true/false, 默认 true % Window - 平滑窗口长度, 默认 11 % Order - 平滑多项式阶数, 默认 3 % 输出: % kappa - 各点的曲率值 % s - 各点的累计弧长归一化 [0,1] p inputParser; addParameter(p, Smooth, true, islogical); addParameter(p, Window, 11, (x) isnumeric(x) x 1); addParameter(p, Order, 3, (x) isnumeric(x) x 1); parse(p, varargin{:}); x x(:); y y(:); % 强制列向量 % 1. 弧长参数化 dx gradient(x); dy gradient(y); s cumsum(hypot(dx, dy)); s s / max(s); % 归一化到 [0,1] % 2. 平滑可选 if p.Results.Smooth x sgolayfilt(x, p.Results.Order, p.Results.Window); y sgolayfilt(y, p.Results.Order, p.Results.Window); end % 3. 用样条求导减少差分噪声 pp_x csape(s, x, variational); pp_y csape(s, y, variational); [dx, ddx] fnval([fnder(pp_x, 1); fnder(pp_x, 2)], s); [dy, ddy] fnval([fnder(pp_y, 1); fnder(pp_y, 2)], s); % 4. 计算曲率 kappa abs(dx .* ddy - dy .* ddx) ./ (dx.^2 dy.^2).^(3/2); kappa(isnan(kappa)) 0; % 处理除零 % 5. 简单验证: 与直接差分结果对比 dxd gradient(x) ./ gradient(s); dyd gradient(y) ./ gradient(s); ddxd gradient(dxd) ./ gradient(s); ddyd gradient(dyd) ./ gradient(s); kappa_diff abs(dxd .* ddyd - dyd .* ddxd) ./ (dxd.^2 dyd.^2).^(3/2); kappa_diff(isnan(kappa_diff)) 0; delta norm(kappa - kappa_diff, inf); fprintf(样条法与差分法偏差峰值: %.6f\n, delta); if delta 0.1 warning(两种方法偏差较大请检查平滑窗口或数据质量。); end end这个函数里的第五步“验证”是我后来加上的作用是把样条结果和纯差分结果做一个交叉验证。如果两者偏差很大基本说明数据质量有问题光看某一个结果很容易被带偏。实际使用的时候平滑窗口默认11可能不够或者过度需要根据数据的采样密度微调。一个经验性的方法是窗口大小大约取“一个完整波形包含点数的五分之一”。比如你的曲线一个周期大概有100个点那窗口取20左右比较合适如果取太大曲率峰值会被削平。7. 不同应用场景下的注意事项7.1 道路轨迹曲率计算车辆轨迹数据通常由GPS或者惯性导航采集坐标很大经纬度或者UTM坐标但局部弯曲幅度很小。这种情况下我在算曲率之前会把整条轨迹做一个去均值处理即减去起点坐标让数值在计算机精度范围内更友好。另外GPS数据噪声很严重尤其是低速状态下位置漂移明显。我用过的方案里卡尔曼平滑和RTS平滑的效果最好但实现成本高如果不想引入太多依赖sgolayfilt配movmedian离群点剔除可以解决大部分问题。经纬度坐标还有个特殊问题不能直接在经纬度上算曲率因为一度经度和一度纬度对应的物理距离不同纬度越高差异越大。可以先通过lldistkm之类的坐标转换方法转成平面坐标或者做墨卡托投影再计算曲率。7.2 图像轮廓曲率计算做图像轮廓曲率时轮廓点往往是像素级别坐标离散化严重。直接用像素坐标算曲率会得到锯齿状结果。我建议先用多边形近似或者B样条拟合轮廓再做均匀采样最后计算曲率。常用的bwboundaries提取轮廓后用Douglas-Peucker算法减点再用样条插值加密。图像轮廓还有一个高频需求是角点检测。实际上曲率的局部极大值位置通常就是角点位置。你算完曲率后用findpeaks提取曲率峰再设置一个阈值筛选就能轻松定位角点。这个方法比直接用corner函数更可控一点阈值调起来也更直观。7.3 三维曲线的曲率计算三维空间曲线 \((x(t), y(t), z(t))\) 的曲率公式是\[ \kappa \frac{|\mathbf{r} \times \mathbf{r}|}{|\mathbf{r}|^3} \]MATLAB里可以用cross函数实现。三维比二维麻烦的地方在于没有“正负曲率”的概念了弯曲方向由法向量表达。如果还需要副法向量和挠率那就要继续算三阶导噪声问题会被进一步放大这时候预处理就显得更加重要。% 三维曲线曲率 r_prime [gradient(x); gradient(y); gradient(z)]; r_dprime [gradient(gradient(x)); gradient(gradient(y)); gradient(gradient(z))]; cross_prod cross(r_prime, r_dprime, 1); kappa3d vecnorm(cross_prod, 2, 1) ./ (vecnorm(r_prime, 2, 1).^3);这个方法比较初级。真正在三维点云或轨迹上用的时候经常会先用多项式拟合局部点云再算曲率效果比全局差分稳定得多。8. 我对曲率计算的几点心得多算、多画、多验证是我能给的最实用的建议。曲率算法本身不复杂真正的复杂度在数据质量控制和边界情况处理上。拿到一组数据不要急着套公式先花点时间看数据、清洗数据、做弧长参数化最后再选合适的差分或样条方法。每写一个版本都用圆和椭圆做基准测试。这样虽然看起来慢但实际是最快的——因为排查错误的时间会大幅减少。最后分享一个小技巧保存曲率计算结果的时候记得同时保存弧长坐标和原始索引的映射关系。否则你算出来哪一段曲率最大却反查不到对应原始点后面的工作就得重来。这个坑我踩过不止一次希望你不用再踩。本文还有配套的精品资源点击获取