General matrix representations for B-splines 论文理解

这篇论文 [1] 比较基础，在很多与 B 样条有关的论文中都能找到对它的引用。

B 样条曲线是参数方程，使用参数 (t) 的多项式表达连续曲线。

多项式的计算是递归计算的，k 阶 B 样条的多项式的参数是 k 阶的，需要从 k-1 阶参数计算。可以看回形针视频 [2] 有一个简单的认识。

注意一下 spline 的 degree 与 order。对于多项式而言 degree 与 order 是相等的，是同义词。而对于 spline，order = degree + 1。[3]

B-splines of order (p+1) are connected piece-wise polynomial functions of degree (p).

所以在 [2] 中提到的 “4 个控制点 3 次 B 样条”的 order 为 4，degree 为 3。B 样条上的任意一点受到 order 个控制点的影响，可以看做是 3 次多项式。

作者在这篇论文提出了一种 B 样条曲线使用 Matrix 表达的方式，这种方式在计算上有优势。

回想到斐波那契数列，也是递归计算，其中有一种计算方式也是用矩阵表达[4]。甚至可以从这个矩阵的表达中得出斐波那契数列的通项公式。

以下正文中的公式引用使用两种形式：本文的公式，如(1)；原文 [1] 的公式，如(1*)。

以下推导过程应该会很啰嗦，我的目的是我自己能够看得懂，通过此文能够毫不费力地了解到论文中所有没有明说的细节。

1. Toeplitz Matrix 表达多项式

Toeplitz Matrix 的定义在 [1] 中有，2.1 的第一句话。

The Toeplitz matrix is one whose elements on any line parallel to the main diagonal are all equal.

如果 Toeplitz Matrix 是一个 lower triangular matrix 的话，它就能表达一个多项式。我尝试用二次型的方式去解释这个矩阵，但是失败了。反正 [2] 2.1 的第二个矩阵 (mathbf{T}) 可以表示一个多项式。

[egin{align} f(x) = a_0 + a_1 x + a_x x^2 + dots + a_{n-1} x^{n-1} (a_{n-1} eq 0) end{align}]

你问我这个 (mathbf{T}) 是多少维的方阵，我也不知道，按需构建。这不重要。重点在 [1] 2.2 对它的应用中。

假设有两个多项式。

[egin{align} g(x) &= c_0 + c_1 x + c_x x^2 + dots + c_{m-1} x^{m-1} (c_{m-1} eq 0) \ q(x) &= d_0 + d_1 x + d_x x^2 + dots + d_{n-1} x^{n-1} (d_{n-1} eq 0) end{align}]

(f(x) = g(x) q(x))，(f(x)) 的最高次应当是 (m + n - 2)。

所以 [1] 2.2 中的矩阵具体的维数应当如下。

[egin{align} f(x) &= g(x) q(x) \ &= mathbf{X} egin{bmatrix} c_0 & phantom{} & phantom{} & phantom{} & phantom{} & phantom{} & 0 \ c_1 & c_0 \ vdots & ddots & ddots \ c_{m - 1} & cdots & ddots \ phantom{} & ddots & cdots & ddots & ddots \ phantom{} & phantom{} & c_{m-1} & cdots & ddots & c_0 \ 0 & phantom{} & phantom{} & c_{m-1} & cdots & c_1 & c_0 end{bmatrix}_{(m + n - 1) imes (m + n - 1)} egin{bmatrix} d_0 \ d_1 \ vdots \ d_{n - 1} \ 0 \ vdots \ 0 end{bmatrix}_{(m + n - 1) imes 1} \ &= mathbf{X} egin{bmatrix} c_0 & phantom{} & phantom{} & 0 \ c_1 & c_0 \ vdots & ddots & ddots \ vdots & ddots & ddots & c_0\ vdots & ddots & ddots & vdots\ c_{m - 1} & c_{m-2} & cdots & c_{m-n} \ phantom{} & ddots & cdots & vdots \ phantom{} & phantom{} & c_{m-1} & c_{m-2} \ 0 & phantom{} & phantom{} & c_{m-1} end{bmatrix}_{(m + n - 1) imes n} egin{bmatrix} d_0 \ d_1 \ vdots \ d_{n - 1} end{bmatrix}_{n imes 1} end{align}]

