曲线曲面的插值拟合

曲线参数化

如何选择参数

B样条插值/拟合的输入通常一组已知的数据点，曲线参数化第一步是找到一组参数，能够将这些点"固定"在某些特定值上。

例如，如果数据点是$\mathbf{D}_0, \ldots, \mathbf{D}_n$，那么需要在曲线的定义域中找到$n+1$个参数 $t_0, \ldots, t_n$，使得数据点$\mathbf{D}_k$对应于参数$t_k$，其中$0 \leq k \leq n$。如果$\mathbf{C}(u)$是一条通过所有数据点的曲线，那么对于所有$0 \leq k \leq n$，我们都有$\mathbf{D}_k = \mathbf{C}(t_k)$。

在下面的图中，我们有7个数据点$(n=6)$，因此需要7个参数来建立数据点和参数的对应关系。

参数的选择有无限多可能，我们可以均匀地划分定义域，或者从定义域中随机选择$n+1$个值。但是，选择不当的参数会导致不可预料的结果。

下图显示了4个数据点和三条插值曲线，每条曲线都使用了不同的参数。其中一条曲线向外侧弯曲，产生了凸起。黑色曲线有一个小的凸起，只有这两者之间的曲线比较符合数据点的分布趋势。因此，参数的选择会影响曲线的形状，从而相应地影响曲线的参数化。

下面，我们将讨论一些参数选择的方法，包括均匀分布法、弦长法和向心法。在确定好一组参数之后，就可以计算相应的节点向量。

关于曲线插值、曲线逼近和曲线拟合的区别

曲线插值、曲线逼近和曲线拟合都是用于构造一条平滑曲线来描述一组离散数据点的方法,但它们的目的和原理有所不同。

1. 曲线插值(Curve Interpolation)
   • 目的是找到一条平滑曲线精确地通过所有给定的数据点。
   • 常用的插值方法包括拉格朗日插值、牛顿插值、三次样条插值等。
   • 例如,给定一组坐标点$(x_1,y_1), (x_2,y_2), ..., (x_n,y_n)$,我们可以使用拉格朗日插值多项式来构造一条曲线精确地通过所有这些点。
2. 曲线逼近(Curve Approximation)
   • 目的是在某种意义下(如最小二乘法)找到最佳拟合曲线,使曲线和数据点之间的偏差最小。
   • 常用的逼近方法包括最小二乘法拟合多项式、指数函数、三角函数等。
   • 例如,给定一组测量数据点,我们可以用最小二乘法拟合一条直线或抛物线等简单函数来近似描述这些数据。
3. 曲线拟合(Curve Fitting)
   • 目的是寻找一个适当的函数形式来很好地拟合给定的数据点，使得总误差（根据不同问题有不同的误差距离的衡量方式，一般是均方误差）最小的那个简单函数。
   • 常用的函数形式包括多项式、高斯函数、指数函数、对数函数等参数化模型。
   • 例如,给定一组有规律的数据点,我们可以尝试用一个高斯函数或对数函数等非线性模型来很好地拟合这些数据。

插值和逼近更注重于构造一条平滑曲线来精确地通过或逼近给定的离散数据点，而拟合则更强调寻找一个合适的参数化模型函数来描述潜在的数据生成规律，从广义上来讲，曲线拟合包含曲线插值和曲线逼近两种情况，当需要拟合曲线精确经过所有给定数据点时,可以采用曲线插值的方法；当允许拟合曲线与数据点之间存在一定偏差时,可以采用曲线逼近的方法；还可以采用其他参数化模型函数(如高斯函数、对数函数等)进行非线性曲线拟合,来描述数据的内在规律。

线性回归

线性回归是曲线拟合的一种特殊情况,其拟合函数形式是一条直线(一次函数)。在线性回归中,我们试图找到一条最佳拟合直线,使得数据点到直线的残差平方和最小，这个过程就是对线性函数进行曲线拟合。

当使用线性回归模型时, 我们实际上是在用一条直线来逼近数据分布，这种用最小二乘法找到最佳拟合直线的过程,本质上属于曲线逼近的一种方法。由于线性回归模型并不要求拟合直线精确通过所有数据点，因此不属于曲线插值的范畴，但是,如果数据点刚好在一条直线上，那么线性回归得到的结果就是一个插值问题的解。

曲线均匀参数化

最简单的参数选择方法是均匀参数化。假设定义域为$[0,1]$，需要$n+1$个均匀分布的参数，那么第一个和最后一个参数必须是0和1，因为我们希望曲线通过第一个和最后一个数据点。因此，我们有$t_0 = 0$和$t_n = 1$。

$n+1$个参数将区间$[0,1]$均匀地分成$n$个子区间，每个子区间的长度是$1/n$，划分后的参数是$0$, $1/n$, $2/n$, $3/n$, …, $(n-1)/n$ 和 $1$。因此，我们有：

$$t_i = \frac{i}{n} \quad i = 0, 1, \ldots, n$$

例如，如果我们需要5个参数，$n = 4$，那么均匀分布的参数是0, $1/4$, $1/2$, $3/4$ 和 1。如果我们需要8个参数，$n = 7$，那么均匀分布的参数是0, $1/7$, $2/7$, $3/7$, $4/7$, $5/7$, $6/7$ 和 1。

如果定义域是$[a, b]$而不是$[0,1]$呢？在这种情况下，$[a, b]$被$n+1$个划分点（包括$a$和$b$）分成$n$个区间。由于这个区间的长度是$b - a$，每个子区间的长度是$(b - a)/n$。因此，划分点（即参数）是：

$$t_i = a + \frac{b - a}{n} \cdot i \quad i = 0, 1, \ldots, n$$

尽管均匀参数化很简单，但它会产生一些不稳定的结果。

例如，当数据点不是均匀分布时，使用均匀分布的参数可能会产生异常形状，如凸起、尖锐点和自相交形成的环。

在下面的左图中，在数据点3处有一个环。在右图中，曲线在数据点1、2和3之间出现剧烈摆动。虽然不能说这些问题是均匀参数化引起的，但它确实比使用其他方法时更容易出现问题。

曲线弦长参数化

如果插值曲线和数据点连成的多边形非常接近，那么两个相邻数据点之间的曲线段长度将非常接近这两个数据点相连的弦长，并且插值曲线的总长度也非常接近数据点连成的多边形的总长度。

在下图中，插值曲线的每一段长度都非常接近其对应的弦的长度，曲线的总长度接近数据点连成的多边形的总长度。因此，如果根据弦长对定义域进行划分，划分后的节点参数与弧长参数化产生的值近似，这就是弦长参数化的优点。

假设数据点是 $\mathbf{D}_0, \mathbf{D}_1, \ldots, \mathbf{D}_n$。$\mathbf{D}_{i-1}$ 和 $\mathbf{D}_i$ 之间的长度为 $|\mathbf{D}_i - \mathbf{D}_{i-1}|$，数据点连成的多边形长度是这些弦的长度之和：

$$L = \sum_{i=1}^{n} |\mathbf{D}_i - \mathbf{D}_{i-1}|$$

因此，从点 $\mathbf{D}_0$ 到点 $\mathbf{D}_k$ 的弦长，记作 $L_k$，与多边形总长度 $L$ 的比率是：

$$L_k=\frac{ \sum_{i=1}^{k} |\mathbf{D}_i - \mathbf{D}_{i-1}|}{L}$$

如果我们要对插值曲线进行弦长参数化，那么定义域必须根据比率 $L_k$ 进行划分。如果定义域是 $[0,1]$，那么参数 $t_k$ 值为：

$$\begin{aligned} &t_{0}=0\\ &t_{k}=\frac{1}{L}\left(\sum_{i=1}^{k}|\mathbf{D}_{i}-\mathbf{D}_{i-1}|\right)\\ &t_{n}=1 \end{aligned}$$

其中 $L$ 是数据点连成的多边形的总长度。这样，就根据弦长比率划分定义域，从而分配参数值

看一个例子：假设我们有4个数据点 $(n = 3)$：$\mathbf{D}_0 = (0,0)$，$\mathbf{D}_1 = (1,2)$，$\mathbf{D}_2 = (3,4)$ 和 $\mathbf{D}_3 = (4,0)$，每条弦的长度是：

$$|\mathbf{D}_1 - \mathbf{D}_0| = \sqrt{5} = 2.236\\ |\mathbf{D}_2 - \mathbf{D}_1| = \sqrt{8} = 2.828\\ |\mathbf{D}_3 - \mathbf{D}_2| = \sqrt{17} = 4.123$$

总长度是：

$$L = 2.236 + 2.828 + 4.123 = 9.187$$

最后，我们得到相应的参数：

$$\begin{array}{rcl} t_0&=&0\\ t_1&=&\frac{|\mathbf{D}_1-\mathbf{D}_0|}{L}=0.2434\\ t_2&=&\frac{|\mathbf{D}_1-\mathbf{D}_0|+|\mathbf{D}_2-\mathbf{D}_1|}{L}=0.5512\\ t_3&=&1 \end{array}$$

下图显示了使用均匀参数化和弦长参数化得到的参数分布

如果定义域是 $[a, b]$ 而不是 $[0,1]$ 怎么办？上式中 $L_k$ 的值在0 到1 之间。由于 $[a, b]$ 的长度是 $b - a$，所以 $L_k(b - a)$（其中$0 \leq k \leq n$) 将 $[0, b - a]$ 按照累计的弦长比例进行划分。因此，在 $[a, b]$上划分后得到的参数值为：

$$\begin{aligned} &t_{0}=a\\ &t_{k}=a+L_{k}(b-a)\\ &t_{n}=b \end{aligned}$$

对多项式曲线使用弦长参数化效果并不完美,弦长只是一个近似值。有时，较长的弦可能会导致其曲线段出现较大隆起。在下图中，黑色和蓝色曲线都使用7个数据点进行插值计算。

两条曲线的形状非常相似，除了最后一段，使用弦长参数化插值的曲线具有较大的摆动，最后几段曲线段截然不同，使用弦长法插值的蓝色曲线与使用均匀法插值的红色曲线相比，有很大的凸起和扭曲。这是弦长法经常出现的问题。

曲线向心参数化

向心参数化由E. T. Y. Lee提出，假设我们驾驶一辆汽车穿过一个弯道绕行的赛道。在急转弯处,我们必须非常小心,以免产生过大的法向加速度(即离心力)，否则我们的车辆可能会失控。为了安全驾驶，E. T. Y. Lee建议沿路径的法向力应与角度变化成正比。向心参数化方法实际上是对这一模型的近似。我们可以将向心参数化当作是对弦长参数化方法的一种扩展。

