目录

参数曲线

网页资料

image

相关概念

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

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

CkC^k连续: \(\forall i \le k\),f(b)f(b)g(m)g(m)ii阶导数在连接点处相等

几何连续性GkG^k连续有如下等价定义:

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

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

注意:C2C^2连续一定曲率连续,但曲率连续不一定C2C^2连续

贝赛尔曲线

给定空间中n+1n+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)

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

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

  • 贝塞尔曲线经过P0\mathbf{P}_0Pn\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]\)

image

原控制点Pk\mathbf{P}_k变为 \(\mathbf{P}_k + \mathbf{v}\)后,新的贝塞尔曲线

\(\begin{aligned} \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{aligned}\)

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

image

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

image

定义第ii次迭代得到的点的编号为i0,,i(ni)i0,\cdots,i(n-i),则控制点(第0次迭代)编号为00,,0n00,\cdots,0n,步骤如下:

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

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

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

image

计算步骤如下:

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

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

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

image

\(\mathbf{P}_{i, j}=(1-u) \mathbf{P}_{i-1, j}+u \mathbf{P}_{i-1, j+1} \quad\left\{\begin{array}{l} i=1,2, \ldots, n \\ j=0,1, \ldots, n-i \end{array}\right.\)

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

image

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

image

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

image

image        image

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

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

  3. 考虑到一条路径传播的贡献只和起点与终点的相对位置有关,和路径形状无关,由0i0in0n0的每条路径的贡献均为ui(1u)niu^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} \mathbf{0} \mathbf{i}=\frac{n !}{i !(n-i) !} u^i(1-u)^{n-i} \mathbf{P}_i=B_{n, i}(u) \mathbf{P}_i\)

由于控制点是与uu无关的常数,因此计算导数C(u)\mathbf{C}^{\prime}(u)可以简化为计算Bn,i(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]\),则C(u)\mathbf{C}^{\prime}(u)满足

\(\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\)

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

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

image

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

C(1)=m(PmPm1)=D(0)=n(Q1Q0)\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)

贝塞尔曲线的导数可写为

C(u)=i=0n1Bn1,i(u){n(Pi+1Pi)}=n[(i=0n1Bn,i(u)Pi+1)(i=0n1Bn,i(u)Pi)]\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}

定义

\(\mathbf{C}_1(u)=\sum_{i=0}^{n-1} B_{n, i}(u) \mathbf{P}_{i+1} \qquad  \mathbf{C}_2(u)=\sum_{i=0}^{n-1} B_{n, i}(u) \mathbf{P}_i\)

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

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

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

image

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

\(\mathbf{D}_i^k=\mathbf{D}_{i+1}^{k-1}-\mathbf{D}_i^{k-1} \quad 0 \leq i \leq n-k\)

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

\(\begin{aligned} \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{aligned}\)

利用de Casteljau算法计算特定uu下的kk阶导数值的流程如下:

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

image        image

  1. 把de Casteljau算法应用在第kk列的点上得到导数值

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

  • 将给定的贝塞尔曲线C(u)\mathbf{C}(u)在某一点切割成[0,u][0,u][u,1][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}\)组成

image        image

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

image

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

  1. 设原曲线为

\(\mathbf{C}(u)=\sum_{i=0}^n B_{n, i}(u) \mathbf{P}_i\)

  1. 由de Casteljau算法的特殊性质可得, \(\mathbf{P}_{k 0}\)由 \(\mathbf{P}_0, \mathbf{P}_1, \cdots, \mathbf{P}_k\)计算得到如下

\(\mathbf{P}_{k 0}=\sum_{i=0}^k B_{k, i}(u) \mathbf{P}_i\)

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

D(t)=k=0nBn,k(t)Pk0=k=0nBn,k(t)[i=0kBk,i(u)Pi]\mathbf{D}(t)=\sum_{k=0}^n B_{n, k}(t) \mathbf{P}_{k 0}=\sum_{k=0}^n B_{n, k}(t)\left[\sum_{i=0}^k B_{k, i}(u) \mathbf{P}_i\right]

  1. 考虑点Ph\mathbf{P}_{h}的系数,只有 \(k \ge h\)时,该系数才非0,计算如下

\(\begin{aligned} \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]\\ & \xlongequal[(tu)^h,n!,j!]{提出} \frac{n !}{h !}(t u)^h \sum_{j=0}^{n-h} \frac{1}{((n-h)-j) ! j !}[t(1-u)]^j\left[(1-t)^{(n-h)-j}\right]\\ &\xlongequal[凑二项展开]{增加(n-h)!项}\frac{n !}{h !(n-h) !}(t u)^h \sum_{j=0}^{n-h} \frac{(n-h) !}{((n-h)-j) ! j !}[t(1-u)]^j\left[(1-t)^{(n-h)-j}\right]\\ &\xlongequal[二项展开式]{合并} \frac{n !}{h !(n-h) !}(t u)^h(1-(t u))^{n-h}\\ &=B_{n, h}(t u) \end{aligned}\)

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