[5] 给了一个具体的例子，可以看一下。

2. Toeplitz Matrix 表达 B 样条的 basis functions

我对 B 样条不熟，各种乱七八糟的多项式，不清晰。希望看完这个能够清晰一些。

(1*) 给出了 B 样条各阶 basis function 的递归公式。

作者说 By means of basis translation from B-spline to power basis, ...，所以 basis function 就从以参数 (t) 作为变量，转换到以 (u) 作为变量。作者说参考的两个参考文献都是书，而且其中的 Practical Guide to Splines 还只能找到扫描版。不容易搜索关键字 power basis。所以我另外找了一篇参考文献 [6]，也是 De Boor 写的。看 [6] 的 8. Condition，是精度的问题，用 power basis 转换变量的取值空间能够更容易控制精度。注意 Condition 这个词，联想到了 Matrix Condition Number（[7] 2.2.2 与 4.3.2）。Matrix 用于表示一个线性系统，Condition Number 表示 forward error / backward error，forward error 是 x 的误差，而 backward error 是 f(x) 的误差。希望 condition number 尽可能小，即 f(x) 是观察值，用观察值去估计实际关心的状态值 x，希望观察值小的抖动，只会造成状态值小的抖动，这样就容许了观察值的观测误差，对测量条件没有这么苛刻。[6] 的 8. Condition 也是这个意思。

[egin{align} B_{j,k}(t) &= {t - t_j over t_{j + k - 1} - t_j} B_{j,k-1}(t) + {t_{j + k} - t over t_{j + k} - t_{j + 1}} B_{j + 1, k -1}(t) \ &= ({t_i - t_j over t_{j + k - 1} - t_j} + {t - t_i over t_{i + 1} - t_i}{t_{i + 1} - t_i over t_{j + k - 1} - t_j}) B_{j,k-1}(t) \ &phantom{=} + ({t_{j + k} - t_i over t_{j + k} - t_{j + 1}} - {t - t_i over t_{i + 1} - t_i}{t_{i + 1} - t_i over t_{j + k} - t_{j + 1}}) B_{j + 1, k -1}(t) \ B_{j,k}(u) &= (d_{0,j} + ud_{1,j}) B_{j,k-1}(u) + (h_{0,j} + uh_{1,j}) B_{j + 1, k -1}(u) \ u &= {t - t_i over t_{i+1} - t_i}, uin [0, 1) end{align}]

上式中出现的 index (i) 是什么意思？公式 (3*) 上方的文字做了解释，(i) 表示当前关注的 (tin[t_i,t_{i+1}))，这个区间内的 B 样条与周围 (k) 个 basis functions 以及它们对应的 control points 有关，(j) 是这 (k) 个部分的 index。

作者给出 (B_{j,k-1}(u)) 的定义。

[egin{align} B_{j,k-1}(u) &= egin{bmatrix} 1 & u & dots & u^{k-2} end{bmatrix} egin{bmatrix} N_{0,j}^{k-1} \ N_{1,j}^{k-1} \ vdots \ N_{k-2,j}^{k-1} end{bmatrix} \ u &= {t - t_i over t_{i+1} - t_i}, uin [0, 1) end{align}]

用 Toeplitz Matrix 表示这个迭代过程，(u) 从 (k-2) 此上升到 (k-1) 次。

