此项目旨在生成常见的样条曲线，项目inspired by 叶飞影
样条曲线：转述知乎答案解释其意义。
预备知识：已知离散的数据，但不知函数表达式，插值和拟合都是为了寻找函数表达式。区别在于，插值得到的函数能够穿过已知的点（在已知的点的函数表达式的值等于已知数值，但容易出现龙格现象），拟合只求函数图形神似而不求穿过已知点。
那么怎么能既穿过已知点又能让函数图形像呢？就是怎么避免龙格现象呢？
答案是分段插值，就是将全部数据分割成若干部分，每个小部分用插值得到不同的函数，最后用很多不同的函数表达原来的序列。问题又来了，不同函数两端衔接不好怎么办？
答案是高次样条差值，既每个分段函数都采用高次函数形式来构造（三次样条差值 就是用x的三次方形式构造）这就保证了得到的多个函数关系式在先接触具有n-1次的连续可导性质（翻译成人话就是衔接保证光滑）
一句话总结：三次样条插值就是将原始长序列分割成若干段构造多个三次函数（每段一个），使得分段的衔接处具有二阶导数连续的性质（也就是光滑衔接）。
其中“三次”只函数基本形式使用三次函数的形式。“样条”是一种手艺，指加工曲面时使得曲面光滑的手艺。“插值”你肯定知道是啥意思了~~
此处答案虽然针对三次样条插值，但是可以很好解释样条插值的概念。
首先是各种曲线的基类
#ifndef __SPLINE_BASE__ #define __SPLINE_BASE__ #include <cstddef> class SplineBase { public: // 计算样条数值 virtual bool BuildSpline(const void* ctrlValuesPtr, int ctrlStride, int ctrlCount, void* splineValuesPtr, int splineStride) const = NULL; // 计算切线数值 virtual bool BuildTangent(const void* ctrlValuesPtr, int ctrlStride, int ctrlCount, void* tangentValuesPtr, int tangentStride) const { return false; } }; inline float YfGetFloatValue(const void* valuesPtr, int stride, int index) { return *(float*)((char*)valuesPtr + stride * index); } #endif //__SPLINE_BASE__
每种曲线的类都会继承此基类，并实现 BuildSpline或者同时实现BuildSpline和Buildtangent两个函数。三种曲线中，三阶及以上阶数的Bezier曲线只通过起始控制点和终点，B样条曲线不会经过其控制点，CatmullRoom曲线会通过所有控制点。BSpline和CatmullRoom样条公式要求4个控制点，然后给出位于第2和第3点之间的光滑曲线(效果如下图，图片来源游戏编程精粹1)。为了得到第1个控制点和第2个控制点之间的点，可以将第一个控制点输入两次，同理第3和第4个控制点之间的曲线需要第四点输入两次。
N阶Bezier 曲线
1. 原理参考我的另外一篇博客：N 阶贝塞尔曲线生成，此篇代码每次都要计算N次方，比较耗时，可以在绘制图形之前，先计算好各阶数值即可。如下代码分别是bezier.h 和 bezier.cpp实现N阶贝塞尔曲线
/**************************************************************** File name : bezier.h Author : Shike Version : 1.0 Create Date : 2020/05/01 Description : Bezier Spline *****************************************************************/ #ifndef __BEZIER__SPLINE__ #define __BEZIER__SPLINE__ #include "spline_base.h" #define YD_MAX_BEZIER_CONTROL_VALUE 33 class Bezier : public SplineBase { public: Bezier(); ~Bezier(); // 设置输出样条值的数目 void SetSplineValuesCount(int count); // 获得输出样条值的数目 int GetSplineValuesCount() const; // 计算样条数值 bool BuildSpline(const void* ctrlValuesPtr, int ctrlStride, int ctrlCount, void* splineValuesPtr, int splineStride) const; protected: void ClearPowT(); void BuildPowT(); double GetValueT(int t, int p) const { return m_pow_t[YD_MAX_BEZIER_CONTROL_VALUE*t + p]; } protected: int m_valuesCount; double *m_pow_t; protected: static void BuildYanghuiTriangle(); //pascal's triangle static int m_yanghuiRowIndex[YD_MAX_BEZIER_CONTROL_VALUE]; static int m_yanghuiTriangle[(YD_MAX_BEZIER_CONTROL_VALUE + 1)*YD_MAX_BEZIER_CONTROL_VALUE / 2]; }; #endif // __BEZIER__SPLINE__
bezier.cpp
#include "bezier.h" #include "assert.h" #include "stdlib.h" #include <cstdio> int Bezier::m_yanghuiRowIndex[YD_MAX_BEZIER_CONTROL_VALUE] = { 0 }; int Bezier::m_yanghuiTriangle[(YD_MAX_BEZIER_CONTROL_VALUE + 1)*YD_MAX_BEZIER_CONTROL_VALUE / 2] = { 0 }; Bezier::Bezier() { if (m_yanghuiTriangle[0] == 0) { BuildYanghuiTriangle(); } m_valuesCount = 0; m_pow_t = NULL; SetSplineValuesCount(100); } Bezier::~Bezier() { ClearPowT(); } void Bezier::BuildYanghuiTriangle() { // 第0行 m_yanghuiRowIndex[0] = 0; m_yanghuiTriangle[0] = 1; int index = 1; int t0, t1; int* lastRow; for (int i = 1; i < YD_MAX_BEZIER_CONTROL_VALUE; i++) { m_yanghuiRowIndex[i] = index; m_yanghuiTriangle[index] = 1; index++; for (int j = 1; j <= i; j++) { lastRow = m_yanghuiTriangle + m_yanghuiRowIndex[i - 1]; t0 = lastRow[j - 1]; t1 = (j < i) ? lastRow[j] : 0; m_yanghuiTriangle[index] = t0 + t1; index++; } } assert(index == (YD_MAX_BEZIER_CONTROL_VALUE + 1)*YD_MAX_BEZIER_CONTROL_VALUE / 2); } //设置输出样条值的数目 void Bezier::SetSplineValuesCount(int count) { if (count < 2) { count = 2; } if (count == m_valuesCount) { return; } m_valuesCount = count; BuildPowT(); } //获得输出样条值的数目 int Bezier::GetSplineValuesCount() const { return m_valuesCount; } void Bezier::ClearPowT() { if (m_pow_t) { free(m_pow_t); m_pow_t = NULL; } } void Bezier::BuildPowT() { ClearPowT(); m_pow_t = (double*)malloc(m_valuesCount*YD_MAX_BEZIER_CONTROL_VALUE * sizeof(double)); double t; for (int i = 0; i < m_valuesCount; i++) { t = i / (m_valuesCount - 1.0f); m_pow_t[i*YD_MAX_BEZIER_CONTROL_VALUE] = 1.0f; for (int j = 1; j < YD_MAX_BEZIER_CONTROL_VALUE; j++) { m_pow_t[i*YD_MAX_BEZIER_CONTROL_VALUE + j] = m_pow_t[i*YD_MAX_BEZIER_CONTROL_VALUE + j - 1] * t; } } } //计算样条数值 bool Bezier::BuildSpline(const void* ctrlValuesPtr, int ctrlStride, int ctrlCount, void* splineValuesPtr, int splineStride) const { if (ctrlCount < 2 || ctrlCount > YD_MAX_BEZIER_CONTROL_VALUE) { return false; } double* destValue; double* srcValue; double v; const int* yanghuiRow = m_yanghuiTriangle + m_yanghuiRowIndex[ctrlCount - 1]; for (int i = 0; i < m_valuesCount; i++) { for (int dim = 0; dim < ctrlStride; dim++) { v = 0.0f; for (int j = 0; j < ctrlCount; j++) { srcValue = (double*)ctrlValuesPtr + ctrlStride * j + dim; v += yanghuiRow[j] * (*srcValue) * GetValueT(i, j) * GetValueT(m_valuesCount - 1 - i, ctrlCount - 1 - j); } destValue = (double*)splineValuesPtr + splineStride * i + dim; *destValue = v; } } return true; }
显示效果如下
BSpline
在数学的子学科数值分析里，B-样条是样条曲线一种特殊的表示形式。它是B-样条基曲线的线性组合。B-样条是贝兹（贝塞尔）曲线的一种一般化，可以进一步推广为非均匀有理B样条(NURBS)，使得我们能给更多一般的几何体建造精确的模型。
常数B样条
常数B样条是最简单的样条。只定义在一个节点距离上，而且不是节点的函数。它只是不同节点段（knot span）的标志函数（indicator function）。
线性B样条
线性B样条定义在两个相邻的节点段上，在节点连续但不可微。
三次B样条
一个片断上的B样条的表达式可以写作：
其中Si是第i个B样条片段，而P是一个控制点集，i和k是局部控制点索引。控制点的集合会是
的集合，其中
是比重，当它增加时曲线会被拉向控制点
，在减小时则把曲线远离该点。片段的整个集合m-2条曲线（
）由m+1个控制点（
）定义，作为t上的一个B样条可以定义为
其中i是控制点数，t是取节点值的全局参数。这个表达式把B样条表示为B样条基函数的线性组合，这也是这个名称的原因。有两类B样条-均匀和非均匀。非均匀B样条相邻控制点间的距离不一定要相等。一个一般的形式是区间随着插入控制点逐步变小到0。
bspline.cpp核心代码
void BSpline::BuildWeights() { ClearWeights(); m_splineWeights = (double **)malloc(4 * sizeof(double*)); m_tangentWeights = (double **)malloc(4 * sizeof(double*)); for (int i = 0; i < 4; i++) { m_splineWeights[i] = (double*)malloc((m_subD) * sizeof(double)); m_tangentWeights[i] = (double*)malloc((m_subD) * sizeof(double)); } double u, u_2, u_3; for (int i = 0; i < m_subD; i++) { u = (float)i / m_subD; u_2 = u * u; u_3 = u_2 * u; // 参见"游戏编程精粹1"P331 m_splineWeights[0][i] = (-1.0f*u_3 + 3.0f*u_2 - 3.0f*u + 1.0f) / 6.0f; m_splineWeights[1][i] = (3.0f*u_3 - 6.0f*u_2 + 0.0f*u + 4.0f) / 6.0f; m_splineWeights[2][i] = (-3.0f*u_3 + 3.0f*u_2 + 3.0f*u + 1.0f) / 6.0f; m_splineWeights[3][i] = (1.0f*u_3 + 0.0f*u_2 + 0.0f*u + 0.0f) / 6.0f; // 参见"游戏编程精粹1"P333 m_tangentWeights[0][i] = (-1.0f*u_2 + 2.0f*u - 1.0f)*0.5f; m_tangentWeights[1][i] = (3.0f*u_2 - 4.0f*u + 0.0f)*0.5f; m_tangentWeights[2][i] = (-3.0f*u_2 + 2.0f*u + 1.0f)*0.5f; m_tangentWeights[3][i] = (1.0f*u_2 + 0.0f*u + 0.0f)*0.5f; } }
显示效果如下
CatmullRoom
所谓样条曲线是指给定一组控制点而得到一条曲线，曲线的大致形状由这些点予以控制，一般可分为插值样条和逼近样条两种，插值样条通常用于数字化绘图或动画的设计，逼近样条一般用来构造物体的表面。CatmullRom样条与上一节所讲的B样条很相似，不同在于CatmullRom样条的曲线会经过其每一个控制点。
catmullroom.cpp核心代码入下
void CatmullRoom::BuildWeights() { ClearWeights(); m_splineWeights = (double **)malloc(4 * sizeof(double*)); m_tangentWeights = (double **)malloc(4 * sizeof(double*)); for (int i = 0; i < 4; i++) { m_splineWeights[i] = (double*)malloc((m_subD) * sizeof(double)); m_tangentWeights[i] = (double*)malloc((m_subD) * sizeof(double)); } double u, u_2, u_3; for (int i = 0; i < m_subD; i++) { u = (float)i / m_subD; u_2 = u * u; u_3 = u_2 * u; // 参见"游戏编程精粹1"P333 m_splineWeights[0][i] = (-1.0f*u_3 + 2.0f*u_2 - 1.0f*u + 0.0f)*0.5f; m_splineWeights[1][i] = (3.0f*u_3 - 5.0f*u_2 + 0.0f*u + 2.0f)*0.5f; m_splineWeights[2][i] = (-3.0f*u_3 + 4.0f*u_2 + 1.0f*u + 0.0f)*0.5f; m_splineWeights[3][i] = (1.0f*u_3 - 1.0f*u_2 + 0.0f*u + 0.0f)*0.5f; m_tangentWeights[0][i] = (-3.0f*u_2 + 4.0f*u - 1.0f)*0.5f; m_tangentWeights[1][i] = (9.0f*u_2 - 10.0f*u + 0.0f)*0.5f; m_tangentWeights[2][i] = (-9.0f*u_2 + 8.0f*u + 1.0f)*0.5f; m_tangentWeights[3][i] = (3.0f*u_2 - 2.0f*u + 0.0f)*0.5f; } }
博主叶飞影博客附带可样条视化工具，可以下载设置控制点检测曲线效果。此博客对应的代码地址 multiple_splines, 代码针对二维数据进行了测试，三维理论上是支持的，可能需要修改一些东西，如果有需求，可在此代码基础上进行修改