\(\mathbf{D}(t)=\sum_{h=0}^n B_{n, h}(t u) \mathbf{P}_h=\mathbf{C}(t u)\)

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

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

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

  3. 新控制点记为 \(\mathbf{Q}_0, \mathbf{Q}_1, \cdots, \mathbf{Q}_{n+1}\),则Q0=P0\mathbf{Q}_0=\mathbf{P}_0Qn+1=Pn\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 \quad 1 \leq i \leq n\)

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

image

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

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

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

B样条(B-spline)曲线

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

  • 定义域(domain)由节点划分

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

首先引入如下定义:

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

iipp次B样条基函数Ni,pN_{i,p}满足Cox-de Boor递归公式(recursion formula):

\(\begin{aligned} & N_{i, 0}(u)= \begin{cases}1 & \text { if } u_i \leq u<u_{i+1} \\ 0 & \text { otherwise }\end{cases} \\ & 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{aligned}\)

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

image

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

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

image        image

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

image

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

  • 非负性: \(\forall i,p,u\),有Ni,p(u)N_{i,p}(u)非负

  • 局部支撑(local support):Ni,p(u)N_{i,p}(u)[ui,ui+p+1)[u_i,u_{i+p+1})上的非零多项式

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

  • 单位分解(partition of unity):区间[ui,ui+1)[u_i,u_{i+1})上所有pp次基函数之和为1

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

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

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

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

    • kk重节点可视为kk个不同的节点移动到同一位置,则原kk个不同的节点组成的k1k-1个非零区间消失
  • kk重节点的内部节点(internal knot)上,非零pp次基函数的数量最多为pk+1p-k+1

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

    • ui1u_{i-1}移动至uiu_i,则原来从ui1u_{i-1}开始的非零基函数变为从uiu_i开始,uiu_i处的非零基函数减少1

    • ui+1u_{i+1}移动至uiu_i,则原来从ui+1u_{i+1}结束的非零基函数变为从uiu_i结束,uiu_i处的非零基函数减少1

给定n+1n+1个控制点 \(\mathbf{P}_0, \mathbf{P}_1, \cdots, \mathbf{P}_n\)和节点向量U={u0,u1,,um}U=\left\{u_0,u_1,\cdots,u_m\right\},则pp次B样条曲线为

\(\mathbf{C}(u)=\sum_{i=0}^n N_{i, p}(u) \mathbf{P}_i\)

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

B样条曲线可分为如下3类:开放(open)B样条曲线,固定(clamped)B样条曲线,闭合(closed)B样条曲线

image

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

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

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

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

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

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

image

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

  • 包裹控制点(wrapping control points)

    • 构造步骤:

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

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

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

image

image

image

  • 包裹节点(wrapping knots)

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

    • 构造步骤:

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

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

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

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

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

    • 曲线定义域为[u0,un][u_0,u_{n}]

image

  • B样条曲线是分段曲线,每一段都是pp次曲线

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

image

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

  • 固定B样条曲线经过P0\mathbf{P}_0Pn\mathbf{P}_n   

    • 由于u0=u1==up=0u_0=u_1=\cdots=u_p=0,有 \(N_{0,0}(u), N_{1,0}(u), \ldots ., N_{p-1,0}(u)\)为零,Np,0(u)N_{p,0}(u)非零,故Np,0(0)=1N_{p,0}(0)=1,则C(0)=P0\mathbf{C}(0)=\mathbf{P}_0,类似地有C(1)=Pn\mathbf{C}(1)=\mathbf{P}_n
  • 强凸包性质:如果 \(u \in [u_i,u_{i+1}]\),则C(u)\mathbf{C}(u)在控制点Pip,Pip+1,,Pi\mathbf{P}_{i-p},\mathbf{P}_{i-p+1},\cdots,\mathbf{P}_{i}形成的凸包中

    • [ui,ui+1)[u_i,u_{i+1})上只有p+1p+1个基函数 \(N_{i, p}(u), \ldots, N_{i-p+1, p}(u), N_{i-p, p}(u)\)非零且和为1
  • 局部修改性质(local modification scheme):改变控制点Pi\mathbf{P}_i 只影响曲线C(u)\mathbf{C}(u)[ui,ui+p+1][u_i,u_{i+p+1}]上的部分

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

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

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

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

优点:

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

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

  • B样条曲线控制更灵活

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

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

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

缺点:

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

节点插入:

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

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

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

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

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

  2. 需要找到pp个新控制点Qi\mathbf{Q}_{i}位于legPi1Pi\mathbf{P}_{i-1}\mathbf{P}_{i}上, \(i \in [k-p+1,k]\),替代p1p-1个原控制点Pk+1,Pk2,,Pkp1\mathbf{P}_{k+1},\mathbf{P}_{k-2},\cdots,\mathbf{P}_{k-p-1}