[egin{align} B_{j,k}(u) &= (d_{0,j} + ud_{1,j}) B_{j,k-1}(u) + (h_{0,j} + uh_{1,j}) B_{j + 1, k -1}(u) \ &= egin{bmatrix} 1 & u & dots & u^{k-1} end{bmatrix} left{ egin{bmatrix} N_{0,j}^{k-1} & 0 \ N_{1,j}^{k-1} & N_{0,j}^{k-1} \ vdots & N_{1,j}^{k-1} \ N_{k-2,j}^{k-1} & vdots \ 0 & N_{k-2,j}^{k-1} end{bmatrix} egin{bmatrix} d_{0,j} \ d_{1,j} end{bmatrix} + egin{bmatrix} N_{0,j+1}^{k-1} & 0 \ N_{1,j+1}^{k-1} & N_{0,j+1}^{k-1} \ vdots & N_{1,j+1}^{k-1} \ N_{k-2,j+1}^{k-1} & vdots \ 0 & N_{k-2,j+1}^{k-1} end{bmatrix} egin{bmatrix} h_{0,j} \ h_{1,j} end{bmatrix} ight} end{align}]

如公式 (4*)，在 (tin[t_i,t_{i+1})) 这个区间内，点的坐标可以写作。

[egin{align} mathbf{c}_{i-k+1} &= egin{bmatrix} B_{i-k+1,k}(u) & B_{i-k+2,k}(u) & cdots & B_{i,k}(u) end{bmatrix} imes egin{bmatrix} mathbf{V}_{i-k+1} \ mathbf{V}_{i-k+2} \ vdots \ mathbf{V}_{i} end{bmatrix} end{align}]

这 (k) 个 basis function 都是 (k-1) 次的。扩展开。

[egin{align} &phantom{=}egin{bmatrix} B_{i-k+1,k}(u) & B_{i-k+2,k}(u) & cdots & B_{i,k}(u) end{bmatrix} \ &= egin{bmatrix} 1 & u & dots & u^{k-1} end{bmatrix} mathbf{M}^k(i) \ mathbf{M}^k(i) &= egin{bmatrix} N_{0,i-k+1}^{k} & N_{0,i-k+2}^{k} & dots & N_{0,i}^{k} \ N_{1,i-k+1}^{k} & N_{1,i-k+2}^{k} & dots & N_{1,i}^{k} \ vdots & vdots & & vdots \ N_{k-1,i-k+1}^{k} & N_{k-1,i-k+2}^{k} & dots & N_{k-1,i}^{k} end{bmatrix} end{align}]

所以，B 样条公式的关注重点在 (mathbf{M}^k(i))，(mathbf{M}^k(i)) 是变化的部分，其他部分相对固定。

3. (mathbf{M}^k(i)) 的迭代计算

(mathbf{M}^k(i)) 第 1 列。

[egin{align} egin{bmatrix} N_{0,i-k+1}^{k} \ N_{1,i-k+1}^{k} \ vdots \ vdots \ N_{k-1,i-k+1}^{k} end{bmatrix} &= egin{bmatrix} N_{0,i-k+1}^{k-1} & 0 \ N_{1,i-k+1}^{k-1} & N_{0,i-k+1}^{k-1} \ vdots & N_{1,i-k+1}^{k-1} \ N_{k-2,i-k+1}^{k-1} & vdots \ 0 & N_{k-2,i-k+1}^{k-1} end{bmatrix} egin{bmatrix} d_{0,i-k+1} \ d_{1,i-k+1} end{bmatrix} + egin{bmatrix} N_{0,i-k+2}^{k-1} & 0 \ N_{1,i-k+2}^{k-1} & N_{0,i-k+2}^{k-1} \ vdots & N_{1,i-k+2}^{k-1} \ N_{k-2,i-k+2}^{k-1} & vdots \ 0 & N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} h_{0,i-k+1} \ h_{1,i-k+1} end{bmatrix} end{align}]

(mathbf{M}^k(i)) 可以每一列分别计算再求和，记得补 (mathbf{0}) 列。

