此项目旨在生成常见的样条曲线,项目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, 代码针对二维数据进行了测试,三维理论上是支持的,可能需要修改一些东西,如果有需求,可在此代码基础上进行修改