假设数据点为 $\mathbf{D}_0, \mathbf{D}_1, \ldots, \mathbf{D}_n$，我们取一个正的"幂"值 $a$，通常$a = 1/2$ 。两个相邻数据点之间的距离通过 $|\mathbf{D}_k - \mathbf{D}_{k-1}|^a$ 而不是 $|\mathbf{D}_k - \mathbf{D}_{k-1}|$ 来计算

因此，数据点连接成的多边形的总长度是：

$$L = \sum_{i=1}^{n} |\mathbf{D}_i - \mathbf{D}_{i-1}|^a$$

多边形上从 $\mathbf{D}_0$ 到 $\mathbf{D}_k$ 的长度与总长度的比是：

$$L_k=\frac{\sum_{i=1}^k|\mathbf{D}_i-\mathbf{D}_{i-1}|^a}{L}$$

因此，使用新的弦长计算方法，得到的累计弦长比例为 $L_0 = 0, L_1, \ldots, L_n = 1$ ，假设定义域为$[0,1]$，划分后的曲线参数为：

$$\begin{array}{rcl} t_0&=&0\\ t_k&=&\frac{1}{L}\left(\sum_{i=1}^k|\mathrm D_i-\mathrm D_{i-1}|^a\right)\\ t_n&=&1 \end{array}$$

如果 $a = 1$，那么向心参数化退化为弦长法。如果 $a < 1$，例如 $a = 1/2$，那么$|\mathbf{D}_k - \mathbf{D}_{k-1}|^a$ 的值会小于 $|\mathbf{D}_k - \mathbf{D}_{k-1}|$。那么，较长的弦（长度大于1的弦）对数据点连成的多边形的总长度的贡献会减小，较短的弦（长度小于1的弦）对数据点连成的多边形的总长度的贡献会增加。因此，向心参数化比弦长法能够更好地处理具有较大曲率的曲线段。

下面重新计算弦长参数化中的例子，我们有四个数据点 $(n = 3)$：$\mathbf{D}_0 = (0,0)$，$\mathbf{D}_1 = (1,2)$，$\mathbf{D}_2 = (3,4)$ 和 $\mathbf{D}_3 = (4,0)$。取 $a = 1/2$，那么每条弦的长度是：

$$|\mathbf{D}_1-\mathbf{D}_0|^{1/2}=\sqrt{\sqrt{5}}=1.495\\|\mathbf{D}_2-\mathbf{D}_1|^{1/2}=\sqrt{2\sqrt{2}}=1.682\\|\mathbf{D}_3-\mathbf{D}_2|^{1/2}=\sqrt{\sqrt{17}}=2.031$$

弦的总长度是：

$$L=\sqrt{\sqrt{5}}+\sqrt{2\sqrt{2}}+\sqrt{\sqrt{17}}=5.208$$

因此，参数是：

$$\begin{array}{rcl} t_0&=&0\\ t_1&=&\frac{|\mathbf{D}_1-\mathbf{D}_0|^{1/2}}{L}=0.2871\\ t_2&=&\frac{|\mathbf{D}_1-\mathbf{D}_0|^{1/2}+|\mathbf{D}_2-\mathbf{D}_1|^{1/2}}{L}=0.6101\\ t_3&=&1 \end{array}$$

下面给出了使用均匀参数化、弦长参数化和向心参数化计算的三组参数的分布情况：

下面看一个极端的例子，下图显示了使用均匀参数化（黑色）、弦长参数化（蓝色）和向心参数化（红色）对4个数据点进行插值得到的B样条曲线。均匀参数化得到的曲线有一个峰值，弦长参数化得到的曲线有两个大的凸起，向心参数化则很好地处理了两个距离非常接近的数据点（1号点和2号点）。

那么可以说向心参数化比其他两种方法更好吗？下面是一个反例。我们有7个数据点，黑色、蓝色和红色曲线是使用均匀参数化、弦长参数化和向心参数化插值得到的曲线。如图所示，均匀参数化产生了一个非常稳定的插值曲线，向心参数化得到的曲线波动大于均匀参数化得到的曲线，而弦长参数化产生的曲线波动幅度最大。

曲线通用方法参数化

1999年，Choong-Gyoo Lim提出了一种有趣的参数化方法。在之前讨论的方法中，我们先确定参数，然后计算一个节点向量。Lim提出的方法称为通用方法，通过使用均匀分布的节点来计算参数，正好相反。

假设我们需要$n+1$个参数，每个数据点对应一个参数，插值的B样条曲线次数是$p$。那么，节点的数量是$m+1$，满足$m = n + p + 1$。Lim认为这些节点应该是均匀分布的。前$p+1$个节点被设置为0，后$p+1$个节点被设置为1，其余的$n-p$个节点均匀的分布在定义域$[0,1]$之间。因此，节点为：

$$\begin{aligned} u_{0}&=\quad u_1=\cdots=u_p=0\\ u_{p+i}&=\quad\frac i{n-p+1}\quad\mathrm{for~}i=1,2,\ldots,n-p\\ u_{m-p}&=\quad u_{m-p+1}=\cdots=u_m=1 \end{aligned}$$

这样我们就有了 $n+1$个 B样条基函数。然后，选择对应基函数达到最大值时的参数作为数据点的参数。如下图所示，其中$n = 6$（7个数据点），$p = 4$，和$m = 11$（12个节点）。由于使用的是clamped类型节点，0和1是重数为5的两个节点，并且只有两个内部节点，分别为1/3和2/3。有$n+1$个B样条基函数，第一个和最后一个的最大值分别为0和1，其他基函数的最大值用竖直蓝色线段标记，对应的参数用黄色点标记出。

看一个例子。假设我们有4个数据点（$n=3$），次数$p=2$。因此，节点的数量是7（即$m = n + p + 1 = 6$）。由于节点是均匀分布的，它们是：

$$u_0 = u_1 = u_2 = 0, \quad u_3 = 0.5, \quad u_4 = u_5 = u_6 = 1$$

然后，我们可以计算B样条基函数。从第0阶开始：

接下来，我们计算次数为1的基函数。由于$N_{0,0}(u)$和$N_{1,0}(u)$都是0，因此$N_{0,1}(u)$在任何区间都是零。同样，$N_{4,1}(u)$在任何地方也都是零。因此次数为1的基函数，只有 $N_{1,1}(u)$、$N_{2,1}(u)$和$N_{3,1}(u)$不为0，如下所示：

次数为2的基函数如下：

下面的图显示了所有四个次数为2的B样条基函数。

不难算出，$N_{0,2}(u)$、$N_{1,2}(u)$、$N_{2,2}(u)$和$N_{3,2}(u)$的最大值分别为1（在$u = 0$处）、2/3（在$u = 1/3$处）、2/3（在$u = 2/3$处）和1（在$u = 1$处）。

因此，使用通用方法，得到的节点向量是${0, 0, 0, 0.5, 1, 1, 1}$，数据点对应的参数是${0, 1/3, 2/3, 1}$。在这种情况下，我们的参数是均匀分布的。

使用通用参数化方法插值得到的曲线有一个非常有用的性质。它是仿射不变的。这意味着，可以通过变换数据点来获得变换后的插值B样条曲线。这与B样条曲线的仿射不变性质类似。

如果在原始插值曲线和变换后的插值B样条曲线中使用相同的节点和参数集，那么通过变换数据点来实现曲线的变换。根据这一点，我们知道均匀参数化也是仿射不变的，因为节点向量是从一组均匀分布的参数计算出来的，这些参数在变换前后没有改变。然而，弦长参数化和向心参数化不是仿射不变的，因为变换后的曲线中，弦长分布可能与原始的不一样，我们会得到一组新的弦长，必须重新计算新的参数。

曲面参数化

由$e+1$行$f+1$列控制点定义的次数为$(p, q)$的B样条曲面的方程如下：

$$\mathbf{S}(u,v) = \sum_{i=0}^{e} \sum_{j=0}^{f} N_{i,p}(u)N_{j,q}(v)\mathbf{P}_{ij}$$

它需要两组参数（参数点）来进行曲面插值和逼近。

假设我们有$m+1$行$n+1$列数据点$\mathbf{D}_{ij}$，其中$0 \leq i \leq m,0 \leq j \leq n$。因此，在$u$方向需要$m+1$个参数$s_0, ..., s_m$，$v$方向需要$n+1$个参数$t_0, ..., t_n$，这样，曲面参数域中的点$(s_c, t_d)$与曲面上的点$\mathbf{S}(s_c, t_d)$对应，用方程表示如下：

$$\mathbf{S}(s_c, t_d) = \sum_{i=0}^{e} \sum_{j=0}^{f} N_{i,p}(s_c)N_{j,q}(t_d)\mathbf{P}_{ij}$$

而$\mathbf{S}(s_c, t_d)$又与数据点$\mathbf{D}_{cd}$相对应，其中 $s_c$ 和 $t_d$ 分别是 $u$ 和 $v$ 方向的参数。

在B样条曲面方程中，$u$方向对应于$N_{i,p}(u)$和$\mathbf{P}_{ij}$中的索引$i$。 $i$ 的范围为从0到$m$，$N_{0,p}(u), N_{1,p}(u), ..., N_{m,p}(u)$是控制点的基函数。因此，在$u$方向，我们需要$m+1$个参数，已知次数 $p$ 和第 $j$ 列数据点，我们可以计算$m+1$个参数$u_{0,j}, u_{1,j}, ..., u_{m,j}$，如下图所示。参数$s_0, s_1, ..., s_m$是每行的平均值。即参数 $s_i$ 是第 $i$ 行参数的平均值，即$s_i = (u_{i,0} + u_{i,1} + ... + u_{i,n})/(n+1)$。
在B样条曲面方程中，$u$方向对应于$N_{i,p}(u)$和$\mathbf{P}_{ij}$中的索引$i$。 $i$ 的范围为从0到$m$，$N_{0,p}(u), N_{1,p}(u), ..., N_{m,p}(u)$是控制点的基函数。因此，在$u$方向，我们需要$m+1$个参数，已知次数 $p$ 和第 $j$ 列数据点，我们可以计算$m+1$个参数$u_{0,j}, u_{1,j}, ..., u_{m,j}$，如下图所示。参数$s_0, s_1, ..., s_m$是每行的平均值。即参数 $s_i$ 是第 $i$ 行参数的平均值，即$s_i = (u_{i,0} + u_{i,1} + ... + u_{i,n})/(n+1)$。