[egin{align} egin{bmatrix} N_{0,i-k+1}^{k} & overbrace{0 cdots 0}^{k-1} \ N_{1,i-k+1}^{k} & 0 cdots 0 \ vdots & vdots \ vdots & vdots \ N_{k-1,i-k+1}^{k} & 0 cdots 0 end{bmatrix} &= egin{bmatrix} N_{0,i-k+1}^{k-1} & 0 \ N_{1,i-k+1}^{k-1} & N_{0,i-k+1}^{k-1} \ vdots & N_{1,i-k+1}^{k-1} \ N_{k-2,i-k+1}^{k-1} & vdots \ 0 & N_{k-2,i-k+1}^{k-1} end{bmatrix} egin{bmatrix} d_{0,i-k+1} & overbrace{0 cdots 0}^{k-1} \ d_{1,i-k+1} & 0 cdots 0 end{bmatrix} + egin{bmatrix} N_{0,i-k+2}^{k-1} & 0 \ N_{1,i-k+2}^{k-1} & N_{0,i-k+2}^{k-1} \ vdots & N_{1,i-k+2}^{k-1} \ N_{k-2,i-k+2}^{k-1} & vdots \ 0 & N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} h_{0,i-k+1} & overbrace{0 cdots 0}^{k-1} \ h_{1,i-k+1} & 0 cdots 0 end{bmatrix} end{align}]

下面使用了矩阵分块。

