参数曲线（贝塞尔曲线、B样条曲线和NURBS曲线）

参数曲线

网页资料

相关概念

齐次坐标：齐次坐标到笛卡尔坐标的变换$(x,y,z,w)\Rightarrow (x/w,y/w,z/w)$唯一，反之则不唯一

有理曲线：齐次形式的参数曲线称为有理曲线

连续性

$C^k$连续：$\forall i\le k$,$f(b)$和$g(m)$的$i$阶导数在连接点处相等

几何连续性$G^k$连续有如下等价定义：

• 以弧长为参数时，$\forall i\le k$,$f(b)$和$g(m)$的$i$阶导数在连接点处相等；

• 存在两个参数化方式，使得$\forall i\le k$,$f(b)$和$g(m)$的$i$阶导数在连接点处相等；

注意：$C^2$连续一定曲率连续，但曲率连续不一定$C^2$连续

贝赛尔曲线

构建贝赛尔曲线

给定空间中$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_n$，由这些控制点定义的贝赛尔曲线为

$\mathbf{C}(u)={\sum }_{i=0}^n{B}_{n,i}(u){\mathbf{P}}_i$

其中系数定义如下

$B_{n,i}(u)=\frac{n!}{i!(n-i)!}{u}^i(1-u{)}^{n-i}=C_n^i{u}^i(1-u{)}^{n-i}$

• 贝塞尔曲线可视为对所有控制点的加权平均

• 线段${\mathbf{P}}_0{\mathbf{P}}_1,{\mathbf{P}}_1{\mathbf{P}}_2,\cdots ,{\mathbf{P}}_{n-1}{\mathbf{P}}_n$称为legs，按此顺序连接形成控制折线（control polyline）/控制多边形（control polygon）

• $B_{n,i}(u)$称为贝塞尔基函数（Bézier basis functions）或伯恩斯坦多项式（Bernstein polynomials）

贝塞尔曲线具有以下性质：

• $n$阶贝赛尔曲线由$n+1$个控制点定义

• 贝塞尔曲线经过${\mathbf{P}}_0$和${\mathbf{P}}_n$，且在这两点处与控制折线相切

• 非负性：所有基函数均非负

• 单位分解（partition of unity）：所有基函数之和为1

• 凸包性质：贝塞尔曲线完全位于给定控制点的凸包内

• 变分递减性质（variation diminishing property）：没有一条直线与贝塞尔曲线相交的次数多于与曲线控制折线相交的次数

• 仿射不变性（affine invariance）：对贝塞尔曲线的仿射变换等价于对控制点的仿射变换

如果$u\in [a,b]$，则可通过变换$\underline{u}=\frac{u-a}{b-a}\in [0,1]$

移动控制点

原控制点${\mathbf{P}}_k$变为${\mathbf{P}}_k+\mathbf{v}$后，新的贝塞尔曲线

$\begin{array}{rl}\mathbf{D}(u)& =\sum _{i=0}^{k-1}{B}_{n,i}(u){\mathbf{P}}_i+B_{n,k}(u)\left({\mathbf{P}}_k+\mathbf{v}\right)+\sum _{i=k+1}^n{B}_{n,i}(u){\mathbf{P}}_i\\ & =\sum _{i=0}^n{B}_{n,i}(u){\mathbf{P}}_i+B_{n,k}(u)\mathbf{v}\\ & =\mathbf{C}(u)+B_{n,k}(u)\mathbf{v}\end{array}$

更改单个控制点，贝塞尔曲线上所有的点都将移动到新位置，曲线形状会发生整体的变化

德卡斯特里奥算法（de Casteljau’s algorithm）

直接法（直接把$u$带入基函数，计算每个基函数与其对应控制点的乘积，最后将它们相加）计算贝塞尔曲线上点的坐标虽然速度快，但是数值不稳定

几何解释

定义第$i$次迭代得到的点的编号为$i0,\cdots ,i(n-i)$，则控制点（第0次迭代）编号为$00,\cdots ,0n$，步骤如下：

1. 先对点$00,\cdots ,0n$依次形成的线段进行操作，对于所有$j\in [0,n-1]$，在线段$0j\sim 0(j+1)$找一点，该点把该线段分成$u:(1-u)$的两截，即得到如图中的点$10,\cdots ,1(n-1)$

2. 对点$10,\cdots ,1(n-1)$依次形成的线段进行同样的操作，得到如图中的点$20,\cdots ,2(n-2)$

3. 该操作重复$n$次之后会产生点$n0$，即得到贝塞尔曲线上的点$\mathbf{C}(u)$

实际计算

计算步骤如下：

1. 把所有控制点$00,\cdots ,0n$排成如图最左边的一列，每相邻的两个控制点之间得到一个新点

2. 新点为左边两个点的加权和，左上方的点权重为$1-u$，左下方的点权重为$u$

3. 按此规则迭代至得到$n0$，即得到贝塞尔曲线上的点$\mathbf{C}(u)$

递归关系