$v$方向的参数计算类似。每行有$n+1$个数据点，因此需要$n+1$个参数。因此，对于第 $i$ 行的数据点，我们可以计算出$n+1$个参数值$v_{i,0}, v_{i,1}, ..., v_{i,n}$，有$m+1$行，这些值可以组织成一个$(m+1) \times (n+1)$矩阵，参数 $t_j$ 是第 $j$ 列参数的平均值，即 $t_j = (v_{0,j} + v_{1,j} + ... + v_{m,j})/(m+1)$。这样，我们可以计算出 $n+1$ 个 $v$ 方向的参数。

算法总结如下：

Input： $(m+1) \times (n+1)$ 个数据点 $\mathbf{D}_{ij}$；
Output： $u$方向的参数 $s_0, ..., s_m$ 和 $v$方向的参数 $t_0, ..., t_n$；

算法：

​ // 计算 $s_0, ..., s_m$
​ for $j$ = $0$ to $n$ do
​ 计算一组 $m+1$ 个参数 $u_{0,j}, u_{1,j}, ..., u_{m,j}$；
​ for $i$ = $0$ to $m$ do
​ $s_i = (u_{i,0} + u_{i,1} + ... + u_{i,n})/(n+1)$；
​ （获得$u$方向的参数）

​ // 计算 $t_0, ..., t_n$
​ for $i$ = $0$ to $m$ do
​ 计算一组 $n+1$ 个参数 $v_{i,0}, v_{i,1}, ..., v_{i,n}$；
​ for $j$ = $0$ to $n$ do
​ $t_j = (v_{0,j} + v_{1,j} + ... + v_{m,j})/(m+1)$；
​ （获得$v$方向的参数）

根据参数 $s_0, s_1, ..., s_m$ 和次数 $p$ ，我们可以计算出$u$方向的节点向量 $U$，根据参数 $t_0, t_1, ..., t_n$ 和次数 $q$ ，我们可以计算出$v$方向的节点向量 $V$。

注意，上述只是一个概念性算法，效率不高。这个算法只适用于均匀参数化、弦长参数化和向心参数化。对于通用参数化方法，因为不涉及数据点，我们可以对每一行和每一列的数据点应用均匀分布的节点来计算参数。

求解线性方程组

插值和逼近的计算过程中都会涉及到求解线性方程组，下面我们讨论一个常用的求解方法。

假设我们有一个 $n \times n$ 的矩阵 $\mathbf{A}$，一个 $n \times h$ 的"常数"项矩阵 $\mathbf{B}$，以及一个 $n \times h$ 的未知矩阵 $\mathbf{X}$，定义如下：

$$\mathbf{A} = \left[\begin{array}{cccc}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{array}\right]_{n\times n}$$

$$\mathbf{B} = \left[\begin{array}{cccc}b_{11}&b_{12}&\cdots&b_{1h}\\b_{21}&b_{22}&\cdots&b_{2h}\\\vdots&\vdots&\ddots&\vdots\\b_{n1}&b_{n2}&\cdots&b_{nh}\end{array}\right]_{n\times h}$$

$$\mathbf{X} = \left[\begin{array}{cccc}x_{11}&x_{12}&\cdots&x_{1h}\\x_{21}&x_{22}&\cdots&x_{2h}\\\vdots&\vdots&\ddots&\vdots\\x_{n1}&x_{n2}&\cdots&x_{nh}\end{array}\right]_{n\times h}$$

它们满足以下关系：

$$\mathbf{B} = \mathbf{A}\cdot\mathbf{X}$$

如果 $\mathbf{A}$ 和 $\mathbf{B}$ 已知，我们需要一种快速的方法来求解 $\mathbf{X}$。有人可能会说：计算矩阵 $\mathbf{A}$ 的逆矩阵 $\mathbf{A}^{-1}$，然后解就是 $\mathbf{X}$ = $\mathbf{A}^{-1}\mathbf{B}$，虽然这是正确的思路，但是这样计算量很大。

LU分解

一个有效的求解 $\mathbf{B}$ = $\mathbf{A}\mathbf{X}$ 的方法是LU分解。虽然高斯消元法和Cholesky等其他方法也能计算，但LU分解是一种加速计算的方法。

LU分解首先将矩阵 $\mathbf{A}$ “分解” 为 $\mathbf{A} = \mathbf{L} \mathbf{U}$的形式，其中 $\mathbf{L}$ 和 $\mathbf{U}$ 分别是下三角和上三角矩阵。如果 $\mathbf{A}$ 是一个 $n \times n$ 矩阵，那么 $\mathbf{L}$ 和 $\mathbf{U}$ 也是 $n \times n$ 矩阵，其形式如下：