[egin{align} mathbf{M}^k(i) &= egin{bmatrix} N_{0,i-k+1}^{k} & N_{0,i-k+2}^{k} & dots & N_{0,i}^{k} \ N_{1,i-k+1}^{k} & N_{1,i-k+2}^{k} & dots & N_{1,i}^{k} \ vdots & vdots & & vdots \ N_{k-1,i-k+1}^{k} & N_{k-1,i-k+2}^{k} & dots & N_{k-1,i}^{k} end{bmatrix} \ &= egin{bmatrix} N_{0,i-k+1}^{k} & 0 & dots & 0 \ N_{1,i-k+1}^{k} & 0 & dots & 0 \ vdots & vdots & & vdots \ N_{k-1,i-k+1}^{k} & 0 & dots & 0 end{bmatrix} + egin{bmatrix} 0 & N_{0,i-k+2}^{k} & dots & 0 \ 0 & N_{1,i-k+2}^{k} & dots & 0 \ vdots & vdots & & vdots \ 0 & N_{k-1,i-k+2}^{k} & dots & 0 end{bmatrix} + cdots + egin{bmatrix} 0 & 0 & dots & N_{0,i}^{k} \ 0 & 0 & dots & N_{1,i}^{k} \ vdots & vdots & & vdots \ 0 & 0 & dots & N_{k-1,i}^{k} end{bmatrix} \ &= egin{bmatrix} N_{0,i-k+1}^{k-1} & 0 \ N_{1,i-k+1}^{k-1} & N_{0,i-k+1}^{k-1} \ vdots & N_{1,i-k+1}^{k-1} \ N_{k-2,i-k+1}^{k-1} & vdots \ 0 & N_{k-2,i-k+1}^{k-1} end{bmatrix} egin{bmatrix} d_{0,i-k+1} & overbrace{0 cdots 0}^{k-1} \ d_{1,i-k+1} & 0 cdots 0 end{bmatrix} + egin{bmatrix} N_{0,i-k+2}^{k-1} & 0 \ N_{1,i-k+2}^{k-1} & N_{0,i-k+2}^{k-1} \ vdots & N_{1,i-k+2}^{k-1} \ N_{k-2,i-k+2}^{k-1} & vdots \ 0 & N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} 0 & d_{0,i-k+2} & overbrace{0 cdots 0}^{k-2} \ 0 & d_{1,i-k+2} & 0 cdots 0 end{bmatrix} \ &+ cdots + egin{bmatrix} N_{0,i}^{k-1} & 0 \ N_{1,i}^{k-1} & N_{0,i}^{k-1} \ vdots & N_{1,i}^{k-1} \ N_{k-2,i}^{k-1} & vdots \ 0 & N_{k-2,i}^{k-1} end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-1} & d_{0,i} \ 0 cdots 0 & d_{1,i} end{bmatrix} \ &+ egin{bmatrix} N_{0,i-k+2}^{k-1} & 0 \ N_{1,i-k+2}^{k-1} & N_{0,i-k+2}^{k-1} \ vdots & N_{1,i-k+2}^{k-1} \ N_{k-2,i-k+2}^{k-1} & vdots \ 0 & N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} h_{0,i-k+1} & overbrace{0 cdots 0}^{k-1} \ h_{1,i-k+1} & 0 cdots 0 end{bmatrix} + egin{bmatrix} N_{0,i-k+3}^{k-1} & 0 \ N_{1,i-k+3}^{k-1} & N_{0,i-k+3}^{k-1} \ vdots & N_{1,i-k+3}^{k-1} \ N_{k-2,i-k+3}^{k-1} & vdots \ 0 & N_{k-2,i-k+3}^{k-1} end{bmatrix} egin{bmatrix} 0 & h_{0,i-k+2} & overbrace{0 cdots 0}^{k-2} \ 0 & h_{1,i-k+2} & 0 cdots 0 end{bmatrix} \ &+ cdots + egin{bmatrix} N_{0,i+1}^{k-1} & 0 \ N_{1,i+1}^{k-1} & N_{0,i+1}^{k-1} \ vdots & N_{1,i+1}^{k-1} \ N_{k-2,i+1}^{k-1} & vdots \ 0 & N_{k-2,i+1}^{k-1} end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-1} & h_{0,i} \ 0 cdots 0 & h_{1,i} end{bmatrix} \ &= egin{bmatrix} N_{0,i-k+1}^{k-1} & 0 \ N_{1,i-k+1}^{k-1} & N_{0,i-k+1}^{k-1} \ vdots & N_{1,i-k+1}^{k-1} \ N_{k-2,i-k+1}^{k-1} & vdots \ 0 & N_{k-2,i-k+1}^{k-1} end{bmatrix} egin{bmatrix} d_{0,i-k+1} & overbrace{0 cdots 0}^{k-1} \ d_{1,i-k+1} & 0 cdots 0 end{bmatrix} \ &+ egin{bmatrix} N_{0,i-k+2}^{k-1} & 0 \ N_{1,i-k+2}^{k-1} & N_{0,i-k+2}^{k-1} \ vdots & N_{1,i-k+2}^{k-1} \ N_{k-2,i-k+2}^{k-1} & vdots \ 0 & N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} h_{0,i-k+1} & d_{0,i-k+2} & overbrace{0 cdots 0}^{k-2} \ h_{1,i-k+1} & d_{1,i-k+2} & 0 cdots 0 end{bmatrix} \ &+ egin{bmatrix} N_{0,i-k+3}^{k-1} & 0 \ N_{1,i-k+3}^{k-1} & N_{0,i-k+3}^{k-1} \ vdots & N_{1,i-k+3}^{k-1} \ N_{k-2,i-k+3}^{k-1} & vdots \ 0 & N_{k-2,i-k+3}^{k-1} end{bmatrix} egin{bmatrix} 0 & h_{0,i-k+2} & d_{0,i-k+3} & overbrace{0 cdots 0}^{k-3} \ 0 & h_{1,i-k+2} & d_{1,i-k+3} & 0 cdots 0 end{bmatrix} \ &+ cdots \ &+ egin{bmatrix} N_{0,i}^{k-1} & 0 \ N_{1,i}^{k-1} & N_{0,i}^{k-1} \ vdots & N_{1,i}^{k-1} \ N_{k-2,i}^{k-1} & vdots \ 0 & N_{k-2,i}^{k-1} end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-2} & h_{0,i-1} & d_{0,i} \ 0 cdots 0 & h_{1,i-1} & d_{1,i} end{bmatrix} \ &+ egin{bmatrix} N_{0,i+1}^{k-1} & 0 \ N_{1,i+1}^{k-1} & N_{0,i+1}^{k-1} \ vdots & N_{1,i+1}^{k-1} \ N_{k-2,i+1}^{k-1} & vdots \ 0 & N_{k-2,i+1}^{k-1} end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-1} & h_{0,i} \ 0 cdots 0 & h_{1,i} end{bmatrix} \ end{align}]