${\mathbf{P}}_{i,j}=(1-u){\mathbf{P}}_{i-1,j}+u{\mathbf{P}}_{i-1,j+1}\{\begin{array}{l}i=1,2,\dots ,n\\ j=0,1,\dots ,n-i\end{array}$

直接按递归关系计算的函数为

上述函数会有大量的重复计算，如下图中的${\mathbf{P}}_{n-2,1}$，导致计算效率低下

特殊性质

如果把某一列中一组连续的点作为贝塞尔曲线的控制点，那利用de Casteljau算法计算时，曲线上的点就是等边三角形内所选点形成的边的对面的顶点，如下图所示

算法的正确性

1. 每个控制点对$n0$的贡献完全通过图示的平行四边形区域传播过去

2. 控制点经过平行四边形中的每条路径都对$n0$有独立的贡献

3. 考虑到一条路径传播的贡献只和起点与终点的相对位置有关，和路径形状无关，由$0i$到$n0$的每条路径的贡献均为$u^i(1-u{)}^{n-i}$

4. 平行四边形内路径的总数恰好为组合数$C_n^i=\frac{n!}{i!(n-i)!}$

5. 易得，所有路径的贡献总和满足贝塞尔曲线关系式

$\frac{n!}{i!(n-i)!}{u}^i(1-u{)}^{n-i}0i=\frac{n!}{i!(n-i)!}{u}^i(1-u{)}^{n-i}{\mathbf{P}}_i=B_{n,i}(u){\mathbf{P}}_i$

贝塞尔曲线的导数

由于控制点是与$u$无关的常数，因此计算导数${\mathbf{C}}^{\prime }(u)$可以简化为计算$B_{n,i}(u)$的导数，满足

$\frac{\mathrm{d}}{\mathrm{d}u}{B}_{n,i}(u)=B_{n,i}^{\prime }(u)=n\left(B_{n-1,i-1}(u)-B_{n-1,i}(u)\right)$

定义${\mathbf{Q}}_i=n\left({\mathbf{P}}_{i+1}-{\mathbf{P}}_i\right),\forall i\in [0,n-1]$，则${\mathbf{C}}^{\prime }(u)$满足

d d u C ( u ) = C ′ ( u ) = ∑ i = 0 n − 1 B n − 1 , i ( u ) { n ( P i + 1 − P i ) } = ∑ i = 0 n − 1 B n − 1 , i ( u ) Q i \frac{\mathrm{d}}{\mathrm{d} u} \mathbf{C}(u)=\mathbf{C}^{\prime}(u)=\sum_{i=0}^{n-1} B_{n-1, i}(u)\left\{n\left(\mathbf{P}_{i+1}-\mathbf{P}_i\right)\right\}=\sum_{i=0}^{n-1} B_{n-1, i}(u)\mathbf{Q}_i dud​C(u)=C′(u)=∑i=0n−1​Bn−1,i​(u){n(Pi+1​−Pi​)}=∑i=0n−1​Bn−1,i​(u)Qi​

${\mathbf{C}}^{\prime }(u)$也被称为原贝塞尔曲线的速度图（hodograph）

连接贝塞尔曲线

• 满足$G^1$连续：前一条曲线的第一条leg和后一条曲线的最后一条leg满足如下图的关系

• 满足$C^1$连续：前一条曲线在$u=1$处的切向量和后一条曲线在$u=0$处的切向量方向和大小都一样

${\mathbf{C}}^{\prime }(1)=m\left({\mathbf{P}}_m-{\mathbf{P}}_{m-1}\right)={\mathbf{D}}^{\prime }(0)=n\left({\mathbf{Q}}_1-{\mathbf{Q}}_0\right)$

导数与de Casteljau算法的关系

贝塞尔曲线的导数可写为

C ′ ( u ) = ∑ i = 0 n − 1 B n − 1 , i ( u ) { n ( P i + 1 − P i ) } = n [ ( ∑ i = 0 n − 1 B n , i ( u ) P i + 1 ) − ( ∑ i = 0 n − 1 B n , i ( u ) P i ) ] \begin{aligned} \mathbf{C}^{\prime}(u) & =\sum_{i=0}^{n-1} B_{n-1, i}(u)\left\{n\left(\mathbf{P}_{i+1}-\mathbf{P}_i\right)\right\} \\ & =n\left[\left(\sum_{i=0}^{n-1} B_{n, i}(u) \mathbf{P}_{i+1}\right)-\left(\sum_{i=0}^{n-1} B_{n, i}(u) \mathbf{P}_i\right)\right] \end{aligned} C′(u)​=i=0∑n−1​Bn−1,i​(u){n(Pi+1​−Pi​)}=n[(i=0∑n−1​Bn,i​(u)Pi+1​)−(i=0∑n−1​Bn,i​(u)Pi​)]​

定义

${\mathbf{C}}_1(u)={\sum }_{i=0}^{n-1}{B}_{n,i}(u){\mathbf{P}}_{i+1}{\mathbf{C}}_2(u)={\sum }_{i=0}^{n-1}{B}_{n,i}(u){\mathbf{P}}_i$

则贝塞尔曲线的导数可表示为

${\mathbf{C}}^{\prime }(u)=n\left({\mathbf{C}}_1(u)-{\mathbf{C}}_2(u)\right)$

因此，理论上，可以使用de Casteljau算法来计算${\mathbf{C}}_1(u)$和${\mathbf{C}}_2(u)$，再计算差值，最后乘以$n$ ，得到${\mathbf{C}}^{\prime }(u)$。故利用de Casteljau算法可以同时计算出曲线$\mathbf{C}(u)$上的点及其切向量${\mathbf{C}}^{\prime }(u)$。下图显示了计算过程。

高阶导数

定义${\mathbf{D}}_i^0={\mathbf{P}}_i,\forall i\in [0,n]$，类似一阶导数的递推关系，贝塞尔曲线$k$阶导数的控制点满足有限差分关系

${\mathbf{D}}_i^k={\mathbf{D}}_{i+1}^{k-1}-{\mathbf{D}}_i^{k-1}0\le i\le n-k$

则贝塞尔曲线的$k$阶导数可表示为

$\begin{array}{rl}{\mathbf{C}}^{[k]}(u)& =n(n-1)(n-2)\cdots (n-k+1)\sum _{i=0}^{n-k}{B}_{n-k,i}(u)\left({\mathbf{D}}_{i+1}^{k-1}-{\mathbf{D}}_i^{k-1}\right)\\ & =n(n-1)(n-2)\cdots (n-k+1)\sum _{i=0}^{n-k}{B}_{n-k,i}(u){\mathbf{D}}_i^k\end{array}$

利用de Casteljau算法计算特定$u$下的$k$阶导数值的流程如下：

1. $k$阶导数的控制点满足如下左图的关系，根据如下右图的思路计算第$k$列的点${\mathbf{D}}_i^k={\sum }_{j=0}^k(-1{)}^{k-j}{C}_k^j{\mathbf{P}}_{i+j}$

2. 把de Casteljau算法应用在第$k$列的点上得到导数值

贝塞尔曲线的分段

贝塞尔曲线的分段需要满足以下条件：

• 将给定的贝塞尔曲线$\mathbf{C}(u)$在某一点切割成$[0,u]$和$[u,1]$两段贝塞尔曲线，且有各自的新控制点

• 新生成的贝塞尔曲线必须和原贝塞尔曲线同阶

贝塞尔曲线分段的新控制点可借助de Casteljau算法得到，以6阶贝塞尔曲线为例：

• 左侧折线由点${\mathbf{P}}_{00}={\mathbf{P}}_0,{\mathbf{P}}_{10},{\mathbf{P}}_{20},{\mathbf{P}}_{30},{\mathbf{P}}_{40},{\mathbf{P}}_{50},{\mathbf{P}}_{60}=\mathbf{C}(u)$组成

• 右侧折线由点${\mathbf{P}}_{60}=\mathbf{C}(u),{\mathbf{P}}_{51},{\mathbf{P}}_{42},{\mathbf{P}}_{33},{\mathbf{P}}_{24},{\mathbf{P}}_{15},{\mathbf{P}}_{06}={\mathbf{P}}_6$组成

• 以上两组控制点也对应如下de Casteljau算法中红圈内顺箭头方向的一组点和蓝圈内逆箭头方向的一组点

算法的正确性

以原贝塞尔曲线$[0,u]$的部分为例，证明步骤如下：

1. 设原曲线为

$\mathbf{C}(u)={\sum }_{i=0}^n{B}_{n,i}(u){\mathbf{P}}_i$

2. 由de Casteljau算法的特殊性质可得，${\mathbf{P}}_{k0}$由${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_k$计算得到如下

${\mathbf{P}}_{k0}={\sum }_{i=0}^k{B}_{k,i}(u){\mathbf{P}}_i$

3. 将由${\mathbf{P}}_{00}={\mathbf{P}}_0,{\mathbf{P}}_{10},\cdots ,{\mathbf{P}}_{(n-1)0},{\mathbf{P}}_{n0}=\mathbf{C}(u)$定义的新贝塞尔曲线$\mathbf{D}(t)$仍由原控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_n$表示，有

$\mathbf{D}(t)={\sum }_{k=0}^n{B}_{n,k}(t){\mathbf{P}}_{k0}={\sum }_{k=0}^n{B}_{n,k}(t)\left[{\sum }_{i=0}^k{B}_{k,i}(u){\mathbf{P}}_i\right]$

4. 考虑点${\mathbf{P}}_h$的系数，只有$k\ge h$时，该系数才非0，计算如下

$\begin{array}{rl}{\mathbf{P}}_h\text{的系数 }& =B_{n,h}(t)B_{h,h}(u)+B_{n,h+1}(t)B_{h+1,h}(u)+\cdots +B_{n,n}(t)B_{n,h}(u)\\ & =\sum _{j=0}^{n-h}{B}_{n,h+j}(t)B_{h+j,h}(u)\\ & =\sum _{j=0}^{n-h}\left[\frac{n!}{(h+j)!(n-(h+j))!}{t}^{h+j}(1-t{)}^{n-(h+j)}\right]\left[\frac{(h+j)!}{h!j!}{u}^h(1-u{)}^j\right]\\ & \underset{(tu{)}^h,n!,j!}{\overset{提出}{=}}\frac{n!}{h!}(tu{)}^h\sum _{j=0}^{n-h}\frac{1}{((n-h)-j)!j!}[t(1-u){]}^j\left[(1-t{)}^{(n-h)-j}\right]\\ & \underset{凑二项展开}{\overset{增加(n-h)!项}{=}}\frac{n!}{h!(n-h)!}(tu{)}^h\sum _{j=0}^{n-h}\frac{(n-h)!}{((n-h)-j)!j!}[t(1-u){]}^j\left[(1-t{)}^{(n-h)-j}\right]\\ & \underset{二项展开式}{\overset{合并}{=}}\frac{n!}{h!(n-h)!}(tu{)}^h(1-(tu){)}^{n-h}\\ & =B_{n,h}(tu)\end{array}$

5. 由于$t\in [0,1]$，且原贝塞尔曲线只取$[0,u]$的部分，则新曲线符合原曲线如下

$\mathbf{D}(t)={\sum }_{h=0}^n{B}_{n,h}(tu){\mathbf{P}}_h=\mathbf{C}(tu)$

贝塞尔曲线的阶数提升（degree elevation）

阶数的提升要求不改变原贝塞尔曲线的形状，步骤如下

1. 假设$n$阶贝塞尔曲线由$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_n$定义，需要增加阶数到$n+1$而不改变形状

2. 由于升阶的曲线也经过${\mathbf{P}}_0$和${\mathbf{P}}_n$，故这两个点一定在新的控制点集合中

3. 新控制点记为${\mathbf{Q}}_0,{\mathbf{Q}}_1,\cdots ,{\mathbf{Q}}_{n+1}$，则${\mathbf{Q}}_0={\mathbf{P}}_0$且${\mathbf{Q}}_{n+1}={\mathbf{P}}_n$，其他点计算如下

${\mathbf{Q}}_i=\frac{i}{n+1}{\mathbf{P}}_{i-1}+\left(1-\frac{i}{n+1}\right){\mathbf{P}}_{i}1\le i\le n$

阶数提升算法类似de Casteljau算法，但比例不是常数

随着阶数提升，控制折线逐渐收敛至贝塞尔曲线

算法的正确性

整体推导思路类似数学归纳法，利用两条曲线的任意阶导数处处相等的性质，通过使$u=0$处的各阶导数得到${\mathbf{Q}}_0,{\mathbf{Q}}_1,\cdots ,{\mathbf{Q}}_{n+1}$

该推导只说明了必要性，推导过程见网页

B样条（B-spline）曲线

B样条基函数

定义

不同于贝塞尔基函数，B样条基函数具有以下特点：

• 定义域（domain）由节点划分

• 只在几个相邻的子区间非零，而非在整个区间内都非零

首先引入如下定义：

概念	定义
节点（knots）	非递减序列$u_0\le u_1\le \cdots \le u_m$
节点向量（knot vector）	U = { u 0 , u 1 , ⋯   , u m } U=\left\{u_0,u_1,\cdots,u_m\right\} U={u0​,u1​,⋯,um​}
第$i$个节点区间（knot span）	半开半闭区间$[u_i,u_{i+1})$
均匀（uniform）节点向量/节点序列	所有节点区间的跨度都相等
非均匀（non-uniform） 节点向量/节点序列	存在节点区间的跨度不相等
单节点（simple knot）	只出现1次的节点$u_i$
重数(multiplicity)为$k$的重节点(multiple knot)	出现$k$次的节点$u_i$，如$u_i=u_{i+1}=\cdots =u_{i+k-1}$，计节点数量时1个$k$重节点一般算$k$个节点

第$i$个$p$次B样条基函数$N_{i,p}$满足Cox-de Boor递归公式（recursion formula）：

$\begin{array}{rl}& N_{i,0}(u)=\{\begin{array}{ll}1& \text{ if }{u}_i\le u<u_{i+1}\\ 0& \text{ otherwise }\end{array}\\ & N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u)\end{array}$

该递归公式可通过三角形的方案计算

特殊性质

• 基函数$N_{i,p}(u)$在$p+1$个区间$\left[u_i,u_{i+1}\right),\left[u_{i+1},u_{i+2}\right),\dots ,\left[u_{i+p},u_{i+p+1}\right)$上非零，或者说在区间$[u_i,u_{i+p+1})$上非零

• 任意节点区间$[u_i,u_{i+1})$上至多有$p+1$个$p$次基函数非零，依次为$N_{i-p,p}(u),N_{i-p+1,p}(u),\cdots ,N_{i,p}(u)$

• $N_{i,p}(u)$是$N_{i,p-1}(u)$和$N_{i+1,p-1}(u)$的线性组合，二者的系数为如下图的距离之比，均为$u$的线性函数且在$[0,1)$内

重要性质

• $N_{i,p}(u)$是关于$u$的$p$次多项式

• 非负性：$\forall i,p,u$，有$N_{i,p}(u)$非负

• 局部支撑（local support）：$N_{i,p}(u)$是$[u_i,u_{i+p+1})$上的非零多项式

• 任意节点区间$[u_i,u_{i+1})$上至多有$p+1$个$p$次基函数非零，依次为$N_{i-p,p}(u),N_{i-p+1,p}(u),\cdots ,N_{i,p}(u)$

• 单位分解（partition of unity）：区间$[u_i,u_{i+1})$上所有$p$次基函数之和为1

• 节点数为$m+1$，$p$次基函数的数量为$n+1$，则有$m=n+p+1$

• 基函数$N_{i,p}(u)$是由多个曲线组合而成的组合曲线（composite curve），其连接点为$[u_i,u_{i+p+1})$内的节点

• $k$重节点处，基函数$N_{i,p}(u)$满足$C^{p-k}$连续

多重节点的影响

• 每个$k$重节点最多减少基函数的$k-1$个非零区间

  • $k$重节点可视为$k$个不同的节点移动到同一位置，则原$k$个不同的节点组成的$k-1$个非零区间消失

• 在$k$重节点的内部节点（internal knot）上，非零$p$次基函数的数量最多为$p-k+1$

  • 单节点$u_i$上至多有$p+1$个非零$p$次基函数

  • 将$u_{i-1}$移动至$u_i$，则原来从$u_{i-1}$开始的非零基函数变为从$u_i$开始，$u_i$处的非零基函数减少1

  • 将$u_{i+1}$移动至$u_i$，则原来从$u_{i+1}$结束的非零基函数变为从$u_i$结束，$u_i$处的非零基函数减少1

B样条曲线

定义

给定$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_n$和节点向量 U = { u 0 , u 1 , ⋯   , u m } U=\left\{u_0,u_1,\cdots,u_m\right\} U={u0​,u1​,⋯,um​}，则$p$次B样条曲线为

$\mathbf{C}(u)={\sum }_{i=0}^n{N}_{i,p}(u){\mathbf{P}}_i$

其中节点$u_i$对应的点$\mathbf{C}(u_i)$为knot point

B样条曲线可分为如下3类：开放（open）B样条曲线，固定（clamped）B样条曲线，闭合（closed）B样条曲线

固定曲线

固定曲线的前$p+1$和后$p+1$个节点必须相同，即首尾为$p+1$重节点

除首尾节点外其余节点分布均匀，故也称准均匀B样条曲线

开放曲线

如果首末节点的重数小于$p+1$，则B样条曲线不经过首末控制点，也不和首末legs相切

“完全支撑”（full support）：节点区间$[u_i,u_{i+1})$上$p$次非零基函数的个数达到最大值$p+1$

在节点区间$[u_0,u_p)$和$[u_{m-p},u_m]$中没有基函数的“完全支撑”，故开放B样条曲线的定义域为$[u_p,u_{m-p}]$，且在部分节点区间未使用的情况下仍由所有控制点定义

考虑$n=13,p=6,m=n+p+i=13$的均匀开放B样条曲线，其定义域为$[u_6,u_{14}]$如下

闭合曲线

构造由$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_n$定义的$p$次闭合曲线的方法共有两种：

• 包裹控制点（wrapping control points）

  • 构造步骤：

    1. 设计均匀节点序列$u_0=0,u_1=1/m,u_1=2/m,\dots ,u_m=1$，此时曲线的定义域为$[u_p,u_{n-p}]$

    2. 让首尾$p$个节点相等，即${\mathbf{P}}_0={\mathbf{P}}_{n-p+1},{\mathbf{P}}_1={\mathbf{P}}_{n-p+2},\cdots ,{\mathbf{P}}_{p-1}={\mathbf{P}}_n$，如下图的静态结构和动态过程

  • 曲线在$\mathbf{C}\left(u_p\right)=\mathbf{C}\left(u_{n-p}\right)$处满足$C^{p-1}$连续

• 包裹节点（wrapping knots）

  • 网页中内容有误，思路可能是利用给定控制点，和与原节点向量部分相同的新节点向量，新构造一段B样条曲线使得原开放曲线闭合

  • 构造步骤：

    1. 添加新的控制点${\mathbf{P}}_{n+1}={\mathbf{P}}_0$，由于首尾控制点重合，此时控制点个数可认为仍是$n+1$

    2. 找到合适的节点序列$u_0,u_1,\cdots ,u_n$，该节点序列的优势在于不必为均匀的

    3. 添加和前$p+2$个节点相等的$p+2$个节点，即$u_{n+1}=u_0,u_{n+2}=u_1,\cdots ,u_{n+p+2}=u_{p+1}$，以引导受${\mathbf{P}}_{n+1}$影响的曲线末端和受${\mathbf{P}}_0$影响的曲线开端重合，此时节点数量为$n+p+2$

    4. 由上述$n+1$个控制点和$n+p+2$个节点定义的$p$次开放曲线即为闭合曲线

  • 曲线在$\mathbf{C}\left(u_0\right)=\mathbf{C}\left(u_{n+1}\right)$处满足$C^{p-1}$连续

  • 曲线定义域为$[u_0,u_n]$

重要性质

• B样条曲线是分段曲线，每一段都是$p$次曲线

  • 一般来说，次数越低，B样条曲线越接近其控制折线（以下三图的次数分别为7、5、3）

• 必满足$m=n+p+1$

• 固定B样条曲线经过${\mathbf{P}}_0$和${\mathbf{P}}_n$

  • 由于$u_0=u_1=\cdots =u_p=0$，有$N_{0,0}(u),N_{1,0}(u),\dots .,N_{p-1,0}(u)$为零，$N_{p,0}(u)$非零，故$N_{p,0}(0)=1$，则$\mathbf{C}(0)={\mathbf{P}}_0$，类似地有$\mathbf{C}(1)={\mathbf{P}}_n$

• 强凸包性质：如果$u\in [u_i,u_{i+1}]$，则$\mathbf{C}(u)$在控制点${\mathbf{P}}_{i-p},{\mathbf{P}}_{i-p+1},\cdots ,{\mathbf{P}}_i$形成的凸包中

  • 在$[u_i,u_{i+1})$上只有$p+1$个基函数$N_{i,p}(u),\dots ,N_{i-p+1,p}(u),N_{i-p,p}(u)$非零且和为1

• 局部修改性质（local modification scheme）：改变控制点${\mathbf{P}}_i$ 只影响曲线$\mathbf{C}(u)$在$[u_i,u_{i+p+1}]$上的部分

• $\mathbf{C}(u)$在$k$重节点处满足$C^{p-k}$连续

• 变分递减性质：如果曲线位于平面（或空间）中，则直线（或平面）与B样条曲线相交的次数不超过与曲线的控制折线相交的次数

• 贝塞尔曲线是B样条曲线的特例：当$n=p$（对应节点数量为$2(p+1)$）且分别有$p+1$个节点固定在首尾时，B样条曲线等价于贝塞尔曲线

• 仿射不变性：对B样条曲线的仿射变换等价于对控制点的仿射变换

优点：

• B样条曲线可以转化为贝塞尔曲线

• B样条曲线满足贝塞尔曲线的所有性质

• B样条曲线控制更灵活

  • 曲线的次数与控制点的数量相对独立，使用较低次数的曲线仍可保持大量控制点

  • 可改变控制点的位置进行局部而非全局修改

  • 强凸包性质带来更精细的形状控制

缺点：

• 不能表示部分曲线，例如圆和椭圆，需要用到NURBS曲线

重要算法

节点插入

单次插入

节点插入：

• 在现有的节点向量中添加一个新的节点，而不改变曲线的形状

• 新节点可以和原节点相等，使得其重数加1

• 控制点数量必须加1（曲线次数改变会影响全局的形状，因此须保持不变）

若新节点$t\in [u_k,u_{k+1})$，则插入步骤为：

1. $\mathbf{C}(t)$位于${\mathbf{P}}_k,{\mathbf{P}}_{k-1},\cdots ,{\mathbf{P}}_{k-p}$定义的凸包内，且其他控制点的基函数为零，故只需要考虑${\mathbf{P}}_k,{\mathbf{P}}_{k-1},\cdots ,{\mathbf{P}}_{k-p}$的变化

2. 需要找到$p$个新控制点${\mathbf{Q}}_i$位于leg${\mathbf{P}}_{i-1}{\mathbf{P}}_i$上，$i\in [k-p+1,k]$，替代$p-1$个原控制点${\mathbf{P}}_{k+1},{\mathbf{P}}_{k-2},\cdots ,{\mathbf{P}}_{k-p-1}$

3. 新控制点${\mathbf{Q}}_i$满足如下公式和图示

$\begin{array}{rl}{\mathbf{Q}}_i& =\left(1-a_i\right){\mathbf{P}}_{i-1}+a_i{\mathbf{P}}_i\\ a_i& =\frac{t-u_i}{u_{i+p}-u_i}\end{array}$

多次插入

如在$[u_k,u_{k+1})$中插入$h$次新节点$t$，则插入步骤为：

1. 系数$a_{i,h}$满足

$a_{i,h}=\frac{t-u_i}{u_{i+p-(h-1)}-u_i}$

2. 插入$h$次新节点$t$于$s$重节点$u_k$，则可：

   1. 列出受影响的$p+1$原始控制点作为第0列

   2. 忽略最后$s$个控制点${\mathbf{P}}_{k-s+1},{\mathbf{P}}_{k-s+2},\cdots ,{\mathbf{P}}_k$

   3. 依次计算第1到第$h$列

   4. 新的控制点如下图蓝色多边形内所示

计算步骤总结如下图

计算过程可理解为如下割角（corner cutting）的过程

德布尔算法（de Boor’s algorithm）

de Boor算法用于计算在特定$u$处的B样条曲线上的点，其整体思想为：

• 增加节点重数会减少该节点处非零基函数的数量

• $k$重节点处最多有$p-k+1$个非零基函数，故$p$重节点（记为$u_i$）处最多有唯一非零基函数

• 考虑单位分解性质，该唯一非零基函数在该$p$重节点处的值恰好为1，即$\mathbf{C}(u)=N_{i,p}(u){\mathbf{P}}_i={\mathbf{P}}_i$

• 如果重复插入新节点$\underline{u}$至其重数为$p$，则最后生成的新控制点即为原B样条曲线在$u$处的点$\mathbf{C}(u)$

de Boor算法的流程可总结如下

得到最后生成的新控制点${\mathbf{P}}_{k-s,p-s}$和割角的过程如下图

de Boor算法和de Casteljau算法的不同如下

	de Boor算法	de Casteljau算法
计算新控制点的系数	由列和控制点的序数决定	为定值$u$和$1-u$
参与计算的控制点	只有$p+1$个受影响的控制点	所有控制点

当节点向量中只有两个重数为$n+1$的节点时，de Boor算法等价于de Casteljau算法，有

$a_i=\frac{u-u_i}{u_{i+n}-u_n}=ui\in [1,n]$

B样条曲线的分段

B样条曲线的分段使用de Boor算法，其余和贝塞尔曲线的分段完全相同

把B样条曲线分为$[0,u]$和$[u,1]$上两段B样条曲线的步骤如下：

1. 选择控制点

   1. 使用de Boor算法插入新节点$u$

   2. 从${\mathbf{P}}_0$开始沿控制折线前进，当遇到原控制点或者在de Boor算法中计算得到的控制点时，选中该点并转向，直至到达并选中$\mathbf{C}(u)$

   3. 上述过程中选中的点作为$[0,u]$上的B样条曲线的控制点

   4. $[u,1]$上的B样条曲线的控制点采用相同的方式获得，只是路径改为从$\mathbf{C}(u)$出发直至到达${\mathbf{P}}_n$

   5. 上述操作如下左图所示，选出的控制点满足下右图的三角形计算方案中蓝色多边形显示的关系

2. 选择节点

   1. 增加$p+1$重节点$u$，使得分段后的两段B样条曲线经过$\mathbf{C}(u)$

   2. $[0,u]$上的B样条曲线的节点依次为$[0,u)$上所有的原节点和$p+1$重节点$u$

   3. $[u,1]$上的B样条曲线的节点依次为$p+1$重节点$u$和$(u,1]$上所有的原节点

如果$p$次B样条曲线在其节点处分段，则每段曲线都是$p$次贝塞尔曲线，但此操作有如下缺点：

• 会引入大量新的控制点

• B样条曲线在$k$重节点处满足$C^{p-k}$连续，且不受移动控制点的影响，但分段后在贝塞尔曲线之间保持该连续性较为困难

NURBS曲线

NURBS（Non-Uniform Rational B-Splines）曲线，即非均匀有理B样条曲线，是B样条曲线通过齐次坐标到有理曲线的推广，可表示圆、椭圆和许多其他不能用多项式表示的曲线

定义

给定$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_{}$和$m+1$个节点的节点向量 U = { u 0 , u 1 , ⋯   , u m } U=\left\{u_0,u_1,\cdots,u_m\right\} U={u0​,u1​,⋯,um​}，则由此定义的$p$阶B样条曲线为

$\mathbf{C}(u)={\sum }_{i=0}^n{N}_{i,p}(u){\mathbf{P}}_i$

把控制点${\mathbf{P}}_i$重写为

${\mathbf{P}}_i=\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]$

将上述${\mathbf{P}}_i$视为齐次坐标，则其与非零常数相乘不改变其位置，故与权重$w_i$相乘得到与${\mathbf{P}}_i$位置相同的新形式

${\mathbf{P}}_i^w=\left[\begin{array}{c}{w}_i{x}_i\\ w_i{y}_i\\ w_i{z}_i\\ w_i\end{array}\right]$

把${\mathbf{P}}_i^w$带入B样条曲线的表达式中有

${\mathbf{C}}^w(u)={\sum }_{i=0}^n{N}_{i,p}(u){\mathbf{P}}_i^w={\sum }_{i=0}^n{N}_{i,p}(u)\left[\begin{array}{c}{w}_i{x}_i\\ w_i{y}_i\\ w_i{z}_i\\ w_i\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{x}_i\right)\\ {\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{y}_i\right)\\ {\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{z}_i\right)\\ {\sum }_{i=0}^n{N}_{i,p}(u)w_i\end{array}\right]$