$$\mathbf{L} = \left[\begin{array}{ccccc}l_{11}&0&0&\cdots&0\\l_{21}&l_{22}&0&\cdots&0\\\vdots&\vdots&\ddots&\ddots&\vdots\\l_{n1}&l_{n2}&l_{n3}&\cdots&l_{nn}\end{array}\right]\quad\mathbf{U} = \left[\begin{array}{ccccc}u_{11}&u_{12}&u_{13}&\cdots&u_{1n}\\0&u_{22}&u_{23}&\cdots&u_{2n}\\0&0&u_{33}&\cdots&u_{3n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\cdots&u_{nn}\end{array}\right]$$

下三角矩阵 $\mathbf{L}$ 在对角线上方的所有元素都为零，上三角矩阵 $\mathbf{U}$ 在对角线下方的所有元素都为零。

如果找到了满足 $\mathbf{A} = \mathbf{L} \mathbf{U}$ 的$\mathbf{L} \mathbf{U}$分解，则原始方程变为 $\mathbf{B} = (\mathbf{L} \mathbf{U}) \mathbf{X}$，根据矩阵的结合律，这个方程可以写为 $\mathbf{B} = \mathbf{L} (\mathbf{U} \mathbf{X})$。

由于 $\mathbf{L}$ 和 $\mathbf{B}$ 是已知的，那么上式方程相当于求解 $\mathbf{B} = \mathbf{L} \mathbf{Y}$，其中 $\mathbf{Y} = \mathbf{U} \mathbf{X}$。然后，由于 $\mathbf{U}$ 和 $\mathbf{Y}$ 是已知的，因此可以求解出 $\mathbf{X}$，就得到了最终的结果。

这样，求解方程$\mathbf{B} = \mathbf{A}\mathbf{X}$ 分解为两步：

1. 从 $\mathbf{B} = \mathbf{L} \mathbf{Y}$ 中解出 $\mathbf{Y}$
2. 从 $\mathbf{Y} = \mathbf{U} \mathbf{X}$ 中解出 $\mathbf{X}$

前向代换

第一步，展开 $\mathbf{B} = \mathbf{L} \mathbf{Y}$ 得到：

$$\begin{bmatrix} b_{11} & b_{12} & \cdots & b_{1h} \\ b_{21} & b_{22} & \cdots & b_{2h} \\ \vdots & \vdots & \ddots & \vdots \\ b_{n1} & b_{n2} & \cdots & b_{nh} \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{bmatrix} \begin{bmatrix} y_{11} & y_{12} & \cdots & y_{1h} \\ y_{21} & y_{22} & \cdots & y_{2h} \\ \vdots & \vdots & \ddots & \vdots \\ y_{n1} & y_{n2} & \cdots & y_{nh} \end{bmatrix}$$

不难看出矩阵 $\mathbf{B}$ 的第 $j$ 列是矩阵 $\mathbf{L}$ 和矩阵 $\mathbf{Y}$ 的第 $j$ 列的乘积，因此，我们一次可以求解出 $\mathbf{Y}$ 中的一列，如下所示：

$$\left[\begin{array}{c}b_1\\b_2\\\vdots\\b_n\end{array}\right] = \left[\begin{array}{ccccc}l_{11}&0&0&\cdots&0\\l_{21}&l_{22}&0&\cdots&0\\\vdots&\vdots&\ddots&\ddots&\vdots\\l_{n1}&l_{n2}&l_{n3}&\cdots&l_{nn}\end{array}\right]\cdot\left[\begin{array}{c}y_1\\y_2\\\vdots\\y_n\end{array}\right]$$

上式相当于方程：

$$\begin{array}{rcl} b_1&=&l_{11}y_1\\ b_2&=&l_{21}y_1&+&l_{22}y_2\\ &\vdots&\vdots\\ b_n&=&l_{n1}y_1&+&l_{n2}y_2&+&l_{n3}y_3&+&\cdots&+&l_{nn}y_n \end{array}$$

从上面的等式中，我们可以算出 $y_1 = \frac{b_1}{l_{11}}$，一旦有了 $y_1$，代入第二个方程就可以得出 $y_2 = \frac{b_2 - l_{21}y_1}{l_{22}}$，有了 $y_1$ 和 $y_2$，代入第三个方程中，我们得出 $y_3 = \frac{b_3 - (l_{31}y_1 + l_{32}y_2)}{l_{33}}$。

因此，我们从第一个方程计算 $y_1$ 并将其代入第二个方程，计算出 $y_2$。有了 $y_1$ 和 $y_2$，将它们代入第三个方程求解出 $y_3$。依此类推，当计算第 $i$ 方程时，$y_1, y_2, \ldots, y_{i-1}$ 可以通过前$i-1$个方程计算出，将它们代入第 $i$ 个方程进行求解，下面是求解 $y_i$ 的公式：

$$y_i = \frac{1}{l_{ii}}\left[b_i-\sum_{k=1}^{i-1}l_{i,k}y_k\right]$$

由于 $y_i$ 的值会被代入下一个方程以求解 $y_{i+1}$ 的值，这个过程被称为前向代换。重复以上过程，可以求解出向量 $\mathbf{Y}$。

下面是算法的流程：

Input： 矩阵 $\mathbf{B}_{n \times h}$ 和下三角矩阵 $\mathbf{L}_{n \times n}$
Output： 矩阵 $\mathbf{Y}_{n \times h}$

算法：

    for $j$ = $1$ to $h$ do // 共有 $h$ 列
        begin // 对每一列执行以下操作
            $y_{1,j} = b_{1,j} / l_{1,1}$; // 计算当前列的 $y_1$
            for $i$ = $2$ to $n$ do // 处理该列的元素
                begin
                    $sum = 0$; // 求解当前列的 $y_i$
                    for $k$ = $1$ to $i-1$ do
                                $sum = sum + l_{i,k} \times y_{k,j}$;
                    $y_{i,j} = (b_{i,j} - sum) / l_{i,i}$;
                end
    end

后向代换

求解出$\mathbf{Y}$后， 我们可以根据 $\mathbf{Y} = \mathbf{U} \mathbf{X}$ 求解出 $\mathbf{X}$。将这个等式展开，只看$\mathbf{Y}$中的某一列和 $\mathbf{X}$ 的相应列，得到：

$$\left[\begin{array}{c}y_1\\y_2\\\vdots\\y_n\end{array}\right]=\left[\begin{array}{ccccc}u_{11}&u_{12}&u_{13}&\cdots&u_{1n}\\0&u_{22}&u_{23}&\cdots&u_{2n}\\0&0&u_{33}&\cdots&u_{3n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\cdots&u_{nn}\end{array}\right]\cdot\left[\begin{array}{c}x_1\\x_2\\\vdots\\x_n\end{array}\right]$$

等价于方程：

$$\begin{array}{cccccccccc}y_1&=&u_{11}x_1&+&u_{12}x_2&+&u_{13}x_3&+&\cdots&+&u_{1n}x_n\\y_2&=&&&u_{22}x_2&+&u_{23}x_3&+&&+&u_{2n}x_n\\&\vdots&&&&&\ddots&&&&\\y_n&=&&&&&&&&u_{nn}x_n\end{array}$$

现在，可以从第 $n$ 方程中立即得到 $x_n$，因为 $x_n = \frac{y_n}{u_{nn}}$。

一旦 $x_n$ 可用，将其代入第 $n-1$ 个方程 $y_{n-1} = u_{n-1,n-1}x_{n-1} + u_{n-1,n}x_n$ 并解出 $x_{n-1}$ 得到：

$$x_{n-1} = \frac{y_{n-1} - u_{n-1,n}x_n}{u_{n-1,n-1}}$$

现在，我们有 $x_n$ 和 $x_{n-1}$。将它们代入第 $n-2$ 个方程 $y_{n-2} = u_{n-2,n-2}x_{n-2} + u_{n-2,n-1}x_{n-1} + u_{n-2,n}x_n$ 并解出 $x_{n-2}$ 得到：

$$x_{n-2} = \frac{y_{n-2} - (u_{n-2,n-1}x_{n-1} + u_{n-2,n}x_n)}{u_{n-2,n-2}}$$

从 $x_n, x_{n-1}$ 和 $x_{n-2}$，我们可以从第 $n-3$ 方程中解出 $x_{n-3}$。依此类推，当 $x_n, x_{n-1}, \ldots, x_{i+1}$ 已知时，我们可以利用以下关系从第 $i$ 方程中解出 $x_i$：

$$x_{i} = \frac{1}{u_{ii}}\left[y_{i}-\sum_{k=i+1}^{n}u_{i,k}x_{k}\right]$$

重复这个过程直到计算出 $x_1$。然后，所有的 $x$ 都已经得到，线性方程组就求解出了。

以下算法总结了这个过程：

Input： 矩阵 $\mathbf{Y}_{n \times h}$ 和上三角矩阵 $\mathbf{U}_{n \times n}$
Output： 矩阵 $\mathbf{X}_{n \times h}$

算法：

    for $j$ = $1$ to $h$ do // 共有 $h$ 列
        begin // 对每一列执行以下操作
            $x_{n,j} = y_{n,j} / u_{n,n}$; // 计算当前列的 $x_n$
            for $i$ = $n-1$ downto $1$ do // 处理该列的元素
                begin
                    $sum = 0$; // 求解当前列的 $x_i$
                    for $k$ = $i+1$ to $n$ do
                        $sum = sum + u_{i,k} \times x_{k,j}$;
                    $x_{i,j} = (y_{i,j} - sum) / u_{i,i}$;
                end
        end

这次我们首先计算 $x_n$ ，然后代入前面的方程计算$x_{n-1}$，最后倒推计算得到 $x_1$，因此，这个过程被称为后向代换。

全局曲线插值

要用一组点拟合出B样条曲线，最简单的方法是使用全局插值方法。

假设我们有$n+1$个数据点$\mathbf{D}_0$, $\mathbf{D}_1$, …, $\mathbf{D}_n$，并希望拟合出一条次数为$p$的B样条曲线，其中$p \leq n$是指定值。我们可以选择一组参数值$t_0, t_1, ..., t_n$，其中参数$t_k$对应于数据点$\mathbf{D}_k$，根据这些参数，计算出$m+1$个节点的节点向量，其中$m = n + p + 1$。因此，我们已知节点向量和次数$p$，现在缺少的是一组$n+1$个控制点。全局插值法就是用来找到这些控制点的一种方法。

全局曲线插值：给定一组$n+1$个数据点，$D_0, D_2, ..., D_n$ 和一个次数$p$，找到一条由$n+1$个控制点定义的次数为$p$的B样条曲线，该曲线按给定顺序经过所有数据点。

解决方案

假设次数为$p$的插值B样条曲线定义如下：

$$C(u) = \sum_{i=0}^{n} N_{i,p}(u)P_i$$

这个B样条曲线有$n+1$个未知控制点。由于参数$t_k$对应于数据点$\mathbf{D}_k$，将$t_k$代入上述方程可以得到：

$$\mathbf{D}_k = C(t_k) = \sum_{i=0}^{n} N_{i,p}(t_k)P_i \quad 0 \leq k \leq n$$

上述方程中有$n+1$个B样条基函数($N_{0,p}(u)$, $N_{1,p}(u)$, $N_{2,p}(u)$, …, 和 $N_{n,p}(u)$)和$n+1$个参数($t_0, t_1, t_2, .., 和 t_n$，将这些$t_k$的值代入$N_{i,p}(u)$中即可得到$(n +1)^2$个值。这些值可以写成$(n +1) \times (n +1)$矩阵$\mathbf{N}$的形式，其中第$k$行是在$t_k$处计算得到的基函数$N_{0,p}(u)$, $N_{1,p}(u)$, $N_{2,p}(u)$, …, 和 $N_{i,p}(u)$的值，如下所示：

$$N = \begin{bmatrix} N_{0,p}(t_0) & N_{1,p}(t_0) & \cdots & N_{n,p}(t_0) \\ N_{0,p}(t_1) & N_{1,p}(t_1) & \cdots & N_{n,p}(t_1) \\ \vdots & \vdots & \ddots & \vdots \\ N_{0,p}(t_n) & N_{1,p}(t_n) & \cdots & N_{n,p}(t_n) \end{bmatrix}$$

我们将向量$\mathbf{D}_k$和$\mathbf{P}_i$也写成两个矩阵$\mathbf{D}$和$\mathbf{P}$，如下所示：

$$\left.\mathbf{D}=\left[\begin{array}{ccccc}{d_{01}}&{d_{02}}&{d_{03}}&{\cdots}&{d_{0s}}\\{d_{11}}&{d_{12}}&{d_{13}}&{\cdots}&{d_{1s}}\\{\vdots}&&&{\ddots}&{\vdots}\\{d_{n1}}&{d_{n2}}&{d_{n3}}&{\cdots}&{d_{ns}}\end{array}\right.\right]\quad\mathbf{P}=\left[\begin{array}{ccccc}{p_{01}}&{p_{02}}&{p_{03}}&{\cdots}&{p_{0s}}\\{p_{11}}&{p_{12}}&{p_{13}}&{\cdots}&{p_{1s}}\\{\vdots}&&&{\ddots}&{\vdots}\\{p_{n1}}&{p_{n2}}&{p_{n3}}&{\cdots}&{p_{ns}}\end{array}\right]$$

$\mathbf{D}_k$表示一个$s$维空间中的向量(即，$\mathbf{D}_k=[d_{k1}, \ldots, d_{ks}]$)，它是矩阵$\mathbf{D}$的第$k$行。类似地，$\mathbf{P}_i$也是一个在$s$维空间中的向量(即，$\mathbf{P}_i = [p_{i1}, \ldots, p_{is}]$)。在三维空间中，我们有$s=3$，在平面上我们有$s=2$。$\mathbf{D}$和$\mathbf{P}$都是$(n +1) \times s$矩阵，因此$\mathbf{D}_k$和$t_i$的关系可以写成以下更简单的形式：

$$\mathbf{D} = \mathbf{N} \cdot \mathbf{P}$$

矩阵$\mathbf{D}$为输入的数据点，矩阵$\mathbf{N}$是给定参数计算出的B样条基函数的值，$\mathbf{D}$和$\mathbf{N}$都是已知的，唯一的未知数是矩阵$\mathbf{P}$，因此上面的方程是一个关于$\mathbf{P}$的线性方程组，求解出$\mathbf{P}$即可得到控制点，就能计算出插值的B样条曲线。

算法实现

$\mathbf{D}$的第$i$列和$\mathbf{P}$的第$i$列记为$\mathbf{D}_i$和$\mathbf{P}_i$。根据上述线性方程组，我们有：

$$\mathbf{d}^i = \mathbf{N} \cdot \mathbf{p}^i$$

根据$\mathbf{N}$和$\mathbf{d}^i$求解出$\mathbf{p}^i$，其中$0\leq i \leq h$，对每个$i$重复此操作，计算出$\mathbf{P}$的每一列。这样，所有控制点都已计算出来。

但是，这是非常低效的。不过，许多数值计算库为我们提供了现成的线性系统求解器，能够有效地求解方程组$\mathbf{D} = \mathbf{N}\mathbf{P}$。

下面是相关算法的步骤：

Input： $n+1$个数据点$\mathbf{D}_0$, … $\mathbf{D}_n$和次数$p$
Output： 次数为$p$的B样条曲线，按照给定的顺序经过所有数据点

算法：

​ 选择一种方式计算$n+1$个参数$t_0, ..., t_n$，得到节点向量$U$
​ for $i$ = $0$ to $n$ do
​ for $j$ = $0$ to $n$ do
​ 计算矩阵$\mathbf{N}$中第$i$行和第$j$列的值$N_{j,p}(t_i)$；

​ for $i$ = $0$ to $n$ do
​ 将数据点$\mathbf{D}_i$存储在矩阵$\mathbf{D}$的第$i$行；

​ 使用线性方程求解器求解$\mathbf{D} = \mathbf{N}\mathbf{P}$得到$\mathbf{P}$；
​ $\mathbf{P}$的第$i$行是控制点$\mathbf{P}_i$；
​ 控制点$\mathbf{P}_0$, …, $\mathbf{P}_n$，节点向量$U$和次数$p$确定了一条插值B样条曲线。

下面左图是一个例子。有9个给定数据点（黑色表示），蓝色点为计算出的控制点。插值曲线上的蓝点小圆点是使用弦长法计算出的节点。可以看出，这些“节点”非常接近数据点，控制多边形也紧贴着数据点连成的多边形。但并不是所有的时候都会如此，下图右图中，数据点连成的多边形和控制多边则形非常不同。

如果通过对连续$p$个参数取平均来计算节点，那么矩阵$\mathbf{N}$中每个元素都是正值，且带状的半带宽小于$p$（即，如果$|i - k| \geq p$，则$N_{i,p}(t_k) = 0$），de Boor在1978年证明了这一点。这意味着通过这样的参数化算法得到的线性方程组可以使用高斯消元法求解。

参数和节点的影响

一般来说，所选参数和节点对结果曲线的影响是无法预测的，**但是，如果弦长分布大致相同，那么四种参数化方法生成的曲线应该相似。**此外，通用参数化方法应该与均匀参数化方法表现相似，因为当节点均匀分布时，B样条基函数的最大值分布是均匀的。但是，当弦长分布剧烈变化时，四种参数选择方式产生的曲线则会大不相同。

下图中的四条曲线是使用四种参数化方法获得的，插值B样条曲线的次数为3。均匀参数化产生了一个尖点，弦长参数化生成的曲线在数据点之间会有较大的波动，向心参数化类似于弦长参数化但表现更好，通用参数化方法紧贴着数据点连成的多边形（比均匀参数化更好），但产生了一个小的环。

(均匀法)

(弦长法)

(向心法)

(通用方法)

参数和节点之间的关系如何？下图显示了所有四种方法的参数和节点分布。通用参数化方法获得的参数和节点比弦长参数化和向心参数化的分布更均匀。弦长参数化获得的较长曲线段在向心参数化中会变短，曲线也不会在数据点之间剧烈摆动。

次数的影响

次数对插值B样条曲线形状的影响也很难预测。可以从下图中观察到，均匀参数化和通用参数化法通常很好地跟踪长弦，但是，它们在短弦上会出现问题。由于参数间距几乎相等，因此对于较短的弦，插值曲线会拉伸得更长一些。因此，我们会看到峰值和自交环。当次数增加时，这种情况会变得更加明显，因为高阶曲线提供了更多的摆动自由度。

	均匀法	弦长法	向心法	通用方法
2阶
3阶
4阶
5阶

对于弦长参数化，上图显示它对长弦的效果不是很好，特别是后面或前面有一些较短的弦，可能会发生大的凸起。次数对上图中显示的插值曲线形状没有显著影响。由于向心参数化是弦长参数化的扩展，因此两者具有相同的特点。然而，向心参数化会降低两个相邻参数之间距离之差的影响，它也具有均匀参数化和通用参数化方法相似的特点。例如，生成的插值曲线会紧贴着较长的弦，当次数增加时，短弦处可能会出现自交环。

为什么这种方法是全局的？

即使使用满足局部修改特性的B样条曲线，这种插值方法也是全局的，因为改变单个数据点的位置会完全改变插值曲线的形状。下图中，黄色点是数据点，其中一个被移动到新位置，以浅蓝色标记，并用红色箭头指示。这九个数据点使用4次B样条曲线采用向心参数化插值得到。

可以看到，移动数据点后的插值曲线（蓝色）和原始曲线有八个数据点是相同的，但八个曲线段都不完全相同，因此，改变单个数据点的位置会全局地改变插值曲线的形状！

全局曲线逼近

在插值中，插值曲线会按照顺序通过所有给定的数据点，曲线可能会在数据点之间疯狂摆动，而不是紧贴着数据点连成的多边形。曲线逼近则是放宽了要求，不严格要求曲线必须经过所有点。在全局逼近中，除了第一个和最后一个数据点之外，曲线不必经过每一个点。

为了衡量一条曲线“逼近”给定数据点连成的多边形的程度，提出了误差距离的概念。

误差距离是数据点与其在曲线上的对应点之间的距离。因此，如果这些误差距离的总和最小，曲线就会紧贴数据点连成的多边形。插值曲线就是一种特殊结果，因为每个数据点的误差距离为零，以这种方式获得的曲线称为逼近曲线。

假设我们有$n+1$个数据点$\mathbf{D}_0, \mathbf{D}_1, ..., \mathbf{D}_n$，并希望找到一条接近多边形形状的B样条曲线，而不用经过每个数据点。为此，我们需要两个额外的条件：控制点的数量（假设为$h+1$）和一个次数（设为$p$），其中$n > h \geq p \geq 1$。逼近比插值更灵活，因为我们不仅可以选择次数，还可以选择控制点的数量。

全局曲线逼近：给定一组$n+1$个数据点，$\mathbf{D}_0, \mathbf{D}_1, ..., \mathbf{D}_n$，一个次数 p，和数量 h，其中 $n > h \geq p \geq 1$，找到一条由 h+1 个控制点定义的次数为 p 的B样条曲线，并满足以下条件：

1. 这条曲线经过第一个点和最后一个点（即$\mathbf{D}_0$ 和 $\mathbf{D}_n$）；
2. 这条曲线上各个点的误差距离之和最小。

有了 $h$ 和 $p$，我们可以确定一组参数和一个节点向量。设参数为 $t_0, t_1, ..., t_n$。参数的数量必须等于数据点的数量。现在，假设次数为 p 的逼近B样条曲线为：

$$C(u) = \sum_{i=0}^{h} N_{i,p}(u)P_i$$

其中 $\mathbf{P}_0, \mathbf{P}_1, ..., \mathbf{P}_h$ 是 $h+1$ 个未知控制点。由于我们希望曲线通过第一个和最后一个数据点，我们有 $\mathbf{D}_0 = \mathbf{C}(0) = \mathbf{P}_0$ 和 $\mathbf{D}_n = \mathbf{C}(1) = \mathbf{P}_h$。因此，只有 $h-1$ 个未知控制点 $\mathbf{P}_1, \mathbf{P}_2, ..., \mathbf{P}_{h-1}$，这样曲线方程变为以下形式：

$$\mathbf{C}(u)=N_{0,p}(u)\mathbf{D}_{0}+\left(\sum_{i=1}^{h-1}N_{i,p}(u)\mathbf{P}_{i}\right)+N_{h,p}(u)\mathbf{D}_{n}$$

最小二乘法

如何测量误差距离？

参数 $t_k$ 对应于数据点 $\mathbf{D}_k$, $\mathbf{D}_k$ 和曲线上$t_k$对应的点之间的距离是 $|\mathbf{D}_k - \mathbf{C}(t_k) |$。由于这个距离的计算需要开根，我们通常选择使用平方距离$|\mathbf{D}_k - \mathbf{C}(t_k) |^2$ 来避免开根计算带来的额外开销。

因此，所有平方误差距离之和是：

$$f(P_1, \ldots, P_{h-1}) = \sum_{k=1}^{n-1} | \mathbf{D}_k - C(t_k) |^2$$

我们的目标是找到一些控制点 $\mathbf{P}_1, ..., \mathbf{P}_{h-1}$，使得函数 $f()$ 最小化。

解决方案

我们将$\mathbf{D}_k - \mathbf{C}(t_k)$改写成另一种形式：

$$\begin{aligned}\mathbf{D}_{k}-\mathbf{C}(t_{k})&=\quad\mathbf{D}_k-\left[N_{0,p}(t_k)\mathbf{D}_0+\left(\sum_{i=1}^{h-1}N_{i,p}(t_k)\mathbf{P}_i\right)+N_{h,p}(t_k)\mathbf{D}_n\right]\\&=\quad(\mathbf{D}_k-N_{0,p}(t_k)\mathbf{D}_0-N_{h,p}(t_k)\mathbf{D}_n)-\left(\sum_{i=1}^{h-1}N_{i,p}(t_k)\mathbf{P}_i\right)\end{aligned}$$

在上面的公式中，$\mathbf{D}_0$, $\mathbf{D}_k$ 和 $\mathbf{D}_n$ 是给定的，而 $N_{0,p}(t_k)$ 和 $N_{h,p}(t_k)$ 是 $N_{0,p}(u)$ 和 $N_{h,p}(u)$ 在 $t_k$ 处的值，为了方便，我们定义一个新的向量 $\mathbf{Q}_k$ 如下：

$$\mathbf{Q}_k = \mathbf{D}_k - N_{0,p}(t_k)\mathbf{D}_0 - N_{h,p}(t_k)\mathbf{D}_h$$

然后，误差距离平方和函数 $f()$ 可以写成以下形式：

$$f(\mathbf{P}_1,\ldots,\mathbf{P}_{h-1})=\sum_{k=1}^{n-1}\left|\mathbf{Q}_k-\left(\sum_{i=1}^{h-1}N_{i,p}(t_k)\mathbf{P}_i\right)\right|^2$$

回想一下恒等式 $\mathbf{x} · \mathbf{x} = | \mathbf{x} |^2$。这意味着向量 $\mathbf{x}$ 与其自身的内积等于 $\mathbf{x}$ 的长度的平方。因此，误差平方项可以重写为：

$$\begin{aligned} &\left|\mathbf{Q}_{k}-\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\right|^{2} \\ &=\left(\mathbf{Q}_{k}-\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\right)\cdot\left(\mathbf{Q}_{k}-\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\right) \\ &=\mathbf{Q}_{k}\cdot\mathbf{Q}_{k}-2\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\cdot\mathbf{Q}_{k}\right)+\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right) \end{aligned}$$

因此，函数$f$变为：

$$\begin{aligned} f(\mathbf{P}_{1},\ldots,\mathbf{P}_{h-1})&=\sum_{k=1}^{n-1}\left[\mathbf{Q}_{k}\cdot\mathbf{Q}_{k}-2\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\cdot\mathbf{Q}_{k}\right)\right.\\ &\quad\left.+\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\right] \end{aligned}$$