[egin{align} &=^{(**)} egin{bmatrix} N_{0,i-k+2}^{k-1} \ N_{1,i-k+2}^{k-1} \ vdots \ N_{k-2,i-k+2}^{k-1} \ 0 end{bmatrix} egin{bmatrix} h_{0,i-k+1} & d_{0,i-k+2} & overbrace{0 cdots 0}^{k-2} end{bmatrix} + egin{bmatrix} 0 \ N_{0,i-k+2}^{k-1} \ N_{1,i-k+2}^{k-1} \ vdots \ N_{k-2,i-k+2}^{k-1} end{bmatrix} egin{bmatrix} h_{1,i-k+1} & d_{1,i-k+2} & overbrace{0 cdots 0}^{k-2} end{bmatrix} \ &+ egin{bmatrix} N_{0,i-k+3}^{k-1} \ N_{1,i-k+3}^{k-1} \ vdots \ N_{k-2,i-k+3}^{k-1} \ 0 end{bmatrix} egin{bmatrix} 0 & h_{0,i-k+2} & d_{0,i-k+3} & overbrace{0 cdots 0}^{k-3} end{bmatrix} + egin{bmatrix} 0 \ N_{0,i-k+3}^{k-1} \ N_{1,i-k+3}^{k-1} \ vdots \ N_{k-2,i-k+3}^{k-1} end{bmatrix} egin{bmatrix} 0 & h_{1,i-k+2} & d_{1,i-k+3} & overbrace{0 cdots 0}^{k-3} end{bmatrix} \ &+ cdots \ &+ egin{bmatrix} N_{0,i}^{k-1} \ N_{1,i}^{k-1} \ vdots \ N_{k-2,i}^{k-1} \ 0 end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-2} & h_{0,i-1} & d_{0,i} end{bmatrix} + egin{bmatrix} 0 \ N_{0,i}^{k-1} \ N_{1,i}^{k-1} \ vdots \ N_{k-2,i}^{k-1} end{bmatrix} egin{bmatrix} overbrace{0 cdots 0}^{k-2} & h_{1,i-1} & d_{1,i} end{bmatrix} \ &= egin{bmatrix} N_{0,i-k+2}^{k-1} & N_{0,i-k+3}^{k-1} & cdots & N_{0,i}^{k-1} \ N_{1,i-k+2}^{k-1} & N_{1,i-k+3}^{k-1} & cdots & N_{1,i}^{k-1} \ vdots & vdots & & vdots \ N_{k-2,i-k+2}^{k-1} & N_{k-2,i-k+3}^{k-1} & cdots & N_{k-2,i}^{k-1} \ 0 & 0 & cdots & 0 end{bmatrix} egin{bmatrix} h_{0,i-k+1} & d_{0,i-k+2} & & & 0 \ & h_{0,i-k+2} & d_{0,i-k+3} & & \ & & ddots & ddots & \ 0 & & & h_{0,i-1} & d_{0,i} \ end{bmatrix} \ &+ egin{bmatrix} 0 & 0 & cdots & 0 \ N_{0,i-k+2}^{k-1} & N_{0,i-k+3}^{k-1} & cdots & N_{0,i}^{k-1} \ N_{1,i-k+2}^{k-1} & N_{1,i-k+3}^{k-1} & cdots & N_{1,i}^{k-1} \ vdots & vdots & & vdots \ N_{k-2,i-k+2}^{k-1} & N_{k-2,i-k+3}^{k-1} & cdots & N_{k-2,i}^{k-1} end{bmatrix} egin{bmatrix} h_{1,i-k+1} & d_{1,i-k+2} & & & 0 \ & h_{1,i-k+2} & d_{1,i-k+3} & & \ & & ddots & ddots & \ 0 & & & h_{1,i-1} & d_{1,i} \ end{bmatrix} \ &= egin{bmatrix} mathbf{M}^{k-1}(i) \ mathbf{0}^T end{bmatrix} egin{bmatrix} h_{0,i-k+1} & d_{0,i-k+2} & & & 0 \ & h_{0,i-k+2} & d_{0,i-k+3} & & \ & & ddots & ddots & \ 0 & & & h_{0,i-1} & d_{0,i} \ end{bmatrix} \ &+ egin{bmatrix} mathbf{0}^T \ mathbf{M}^{k-1}(i) end{bmatrix} egin{bmatrix} h_{1,i-k+1} & d_{1,i-k+2} & & & 0 \ & h_{1,i-k+2} & d_{1,i-k+3} & & \ & & ddots & ddots & \ 0 & & & h_{1,i-1} & d_{1,i} \ end{bmatrix} \ &= egin{bmatrix} mathbf{M}^{k-1}(i) \ mathbf{0}^T end{bmatrix} egin{bmatrix} 1-d_{0,i-k+2} & d_{0,i-k+2} & & & 0 \ & 1-d_{0,i-k+3} & d_{0,i-k+3} & & \ & & ddots & ddots & \ 0 & & & 1-d_{0,i} & d_{0,i} \ end{bmatrix} \ &+ egin{bmatrix} mathbf{0}^T \ mathbf{M}^{k-1}(i) end{bmatrix} egin{bmatrix} -d_{1,i-k+2} & d_{1,i-k+2} & & & 0 \ & -d_{1,i-k+3} & d_{1,i-k+3} & & \ & & ddots & ddots & \ 0 & & & -d_{1,i} & d_{1,i} \ end{bmatrix} end{align}]