${\mathbf{C}}^w(u)$是其次坐标形式的原始B样条曲线，其除以第四个坐标后可转换回笛卡尔坐标

$\mathbf{C}(u)=\left[\begin{array}{c}\frac{{\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{x}_i\right)}{{\sum }_{i=0}^n{N}_{i,p}(u)w_i}\\ \frac{{\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{y}_i\right)}{{\sum }_{i=0}^n{N}_{i,p}(u)w_i}\\ \frac{{\sum }_{i=0}^n{N}_{i,p}(u)\left(w_i{z}_i\right)}{{\sum }_{i=0}^n{N}_{i,p}(u)w_i}\\ 1\end{array}\right]={\sum }_{i=0}^n\frac{N_{i,p}(u)w_i}{{\sum }_{j=0}^n{N}_{j,p}(u)w_j}\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]$

最后得到的表达式为

$\mathbf{C}(u)=\frac{1}{{\sum }_{i=0}^n{N}_{i,p}(u)w_i}{\sum }_{i=0}^n{N}_{i,p}(u)w_i{\mathbf{P}}_i$

NURBS曲线有如下易得的结论：

• 若所有权中均为1，则NURBS曲线退化为B样条曲线

• NURBS曲线是有理的

• 三维空间中的NURBS曲线仅仅是四维空间中B样条曲线的射影