让函数$f$对每个$\mathbf{P}_{g}$求偏导，并找到这些偏导数的零点。在这些共同零点处，函数$f$达到最小值。

在计算函数$f$对每个$\mathbf{P}_{g}$的偏导时，$\mathbf{Q}_k$和$N_{i,p}(t_k)$都是常数，它们对$\mathbf{P}_{g}$的偏导数是零。因此，我们有：

$$\frac{\partial}{\partial\mathbf{P}_{g}}\left(\mathbf{Q}_{k}\cdot\mathbf{Q}_{k}\right)=0$$

考虑求和中的第二项，即 $N_{i,p}(t_k) \mathbf{P}_i \mathbf{Q}_k$ 的总和。每个子项的偏导数计算如下：

$$\frac{\partial}{\partial\mathbf{P}_{g}}\left(N_{i,p}(t_{k})\mathbf{P}_{i}\cdot\mathbf{Q}_{k}\right)=N_{i,p}(t_{k})\frac{\partial\mathbf{P}_{i}}{\partial\mathbf{P}_{g}}\cdot\mathbf{Q}_{k}$$

$\mathbf{P}_i$ 对$\mathbf{P}_{g}$ 的偏导数仅在 $i = g$ 时非零。因此，第二项的偏导数如下：

$$\frac{\partial}{\partial\mathbf{P}_{g}}\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\cdot\mathbf{Q}_{k}\right)=N_{g,p}(t_{k})\mathbf{Q}_{k}$$