image

  1. 新控制点Qi\mathbf{Q}_{i}满足如下公式和图示

\(\begin{aligned} \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{aligned}\)

image

如在[uk,uk+1)[u_k,u_{k+1})中插入hh次新节点tt,则插入步骤为:

  1. 系数ai,ha_{i,h}满足

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

  1. 插入hh次新节点ttss重节点uku_k,则可:

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

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

    3. 依次计算第1到第hh

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

image

计算步骤总结如下图

image

计算过程可理解为如下割角(corner cutting)的过程

image

de Boor算法用于计算在特定uu处的B样条曲线上的点,其整体思想为:

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

  • kk重节点处最多有pk+1p-k+1个非零基函数,故pp重节点(记为uiu_i)处最多有唯一非零基函数

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

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

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

image

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

image

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

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

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

\(a_i=\frac{u-u_i}{u_{i+n}-u_n}=u \quad i \in[1,n]\)

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

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

  1. 选择控制点

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

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

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

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

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

image        image

  1. 选择节点

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

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

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

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

  • 会引入大量新的控制点

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

NURBS曲线

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

给定n+1n+1个控制点P0,P1,,P\mathbf{P}_{0},\mathbf{P}_{1},\cdots,\mathbf{P}_{}m+1m+1个节点的节点向量U={u0,u1,,um}U=\left\{u_0,u_1,\cdots,u_m\right\},则由此定义的pp阶B样条曲线为

\(\mathbf{C}(u)=\sum_{i=0}^n N_{i, p}(u) \mathbf{P}_i\)

把控制点Pi\mathbf{P}_{i}重写为

Pi=[xiyizi1]\mathbf{P}_i=\left[\begin{array}{c} x_i \\ y_i \\ z_i \\ 1 \end{array}\right]

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

Piw=[wixiwiyiwiziwi]\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]

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

Cw(u)=i=0nNi,p(u)Piw=i=0nNi,p(u)[wixiwiyiwiziwi]=[i=0nNi,p(u)(wixi)i=0nNi,p(u)(wiyi)i=0nNi,p(u)(wizi)i=0nNi,p(u)wi]\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]

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

C(u)=[i=0nNi,p(u)(wixi)i=0nNi,p(u)wii=0nNi,p(u)(wiyi)i=0nNi,p(u)wii=0nNi,p(u)(wizi)i=0nNi,p(u)wi1]=i=0nNi,p(u)wij=0nNj,p(u)wj[xiyizi1]\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+1n+1个控制点P0,P1,,P\mathbf{P}_{0},\mathbf{P}_{1},\cdots,\mathbf{P}_{}和对应的非负权重wiw_i,以及m+1m+1个节点的节点向量U={u0,u1,,um}U=\left\{u_0,u_1,\cdots,u_m\right\},则由此定义的pp次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}\)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  • 射影不变性(projective invariance):对NURBS曲线的射影变换等价于对控制点的射影変換

节点插入的步骤如下(例子见网页):

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

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

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

把四维贝塞尔曲线射影到超平面w=1w=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阶贝塞尔曲线为抛物线

  • 假设圆锥曲线经过P0,P2\mathbf{P}_0,\mathbf{P}_2,且在P0,P2\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+2 b x y+c y^2+2 d x+2 e y+f=0\)

  • 如果 \(f \neq 0\),则有只包含五个未知数的方程如下

\(a x^2+2 b x y+c y^2+2 d x+2 e y+1=0\)

  • 方程的梯度如下

\(\nabla_{p(x, y)}=\langle 2 a x+2 b y+2 d, 2 b x+2 c y+2 e\rangle\)

  • 目前可以得到四个方程

    1. P0\mathbf{P}_0在曲线上: \(a U_0^2+2 b U_0 V_0+c V_0^2+2 d U_0+2 e V_0+1=0\)

    2. P2\mathbf{P}_2在曲线上: \(a U_2^2+2 b U_2 V_2+c V_2^2+2 d U_2+2 e V_2+1=0\)

    3. P0\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. P2\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}\)

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

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

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

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

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

  • 通过把P0,P2\mathbf{P}_0,\mathbf{P}_2置于原点两侧,且P0P2\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}|}{\left|\mathrm{MP}_1\right|} \mathbf{P}_1\)

image

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

    • \(w > 1\)时,为双曲线(hyperbola)

    • \(w = 1\)时,为抛物线(parabola)

    • \(w < 1\)时,为椭圆(ellipse)

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

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

image

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

  • 用三角形定义整圆,则节点为0,0,0,1/3,1/3,2/3,2/3,1,1,10,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,10,0,0,1/4,1/4,1/2,1/2,3/4,3/4,1,1,1

image