0 引言

拟合的方式分为插值和逼近两种。拟合问题的输入通常包括一些几何数据，例如点和导矢。输出时Nurbs曲线曲面，也就是它的控制点、节点矢量和权因子。此外，次数 $p$ (曲面为$(p,q)$)必须由适用者输入或由算法选择一个合适的值。如果要求曲线是$C^r$连续的，则选定的次数$p$必须满足$p>=r+1$。(假定内节点的重复度均为1)。如果没有其他要求，对于插值来说选择$p=r+1$一般就足够了。对于逼近来说$p>r+1$可能产生更好的结果。

拟合算法分为两类

• 全局拟合算法——通常需要求解一个方程组或优化问题。若给定的数据只有点和导矢，并且只有控制点是未知的(次数、节点矢量和权因子都已预先选定)，那么这个系统是线性的。如果是更复杂的数据如曲率，被指定，或者节点和、或权重也是系统的未知量，那么得到的系统是非线性的。
• 局部拟合算法——在本质上更几何化一些：分段地创建曲线或曲面，每一步只用到局部的数据。

以下仅介绍曲线曲面的全局拟合

首先给出非有理B样条曲线曲面的公式

B样条基函数

令$U=U_0,u_1,...u_m$是一个单调不减的实例序列,即$u_i\le u_{i+1},i=0,1,...,m-1$.其中,$u_i$称为节点，$U$称为节点矢量，用$N_{i,p}(u)$表示第$i$个$p$次$(p+1$阶$)B$样条基函数，其定义为
$$\begin{gathered} {N}_{i,0}(u)=\{\begin{array}{l}1，u_i\le u\le u_{i+1}\\ 0，其他\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{gathered}$$
其中，

1.$N_{i,0}(u)$是一个阶梯函数，他在半开区间$u\in [u_i,u_{i+1})$外都为零；

2.当$p>0$时，$N_{i,p(u)}$是两个$p-1$次基函数的线性组合；

3.计算一组基函数时需要事先指定节点矢量$U$和次数$p$；

4.规定$0/0=0$；

B样条基函数的性质

1.局部支撑性：若$u\notin [u_i,u_{i+p+1})$，则$N_{i,p}(u)=0$；

2.在任意给定的节点区间$[u_j,u_{j+1})$内，最多$p+1$个$N_{i,p}$是非零的，他们是$N_{j-p,p},...,N_{j,p}$；

3.非负性：对于所有的$i,p$和$u$,有$N_{i,p}(u)\ge 0$

4.规范性：对于任意的节点区间$[u_i,u_{i+1})$,当$u\in [u_i,ui+1)$时,${\sum }_{j=i-p}^i(u)=1$.

5.等其他性质

1 全局插值

1 给定数据点的全局曲线插值

输入：

• 一组数据点，$Q_k,k=0,1,...,n$
• 每一个数据的参数值${\overline{u}}_k$
• 节点矢量$U=u_0,u_1,...,u_m$

输出：

• $n+1$个控制点$P_i$

求解过程：

由输入，可以建立一个系数矩阵为$(n+1)\times (n+1)$线性方程组
$$Q_k=C({\overline{u}}_k)=\sum _{i=0}^n{N}_{i,p}({\overline{u}}_k)P_i$$

• 确定参数值${\overline{u}}_k$
• 确定节点矢量$U$

确定参数值${\overline{u}}_k$,假定$u\in [0,1]$

• 均匀参数化：不推荐

• 弦长参数化：

  令总弦长$d={\sum }_{k=1}^n|Q_k-Q_{k-1}|$,则${\overline{u}}_0=0,{\overline{u}}_n=1$,
  $${\overline{u}}_k={\overline{u}}_{k-1}+\frac{|Q_k-Q_{k-1}|}{d},k=1,2,...,n-1$$

• 向心参数化:

  令总弦长$d={\sum }_{k=1}^n|Q_k-Q_{k-1}|$,则${\overline{u}}_0=0,{\overline{u}}_n=1$,
  $${\overline{u}}_k={\overline{u}}_{k-1}+\frac{\sqrt{|Q_k-Q_{k-1}|}}{d},k=1,2,...,n-1$$

确定节点矢量U

• 等距分布(不推荐，弦长参数和向心参数联立可能会产生奇异方程组)：
  $$\begin{gathered} u_0=...=u_p=0,u_{m-p}=...=u_m=1, \\ {u}_{j+p}=\frac{j}{n-p+1},j=1,2,...,n-p \end{gathered}$$

• 取平均值法：
  $$\begin{gathered} u_0=...=u_p=0,u_{m-p}=...=u_m=1, \\ {u}_{j+p}=\frac{1}{p}\sum _{i=j}^{j+p+1}{\overline{u}}_i,j=1,2,...,n-p \end{gathered}$$

2 端点导矢指定的全局曲线插值

输入：

• 一组数据点，$Q_k,k=0,1,...,n$
• 每一个数据的参数值${\overline{u}}_k$
• 节点矢量$U=u_0,u_1,...,u_m$
• 首末端点的一阶导矢$D_0$和$D_n$

输出：

• $n+1$个控制点$P_i$

求解过程：

由于每个导矢都会增加一个节点和一个控制点，即增加一个线性方程。可以建立一个系数矩阵为$(n+3)\times (n+3)$线性方程组
$$C({\overline{u}}_k)=\sum _{i=0}^{n+2}{N}_{i,p}({\overline{u}}_k)P_i$$

• 确定参数值${\overline{u}}_k$，这里只是多了两个节点。
• 确定节点矢量$U$

3 三次样条曲线插值

上述方法使用于次数$p>1$的情况，当p=3时有一个效率更高的方法生成常用的C^2连续的三次样条曲线。${\overline{u}}_k$的计算和前面一样。节点矢量如下确定
$$\begin{gathered} u_0=...=u_3=0,u_{n+3}=...=u_{n+6}=1, \\ {u}_{j+3}={\overline{u}}_j,j=1,2,...,n-1 \end{gathered}$$
曲线刚好在节点处插值于$Q_k$，开始的两个方程和最后两个方程分别为
$$\begin{gathered} P_0=Q_0 \\ -P_0+P_1=\frac{u_4}{3}{D}_0 \\ -P_{n+1}+P_{n+2}=\frac{1-u_{n+2}}{3}{D}_n,P_{n+2}=Q_n \end{gathered}$$

在每一个内节点处，只有三个非零的三次基函数，所以剩下的$n-1$个方程可以写为如下形式

$$Q_k=C({\overline{u}}_k)=N_{k,3}({\overline{u}}_k)P_k+N_{k+1,3}({\overline{u}}_k)P_{k+1}+N_{k+2,3}({\overline{u}}_k)P_{k+2}$$

其中$k=1,\mathrm{2...},n-1$，令

$$a_k=N_{k,3}({\overline{u}}_k),b_k=N_{k+1,3}({\overline{u}}_k),c_k=N_{k+2,3}({\overline{u}}_k)$$

则得到三对角方程组系统

$$\left[\begin{array}{c}{R}_2\\ R_3\\ ...\\ R_{n-1}\\ R_n\end{array}\right]=\left[\begin{array}{c}{Q}_1-a_1{P}_1\\ Q_2\\ ...\\ Q_{n-2}\\ Q_{n-1}-c_{n-1}{P}_{n=1}\end{array}\right]=\left[\begin{array}{ccccccc}{b}_1& c_1& 0& ...& 0& 0& 0\\ a_2& b_2& c_2& ...& 0& 0& 0\\ ...& & & & & & \\ 0& 0& 0& ...& a_{n-2}& b_{n-2}& c_{n-2}\\ 0& 0& 0& ...& 0& a_{n-1}& b_{n-1}\end{array}\right]\left[\begin{array}{c}{P}_2\\ P_3\\ ...\\ P_{n-1}\\ P_n\end{array}\right]$$

4 全局曲面插值

给定$(n+1)\times (m+1)$个数据点$Q_{k,l},k=0,1,...,n$和$l=0,1,...,m$，创建一个非有理的$(p,q)$次$B$样条曲面使其插值与这些点，即
$$Q_{k,l}=S({\overline{u}}_k,{\overline{v}}_l)=\sum _{i=0}^n\sum _{i=0}^m{N}_{i,p}({\overline{u}}_k)N_{j,q}({\overline{v}}_l)P_{i,j}$$
计算合理的$({\overline{u}}_k,{\overline{v}}_l)$值以及节点矢量$U$和$V$，一个常用的方法是对每个$l$，按照曲线插值的方法计算参数${\overline{u}}_0^l,{\overline{u}}_1^l...,{\overline{u}}_n^l$，然后通过对所有的${\overline{u}}_k^l(l=0,1,...,m)$取平均值得到${\overline{u}}_k^l$,即
$${\overline{u}}_k=\frac{1}{m+1}\sum _{l=0}^m{\overline{u}}_k^l,k=0,1,...,n$$

考虑控制点的计算，$S(u,v)$是张量积曲面，$P_{i,j}$可以通过一系列曲线插值来更简单和更高效地得到。对于固定的$l$
$$Q_{k,l}=\sum _{i=0}^n{N}_{i,p}(\sum _{j=0}^m{N}_{j,q}{P}_{i,j})=\sum _{i=0}^n{N}_{i,p}({\overline{u}}_k)R_{i,l}$$
其中
$$R_{i,l}=\sum _{j=0}^m{N}_{j,q}{P}_{i,j}$$
$Q_{k,l}(k=0,...,n)$恰好是$R_{i,l}$的数据点，$R_{i,l}$恰好是$P_{i,j}$的数据点。

算法步骤：

1.用节点矢量$U$和参数${\overline{u}}_k$做$m+1$次曲线插值：对于$l=0,...,m$，分别构造插值于点$Q_{0,l},...,Q_{n,l}$的曲线，由此得到$R_{i,l}$;

2.用节点矢量$V$和参数${\overline{v}}_l$做$n+1$次曲线插值：对于$i=0,...,n$，分别构造插值于$R_{i,0},...,R_{i,m}$的曲线，得到曲面的控制点$P_{i,j}$;

显然，此算法关于参数$u$和$v$是对称的，因而同样的曲面还可以按照如下方式得到：

1.做$n+1$次曲线插值(对于$k=0,...,n$，分别构造插值于$Q_{k,0},...,Q_{k,m}$的曲线)得到点$R_{k,j}$(等参曲线$S({\overline{u}}_k,v)$的控制点)；

2.然后做$m+1$次曲线插值(对于$j=0,...,m$，分别构造插值于$R_{0,j},...,R_{n,j}$的曲线)得到曲面的控制点$P_{i,j}$;

3 全局逼近

全局逼近的两种方式：

1 从最小的开始

• 从最少的或较少数量的控制点开始
• 使用一种全局拟合方法，构造一条逼近于所给数据点的曲线或曲面；
• 检查曲线或曲面与给定点数据的偏差；
• 如果每点处的偏差都小于E，返回；否则，增加控制点的数目然后转到步骤1b；

2 从最多的开始

1 最小二乘曲线逼近

输入:

• $Q_0,...Q_m(m>n,n\ge p,p\ge 1)$

试图找到一条$p$次非有理曲线
$$
C(u)=\sum^n_{i=0}N_{i,p}(u)P_i,u\in[0,1]

$$
满足条件：

• $Q_0=C(0),Q_m=C(1)$

• 其余的数据点$Q_k$在最小二乘的意义下被逼近，即
  $$\sum _{k=1}^{m-1}|Q_k-C({\overline{u}}_k){|}^2$$
  关于n+1个变量$P_i$达到最小；${\overline{u}}_k$是预先计算好的参数值。注意生成的曲线一般不精确地通过数据点$Q_k$，并且，一般$C({\overline{u}}_k)$不是曲面上与$Q_k$最接近的点。令
  $$R_k=Q_k-N_{0,p}({\overline{u}}_k)Q_0-N_{n,p}({\overline{u}}_k)Q_m,k=1,...,m-1$$
  然后令
  $$\begin{gathered} f=\sum _{k=1}^{m-1}|Q_k-C({\overline{u}}_k){|}^2 \\ =\sum _{k=1}^{m-1}|R_k-\sum _{i=1}^{n-1}{N}_{i,p}({\overline{u}}_k)P_i{|}^2 \\ =\sum _{k=1}^{m-1}[R_k\cdot R_k-2\sum _{i=1}^{n-1}{N}_{i,p}({\overline{u}}_k)(R_k\cdot P_i)+(\sum _{i=1}^{n-1}{N}_{i,p}({\overline{u}}_k)P_i)(\sum _{i=1}^{n-1}{N}_{i,p}({\overline{u}}_k)P_i)] \end{gathered}$$
  $f$是关于$n-1$个变量$P_1,...,P_{n-1}$的标量值函数，应用标准的线性最小二乘拟合，为使目标函数$f$最小，令$f$关于$n-1$个未知控制点$P_l$的偏导都等于零。它的第$l$个偏导为
  $$\frac{\partial f}{\partial P_l}=\sum _{k=1}^{m-1}(-2N_{l,p}({\overline{u}}_k)R_k+2N_{l,p}({\overline{u}}_k)\sum _{i=1}^{n-1}{N}_{i,p}({\overline{u}}_k)P_i)$$

  $$\sum _{i=1}^{n-1}(\sum _{k=1}^{m-1}{N}_{l,p}({\overline{u}}_k)N_{i,p}({\overline{u}}_k))P_i=\sum _{k=1}^{m-1}{N}_{l,p}({\overline{u}}_k)R_k$$

  式中线性方程组以控制点$P_1,...,P_{n-1}$为未知量，$l=1,2,...,n-1$,则得到一个含$n-1$个未知量和$n-1$个方程的线性方程组
  $$(N^{T}N)P=R$$
  这里，$N$是由标量组成的$(m-1)\times (n-1)$的矩阵
  $$N=\left[\begin{array}{ccc}{N}_{1,p}({\overline{u}}_1)& ...& N_{n-1,p}({\overline{u}}_1)\\ & ...& \\ N_{1,p}({\overline{u}}_{m-1})& ...& N_{n-1,p}({\overline{u}}_{m-1})\end{array}\right]$$
  $R$是由$n-1$个点组成的列向量：
  $$R=\left[\begin{array}{cc}{N}_{1,p}({\overline{u}}_1)R_1+...+N_{1,p}({\overline{u}}_1)R_{m-1}& \\ & ...\\ N_{n-1,p}({\overline{u}}_{m-1})R_1+...+N_{n-1,p}({\overline{u}}_{m-1})R_{m-1}& \end{array}\right]$$

2 带权的约束最小二乘曲线拟合

3 最小二乘曲面逼近

4 在规定精度内的逼近