第三项的偏导数比较复杂，需要用到求导规则 $(f.g)' = f'.g + f.g'$ ：

$$\begin{gathered} {\frac{\partial}{\partial\mathbf{P}_{g}}}\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right) \\ =\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\frac{\partial\mathbf{P}_{i}}{\partial\mathbf{P}_{g}}\right)+\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\frac{\partial\mathbf{P}_{i}}{\partial\mathbf{P}_{g}}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right) \\ =2\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\frac{\partial\mathbf{P}_{i}}{\partial\mathbf{P}_{g}}\right) \end{gathered}$$

由于$\mathbf{P}_i$ 对 $\mathbf{P}_{g}$ 的偏导数在 $i$ 不等于 $g$ 时为零，求和中第三项对 $\mathbf{P}_{g}$ 的偏导数是：

$$\frac{\partial}{\partial\mathbf{P}_{g}}\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)\cdot\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)=2N_{g,p}(t_{k})\left(\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}\right)$$

综上，函数 $f()$ 对 $\mathbf{P}_{g}$ 的偏导数是：

$$\frac{\partial f}{\partial\mathbf{P}_{g}}=-2N_{g,p}(t_{k})\mathbf{Q}_{k}+2N_{g,p}(t_{k})\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}$$

令其等于零，我们得到：

$$\sum_{k=1}^{n-1}N_{g,p}(t_k)\sum_{i=1}^{h-1}N_{i,p}(t_k)\mathbf{P}_i=\sum_{k=1}^{n-1}N_{g,p}(t_k)\mathbf{Q}_k$$

由于我们有 $h-1$ 个变量，$g$ 的取值从 1 到 $h-1$ ，因此有 $h-1$ 个这样的方程。

我们定义以下三个矩阵：

$$\left.\mathbf{P}=\left[\begin{array}{c}\mathbf{P}_1\\\mathbf{P}_2\\\vdots\\\mathbf{P}_{h-1}\end{array}\right.\right]\quad\mathbf{Q}=\left[\begin{array}{c}\sum_{k=1}^{n-1}N_{1,p}(t_k)\mathbf{Q}_k\\\sum_{k=1}^{n-1}N_{2,p}(t_k)\mathbf{Q}_k\\\vdots\\\sum_{k=1}^{n-1}N_{h-1,p}(t_k)\mathbf{Q}_k\end{array}\right]$$

$$\mathbf{N}=\left[\begin{array}{cccc}N_{1,p}(t_1)&N_{2,p}(t_1)&\cdots&N_{h-1,p}(t_1)\\N_{1,p}(t_2)&N_{2,p}(t_2)&\cdots&N_{h-1,p}(t_2)\\\vdots&&\ddots&\vdots\\N_{1,p}(t_{n-1})&N_{2,p}(t_{n-1})&\cdots&N_{h-1,p}(t_{n-1})\end{array}\right]$$

在这里，$\mathbf{P}$ 的第 $k$ 行是向量 $\mathbf{P}_k$，$\mathbf{Q}$ 的第 $k$ 行是上述第 $k$ 个方程的右侧，$\mathbf{N}$ 的第 $k$ 行是 $N_{1,p}(u), N_{2,p}(u), ..., N_{h-1,p}(u)$ 在$t_k$处的值。因此，如果输入的数据点是 $s$ 维向量，$\mathbf{P}$, $\mathbf{N}$ 和 $\mathbf{Q}$ 分别是 $(h-1)×s, (n-1)×(h-1) $和 $(h-1)×s$ 矩阵。

现在，让我们重写第 $g$ 个线性方程：

$$\sum_{k=1}^{n-1}N_{g,p}(t_{k})\sum_{i=1}^{h-1}N_{i,p}(t_{k})\mathbf{P}_{i}=\sum_{k=1}^{n-1}N_{g,p}(t_{k})\mathbf{Q}_{k}$$

换一种形式，以便可以看出 $\mathbf{P}_i$ 的系数：

$$\sum_{i=1}^{h-1}\left(\sum_{k=1}^{n-1}N_{g,p}(t_k)N_{i,p}(t_k)\right)\mathbf{P}_i=\sum_{k=1}^{n-1}N_{g,p}(t_k)\mathbf{Q}_k$$

最后，$\mathbf{P}_i$ 的系数是：

$$\begin{aligned}\sum_{k=1}^{n-1}N_{g,p}(t_k)N_{i,p}(t_k)\end{aligned}$$

如果观察矩阵 $\mathbf{N}$，就会发现 $N_{g,p}(t_1), N_{g,p}(t_2), ..., N_{g,p}(t_{n-1})$ 是 $\mathbf{N}$ 的第 $g$ 列，而 $N_{i,p}(t_1), N_{i,p}(t_2), ..., N_{i,p}(t_{n-1})$ 是 $\mathbf{N}$ 的第 $i$ 列。$\mathbf{N}$ 的第 $g$ 列是 $\mathbf{N}$ 转置矩阵 $\mathbf{N}^T$ 的第 $g$ 行，而 $\mathbf{P}_i$ 的系数是 $\mathbf{N}$^T 的第 $g$ 行和 $\mathbf{N}$ 的第 $i$ 列的内积。有了这个观察，线性方程组可以重写为：

$$\begin{pmatrix}\mathrm{N}^T\mathrm{N}\end{pmatrix}\mathrm{P}=\mathrm{Q}$$

由于 $\mathbf{N}$ 和 $\mathbf{Q}$ 是已知的，求解关于 $\mathbf{P}$的线性方程组就得到了所需的控制点。

算法

Input： $n+1$个数据点$\mathbf{D}_0$, $\mathbf{D}_1$, …, $\mathbf{D}_n$，次数$p$，和所需的控制点数量$h+1$；
Output： 一条由$h+1$个控制点定义的次数为$p$的B样条曲线，由给定的数据点逼近而成；

算法：

    计算一组参数$t_0, ..., t_n$和节点向量$U$
    让$\mathbf{P}_0=\mathbf{D}_0$且$\mathbf{P}_h=\mathbf{D}_n$
    for $k$ = $1$ to $n-1$ do
        通过下式计算$\mathbf{Q}_k$：
            $\mathbf{Q}_{k}=\mathbf{D}_{k}-N_{0,p}(t_{k})\mathbf{D}_{0}-N_{h,p}(t_{k})\mathbf{D}_{n}$
    for $i$ = $1$ to $h-1$ do
        计算下式并存储到矩阵$\mathbf{Q}$的第$i$行中：
            $\sum_{k=1}^{n-1}N_{i,p}(t_k)\mathbf{Q}_k$
    得到矩阵$\mathbf{Q}$

    for $k$ = $1$ to $n-1$ do
        for $i$ = $1$ to $h-1$ do
            计算$N_{i,p}(t_k)$并存储到矩阵$\mathbf{N}$的第$k$行第$i$列
    得到矩阵$\mathbf{N}$
    计算$\mathbf{M} = \mathbf{N}^T \mathbf{N}$
    根据方程$\mathbf{M} \cdot \mathbf{P} = \mathbf{Q}$求解出$\mathbf{P}$
    矩阵$\mathbf{P}$的第$i$行$\mathbf{P}_i$就是控制点集
    控制点$\mathbf{P}_0$, …, $\mathbf{P}_h$，节点向量$U$和次数$p$已知，可以得到逼近的B样条曲线

控制点数量和次数的影响

显然，数据点数量会影响逼近曲线的形状。那么次数 $p$ 和控制点数量对曲线形状有什么影响？下图显示了10个数据点（n = 9）在不同次数和控制点数量下的逼近曲线。每一行是次数相同但控制点数量不同的逼近曲线，每一列是控制点数量相同但次数不同的逼近曲线，所有曲线都使用向心参数化计算其参数。

	控制点数量=4	控制点数量=5	控制点数量=6	控制点数量=7
次数=2
次数=3
次数=4
次数=5

可以看出，通常次数较低的曲线无法很好地逼近数据点连成的多边形，而较高次数的曲线结果更好（更接近数据多边形）。类似地，控制点越多，逼近曲线的灵活性就越高。因此，在每一行中，随着控制点数量的增加，曲线更接近数据点连成的多边形。

那是否应该使用高阶曲线和更多的控制点？答案是不，因为全局逼近只需要比全局插值更少的控制点。如果控制点的数量等于数据点的数量，全局逼近就变成了全局插值，我们可以直接使用全局插值！至于次数，只要生成的曲线能够接近数据多边形的形状，我们希望次数尽可能的小。

为什么这种方法是全局的？

这种逼近方法是全局的，因为改变一个数据点的位置会导致整个曲线发生变化。下图中的黄点是给定的数据点，用一条由5个控制点定义的3次B样条曲线（$n = 7, p=3$ 和 $h = 4$）来逼近数据点。参数是使用向心参数化计算获得的。假设第四个数据点被移动到蓝色点标记的新位置。蓝色曲线是新生成的插值B样条曲线。

如图所示，除了经过第一个和最后一个数据点外，原始曲线和新曲线的形状截然不同，因此，修改数据点引起的变化是全局的！

计算节点向量

得到一组参数后，我们就可以计算节点向量，参数点相当于控制点。假设我们有$n+1$个参数$t_0, t_1, ..., t_n$ 和 次数 $p$。对于一个次数为$p$的B样条曲线，我们需要$m+1$个节点，$m$满足条件$m = n + p + 1$。如果B样条曲线是clamped的，那么这些节点的分布为：

$u_0 = u_1 = ... = u_p = 0, u_{p+1}, ..., u_{m-p-1}, u_{m-p} = u_{m-p+1} = ... = u_m = 1.$

前 $p+1$ 个和最后 $p+1$ 个节点分别是0和1。其余$n-p$个节点可以均匀分布，也可以是其他非均匀分布。假设其余的$n-p$个内部节点均匀分布。那么有，$u_p = 0, u_{p+1}, ..., u_{m-p-1}, u_{m-p} = 1$将$[0,1]$ 分成 $n-p+1$个子区间。因此，完整的节点为：

$$\begin{aligned}u_{0}&=\quad u_1=\cdots=u_p=0\\u_{j+p}&=\quad\frac j{n-p+1}\quad\mathrm{for~}j=1,2,\ldots,n-p\\u_{m-p}&=\quad u_{m-p+1}=\cdots=u_m=1\end{aligned}$$