(**): 为什么第一项与最后一项都为 0？这个等号前方的式子第一项与最后一项的“(N) 的部分”分别是 (B_{i-k+1,k-1}(u), B_{i+1,k-1}(u)) 的系数部分，是 (k-1) 阶的系数。结合公式 (3*)，知道 (k) 阶的区间涉及范围是 ([i-k+1, i])。所以 (k-1) 阶的区间涉及范围是 ([i-k+2, i])，没有包含 (i-k+1,i+1)，所以在 (k-1) 阶时这两项为 0。

(mathbf{M}^k(i)) 的初始值，一阶只含有一个 (N^1_{0,i})，(mathbf{M}^1(i)=egin{bmatrix}1end{bmatrix})，即常值函数。

4. Uniform B-splines 的 (mathbf{M}^k(i)) 通项公式

Uniform B-splines 的每一个区间的长度相等，假设为 1。可以将 (mathbf{M}^k(i)) 的 (d_{0,*}, d_{1,*}) 写成分数，并且与 (i) 无关系，所以不管有多少结点，都只用算一次 (mathbf{M}^k(i))，也就可以直接写作 (mathbf{M}^k)。得到的结果如公式 (7*)。

[egin{align} mathbf{M}^k &= {1 over k-1} left{ egin{bmatrix} mathbf{M}^{k-1} \ mathbf{0}^T end{bmatrix} egin{bmatrix} 1 & k-2 & & & 0 \ & 2 & k-3 & & \ & & ddots & ddots & \ 0 & & & k-1 & 0 \ end{bmatrix} + egin{bmatrix} mathbf{0}^T \ mathbf{M}^{k-1} end{bmatrix} egin{bmatrix} -1 & 1 & & & 0 \ & -1 & 1 & & \ & & ddots & ddots & \ 0 & & & -1 & 1 \ end{bmatrix} ight} end{align}]

通项公式怎么计算，我不知道。

论文之后计算了一些偏导，但是符号系统实在让人看不懂。需要用的时候自然能推导出。

Reference

[1] Qin, Kaihuai. "General matrix representations for B-splines." The Visual Computer 16.3-4 (2000): 177-186.

[2] 【回形针PaperClip】如何设计一个逼真的三维模型.

[3] Wikipedia: De Boor's algorithm.

[4] Wikipedia: Fibonacci number 3.3 Matrix form.

[5] Unknown: Matrixces and Polynomials.

[6] De Boor, Carl. B (asic)-Spline Basics. No. MRC-TSR-2952. WISCONSIN UNIV-MADISON MATHEMATICS RESEARCH CENTER, 1986.

[7] Solomon, Justin. Numerical algorithms: methods for computer vision, machine learning, and graphics. CRC press, 2015.