重要性质

给定$n+1$个控制点${\mathbf{P}}_0,{\mathbf{P}}_1,\cdots ,{\mathbf{P}}_{}$和对应的非负权重$w_i$，以及$m+1$个节点的节点向量 U = { u 0 , u 1 , ⋯   , u m } U=\left\{u_0,u_1,\cdots,u_m\right\} U={u0​,u1​,⋯,um​}，则由此定义的$p$次NURBS曲线为

$\mathbf{C}(u)={\sum }_{i=0}^n{R}_{i,p}(u){\mathbf{P}}_i$

其中

$R_{i,p}(u)=\frac{N_{i,p}(u)w_i}{{\sum }_{j=0}^n{N}_{j,p}(u)w_j}$

NURBS基函数的重要性质

• $R_{i,p}(u)$是关于$u$的$p$次有理函数

• 非负性：$\forall i,p$，有$R_{i,p}(u)\ge 0$

• 局部支撑：假设$w_i>0$，则和$N_{i,p}(u)$一样，$R_{i,p}(u)$在$[u_i,u_{i+p+1})$上非零

• 在任意节点区间$[u_i,u_{i+p+1})$上，至多有$p+1$个非零$p$次基函数$R_{i-p,p}(u),R_{i-p+1,p}(u),,\cdots ,R_{i,p}(u)$