例如，如果我们有6个（$n = 5$）参数点和已知B样条次数 $p = 3$，那么我们应该找到 $(n + p + 1) + 1 = (5 + 3 + 1) + 1 = 10$个节点（即，$m = 9$）。由于我们使用的是clamped曲线，前四个和最后四个（即，$p + 1$）节点应该是0和1，节点为0, 0, 0, 0, $u_4$, $u_5$, 1, 1, 1, 1, 两个内部节点$u_4,u_5$ 将$[0,1]$分成三个子区间，每个区间的长度为1/3。因此，完整的节点向量为：

$$\{0, 0, 0, 0, \frac{1}{3}, \frac{2}{3}, 1, 1, 1, 1\}.$$

均匀分布的节点向量计算起来非常简单，不过，不推荐使用这种方法，因为如果采用弦长参数化得到参数，并使用均匀分布法计算节点，那么最后得到的线性方程组是奇异的。

另一种计算节点向量的方法是由de Boor提出，即平均参数法，计算公式如下：

$$\begin{aligned}u_{0}&=\quad u_1=\cdots=u_p=0\\u_{j+p}&=\quad\frac1p\sum_{i=j}^{j+p-1}t_i\quad\text{for }j=1,2,\ldots,n-p\\u_{m-p}&=\quad u_{m-p+1}=\cdots=u_m=1\end{aligned}$$

因此，第一个内部节点是 $p$ 个参数 $t_1, t_2, ..., t_p$ 的平均值；第二个内部节点是后面 $p$ 个参数$t_2, t_3, ..., t_{p+1}$的平均值。假设我们有6个（$n = 5$）参数，如下所示：

$$\begin{array}{|c|c|c|c|c|c|c|} \hline t_0 & t_1 & t_2 & t_3 & t_4 & t_5 \\ \hline 0 & \frac{1}{4} & \frac{1}{3} & \frac{2}{3} & \frac{3}{4} & 1 \\ \hline \end{array}$$

假设我们要计算一个次数为$3$的B样条曲线的节点向量。需要10个节点（$m = 9$）。这10个节点中，前四个和后四个节点分别是0和1。因此，第一个内部节点是参数$\frac{1}{4}, \frac{1}{3}, \frac{2}{3}$的平均值，第二个内部节点是接下来的三个参数$\frac{1}{3}, \frac{2}{3}, \frac{3}{4}$的平均值。得到节点向量是：

$$\begin{array}{|c|c|c|c|c|} \hline u_0 = u_1 = ... = u_2 = u_3 & u_4 & u_5 & u_6 = u_7 = ... = u_8 = u_9 \\ \hline 0 & \frac{1}{4} + \frac{1}{3} + \frac{2}{3} / 3 = \frac{5}{12} & \frac{1}{3} + \frac{2}{3} + \frac{3}{4} / 3 = \frac{7}{12} & 1 \\ \hline \end{array}$$

下图说明了参数的位置和得到的节点位置。注意$0(4)$和$1(4)$表示$0$和$1$是四重节点（即重数=4）。可以看到，第一个非零节点区间$[0, \frac{5}{12})$包含两个参数，第二个非零节点区间 $[\frac{5}{12}, \frac{7}{12})$ 不包含参数，第三个非零节点区间 $[\frac{7}{12}, 1)$ 包含两个参数。

全局曲面插值

假设我们有$m+1$行$n+1$列个数据点$\mathbf{D}_{ij}$（$0 \leq i \leq m$ ， $0 \leq j \leq n$），希望找到一个包含所有这些点的次数为$(p, q)$的B样条曲面。与曲线类似，我们已知数据点和次数 $p$ 和 $q$ ，要定义一个B样条插值曲面，还需要两个方向的节点向量$\mathbf{U}$和$\mathbf{V}$，以及一组控制点。控制点的数量要和数据点的数量相等（即有$(m+1) \times (n+1)$个控制点）。

之前在 曲面参数化和节点向量计算 中讨论过，我们可以计算出 $u$ 和 $v$ 两个方向的参数$s_c$（$0 \leq c \leq m$）和$t_d$（$0 \leq d \leq n$），同时还可以计算出 $u$ 和 $v$ 方向的节点向量$\mathbf{U}$和$\mathbf{V}$。因此，剩下的就是要找到所需的控制点！

全局曲面插值：给定一个由$(m+1) \times (n+1)$个数据点组成的网格 $\mathbf{D}_{ij}$（$0 \leq i \leq m$ 和 $0 \leq j \leq n$）和次数 $(p, q)$，找到一个由 $(m+1) \times (n+1)$ 个控制点定义的次数为 $(p, q)$ 的B样条曲面，该曲面包含所有的数据点。

解决方案

设B样条曲面定义如下：

$$S(u,v) = \sum_{i=0}^{m} \sum_{j=0}^{n} N_{i,p}(u)N_{j,q}(v)P_{ij}$$

由于它经过所有数据点，并且参数 $s_c$ 和 $t_d$ 对应于数据点$\mathbf{D}_{cd}$，将$u = s_c$和$v = t_d$代入曲面方程得到：

$$D_{cd} = S(s_c, t_d) = \sum_{i=0}^{m} \sum_{j=0}^{n} N_{i,p}(s_c)N_{j,q}(t_d)P_{ij}$$

由于$N_{i,p}(s_c)$ 与索引 $j$ 无关，因此可以将它从对 $j$ 的求和中提出来：

$$\begin{aligned}\mathbf{D}_{cd}&=\mathbf{S}(s_{c},t_{d})=\sum_{i=0}^{m}\sum_{j=0}^{n}N_{i,p}(s_{c})N_{j,q}(t_{d})\mathbf{P}_{ij}\\&=\sum_{i=0}^{m}N_{i,p}(s_{c})\left(\sum_{j=0}^{n}N_{j,q}(t_{d})\mathbf{P}_{ij}\right)\end{aligned}$$

我们发现索引 $i$ 只出现在$\mathbf{P}_{ij}$中。因此，我们可以将括号中的表达式定义成一个新项，如下所示：

$$\mathbf{Q}_{id}=\sum_{j=0}^nN_{j,q}(t_d)\mathbf{P}_{ij}$$

准确地说，如果 $i$ 为某一固定值，$\mathbf{Q}_{id}$ 就是 $q$ 次B样条曲线（等式右侧为 $q$ 次B样条曲线的定义）在 $t_d$ 处的值，该曲线由 $\mathbf{P}$的第 $i$ 行上的$n+1$个未知控制点（即，$\mathbf{P}_{i0}$, $\mathbf{P}_{i1}$, …, $\mathbf{P}_{in}$）定义 。将$\mathbf{Q}_{id}$ 代入上面的$\mathbf{D}_{cd}$方程得到：

$$D_{cd} = \sum_{i=0}^{m} N_{i,p}(s_c) Q_{id}$$

因此，数据点$\mathbf{D}_{cd}$是 $p$ 次B样条曲线（等式右侧为$p$次B样条曲线的定义）在$s_c$处求值的结果，该曲线由$\mathbf{Q}$的第$d$列的$m+1$个未知控制点（即，$\mathbf{Q}_{0d}$, $\mathbf{Q}_{1d}$, …, $\mathbf{Q}_{md}$）定义。

对每个$c$（$0 \leq c \leq m$）重复上述步骤，数据点的第$d$列（即，$\mathbf{D}_{0d}$, $\mathbf{D}_{1d}$, …, $\mathbf{D}_{md}$）是根据$\mathbf{Q}$的第$d$列（即，$\mathbf{Q}_{0d}$, $\mathbf{Q}_{1d}$, …, $\mathbf{Q}_{md}$）和参数$s_0, s_1, ..., s_m$ 计算获得。由于数据点$\mathbf{D}_{0d}$, $\mathbf{D}_{1d}$, …, $\mathbf{D}_{md}$，次数 $p$ 和参数$s_0, s_1, ..., s_m$是已知的，问题总结如下：

给定次数$p$，参数$s_0, s_1, ..., s_m$和 第$d$列数据点$\mathbf{D}_{0d}$, $\mathbf{D}_{1d}$, …, $\mathbf{D}_{md}$，求第 $d$ 列的控制点 $\mathbf{Q}_{0d}$, $\mathbf{Q}_{1d}$, …, $\mathbf{Q}_{md}$

因此，这其实就是一个曲线插值问题！全局曲线插值方法可以应用到数据点的每一列，从而求得每一列控制点$\mathbf{Q}_{cd}$。我们有$n+1$列数据点，我们将得带$n+1$列控制点$\mathbf{Q}$。

现在，我们使用相同的思路，应用到$\mathbf{Q}_{id}$的等式中，该等式如下：

$$\mathbf{Q}_{id}=\sum_{j=0}^{n}N_{j,q}(t_{d})\mathbf{P}_{ij}$$

在这个等式中，$\mathbf{Q}$的第$i$行数据点（即$\mathbf{Q}_{i0}$, $\mathbf{Q}_{i1}$, …, $\mathbf{Q}_{in}$）是由$n+1$个未知控制点$\mathbf{P}_{i0}$, $\mathbf{P}_{i1}$, …, $\mathbf{P}_{in}$定义的$q$次B样条曲线 在参数 $t_0, t_1, ..., t_n$处求值得到的点。因此，已知次数 $q$ 和参数 $t_0, t_1, ..., t_n$ 对$\mathbf{Q}$的第$i$行进行曲线插值可以得到 第$i$行 的控制点。

