Interpolation and Approximation

1. 参数选择

参数选择的目的是：对于给定数据点$D_0,D_1,...,D_n$，在曲线的定义域内找到满足条件$D_k=C(t_k),0<=k<=n$的$n+1$个参数点。
参数的选择会影响曲线的形状，进而影响曲线的参数化。

1.1 The Uniformly Spaced Method

区间为[0,1]时:

区间为[a,b]时:

当数据点不是均匀分布时，使用等间距参数可能会产生不规则的形状，例如大凸起、尖峰和环路。

1.2 The Chord Length Method

如果一条插值曲线非常接近数据多边形，则相邻两个数据点之间的曲线段长度将非常接近这两个数据点的弦长，并且插值曲线的长度也将非常接近数据多边形的总长度.
计算公式：
区间为[0,1]时:


区间为[a,b]时:

弦长方法通常来说表现不错，但有时候，当弦长较长时，会导致曲线的凸出量较大，如下图所示：

1.3 The Centripetal Method

向心法是弦长法的扩展，相邻点的距离通过$|D_k-D_{k-1}{|}^a$衡量，其中$a>0$.
计算公式如下：



讨论:
当$a=1$时，向心法退化为弦长法；

当$a<1$时，$|D_k-D_{k-1}{|}^a<|D_k-D_{k-1}|$，较长弦（即长度>1）对数据多边形长度的影响减小，而较短弦（即长度<1）对数据多边形长度的影响增加。由于这个特点，向心法比弦长法更能处理尖峰。
效果图如下：

但向心法也不一定效果比其他两种方法好，如下图所示，当数据点分布较为对称时，向心法的效果不如等距法好：

2. 节点向量生成

2.1 Knot Vector Generation

得到一组参数后，就可以生成一个节点向量.
假设得到的$n+1$个参数分别为：$t_0,t_1,...,t_n$，对于$p$次B样条曲线，需要$m+1$个节点，其中$m=n+p+1$. 如果曲线为Clamped，则$u_0=u_1=...=u_p=0，u_{m-p}=u_{m-p+1}=...=u_m=1$，剩余的$n-p$个节点可采用均匀分布或者其他方法。

均匀分布：

等间距的节点向量并不需要知道参数$t$，而且生成起来非常简单。但是，不建议使用这种方法，因为如果将其与弦长法一起用于全局插值，则线性方程组将是奇异的。

de Boor提出的参数平均法：

2.2 The Universal Method

在前面讨论的方法中，都是先确定参数，然后生成节点向量。
Choong-Gyoo Lim提出了一种用等间距节点来计算参数的通用方法。
计算方法如下：
首先定义均匀分布的节点向量：
由节点向量可以得到n+1个B样条基函数；
接下来选择使每个基函数达到最大值时，所对应的点作为参数。
如下图所示，黄色点即为选择的参数点：

3. 曲面的参数与节点向量

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

假设存在$m+1$行，$n+1$列数据点，$D_{ij}$，其中，$0<=i<=m$，$0<=j<=n$，从而$u$方向上需要$m+1$个参数，$s_0,...,s_m$，$v$方向上需要$n+1$个参数，$t_0,...,t_n$，使得曲面上的点$S(s_c,t_d)$与数据点$D_{cd}$相对应。

参数计算方法：

1. 分别对每一列数据点，如第$j$列，$D_{0,j},D_{1,j},...,D_{m,j}$计算$m+1$个参数，得到参数$u_{0,j},u_{1,j},...,u_{m,j}$，如下图所示，
   
   得到$u$方向的参数矩阵后，对每一行求平均，得到$s_0,s_1,...,s_m$.
2. 同理，分别对每一行数据点，如第$i$行，$D_{i,0},D_{i,1},...,D_{i,n}$，计算$n+1$个参数，得到参数$v_{i,0},u_{i,1},...,u_{i,n}$;
   得到$v$方向的参数矩阵后，对每一列求平均，得到$t_0,t_1,...,t_n$.

算法如下：