• 单位分解：节点区间$[u_i,u_{i+p+1})$上的所有非零$p$次基函数之和为1

• 节点数为$m+1$，$p$次基函数的数量为$n+1$，则有$m=n+p+1$

• 基函数$R_{i,p}(u)$是由多个曲线组合而成的组合曲线，其连接点为$[u_i,u_{i+p+1})$内的节点

• $k$重节点处，基函数$R_{i,p}(u)$满足$C^{p-k}$连续

• 如果$\forall i$，有$w_i=c\ne 0$，则$R_{i,p}(u)=N_{i,p}(u)$

NURBS曲线的重要性质

• NURBS曲线$\mathbf{C}(u)$为分段曲线，每一段为$p$次有理贝塞尔曲线

• 必有$m=n+p+1$

• 强凸包性质：如果$u\in [u_i,u_{i+1}]$且所有权重非负，则$\mathbf{C}(u)$在控制点${\mathbf{P}}_{i-p},{\mathbf{P}}_{i-p+1},\cdots ,{\mathbf{P}}_i$形成的凸包中

• 局部修改性质：改变控制点${\mathbf{P}}_i$ 只影响曲线$\mathbf{C}(u)$在$[u_i,u_{i+p+1}]$上的部分

• $\mathbf{C}(u)$在$k$重节点处满足$C^{p-k}$连续

