B样条曲线介绍和实现(等值线平滑)

最近得到个任务是把之前的等值线光滑一下，思考良久决定用B样条去做平滑。虽然之前经常看到这个词条，但一直有正面接触，就趁着这次撸下来，做个总结吧。

B样条曲线

B样条是通过逼近一组控制点生成的曲线，它有如下计算表示：

$P(u)=\sum_{k=0}^{n}p_{k}N_{k,d}(u)$

其中 $p_{k}$ 是输入的一组n+1个控制点中第k个控制点。 $N_{k,d}$ 为B样条混合函数，是一种调和函数。 $N_{k,d}$ 中的k表示第k个混合函数，有多少个控制点就有多少个混合函数。 $N_{k,d}$ 中的d表示次数，如三次B样条曲线d为3。有的资料会用阶数表示d(阶数=次数+1)。我们这里就用次数表示d。d可以赋值为1到n之间的任意整数。虽然也可以设置d=0，但是这样就是折线本身了。B样条的混合函数 $N_{k,d}$ 由Cox-deBoor递归公式定义：

$N_{k,0}=\left\{\begin{matrix} 1, & u_{k}\leq u\leq u_{k+1}\\ 0,& otherwise \end{matrix}\right$

$N_{k,d}=\frac{u-u_{k}}{u_{k+d}-u_{k}}N_{k,d-1}(u)+\frac{u_{k+d+1}-u}{u_{k+d+1}-u_{k+1}}N_{k+1,d-1}(u)$

其中 $u_{k}$ 是节点向量 $U=\left \{ u_{0},u_{1},...,u_{n+d+1}\right \}$ 的索引为k的(总共n+d+2个节点)节点，均匀节点区间通常使用[0,1]或者[0，n+d+1]。

节点向量

理解B样条曲线最重要就是理解什么是结点向量。区别于控制点，控制点表示用n+1个点去控制曲线，而节点表示将曲线划分为n+d+1段，如下图所示：

直线上的6个点投影到曲线上对应的6个点就是节点向量U。我们在实践中往往U是均匀分布的即 $u_{k+1}-u_{k}=c$ ，c为常数。下图展示了一个常见的2次B样条曲线：

折线上的点就是我们的控制点，曲线上红色的点即为我们的节点 $u_{k}$ 。这里会发现一个奇怪的现象，明明节点数为n+d+2，控制点为n+1，为什么图中节点数少于控制点。这时由于混合函数的性质。

混合函数性质

我们把混合函数的计算进行排列：

每个混合函数 $N_{k,d}(u)$ 都定义在 $u_{k}$ 为起点，取值范围为d+1的子区间上。

比如 $N_{1,3}$ ，其定义域为 $[u_{1},u_{5})$ ，因此其他情况下 $N_{1,3}$ 取值为0。这表明了在k=1，即第一个控制点影响范围为 $[u_{1},u_{5})$ 。如下图所示：

每个样条曲线段(两个相邻节点间)受d+1个控制点影响

比如我们取节点范围为 $[u_{3},u_{4})$ ，根据混合函数，发现它影响 $N_{0,3}$ ， $N_{1,3}$ ， $N_{2,3}$ ， $N_{3,3}$ ，表明其只受第0到第3个控制点的影响。如下图所示。

我们发现从 $[u_{0},u_{3})$ 的线段无法满足4个影响点，同样的 $[u_{n+1},u_{n+d+1})$ 范围也无法取到4个点，因此这些范围被定义为无效范围。即：

节点记为 $\left \{ u_{0},u_{1},...,u_{n+d+1}\right \}$ ，所生成的B样条线曲线仅定义在节点范围 $[u_{d},u_{n+1})$ 上。

B样条代码实现

//边界情况(闭合或者相切)
enum class B_Spline_Border{TANGENT,CLOSED};class B_SPLINE{
public:B_SPLINE(){}B_SPLINE(int d,const std::vector<float>& p);//获取P(u)void getPmiu(float u,float& px,float& py);//计算基函数递归式float Cox(float u,int i,int p);//获取B样条线序列void getUniPos(int numU,std::vector<float>& p);
private://次数(阶数-1)int m_D;//结点向量(uk)std::vector<float> m_VecKnots;//控制点std::vector<float> m_Pos;//控制段数n(控制点=n+1)int m_N;//边界情况B_Spline_Border m_Bor;
};B_SPLINE::B_SPLINE(int d,const std::vector<float>& p){m_D=d;m_N=int(p.size()/2)-1;if(p[0]==p[2*m_N]&&p[1]==p[2*m_N+1]){//闭合情况m_Bor=B_Spline_Border::CLOSED;for (int i=0; i<(m_D-1)/2; i++) {//首控制点前添加控制点m_Pos.push_back(p[2*(m_N+i-(m_D-1)/2)]);m_Pos.push_back(p[2*(m_N+i-(m_D-1)/2)+1]);}for (int i=0; i<p.size(); i++) {m_Pos.push_back(p[i]);}for (int i=0; i<m_D/2; i++) {//尾控制点后添加控制点m_Pos.push_back(p[2*(i+1)]);m_Pos.push_back(p[2*(i+1)+1]);}}else{//非闭合采用首尾结点重复度为d+1的方法生成切线m_Bor=B_Spline_Border::TANGENT;for (int i=0; i<m_D; i++) {m_Pos.push_back(p[0]);m_Pos.push_back(p[1]);}for (int i=0; i<p.size(); i++) {m_Pos.push_back(p[i]);}for (int i=0; i<m_D; i++) {m_Pos.push_back(p[2*m_N]);m_Pos.push_back(p[2*m_N+1]);}}m_N=int(m_Pos.size()/2)-1;
}void B_SPLINE::getPmiu(float u,float& px,float& py){px=0;py=0;for (int i=0; i<m_N+1; i++) {px+=m_Pos[2*i]*Cox(u, i, m_D);py+=m_Pos[2*i+1]*Cox(u, i, m_D);}}float B_SPLINE::Cox(float u,int i,int p){if(p==0){if (u>=m_VecKnots[i]&&u<m_VecKnots[i+1]) {return 1;}else{return 0;}}else{float a,b;float dev=(m_VecKnots[i+p]-m_VecKnots[i]);a=dev?(u-m_VecKnots[i])/dev:0;dev=(m_VecKnots[i+p+1]-m_VecKnots[i+1]);b=dev?(m_VecKnots[i+p+1]-u)/dev:0;return a*Cox(u,i,p-1)+b*Cox(u, i+1, p-1);}
}void B_SPLINE::getUniPos(int numU,std::vector<float>& p){//均匀样条m_VecKnots.clear();for (int i=0; i<=m_N+m_D+1; i++) {m_VecKnots.push_back(i);}float x,y;p.clear();for(int i=0;i<=numU;i++){float u=float(i)/numU*(m_VecKnots.size()-1);if (u>=m_VecKnots[m_D]&&u<m_VecKnots[m_N+1]) {getPmiu(u,x,y);p.push_back(x);p.push_back(y);}}//由于右区间是开区间(un+1这个点未定义)，所以最后要闭合得手动连接首尾点(但如果细分够细，也可以选择忽略)if (m_Bor==B_Spline_Border::CLOSED) {p.push_back(p[0]);p.push_back(p[1]);}}

其中Cox是我们的混合函数递归式。这里还可以直接根据解析式写出代码，避免递归造成的时间浪费，但是根据定义来撸代码也是挺爽的。由于是用于等值线的绘制，我们这里使用了两种边界处理方法，开放均匀B样条和闭合B均匀样条。

均匀周期性B样条和开放均匀B样条

均匀周期性B样条表示我们的节点向量是均匀分布的。而周期性表示，当我们的节点向量是均匀分布时，所有的混合函数都是可以通过平移得到的。如下图所示：

开放均匀B样条(也称开放B样条)表示我们可以只改变首尾为d+1次重复，而中间的其他节点向量仍然是均匀的。比如节点向量{0,0,0,0,1,2,3,4,5,6,7,8,9,9,9,9}。我们发现该结点向量首尾各重复了4次(也可称为重复度4)，这种情况生成的B样条即为开放均匀的。如下图的三次B样条所示：

上左图为均匀周期B样条曲线，右图为开放均匀B样条。我们发现右图可以使曲线和首尾两个控制点相切，这里只要保证首尾的重复度为d+1生成的B样条曲线就是开放B样条曲线，开放B样条曲线也可以称为Clamped-B样条曲线。

（PS：许多网上的资料都说开放B样条为不闭合的B样条曲线，个人还是不太赞同的。个人比较认同《计算机图形学》这本书中对开放B样条的定义，即首尾重复d+1次）

闭合B样条曲线

我们如果想获得如下闭合的样条曲线应该怎么做呢：

非常常用的办法为包裹(wrapping)控制点。

包裹控制点：

(1)设计一个均匀 m+1 (m=d+n+1)个节点的节点序列：u0 = 0, u1 = 1/m, u1 = 2/m, ..., um = 1。注意曲线的定义域是 $[u_{d},u_{n+1})$ .

(2)包裹头d个和最后d个控制点。更准确地，设P0 = Pn-d+1, P1 = Pn-d+2, ..., Pd-2 = Pn-1 and Pd-1 = Pn.

如果我们输入的d个控制点首尾是相等的，即连接的情况下：如果d为奇数，则首控制点之前和尾控制点之后各添加int(d/2)个控制点。如果d为偶数则在首控制点之前添加int((d-1)/2)个控制点，尾控制点之后添加int(d/2)个控制点，如上述代码所示：

    if(p[0]==p[2*m_N]&&p[1]==p[2*m_N+1]){//闭合情况m_Bor=B_Spline_Border::CLOSED;for (int i=0; i<(m_D-1)/2; i++) {//首控制点前添加控制点m_Pos.push_back(p[2*(m_N+i-(m_D-1)/2)]);m_Pos.push_back(p[2*(m_N+i-(m_D-1)/2)+1]);}for (int i=0; i<p.size(); i++) {m_Pos.push_back(p[i]);}for (int i=0; i<m_D/2; i++) {//尾控制点后添加控制点m_Pos.push_back(p[2*(i+1)]);m_Pos.push_back(p[2*(i+1)+1]);}}

如下图所示

根据上述代码绘制的两种曲线效果如下：

左图为闭合效果，右图为相切效果。