接下来，
由参数值$s_0,s_1,...,s_m$和次数$p$，可以得到$u$方向上的节点向量$U$;
由参数值$t_0,t_1,...,t_n$和次数$q$，可以得到$v$方向上的节点向量$V$.

4. 线性方程组的求解

假设 $n$×$n$的系数矩阵$A$，$n$×$h$的常数项矩阵$B$，待求解的$n$×$h$矩阵$X$定义如下：

并且满足如下关系：

LU分解：
将矩阵$A$分解为：$A=L.U$，其中$L$为下三角矩阵，$U$为上三角矩阵.
如图所示：

代入原方程，可得：
$$B=(L.U).X$$

$$B=L.(U.X)$$

令$Y=U.X$，得，
$$B=L.Y$$

因而，方程$B=A.X$的求解可分为两步：

1. 由方程$B=L.Y$求解$Y$.
2. 由方程$Y=U.X$求解$X$.

正向代换：
$$B=L.Y$$
计算过程：




算法如下：

反向代换：
$$Y=U.X$$
计算过程：





算法如下：

5. 曲线插值

问题描述：给定一组$n+1$个数据点$D_0,D_2,...,D_n$和次数$p$，求一条由$n+1$个控制点定义的$p$次B样条曲线，该曲线按给定顺序通过所给的数据点。

求解：
假设数据点和控制点均为$s$维空间中的点，则有，





求解得到矩阵$P$后便可得到控制点，从而得到所求的B样条插值曲线。

算法：
一列一列求解$P$


插值效果图：

注意：
如果节点是通过将连续$p$个参数求平均的方式获得的，则矩阵$N$是正定的，并且有：当$|i-k|>=p$时，$N_{i,p}(t_k)=0$（de Boor,1978证明）.这表明上述算法得到的线性方程组可以用高斯消元法求解。

参数和节点对所构建曲线的影响:
如果弦长分布大致相同，那么所有四种参数选择方法的性能应该是相似的。此外，由于具有均匀节点的B样条基函数的最大值分布非常均匀，因此通用方法的性能应该与均匀方法相似。同样，向心法的工作原理应该类似于弦长法，因为前者是后者的延伸。然而，当弦长分布发生剧烈变化时，情况并非总是如此。
如下图所示：

参数和节点的关系：

与弦长法和向心法相比，通用法得到的参数和节点分布更均匀。此外，从向心法得到的参数和结使较短（resp，longer）的弦拉伸得更长（resp，short），因此分布更均匀。因此，用向心法得到的弦长曲线，在曲线上不会产生较大的晃动。

次数的影响：

等间距法和通用法通常能很好地跟踪较长的弦，但对于较短的弦则存在问题。因为参数的间距相等或几乎相等，所以对于较短的弦，插值曲线必须拉伸得稍微长一点。结果，我们看到了峰和环。高次曲线会使这种情况变得更糟，因为高阶曲线提供了更多的自由度。

至于弦长法，上面的图表明，对于较长的弦，尤其是后面或前面有许多短弦的弦，可能会出现大的凸起，因此，弦长法并不适用。阶数对上述插值曲线的形状没有显著影响。

由于向心法是弦长法的扩展，它们具有相同的特点。然而，由于向心法具有使两个相邻参数之间的距离趋于均匀的趋势，它也具有均与法和通用法的特点。例如，生成的插值曲线紧跟长弦，当阶数增加时，短弦可能会出现循环。事实上，上述结果表明向心法和通用法的计算结果是相似的。

插值的全局性：
改变单个数据点的位置会完全改变插值曲线的形状

6. 曲线逼近

问题描述：
给定一组$n+1$个数据点$D_0,D_2,...,D_n$和次数$p$，控制点个数参数$h$，求一条由$h+1$个控制点定义的$p$次B样条曲线，该曲线满足以下两个条件：

1. 曲线经过第一个和最后一个数据点$D_0,D_n$.
2. 该曲线在最小二乘意义下逼近数据多边形.

求解：
1.曲线模型：