• 变分递减性质：如果曲线位于平面（或空间）中，则直线（或平面）与B样条曲线相交的次数不超过与曲线的控制折线相交的次数

• 贝塞尔曲线和B样条曲线是NURBS曲线的特例

  • 当所有权重相等，NURBS曲线等价于B样条曲线

  • 当所有权重相等，$n=p$（对应节点数量为$2(p+1)$）且分别有$p+1$个节点固定在首尾时，B样条曲线等价于贝塞尔曲线

• 射影不变性（projective invariance）：对NURBS曲线的射影变换等价于对控制点的射影変換

节点插入算法

节点插入的步骤如下（例子见网页）：

1. 把给定的三维NURBS曲线转换为四维的B样条曲线

2. 对四维的B样条曲线进行节点插入

3. 把新的控制点集射影回三维

有理贝塞尔曲线

把四维贝塞尔曲线射影到超平面$w=1$可得到三维有理贝塞尔曲线，其表达式为

$\mathbf{C}(u)={\sum }_{i=0}^n{R}_{i,p}(u){\mathbf{P}}_i$

其中

$R_{i,p}(u)=\frac{B_{i,p}(u)w_i}{{\sum }_{j=0}^n{B}_{j,p}(u)w_j}$

有理贝塞尔曲线具有以下特点：