一旦找到所有行的控制点，这些控制点连同两个节点向量和次数$p$和$q$就定义了一个插值B样条曲面。因此，使用B样条的曲面插值等价于 ($m+1$）+（$n+1$）次曲线插值！

算法

Input： $(m+1) \times (n+1)$ 个数据点 $\mathbf{D}_{ij}$ 和次数 ($p, q$)；
Output： 一个经过所有数据点的次数为 ($p, q$) 的B样条曲面；

    计算$u$方向的参数$s_0, s_1, ..., s_m$和节点向量$U$；
    计算$v$方向的参数$t_0, t_1, ..., t_n$和节点向量$V$；
    for $d$ = $0$ to $n$ do // 第$d$列
        begin // 计算第$d$列控制点$\mathbf{Q}_d$
            已知次数$p$、参数$s_0, s_1, ..., s_m$和节点向量$\mathbf{U}$，对$\mathbf{D}$的第$d$列数据点应用曲线插值
            （即$\mathbf{D}_{0d}, \mathbf{D}_{1d}, ..., \mathbf{D}_{md}$）
            结果为第$d$列控制点$\mathbf{Q}_{0d}, \mathbf{Q}_{1d}, ..., \mathbf{Q}_{md}$
        end

    for $c$ = $0$ to $m$ do // 第$c$行
        begin // 计算第$c$行控制点$\mathbf{P}_c$
            已知次数$q$、参数$t_0, t_1, ..., t_n$和节点向量$\mathbf{V}$，对$\mathbf{Q}$的第$c$行应用曲线插值
            （即$\mathbf{Q}_{c0}, \mathbf{Q}_{c1}, ..., \mathbf{Q}_{cn}$）
            结果为第$c$行控制点$\mathbf{P}_{c0}, \mathbf{P}_{c1}, ..., \mathbf{P}_{cn}$
        end

得到 $(m+1) \times (n+1)$ 个控制点 $\mathbf{P}_{ij}$，结合次数 $(p, q)$ 和节点向量 $\mathbf{U}$ 和 $\mathbf{V}$ ，可以定义一个插值B样条曲面，经过给定的数据点。

第一个 for 循环中使用的矩阵 $\mathbf{N}$在循环中计算时不会发生变换，如果我们按照上面的算法步骤计算，最终会求解方程组 $\mathbf{D}=\mathbf{N}\mathbf{Q}$ $n+1$次。因此， 为了加速计算，在第一个 for 循环开始之前，应该先计算 $\mathbf{N}$ 的LU分解，这样每次迭代中的插值都是简单的前向代换和后向代换。类似地，第二个 for 循环中的矩阵 $\mathbf{N}$，它的LU分解 也应该在第二个 for 开始之前计算。如果不这样做，我们将执行 $(m+1 + n+1 = m + n + 2)$ LU分解，而只需要2次就足够了。

为什么这种方法是全局的？

因为插值曲面是通过 $m + n + 2$ 次全局曲线插值得到的，由于曲线的插值是全局的，因此这种曲面插值技术也是全局的。
下图显示了移动一个数据点对插值曲面的影响。这个B样条曲面是通过6行5列$(m=5，n=4)$ 数据点插值生成，曲面次数为(3,3)。第一行的图像显示了节点的等参数曲线，黄圈标记了要移动的数据点。
第二行显示了对应的曲面。显然，改变一个数据点的位置后，蓝色节点曲线的形状发生了剧烈变化，并且向邻近的数据点靠近。影响也会向右侧的节点传播。
从曲面的图中可以看出。数据点移动到右侧后，在左端会产生了一个大的凸起。前侧的边界曲线的形状也稍微发生了改变。

显然，改变数据点的位置会影响到整个曲面，这是一种全局的方法。

	移动前	移动后
节点曲线
曲面

全局曲面逼近

我们希望找到一个B样条曲面来逼近$(m+1) \times (n+1)$个数据点，逼近曲面并不需要经过所有给定的数据点，与曲线逼近一样，我们可以控制逼近曲面的次数和控制点的数量。因此，除了数据点已知以外，还需要确定$u$方向和$v$方向的次数 $p$ 和 $q$ ，以及控制点数量，(假设为$e+1$行$f+1$列)。对于B样条曲面，输入值必须满足 $m > e \geq p \geq 1$ 和 $n > f \geq q \geq 1$ 。

全局曲面逼近：给定$(m+1) \times (n+1)$个数据点组成的网格 $\mathbf{D}_{ij} (0 \leq i \leq m ,0 \leq j \leq n)$，找到一个由$(e+1) \times (f+1)$个控制点 $\mathbf{P}_{ij} (0 \leq i \leq e , 0 \leq j \leq f)$ 定义的次数为 $(p,q)$ 的B样条曲面，该曲面能够逼近给定的数据点网格。

解决方案

设由$(e+1) \times (f+1)$个控制点 $\mathbf{P}_{ij}$ 定义的次数为 $(p,q)$ 的B样条曲面定义如下：

$$S(u,v) = \sum_{i=0}^{e} \sum_{j=0}^{f} N_{i,p}(u)N_{j,q}(v)P_{ij}$$

由于有$m+1$行$n+1$列数据点，因此我们需要$m+1$个$u$方向的参数，$s_0, s_1, ..., s_m$，和$n+1$个$v$方向的参数，$t_0, t_1, ..., t_n$。这些参数的计算在 曲面的参数和节点向量 的讨论中给出了计算方法。

有了这些参数，曲面上与数据点 $\mathbf{D}_{cd}$ 对应的点的计算方式如下：

$$S(s_c, t_d) = \sum_{i=0}^{e} \sum_{j=0}^{f} N_{i,p}(s_c)N_{j,q}(t_d)\mathbf{P}_{ij}$$

$\mathbf{D}_{cd}$ 与曲面上对应的点的误差距离平方为：

$$\left|\mathrm{D}_{cd}-\mathrm{S}(s_{c},t_{d})\right|^{2}$$

因此，所有误差距离平方的总和为：

$$f(P_{00}, P_{01}, ..., P_{ef}) = \sum_{c=0}^{m} \sum_{d=0}^{n} |\mathbf{D}_{cd} - S(s_c, t_d)|^2$$

这是关于$(e+1) \times (f+1)$个未知控制点 $\mathbf{P}_{ij}$ 的函数。为了求$f()$的最小值，我们计算$f()$的偏导数并令其等于0：

$$\frac{\partial f}{\partial\mathbf{P}_{ij}}=0$$

然后，我们可以得到$(e+1) \times (f+1)$个方程，它们的共同零点就是所求的控制点。然而这些方程不是线性的，解决非线性方程组非常耗时，与其追求最优解，不如寻找一个合理的、但不是函数 f() 最小值的解。

为了找到一个非最优解，我们采用与全局曲面插值中类似的方法。对每个数据点列应用曲线逼近以计算一些中间数据点。

这样，每列 $m+1$ 个数据点就能计算出$e+1$个“中间”数据点。由于有 $n+1$ 列，这些“中间”数据点构成了一个$(e+1) \times (n+1)$的网格。然后，对这些中间数据点的每一行应用曲线逼近 计算出所需的控制点。由于每一行有 $n+1$ 个“中间”数据点，并且有$e+1$ 行，因此每行都会产生 $f+1$ 个控制点，最终我们将得到 $(e+1) \times (f+1)$ 个控制点。

算法

Input： $(m+1) \times (n+1)$个数据点，记作$\mathbf{D}_{ij}$，次数$(p,q)$，以及控制点的数量$(e+1) \times (f+1)$；
Output： 逼近的$(p, q)$次B样条曲面；

    计算$u$方向的参数$s_0, s_1, ..., s_m$和节点向量$\mathbf{U}$；
    计算$v$方向的参数$t_0, t_1, ..., t_n$和节点向量$\mathbf{V}$；
    for $d$ = $0$ to $n$ do // $\mathbf{D}$的第$d$列
        begin // 计算"中间数据点"$\mathbf{Q}$
            已知次数$p$、参数$s_0, s_1, ..., s_m$和节点向量$\mathbf{U}$，对给定数据点的第$d$列应用曲线逼近（即$\mathbf{D}_{0d}$, $\mathbf{D}_{1d}$, …, $\mathbf{D}_{md}$）。得到第$d$列"中间数据点"$\mathbf{Q}_{0d}, \mathbf{Q}_{1d}, ..., \mathbf{Q}_{ed}$
            // $\mathbf{Q}$是一个$(e+1) \times (n+1)$矩阵
        end

    for $c$ = $0$ to $e$ do // $\mathbf{Q}$的第$c$行
        begin // 计算所需的控制点$\mathbf{P}$
            已知次数$q$、参数$t_0, t_1, ..., t_n$和节点向量$\mathbf{V}$，对$\mathbf{Q}$的第$c$行应用曲线逼近（即$\mathbf{Q}_{c0}$, $\mathbf{Q}_{c1}$, …, $\mathbf{Q}_{cn}$）
            得到第$c$行控制点$\mathbf{P}_{c0},\mathbf{P}_{c1}, ..., \mathbf{P}_{cf}$
            // $\mathbf{P}$是一个$(e+1) \times (f+1)$矩阵
        end

计算得到$(e+1) \times (f+1)$个控制点，记作$\mathbf{P}_{ij}$，次数$(p, q)$和节点向量$\mathbf{U}$和$\mathbf{V}$已知，因此可以得到逼近的B样条曲面。

在这个算法中，先对$\mathbf{D}$的列进行了$n+1$次曲线逼近，然后对“中间”数据点的行进行了$e+1$次曲线逼近。因此，一共执行了$n + e + 2$次曲线逼近。另外，我们也可以对数据点的每一行进行 $m+1$ 次曲线逼近，创建$(m+1) \times (f+1)$个“中间”数据点。然后，对这些“中间”数据点的每一列进行 $f+1$ 次曲线逼近，最终得到 $(e+1) \times (f+1)$个控制点。通过这种方式，一共执行了$m + f + 2$ 次曲线逼近。

由于这种算法并没有使得误差度量函数 $f()$ 达到最小值，因此它不是最优解，但对于许多应用场景来说，已经满足需要。

注意，上面的算法中，第一个for循环中使用的矩阵 $\mathbf{N}^T·\mathbf{N}$ 不会改变，如果我们按照上述步骤直接计算，最终将会求解方程$\mathbf{Q} =(\mathbf{N}^T·\mathbf{N})\mathbf{P}$ $n+1$ 次。这明显有重复的计算。为了加速计算，在第一个for循环开始之前，应该对$\mathbf{N}^T·\mathbf{N}$ 进行 $\mathbf{L}\mathbf{U}$ 分解，在每次循环中，只会计算一个前向代换和后向代换。类似地，第二个for循环会使用一个新的矩阵$\mathbf{N}^T·\mathbf{N}$。也应该在第二个for开始之前对它进行$\mathbf{L}\mathbf{U}$ 分解。如果不这样，将会执行$(e+1 + (n+1)) = e + n + 2$ 次 $\mathbf{L}\mathbf{U}$ 分解，而实际上只需要2次就足够了。

插值和逼近比较

下面是一个6行5列数据点定义的网格。我们用它来说明插值方法和逼近方法之间的区别。第一列的图像是移动黄色数据点之前的结果，第二列是移动了黄色数据点之后的结果。

第一行和第二行显示了使用弦长参数化确定参数的全局插值结果。插值曲面的次数为(3,3)。红色和蓝色曲线是节点处的等参数曲线。
虽然数据点只是稍微移动了位置，但蓝色节点曲线的形状发生了剧烈变化。两个插值曲面之间的差异更大。将数据点向右移动会在曲面左端产生了一个大的凸起。

最后两行显示了全局逼近的结果。$u$和$v$方向的控制点数量分别为5和4。显然，改变数据点位置对逼近曲面形状的影响没有那么明显。

结果表明，逼近曲面对数据点位置的改变不那么敏感，并且逼近曲面比插值曲面更好地贴近数据点网格的形状。

	移动前	移动后
插值曲线
插值曲面
逼近曲线
逼近曲面