目录1.1 核心原理1.2 解决的问题1.3 目的2. 样条的性质基于数值分析教材3. 代码、函数详细说明3.1 Lei_SplineInterpolator.h代码3.2 Lei_SplineInterpolator.cpp代码3.3 main函数测试代码3.4运行结果3.6 函数说明◆ Lei_SplineInterpolator() / ~Lei_SplineInterpolator()◆ setPoints()◆ setBoundaryConditions()◆ calculateSecondDerivatives()◆ interpolate()◆ interpolateArray()◆ generateDenseData()◆ linearInterpolate()4.应用注意事项4.1边界条件的物理意义选择极其重要4.2严禁外推Extrapolation4.3浮点精度与病态矩阵4.4性能优化策略4.5线程安全性5.参考资料1.1 核心原理三次样条插值Cubic Spline Interpolation是一种分段多项式插值方法。它的数学本质是在给n个数据点之间构造n−1段三次多项式:。为了确保曲线平滑三次样条要求在所有内部节点处满足函数值连续、一阶导数连续、二阶导数连续。在实际求解时为了降低计算复杂度通常引入节点处的二阶导数作为未知量。通过推导最终转化为求解一个严格对角占优的三对角线性方程组三弯矩方程并利用追赶法Thomas Algorithm在O(n)的时间复杂度内解出所有从而确定整条曲线。1.2 解决的问题高次多项式插值的震荡问题龙格现象使用单一高次多项式如拉格朗日插值拟合大量数据点时边缘会出现剧烈的震荡导致结果完全不可用。三次样条通过分段低次多项式完美规避了这一问题。线性插值的折线问题线性插值直线连接虽然稳定但在节点处一阶导数不连续有尖角二阶导数为0。这在物理意义上意味着速度和加速度存在突变用于曲面重构或加工轨迹时会导致机床剧烈震动。1.3 目的在离散的工程测量点之间构建一条不仅穿过所有已知点而且整体极度平滑曲率连续的参数曲线。它是点云预处理、曲面重构、数控加工轨迹光顺的数学基石。2. 样条的性质基于数值分析教材在教材中三次样条函数S(x)通常具备以下四条核心性质性质1插值性曲线严格通过给定的所有数据点即。性质2一阶导数连续性在相邻两段样条的连接点处一阶导数斜率是连续的即。这确保了曲线没有折角。性质3二阶导数连续性在相邻两段样条的连接点处二阶导数曲率也是连续的即。这确保了曲线曲率平滑过渡。性质4边界条件端点条件为了唯一定义样条曲线必须补充两个额外条件。工程中常用的有4a. 自然样条Natural Spline假设两端点不受弯矩S′′(x0)0且 S′′(xn−1)0。适用于毫无端点物理约束的自由曲线。4c. 钳制样条Clamped Spline给定端点的一阶导数斜率S′(x0)v0且 S′(xn−1)vn−1。适用于已知速度、转角等物理约束的场景。3. 代码、函数详细说明3.1 Lei_SplineInterpolator.h代码/** * file Lei_SplineInterpolator.h * brief 基于数值分析的三次样条插值器。 * details 支持自然样条与钳制样条包含数据排序、边界条件处理及安全防御机制。 * author Lei * date 2026-10-11 */ #ifndef LEI_SPLINE_INTERPOLATOR_H #define LEI_SPLINE_INTERPOLATOR_H #include vector #include cmath #include stdexcept #include iostream #include iomanip #include limits #include algorithm /** * class Lei_SplineInterpolator * brief 三次样条插值器 * * 该类提供了一套完整的三次样条插值解决方案。它在内部对输入数据点进行排序 * 基于三对角矩阵和追赶法Thomas Algorithm求解二阶导数并支持自然样条 * 和钳制样条两种边界条件。 */ class Lei_SplineInterpolator { public: /** * brief 构造函数 */ Lei_SplineInterpolator() default; /** * brief 设置原始数据点 (X, Y) * details 内部会对 X 数据进行自动排序。排序后如果发现 X 存在重复值将抛出异常。 * param x 自变量数组 * param y 因变量数组 * throws std::invalid_argument 如果 X 和 Y 长度不一致或数据点少于 2 个或存在重复的 X 值 */ void setPoints(const std::vectordouble x, const std::vectordouble y); /** * brief 设置边界条件 * param initialBoundary 初始边界的一阶导数默认 1.0E30表示自然样条 * param finalBoundary 最终边界的一阶导数默认 1.0E30表示自然样条 * note 若输入 1.0E30内部自动切换为自然样条边界二阶导数为 0 * throws std::invalid_argument 如果输入为 NaN */ void setBoundaryConditions(double initialBoundary 1.0E30, double finalBoundary 1.0E30); /** * brief 计算二阶导数数组 g(x_i) * details 组装三对角矩阵并使用追赶法求解。这是样条插值的核心步骤。 * return 包含 n 个元素的双精度数组即 g(x_0) 到 g(x_{n-1}) * throws std::runtime_error 如果数据点尚未设置 */ std::vectordouble calculateSecondDerivatives(); /** * brief 单点插值 * param xi 待插值的 x 坐标 * return 插值结果 y 坐标 * throws std::runtime_error 如果尚未调用 calculateSecondDerivatives */ double interpolate(double xi) const; /** * brief 批量插值 * param xiArray 待插值的 x 坐标数组 * return 插值结果数组 * throws std::runtime_error 如果尚未调用 calculateSecondDerivatives */ std::vectordouble interpolateArray(const std::vectordouble xiArray) const; /** * brief 一维数组线性插值 * param xi 待插值的 x 坐标 * return 线性插值结果 y 坐标。若 xi 超出范围则返回边界值。 * throws std::runtime_error 如果数据点尚未设置 */ double linearInterpolate(double xi) const; /** * brief 根据插值因子生成更密集的数据点 (对应数组插值因子功能) * param factor 插值因子必须大于等于 1。表示在每个原始区间内插入的细分数。 * param outX 输出的密集 X 数组 * param outY 输出的密集 Y 数组 (按样条或线性插值计算得出) * param useSpline 是否使用三次样条插值 (true)否则使用线性插值 (false) * throws std::invalid_argument 如果插值因子小于 1或数据尚未设置 * throws std::runtime_error 如果未预先计算二阶导数且 useSpline 为 true */ void generateDenseData(int factor, std::vectordouble outX, std::vectordouble outY, bool useSpline true) const; private: std::vectordouble m_x; /// 经过排序的原始 X 数组 std::vectordouble m_y; /// 对应排序后的 Y 数组 std::vectordouble m_M; /// 二阶导数数组 g(x_i) double m_initialBoundary; /// 初始边界条件一阶导数 double m_finalBoundary; /// 最终边界条件一阶导数 bool m_isNaturalSpline; /// 是否为自然样条 bool m_isDataValid; /// 数据是否有效 bool m_isBuilt; /// 是否已计算二阶导数 /** * brief 求解三对角线性方程组 (追赶法 / Thomas Algorithm) * param a 下对角线 * param b 主对角线 * param c 上对角线 * param d 右端向量 * return 解向量 */ std::vectordouble solveThomas(const std::vectordouble a, const std::vectordouble b, const std::vectordouble c, const std::vectordouble d) const; }; #endif // LEI_SPLINE_INTERPOLATOR_H3.2 Lei_SplineInterpolator.cpp代码/** * file Lei_SplineInterpolator.cpp * brief Lei_SplineInterpolator 类的实现 */ #include Lei_SplineInterpolator.h void Lei_SplineInterpolator::setPoints(const std::vectordouble x, const std::vectordouble y) { // 1. 防御性检查长度一致且不少于2个点 if (x.size() ! y.size()) { throw std::invalid_argument([防御性编程] X 和 Y 数组的长度必须一致。); } if (x.size() 2) { throw std::invalid_argument([防御性编程] 至少需要 2 个数据点才能进行插值。); } // 2. 处理乱序数据将 X 和 Y 绑定后排序 std::vectorstd::pairdouble, double points; points.reserve(x.size()); for (size_t i 0; i x.size(); i) { // 检查是否为 NaN if (std::isnan(x[i]) || std::isnan(y[i])) { throw std::invalid_argument([防御性编程] 输入数据包含 NaN。); } points.emplace_back(x[i], y[i]); } // 按 X 坐标进行升序排序 std::sort(points.begin(), points.end(), [](const std::pairdouble, double a, const std::pairdouble, double b) { return a.first b.first; }); // 3. 排序后检查 X 是否严格递增是否存在重复的 X 值 m_x.resize(points.size()); m_y.resize(points.size()); for (size_t i 0; i points.size(); i) { m_x[i] points[i].first; m_y[i] points[i].second; if (i 0 std::abs(m_x[i] - m_x[i - 1]) 1e-12) { throw std::invalid_argument([防御性编程] 排序后发现 X 坐标存在重复值无法进行样条插值。); } } m_isDataValid true;//数据有效标志 m_isBuilt false;//是否计算二阶导数标志 } void Lei_SplineInterpolator::setBoundaryConditions(double initialBoundary, double finalBoundary) { // 防御性检查边界值不能是 NaN if (std::isnan(initialBoundary) || std::isnan(finalBoundary)) { throw std::invalid_argument([防御性编程] 边界条件不能为 NaN。); } // 工程约定大于等于 1.0E30 视为无穷大即自然样条二阶导数为 0 const double INFINITY_BOUNDARY_FLAG 1.0E30; m_isNaturalSpline false; if (initialBoundary INFINITY_BOUNDARY_FLAG finalBoundary INFINITY_BOUNDARY_FLAG) { m_isNaturalSpline true;//自然样条 } m_initialBoundary initialBoundary;//初始边界条件 m_finalBoundary finalBoundary;//最终边界条件 m_isBuilt false; // 边界条件改变需要重新计算 } std::vectordouble Lei_SplineInterpolator::calculateSecondDerivatives() { if (!m_isDataValid) { throw std::runtime_error([防御性编程] 数据点尚未设置或无效无法计算二阶导数。); } if (m_isBuilt) return m_M;//已经有了二阶导数 size_t n m_x.size(); m_M.resize(n, 0.0); // 只有两个点退化为线性插值二阶导数均为 0 if (n 2) { m_M[0] 0.0; m_M[1] 0.0; m_isBuilt true; return m_M; } // 计算步长 h std::vectordouble h(n - 1); for (size_t i 0; i n - 1; i) { h[i] m_x[i 1] - m_x[i];//δi } // 三对角矩阵的系数 std::vectordouble a(n, 0.0); // 下对角线 std::vectordouble b(n, 0.0); // 主对角线 std::vectordouble c(n, 0.0); // 上对角线 std::vectordouble d(n, 0.0); // 右端向量 // 1. 组装内点方程 (i 1 to n-2) for (size_t i 1; i n - 1; i) { a[i] h[i - 1]; b[i] 2.0 * (h[i - 1] h[i]); c[i] h[i]; d[i] 6.0 * ((m_y[i 1] - m_y[i]) / h[i] - (m_y[i] - m_y[i - 1]) / h[i - 1]); } // 2. 处理边界条件 if (m_isNaturalSpline) { // 自然样条M_0 0, M_{n-1} 0 b[0] 1.0; c[0] 0.0; d[0] 0.0; a[n - 1] 0.0; b[n - 1] 1.0; d[n - 1] 0.0; } else { // 钳制样条给定一阶导数 b[0] 2.0 * h[0]; c[0] h[0]; d[0] 6.0 * ((m_y[1] - m_y[0]) / h[0] - m_initialBoundary); a[n - 1] h[n - 2]; b[n - 1] 2.0 * h[n - 2]; d[n - 1] 6.0 * (m_finalBoundary - (m_y[n - 1] - m_y[n - 2]) / h[n - 2]); } // 3. 追赶法求解 m_M solveThomas(a, b, c, d); m_isBuilt true; return m_M; } double Lei_SplineInterpolator::interpolate(double xi) const { if (!m_isBuilt) { throw std::runtime_error([防御性编程] 尚未计算二阶导数请先调用 calculateSecondDerivatives()。); } if (std::isnan(xi)) { return std::numeric_limitsdouble::quiet_NaN(); } // 处理边界外的情况返回 NaN 或抛出异常这里选择返回 NaN 供上层决定 if (xi m_x.front() || xi m_x.back()) { return std::numeric_limitsdouble::quiet_NaN(); } size_t n m_x.size(); if (n 2) { // 两点退化为线性 double t (xi - m_x[0]) / (m_x[1] - m_x[0]); return m_y[0] t * (m_y[1] - m_y[0]); } // 找到 xi 所在的区间 [x_i, x_{i1}] auto it std::upper_bound(m_x.begin(), m_x.end(), xi);//二分查找必须是已经排好序的 size_t i std::distance(m_x.begin(), it) - 1; if (i n - 1) i n - 2; double h_i m_x[i 1] - m_x[i]; double A (m_x[i 1] - xi) / h_i; double B 1.0 - A; // 权重系数 C 和 D double C (1.0 / 6.0) * (A * A * A - A) * (h_i * h_i); double D (1.0 / 6.0) * (B * B * B - B) * (h_i * h_i); return A * m_y[i] B * m_y[i 1] C * m_M[i] D * m_M[i 1]; } std::vectordouble Lei_SplineInterpolator::interpolateArray(const std::vectordouble xiArray) const { if (!m_isBuilt) { throw std::runtime_error([防御性编程] 尚未计算二阶导数无法进行批量插值。); } std::vectordouble result; result.reserve(xiArray.size()); for (double xi : xiArray) { result.push_back(interpolate(xi)); } return result; } double Lei_SplineInterpolator::linearInterpolate(double xi) const { //实现一维数组的线性插值Linear Interpolation也就是把相邻的两个数据点用直线连接起来求中间任意位置的值。 if (!m_isDataValid || m_x.empty()) { throw std::runtime_error([防御性编程] 数据点尚未设置无法进行线性插值。); } if (std::isnan(xi)) { return std::numeric_limitsdouble::quiet_NaN(); } size_t n m_x.size(); if (xi m_x.front()) return m_y.front(); if (xi m_x.back()) return m_y.back(); auto it std::upper_bound(m_x.begin(), m_x.end(), xi);//二分查找必须是已经排好序的 size_t i std::distance(m_x.begin(), it) - 1; if (i n - 1) i n - 2; double h_i m_x[i 1] - m_x[i]; return m_y[i] (xi - m_x[i]) * (m_y[i 1] - m_y[i]) / h_i; } std::vectordouble Lei_SplineInterpolator::solveThomas(const std::vectordouble a, const std::vectordouble b, const std::vectordouble c, const std::vectordouble d) const { size_t n d.size(); std::vectordouble c_prime(n, 0.0); std::vectordouble d_prime(n, 0.0); std::vectordouble x(n, 0.0); // 防御性检查主对角线元素不能为 0一般情况下对角占优不会发生但工程上需要兜底 if (std::abs(b[0]) 1e-15) { throw std::runtime_error([防御性编程] 追赶法求解失败主对角线元素接近 0。); } // 向前消元 c_prime[0] c[0] / b[0]; d_prime[0] d[0] / b[0]; for (size_t i 1; i n; i) { double m b[i] - a[i] * c_prime[i - 1]; if (std::abs(m) 1e-15) { throw std::runtime_error([防御性编程] 追赶法求解失败出现除零错误。); } c_prime[i] (i n - 1) ? c[i] / m : 0.0; d_prime[i] (d[i] - a[i] * d_prime[i - 1]) / m; } // 向后回代 x[n - 1] d_prime[n - 1]; for (int i static_castint(n) - 2; i 0; --i) { x[i] d_prime[i] - c_prime[i] * x[i 1]; } return x; } void Lei_SplineInterpolator::generateDenseData(int factor, std::vectordouble outX, std::vectordouble outY, bool useSpline) const { // 防御性编程校验因子和状态 if (factor 1) { throw std::invalid_argument([防御性编程] 插值因子必须大于等于 1。); } if (!m_isDataValid || m_x.empty()) { throw std::runtime_error([防御性编程] 数据点尚未设置无法生成密集数据。); } if (useSpline !m_isBuilt) { throw std::runtime_error([防御性编程] 使用样条插值前请先调用 calculateSecondDerivatives()。); } size_t n m_x.size(); outX.clear(); outY.clear(); // 预留空间防止频繁内存分配 (假设原始区间数为 n-1) // 总点数 (n-1) * factor 1 (如果包含首尾点) outX.reserve((n - 1) * factor 1); outY.reserve((n - 1) * factor 1); // 如果是因子 1直接复制原始数据 if (factor 1) { outX m_x; outY m_y; return; } // 遍历每个区间 [x_i, x_{i1}] for (size_t i 0; i n - 1; i) { double step (m_x[i 1] - m_x[i]) / static_castdouble(factor); // 在当前区间内生成 factor 个点 (包含起点不包含终点) for (int j 0; j factor; j) { double xi m_x[i] j * step; outX.push_back(xi); if (useSpline) { outY.push_back(interpolate(xi)); } else { outY.push_back(linearInterpolate(xi)); } } } // 补上最后一个原始数据点 outX.push_back(m_x.back()); outY.push_back(m_y.back()); }3.3 main函数测试代码#include Lei_SplineInterpolator.h #include iostream #include vector #include iomanip int main() { // 故意打乱原始数据 X 的顺序验证内部排序功能 std::vectordouble X { 3.00, 0.00, 5.00, 1.20, 2.00, 3.50, 4.00 }; std::vectordouble Y { 150.00, 1.00, 220.00, 2.70, 7.40, 140.00, 80.00 }; try { Lei_SplineInterpolator spline; spline.setPoints(X, Y); spline.setBoundaryConditions(1.0E30, 1.0E30); // 自然样条 spline.calculateSecondDerivatives(); // // 测试 1指定单点插值 // std::cout 1. 单点插值测试: std::endl; std::cout std::fixed std::setprecision(2); std::cout X 1.36, Y spline.interpolate(1.36) (预期 -1.56) std::endl; // // 测试 2使用插值因子生成密集数组 // int factor 4; // 在每个原始数据点之间插入 4 个点 std::vectordouble denseX, denseY; std::cout \n2. 使用插值因子 factor 生成密集数据 (样条插值): std::endl; spline.generateDenseData(factor, denseX, denseY, true); std::cout 生成的密集数据点总数: denseX.size() std::endl; for (size_t i 0; i denseX.size(); i) { std::cout X denseX[i] , Y denseY[i] std::endl; } // 对比 LabVIEW 的示例插值数据 std::cout \n3. 与标准数据对比 (取特定 X 验证): std::endl; std::vectordouble checkX { 0.00, 0.24, 0.48, 0.72, 0.96, 1.20 }; for (double xi : checkX) { double yi spline.interpolate(xi); std::cout X xi , 插值 Y yi std::endl; } // // 测试 3防御性编程测试 // std::cout \n4. 防御性编程测试: std::endl; // 测试非法插值因子 try { spline.generateDenseData(0, denseX, denseY, true); } catch (const std::exception e) { std::cout 捕捉到异常 (因子为0): e.what() std::endl; } // 测试越界插值 double nan_val spline.interpolate(6.0); std::cout 越界插值 X6.0, 返回值是 NaN 吗? (std::isnan(nan_val) ? 是 : 否) std::endl; } catch (const std::exception e) { std::cerr 发生异常: e.what() std::endl; return -1; } return 0; }3.4运行结果3.6 函数说明◆ Lei_SplineInterpolator() / ~Lei_SplineInterpolator()Lei_SplineInterpolator() default;描述构造函数与析构函数。默认构造不分配额外内存内部状态初始为未设置数据m_isDataValid falsem_isBuilt false。◆ setPoints()void Lei_SplineInterpolator::setPoints(const std::vectordouble x,const std::vectordouble y);描述设置原始数据点并对输入数据进行排序。内部将 X 和 Y 绑定为 std::pair使用 std::sort 进行升序排序排序后严格检查是否存在重复的 X 坐标防除零。Parameters[in] x自变量数组顺序可不一致内部会自动排序。[in] y因变量数组必须与 X 长度一致。Exceptionsstd::invalid_argument当 X 与 Y 长度不一致、数据点少于 2 个、存在 NaN、或排序后 X 存在重复值时抛出。◆ setBoundaryConditions()void Lei_SplineInterpolator::setBoundaryConditions(double initialBoundary 1.0E30, double finalBoundary 1.0E30);描述设置样条的边界条件。默认值为 1.0E30代表无穷大当检测到两端点均大于等于该阈值时自动设置为自然样条性质4a否则根据给定的导数值切换为钳制样条性质4c。Parameters[in] initialBoundary初始边界的一阶导数 g(x0)。[in] finalBoundary最终边界的一阶导数 g(xn-1)。Exceptionsstd::invalid_argument输入值为 NaN 时抛出。◆calculateSecondDerivatives()std::vectordouble Lei_SplineInterpolator::calculateSecondDerivatives();描述计算样条曲线的核心引擎。组装三对角矩阵三弯矩方程使用追赶法求解所有节点处的二阶导数 Mig′′(xi)。若只有 2 个数据点则退化为线性插值所有二阶导数为 0。Returns包含n个元素的 std::vectordouble即 g′′(x0) 到 g′′(xn−1)。Exceptionsstd::runtime_error未设置有效数据点时抛出。std::runtime_error追赶法求解出现除零时抛出。◆interpolate()double Lei_SplineInterpolator::interpolate(double xi) const;描述单点插值求值。使用二分查找std::upper_bound定位区间利用局部 A、B、C、D 权重公式计算插值结果。Parameters[in] xi待插值的单点x坐标。Returns插值结果y。若 xi 超出原始数据范围返回NaN静默失败防御性编程。Exceptionsstd::runtime_error尚未调用 calculateSecondDerivatives 时抛出。◆interpolateArray()std::vectordouble Lei_SplineInterpolator::interpolateArray(const std::vectordouble xiArray) const;描述批量插值。循环调用单点求值函数。由于每次调用都会执行二分查找适用于待求点数量适中或乱序的情况。Parameters[in] xiArray待插值的 x 坐标数组。Returns对应的插值结果数组。◆ generateDenseData()void Lei_SplineInterpolator::generateDenseData(int factor,std::vectordouble outX,std::vectordouble outY,bool useSpline true) const;描述利用插值因子Factor生成等间距的密集网格用于重采样和曲线绘制。在每个原始数据区间内插入factor - 1个新点并在循环中采用流式递推非二分查找实现O(n) 的高效计算。Parameters[in] factor插值因子必须 ≥ 1。[out] outX输出的密集 X 坐标数组。[out] outY输出的密集 Y 坐标数组。[in] useSplinetrue 使用样条插值false 使用线性插值。Exceptionsstd::invalid_argument因子小于 1 或数据未设置时抛出。std::runtime_error未预先计算二阶导数且 useSpline 为 true 时抛出◆linearInterpolate()double Lei_SplineInterpolator::linearInterpolate(double xi) const;描述一维线性插值折线连接。作为样条的降级兜底方案。其特点是绝对无“过冲Overshoot”现象。在越界处返回端点值Clamping而非 NaN这是工业现场的保守逻辑。Parameters[in] xi待插值的单点 x 坐标。Returns线性插值结果 y4.应用注意事项4.1边界条件的物理意义选择极其重要默认的 1e30 代表自然样条自由端。如果你的实际工况是悬臂梁或已知端点的确切斜率必须手动填入端点的一阶导数切换为钳制样条否则强制使用自然样条会导致端点附近出现严重的几何偏离。4.2严禁外推Extrapolation代码中 interpolate 在超出原始范围时返回 NaN。这是极其正确的防御机制。三次多项式在数据范围外的行为呈指数级发散机床若使用外推数据进行加工必定发生撞机事故。请务必在调用侧检查 NaN。4.3浮点精度与病态矩阵如果数据点之间的步长hi相差极大例如一个区间间距是 0.001另一个是 1000三对角矩阵可能会接近病态。代码中已加入1e-12的步长校验但在实际工程中如果点云间距差异过大建议先进行“弧长参数化”或对 X 轴进行归一化处理。4.4性能优化策略对于单点查询如鼠标交互使用 interpolate对于批量且已排序的查询建议结合双指针法将 interpolateArray 进一步优化至 O(n)绘制图表或生成NC代码请务必使用 generateDenseData。4.5线程安全性该类不是线程安全的。calculateSecondDerivatives 会修改内部状态且计算量大。在Qt中强烈建议将此类的构建与计算放在QThread或 QtConcurrent的子线程中执行通过信号槽将结果std::vector传回UI线程刷新避免界面卡顿。5.参考资料
目录1.1 核心原理1.2 解决的问题1.3 目的2. 样条的性质基于数值分析教材3. 代码、函数详细说明3.1 Lei_SplineInterpolator.h代码3.2 Lei_SplineInterpolator.cpp代码3.3 main函数测试代码3.4运行结果3.6 函数说明◆ Lei_SplineInterpolator() / ~Lei_SplineInterpolator()◆ setPoints()◆ setBoundaryConditions()◆ calculateSecondDerivatives()◆ interpolate()◆ interpolateArray()◆ generateDenseData()◆ linearInterpolate()4.应用注意事项4.1边界条件的物理意义选择极其重要4.2严禁外推Extrapolation4.3浮点精度与病态矩阵4.4性能优化策略4.5线程安全性5.参考资料1.1 核心原理三次样条插值Cubic Spline Interpolation是一种分段多项式插值方法。它的数学本质是在给n个数据点之间构造n−1段三次多项式:。为了确保曲线平滑三次样条要求在所有内部节点处满足函数值连续、一阶导数连续、二阶导数连续。在实际求解时为了降低计算复杂度通常引入节点处的二阶导数作为未知量。通过推导最终转化为求解一个严格对角占优的三对角线性方程组三弯矩方程并利用追赶法Thomas Algorithm在O(n)的时间复杂度内解出所有从而确定整条曲线。1.2 解决的问题高次多项式插值的震荡问题龙格现象使用单一高次多项式如拉格朗日插值拟合大量数据点时边缘会出现剧烈的震荡导致结果完全不可用。三次样条通过分段低次多项式完美规避了这一问题。线性插值的折线问题线性插值直线连接虽然稳定但在节点处一阶导数不连续有尖角二阶导数为0。这在物理意义上意味着速度和加速度存在突变用于曲面重构或加工轨迹时会导致机床剧烈震动。1.3 目的在离散的工程测量点之间构建一条不仅穿过所有已知点而且整体极度平滑曲率连续的参数曲线。它是点云预处理、曲面重构、数控加工轨迹光顺的数学基石。2. 样条的性质基于数值分析教材在教材中三次样条函数S(x)通常具备以下四条核心性质性质1插值性曲线严格通过给定的所有数据点即。性质2一阶导数连续性在相邻两段样条的连接点处一阶导数斜率是连续的即。这确保了曲线没有折角。性质3二阶导数连续性在相邻两段样条的连接点处二阶导数曲率也是连续的即。这确保了曲线曲率平滑过渡。性质4边界条件端点条件为了唯一定义样条曲线必须补充两个额外条件。工程中常用的有4a. 自然样条Natural Spline假设两端点不受弯矩S′′(x0)0且 S′′(xn−1)0。适用于毫无端点物理约束的自由曲线。4c. 钳制样条Clamped Spline给定端点的一阶导数斜率S′(x0)v0且 S′(xn−1)vn−1。适用于已知速度、转角等物理约束的场景。3. 代码、函数详细说明3.1 Lei_SplineInterpolator.h代码/** * file Lei_SplineInterpolator.h * brief 基于数值分析的三次样条插值器。 * details 支持自然样条与钳制样条包含数据排序、边界条件处理及安全防御机制。 * author Lei * date 2026-10-11 */ #ifndef LEI_SPLINE_INTERPOLATOR_H #define LEI_SPLINE_INTERPOLATOR_H #include vector #include cmath #include stdexcept #include iostream #include iomanip #include limits #include algorithm /** * class Lei_SplineInterpolator * brief 三次样条插值器 * * 该类提供了一套完整的三次样条插值解决方案。它在内部对输入数据点进行排序 * 基于三对角矩阵和追赶法Thomas Algorithm求解二阶导数并支持自然样条 * 和钳制样条两种边界条件。 */ class Lei_SplineInterpolator { public: /** * brief 构造函数 */ Lei_SplineInterpolator() default; /** * brief 设置原始数据点 (X, Y) * details 内部会对 X 数据进行自动排序。排序后如果发现 X 存在重复值将抛出异常。 * param x 自变量数组 * param y 因变量数组 * throws std::invalid_argument 如果 X 和 Y 长度不一致或数据点少于 2 个或存在重复的 X 值 */ void setPoints(const std::vectordouble x, const std::vectordouble y); /** * brief 设置边界条件 * param initialBoundary 初始边界的一阶导数默认 1.0E30表示自然样条 * param finalBoundary 最终边界的一阶导数默认 1.0E30表示自然样条 * note 若输入 1.0E30内部自动切换为自然样条边界二阶导数为 0 * throws std::invalid_argument 如果输入为 NaN */ void setBoundaryConditions(double initialBoundary 1.0E30, double finalBoundary 1.0E30); /** * brief 计算二阶导数数组 g(x_i) * details 组装三对角矩阵并使用追赶法求解。这是样条插值的核心步骤。 * return 包含 n 个元素的双精度数组即 g(x_0) 到 g(x_{n-1}) * throws std::runtime_error 如果数据点尚未设置 */ std::vectordouble calculateSecondDerivatives(); /** * brief 单点插值 * param xi 待插值的 x 坐标 * return 插值结果 y 坐标 * throws std::runtime_error 如果尚未调用 calculateSecondDerivatives */ double interpolate(double xi) const; /** * brief 批量插值 * param xiArray 待插值的 x 坐标数组 * return 插值结果数组 * throws std::runtime_error 如果尚未调用 calculateSecondDerivatives */ std::vectordouble interpolateArray(const std::vectordouble xiArray) const; /** * brief 一维数组线性插值 * param xi 待插值的 x 坐标 * return 线性插值结果 y 坐标。若 xi 超出范围则返回边界值。 * throws std::runtime_error 如果数据点尚未设置 */ double linearInterpolate(double xi) const; /** * brief 根据插值因子生成更密集的数据点 (对应数组插值因子功能) * param factor 插值因子必须大于等于 1。表示在每个原始区间内插入的细分数。 * param outX 输出的密集 X 数组 * param outY 输出的密集 Y 数组 (按样条或线性插值计算得出) * param useSpline 是否使用三次样条插值 (true)否则使用线性插值 (false) * throws std::invalid_argument 如果插值因子小于 1或数据尚未设置 * throws std::runtime_error 如果未预先计算二阶导数且 useSpline 为 true */ void generateDenseData(int factor, std::vectordouble outX, std::vectordouble outY, bool useSpline true) const; private: std::vectordouble m_x; /// 经过排序的原始 X 数组 std::vectordouble m_y; /// 对应排序后的 Y 数组 std::vectordouble m_M; /// 二阶导数数组 g(x_i) double m_initialBoundary; /// 初始边界条件一阶导数 double m_finalBoundary; /// 最终边界条件一阶导数 bool m_isNaturalSpline; /// 是否为自然样条 bool m_isDataValid; /// 数据是否有效 bool m_isBuilt; /// 是否已计算二阶导数 /** * brief 求解三对角线性方程组 (追赶法 / Thomas Algorithm) * param a 下对角线 * param b 主对角线 * param c 上对角线 * param d 右端向量 * return 解向量 */ std::vectordouble solveThomas(const std::vectordouble a, const std::vectordouble b, const std::vectordouble c, const std::vectordouble d) const; }; #endif // LEI_SPLINE_INTERPOLATOR_H3.2 Lei_SplineInterpolator.cpp代码/** * file Lei_SplineInterpolator.cpp * brief Lei_SplineInterpolator 类的实现 */ #include Lei_SplineInterpolator.h void Lei_SplineInterpolator::setPoints(const std::vectordouble x, const std::vectordouble y) { // 1. 防御性检查长度一致且不少于2个点 if (x.size() ! y.size()) { throw std::invalid_argument([防御性编程] X 和 Y 数组的长度必须一致。); } if (x.size() 2) { throw std::invalid_argument([防御性编程] 至少需要 2 个数据点才能进行插值。); } // 2. 处理乱序数据将 X 和 Y 绑定后排序 std::vectorstd::pairdouble, double points; points.reserve(x.size()); for (size_t i 0; i x.size(); i) { // 检查是否为 NaN if (std::isnan(x[i]) || std::isnan(y[i])) { throw std::invalid_argument([防御性编程] 输入数据包含 NaN。); } points.emplace_back(x[i], y[i]); } // 按 X 坐标进行升序排序 std::sort(points.begin(), points.end(), [](const std::pairdouble, double a, const std::pairdouble, double b) { return a.first b.first; }); // 3. 排序后检查 X 是否严格递增是否存在重复的 X 值 m_x.resize(points.size()); m_y.resize(points.size()); for (size_t i 0; i points.size(); i) { m_x[i] points[i].first; m_y[i] points[i].second; if (i 0 std::abs(m_x[i] - m_x[i - 1]) 1e-12) { throw std::invalid_argument([防御性编程] 排序后发现 X 坐标存在重复值无法进行样条插值。); } } m_isDataValid true;//数据有效标志 m_isBuilt false;//是否计算二阶导数标志 } void Lei_SplineInterpolator::setBoundaryConditions(double initialBoundary, double finalBoundary) { // 防御性检查边界值不能是 NaN if (std::isnan(initialBoundary) || std::isnan(finalBoundary)) { throw std::invalid_argument([防御性编程] 边界条件不能为 NaN。); } // 工程约定大于等于 1.0E30 视为无穷大即自然样条二阶导数为 0 const double INFINITY_BOUNDARY_FLAG 1.0E30; m_isNaturalSpline false; if (initialBoundary INFINITY_BOUNDARY_FLAG finalBoundary INFINITY_BOUNDARY_FLAG) { m_isNaturalSpline true;//自然样条 } m_initialBoundary initialBoundary;//初始边界条件 m_finalBoundary finalBoundary;//最终边界条件 m_isBuilt false; // 边界条件改变需要重新计算 } std::vectordouble Lei_SplineInterpolator::calculateSecondDerivatives() { if (!m_isDataValid) { throw std::runtime_error([防御性编程] 数据点尚未设置或无效无法计算二阶导数。); } if (m_isBuilt) return m_M;//已经有了二阶导数 size_t n m_x.size(); m_M.resize(n, 0.0); // 只有两个点退化为线性插值二阶导数均为 0 if (n 2) { m_M[0] 0.0; m_M[1] 0.0; m_isBuilt true; return m_M; } // 计算步长 h std::vectordouble h(n - 1); for (size_t i 0; i n - 1; i) { h[i] m_x[i 1] - m_x[i];//δi } // 三对角矩阵的系数 std::vectordouble a(n, 0.0); // 下对角线 std::vectordouble b(n, 0.0); // 主对角线 std::vectordouble c(n, 0.0); // 上对角线 std::vectordouble d(n, 0.0); // 右端向量 // 1. 组装内点方程 (i 1 to n-2) for (size_t i 1; i n - 1; i) { a[i] h[i - 1]; b[i] 2.0 * (h[i - 1] h[i]); c[i] h[i]; d[i] 6.0 * ((m_y[i 1] - m_y[i]) / h[i] - (m_y[i] - m_y[i - 1]) / h[i - 1]); } // 2. 处理边界条件 if (m_isNaturalSpline) { // 自然样条M_0 0, M_{n-1} 0 b[0] 1.0; c[0] 0.0; d[0] 0.0; a[n - 1] 0.0; b[n - 1] 1.0; d[n - 1] 0.0; } else { // 钳制样条给定一阶导数 b[0] 2.0 * h[0]; c[0] h[0]; d[0] 6.0 * ((m_y[1] - m_y[0]) / h[0] - m_initialBoundary); a[n - 1] h[n - 2]; b[n - 1] 2.0 * h[n - 2]; d[n - 1] 6.0 * (m_finalBoundary - (m_y[n - 1] - m_y[n - 2]) / h[n - 2]); } // 3. 追赶法求解 m_M solveThomas(a, b, c, d); m_isBuilt true; return m_M; } double Lei_SplineInterpolator::interpolate(double xi) const { if (!m_isBuilt) { throw std::runtime_error([防御性编程] 尚未计算二阶导数请先调用 calculateSecondDerivatives()。); } if (std::isnan(xi)) { return std::numeric_limitsdouble::quiet_NaN(); } // 处理边界外的情况返回 NaN 或抛出异常这里选择返回 NaN 供上层决定 if (xi m_x.front() || xi m_x.back()) { return std::numeric_limitsdouble::quiet_NaN(); } size_t n m_x.size(); if (n 2) { // 两点退化为线性 double t (xi - m_x[0]) / (m_x[1] - m_x[0]); return m_y[0] t * (m_y[1] - m_y[0]); } // 找到 xi 所在的区间 [x_i, x_{i1}] auto it std::upper_bound(m_x.begin(), m_x.end(), xi);//二分查找必须是已经排好序的 size_t i std::distance(m_x.begin(), it) - 1; if (i n - 1) i n - 2; double h_i m_x[i 1] - m_x[i]; double A (m_x[i 1] - xi) / h_i; double B 1.0 - A; // 权重系数 C 和 D double C (1.0 / 6.0) * (A * A * A - A) * (h_i * h_i); double D (1.0 / 6.0) * (B * B * B - B) * (h_i * h_i); return A * m_y[i] B * m_y[i 1] C * m_M[i] D * m_M[i 1]; } std::vectordouble Lei_SplineInterpolator::interpolateArray(const std::vectordouble xiArray) const { if (!m_isBuilt) { throw std::runtime_error([防御性编程] 尚未计算二阶导数无法进行批量插值。); } std::vectordouble result; result.reserve(xiArray.size()); for (double xi : xiArray) { result.push_back(interpolate(xi)); } return result; } double Lei_SplineInterpolator::linearInterpolate(double xi) const { //实现一维数组的线性插值Linear Interpolation也就是把相邻的两个数据点用直线连接起来求中间任意位置的值。 if (!m_isDataValid || m_x.empty()) { throw std::runtime_error([防御性编程] 数据点尚未设置无法进行线性插值。); } if (std::isnan(xi)) { return std::numeric_limitsdouble::quiet_NaN(); } size_t n m_x.size(); if (xi m_x.front()) return m_y.front(); if (xi m_x.back()) return m_y.back(); auto it std::upper_bound(m_x.begin(), m_x.end(), xi);//二分查找必须是已经排好序的 size_t i std::distance(m_x.begin(), it) - 1; if (i n - 1) i n - 2; double h_i m_x[i 1] - m_x[i]; return m_y[i] (xi - m_x[i]) * (m_y[i 1] - m_y[i]) / h_i; } std::vectordouble Lei_SplineInterpolator::solveThomas(const std::vectordouble a, const std::vectordouble b, const std::vectordouble c, const std::vectordouble d) const { size_t n d.size(); std::vectordouble c_prime(n, 0.0); std::vectordouble d_prime(n, 0.0); std::vectordouble x(n, 0.0); // 防御性检查主对角线元素不能为 0一般情况下对角占优不会发生但工程上需要兜底 if (std::abs(b[0]) 1e-15) { throw std::runtime_error([防御性编程] 追赶法求解失败主对角线元素接近 0。); } // 向前消元 c_prime[0] c[0] / b[0]; d_prime[0] d[0] / b[0]; for (size_t i 1; i n; i) { double m b[i] - a[i] * c_prime[i - 1]; if (std::abs(m) 1e-15) { throw std::runtime_error([防御性编程] 追赶法求解失败出现除零错误。); } c_prime[i] (i n - 1) ? c[i] / m : 0.0; d_prime[i] (d[i] - a[i] * d_prime[i - 1]) / m; } // 向后回代 x[n - 1] d_prime[n - 1]; for (int i static_castint(n) - 2; i 0; --i) { x[i] d_prime[i] - c_prime[i] * x[i 1]; } return x; } void Lei_SplineInterpolator::generateDenseData(int factor, std::vectordouble outX, std::vectordouble outY, bool useSpline) const { // 防御性编程校验因子和状态 if (factor 1) { throw std::invalid_argument([防御性编程] 插值因子必须大于等于 1。); } if (!m_isDataValid || m_x.empty()) { throw std::runtime_error([防御性编程] 数据点尚未设置无法生成密集数据。); } if (useSpline !m_isBuilt) { throw std::runtime_error([防御性编程] 使用样条插值前请先调用 calculateSecondDerivatives()。); } size_t n m_x.size(); outX.clear(); outY.clear(); // 预留空间防止频繁内存分配 (假设原始区间数为 n-1) // 总点数 (n-1) * factor 1 (如果包含首尾点) outX.reserve((n - 1) * factor 1); outY.reserve((n - 1) * factor 1); // 如果是因子 1直接复制原始数据 if (factor 1) { outX m_x; outY m_y; return; } // 遍历每个区间 [x_i, x_{i1}] for (size_t i 0; i n - 1; i) { double step (m_x[i 1] - m_x[i]) / static_castdouble(factor); // 在当前区间内生成 factor 个点 (包含起点不包含终点) for (int j 0; j factor; j) { double xi m_x[i] j * step; outX.push_back(xi); if (useSpline) { outY.push_back(interpolate(xi)); } else { outY.push_back(linearInterpolate(xi)); } } } // 补上最后一个原始数据点 outX.push_back(m_x.back()); outY.push_back(m_y.back()); }3.3 main函数测试代码#include Lei_SplineInterpolator.h #include iostream #include vector #include iomanip int main() { // 故意打乱原始数据 X 的顺序验证内部排序功能 std::vectordouble X { 3.00, 0.00, 5.00, 1.20, 2.00, 3.50, 4.00 }; std::vectordouble Y { 150.00, 1.00, 220.00, 2.70, 7.40, 140.00, 80.00 }; try { Lei_SplineInterpolator spline; spline.setPoints(X, Y); spline.setBoundaryConditions(1.0E30, 1.0E30); // 自然样条 spline.calculateSecondDerivatives(); // // 测试 1指定单点插值 // std::cout 1. 单点插值测试: std::endl; std::cout std::fixed std::setprecision(2); std::cout X 1.36, Y spline.interpolate(1.36) (预期 -1.56) std::endl; // // 测试 2使用插值因子生成密集数组 // int factor 4; // 在每个原始数据点之间插入 4 个点 std::vectordouble denseX, denseY; std::cout \n2. 使用插值因子 factor 生成密集数据 (样条插值): std::endl; spline.generateDenseData(factor, denseX, denseY, true); std::cout 生成的密集数据点总数: denseX.size() std::endl; for (size_t i 0; i denseX.size(); i) { std::cout X denseX[i] , Y denseY[i] std::endl; } // 对比 LabVIEW 的示例插值数据 std::cout \n3. 与标准数据对比 (取特定 X 验证): std::endl; std::vectordouble checkX { 0.00, 0.24, 0.48, 0.72, 0.96, 1.20 }; for (double xi : checkX) { double yi spline.interpolate(xi); std::cout X xi , 插值 Y yi std::endl; } // // 测试 3防御性编程测试 // std::cout \n4. 防御性编程测试: std::endl; // 测试非法插值因子 try { spline.generateDenseData(0, denseX, denseY, true); } catch (const std::exception e) { std::cout 捕捉到异常 (因子为0): e.what() std::endl; } // 测试越界插值 double nan_val spline.interpolate(6.0); std::cout 越界插值 X6.0, 返回值是 NaN 吗? (std::isnan(nan_val) ? 是 : 否) std::endl; } catch (const std::exception e) { std::cerr 发生异常: e.what() std::endl; return -1; } return 0; }3.4运行结果3.6 函数说明◆ Lei_SplineInterpolator() / ~Lei_SplineInterpolator()Lei_SplineInterpolator() default;描述构造函数与析构函数。默认构造不分配额外内存内部状态初始为未设置数据m_isDataValid falsem_isBuilt false。◆ setPoints()void Lei_SplineInterpolator::setPoints(const std::vectordouble x,const std::vectordouble y);描述设置原始数据点并对输入数据进行排序。内部将 X 和 Y 绑定为 std::pair使用 std::sort 进行升序排序排序后严格检查是否存在重复的 X 坐标防除零。Parameters[in] x自变量数组顺序可不一致内部会自动排序。[in] y因变量数组必须与 X 长度一致。Exceptionsstd::invalid_argument当 X 与 Y 长度不一致、数据点少于 2 个、存在 NaN、或排序后 X 存在重复值时抛出。◆ setBoundaryConditions()void Lei_SplineInterpolator::setBoundaryConditions(double initialBoundary 1.0E30, double finalBoundary 1.0E30);描述设置样条的边界条件。默认值为 1.0E30代表无穷大当检测到两端点均大于等于该阈值时自动设置为自然样条性质4a否则根据给定的导数值切换为钳制样条性质4c。Parameters[in] initialBoundary初始边界的一阶导数 g(x0)。[in] finalBoundary最终边界的一阶导数 g(xn-1)。Exceptionsstd::invalid_argument输入值为 NaN 时抛出。◆calculateSecondDerivatives()std::vectordouble Lei_SplineInterpolator::calculateSecondDerivatives();描述计算样条曲线的核心引擎。组装三对角矩阵三弯矩方程使用追赶法求解所有节点处的二阶导数 Mig′′(xi)。若只有 2 个数据点则退化为线性插值所有二阶导数为 0。Returns包含n个元素的 std::vectordouble即 g′′(x0) 到 g′′(xn−1)。Exceptionsstd::runtime_error未设置有效数据点时抛出。std::runtime_error追赶法求解出现除零时抛出。◆interpolate()double Lei_SplineInterpolator::interpolate(double xi) const;描述单点插值求值。使用二分查找std::upper_bound定位区间利用局部 A、B、C、D 权重公式计算插值结果。Parameters[in] xi待插值的单点x坐标。Returns插值结果y。若 xi 超出原始数据范围返回NaN静默失败防御性编程。Exceptionsstd::runtime_error尚未调用 calculateSecondDerivatives 时抛出。◆interpolateArray()std::vectordouble Lei_SplineInterpolator::interpolateArray(const std::vectordouble xiArray) const;描述批量插值。循环调用单点求值函数。由于每次调用都会执行二分查找适用于待求点数量适中或乱序的情况。Parameters[in] xiArray待插值的 x 坐标数组。Returns对应的插值结果数组。◆ generateDenseData()void Lei_SplineInterpolator::generateDenseData(int factor,std::vectordouble outX,std::vectordouble outY,bool useSpline true) const;描述利用插值因子Factor生成等间距的密集网格用于重采样和曲线绘制。在每个原始数据区间内插入factor - 1个新点并在循环中采用流式递推非二分查找实现O(n) 的高效计算。Parameters[in] factor插值因子必须 ≥ 1。[out] outX输出的密集 X 坐标数组。[out] outY输出的密集 Y 坐标数组。[in] useSplinetrue 使用样条插值false 使用线性插值。Exceptionsstd::invalid_argument因子小于 1 或数据未设置时抛出。std::runtime_error未预先计算二阶导数且 useSpline 为 true 时抛出◆linearInterpolate()double Lei_SplineInterpolator::linearInterpolate(double xi) const;描述一维线性插值折线连接。作为样条的降级兜底方案。其特点是绝对无“过冲Overshoot”现象。在越界处返回端点值Clamping而非 NaN这是工业现场的保守逻辑。Parameters[in] xi待插值的单点 x 坐标。Returns线性插值结果 y4.应用注意事项4.1边界条件的物理意义选择极其重要默认的 1e30 代表自然样条自由端。如果你的实际工况是悬臂梁或已知端点的确切斜率必须手动填入端点的一阶导数切换为钳制样条否则强制使用自然样条会导致端点附近出现严重的几何偏离。4.2严禁外推Extrapolation代码中 interpolate 在超出原始范围时返回 NaN。这是极其正确的防御机制。三次多项式在数据范围外的行为呈指数级发散机床若使用外推数据进行加工必定发生撞机事故。请务必在调用侧检查 NaN。4.3浮点精度与病态矩阵如果数据点之间的步长hi相差极大例如一个区间间距是 0.001另一个是 1000三对角矩阵可能会接近病态。代码中已加入1e-12的步长校验但在实际工程中如果点云间距差异过大建议先进行“弧长参数化”或对 X 轴进行归一化处理。4.4性能优化策略对于单点查询如鼠标交互使用 interpolate对于批量且已排序的查询建议结合双指针法将 interpolateArray 进一步优化至 O(n)绘制图表或生成NC代码请务必使用 generateDenseData。4.5线程安全性该类不是线程安全的。calculateSecondDerivatives 会修改内部状态且计算量大。在Qt中强烈建议将此类的构建与计算放在QThread或 QtConcurrent的子线程中执行通过信号槽将结果std::vector传回UI线程刷新避免界面卡顿。5.参考资料