• 不具有局部修改性质

• 满足射影不变性质（仿射变换是射影变换的子集）

圆锥曲线

唯一确定圆锥曲线的五个条件：

• 由三个非共线点${\mathbf{P}}_0=\left(U_0,V_0\right),{\mathbf{P}}_1=\left(U_1,V_1\right),{\mathbf{P}}_2=\left(U_2,V_2\right)$定义的2阶贝塞尔曲线为抛物线

• 假设圆锥曲线经过${\mathbf{P}}_0,{\mathbf{P}}_2$，且在${\mathbf{P}}_0,{\mathbf{P}}_2$处分别和${\mathbf{P}}_0{\mathbf{P}}_1,{\mathbf{P}}_2{\mathbf{P}}_1$相切，用2次隐式方程表示如下，其中6个系数为未知数

$p(x,y)=a{x}^2+2bxy+c{y}^2+2dx+2ey+f=0$

• 如果$f\ne 0$，则有只包含五个未知数的方程如下

$a{x}^2+2bxy+c{y}^2+2dx+2ey+1=0$

• 方程的梯度如下

${\nabla }_{p(x,y)}=⟨2ax+2by+2d,2bx+2cy+2e⟩$

• 目前可以得到四个方程

  1. ${\mathbf{P}}_0$在曲线上：$a{U}_0^2+2b{U}_0{V}_0+c{V}_0^2+2d{U}_0+2e{V}_0+1=0$

  2. ${\mathbf{P}}_2$在曲线上：$a{U}_2^2+2b{U}_2{V}_2+c{V}_2^2+2d{U}_2+2e{V}_2+1=0$

  3. 在${\mathbf{P}}_0$处和${\mathbf{P}}_0{\mathbf{P}}_1$相切，则切向量和法向量垂直：$\frac{V_1-V_0}{U_1-U_0}=-\frac{a{U}_0+b{V}_0+d}{b{U}_0+c{V}_0+e}$

  4. 在${\mathbf{P}}_2$处和${\mathbf{P}}_2{\mathbf{P}}_1$相切，则切向量和法向量垂直：$\frac{V_2-V_1}{U_2-U_1}=-\frac{a{U}_2+b{V}_2+d}{b{U}_2+c{V}_2+e}$