2.最小二乘意义下的目标函数推导：






3.对单个控制点$P_g$求导:






令导数为0可得，

4.矩阵形式：

又由，


可得：

求解方程得到$P$，即可得到控制点，从而得到所求曲线。

算法：

次数和控制点个数对所构建曲线的影响:

低次曲线通常不能很好地逼近数据多边形,次数越高的曲线产生的结果越好；同样，控制点越多，逼近曲线的灵活性就越高。
但并不是曲线次数越高，控制点个数越多就好，因为全局逼近一般比全局插值使用更少的控制点。如果控制点的个数等于数据点的个数，全局逼近就变成了全局插值，就可以用全局插值了！至于次数，只要生成的曲线能捕捉到数据多边形的形状即可。
逼近的全局性:
更改数据点的位置会导致整个曲线发生变化。

7. 曲面插值

问题描述：
给定一个由$(m+1)\times (n+1)$个数据点$D_{ij}(0<=i<=m,0<=j<=n)$组成的网格，次数$(p,q)$，求一个由$(m+1)\times (n+1)$控制点定义的，按给定顺序经过所有数据点的B样条曲面。

推导过程：

控制点的求解分以下两步：

1. 对于矩阵$D$的每一列数据点，如第$d$列，由次数$p$，参数$s_0,s_1,...s_m$，数据点$D_{0,d},D_{1,d},...,D_{m,d}$，求得第$d$列的控制点$Q_{0,d},Q_{1,d},...,Q_{m,d}$，得到$m+1$行，$n+1$列的矩阵$Q$；
2. 对于矩阵$Q$的每一行数据点，如第$i$行，由次数$q$，参数$t_0,t_1,...t_n$，数据点$Q_{i,0},Q_{i,1},...,Q_{i,n}$，求得第$i$行的控制点$P_{i,0},P_{i,1},...,P_{i,n}$.

算法如下：

由$m+1$行，$n+1$列的控制点矩阵$P_{ij}$，次数$(p,q)$，节点向量$U,V$，可以得到所求的B样条曲面。
插值效果图：

插值的全局性：
由于插值曲面是通过m+n+2条全局曲线插值得到的，因此这种曲面插值是全局的。

8. 曲面逼近

问题描述：
给定一组$(m+1)\times (n+1)$个数据点$D_{ij}(0<=i<=m，0<=j<=n)$，次数$(p,q)$，且满足$m>e>=p>=1$和$n>f>=q>=1$的$e$和$f$，找到由$(e+1)\times (f+1)$控制点$P_{ij}$$(0<=i<=e，0<=j<=f)$定义的，次数为$(p,q)$的B样条曲面，该曲面按给定的顺序逼近于数据点网格。

求解：
$u$方向的$m+1$个参数分别为$s_0,s_1,...,s_m$，$v$方向上的$n+1$个参数分别为$t_0,t_1,...,t_n$.


误差函数：

对每个控制点求偏微分：

从而可以得到$(e+1)\times (f+1)$个非线性方程，求解非线性方程组非常耗时。

我们可以找到一个不使函数$f$最小化的合理的解，而不是寻求最优解。

方法：
采用和曲面插值类似的方法，控制点的求解分为以下两步：

1. 对$D$中的每一列$m+1$数据点，进行曲线逼近，分别得到$e+1$个中间数据点，由于共有$n+1$列，从而可以得到$(e+1)\times (n+1)$的中间数据点网格$Q$；
2. 对$Q$中的每一行$n+1$数据点，进行曲线逼近，分别得到$f+1$个控制点，共有$e+1$行，从而可以得到$(e+1)\times (f+1)$的控制点矩阵$P$.

算法如下：

由计算得到的$(e+1)\times (f+1)$控制点矩阵$P_{ij}$，次数$(p,q)$，节点向量$U$，$V$，可以得到给定数据点的近似B样条曲面.
效果图如下：

尽管该算法不是最小化函数$f$，不是一个最优算法，它仍适合于许多应用。