• 第五个方程则需要再找到一个在曲线上的点，改点应位于控制三角形的内部，以保持凸包性质

• 该点受控制点${\mathbf{P}}_1$的权重的影响，故可使用由${\mathbf{P}}_0,{\mathbf{P}}_1,{\mathbf{P}}_2$定义的有理贝塞尔曲线，其权重依次为1、$w$、1，则控制点的系数如下

$\begin{array}{rl}& B_{2,0}(u)=(1-u{)}^2\\ & B_{2,1}(u)=2(1-u)u\\ & B_{2,2}(u)=u^2\end{array}$

• 2次有理贝塞尔曲线的表达式为

$\mathbf{C}(u)=\frac{1}{(1-u{)}^2+2(1-u)uw+u^2}\left((1-u{)}^2{\mathbf{P}}_0+2(1-u)uw{\mathbf{P}}_1+u^2{\mathbf{P}}_2\right)$

• 通过把${\mathbf{P}}_0,{\mathbf{P}}_2$置于原点两侧，且${\mathbf{P}}_0{\mathbf{P}}_2$的中点位于原点，则有${\mathbf{P}}_0=-{\mathbf{P}}_2$，进一步有$\mathbf{C}(0.5)=\frac{w}{1+w}{\mathbf{P}}_1=\frac{|\mathrm{MX}|}{|{\mathrm{MP}}_1|}{\mathbf{P}}_1$

• 根据射影几何中的定理，由三个非共线控制点${\mathbf{P}}_0,{\mathbf{P}}_1,{\mathbf{P}}_2$和其各自的权重1、$w$、1定义的有理贝塞尔曲线满足：

  • $w>1$时，为双曲线（hyperbola）

  • $w=1$时，为抛物线（parabola）

  • $w<1$时，为椭圆（ellipse）

圆弧和圆

圆弧使用有理贝塞尔曲线表示，整圆使用NURBS曲线表示，详情见网页

圆弧

圆是椭圆的特例，根据下图所示关系和圆锥曲线小节内容，$w=\sin (a)$时可表示弦${\mathbf{P}}_0{\mathbf{P}}_2$对应的圆弧

整圆

多个圆弧拼接在一起即为整圆：

• 用三角形定义整圆，则节点为$0,0,0,1/3,1/3,2/3,2/3,1,1,1$

• 用正方形定义整圆，则节点为$0,0,0,1/4,1/4,1/2,1/2,3/4,3/4,1,1,1$