B样条曲线曲面拟合
0 引言
拟合的方式分为插值和逼近两种。拟合问题的输入通常包括一些几何数据,例如点和导矢。输出时Nurbs曲线曲面,也就是它的控制点、节点矢量和权因子。此外,次数 p p p (曲面为 ( p , q ) (p,q) (p,q))必须由适用者输入或由算法选择一个合适的值。如果要求曲线是 C r C^r Cr连续的,则选定的次数 p p p必须满足 p > = r + 1 p>=r+1 p>=r+1。(假定内节点的重复度均为1)。如果没有其他要求,对于插值来说选择 p = r + 1 p=r+1 p=r+1一般就足够了。对于逼近来说 p > r + 1 p>r+1 p>r+1可能产生更好的结果。
拟合算法分为两类
- 全局拟合算法——通常需要求解一个方程组或优化问题。若给定的数据只有点和导矢,并且只有控制点是未知的(次数、节点矢量和权因子都已预先选定),那么这个系统是线性的。如果是更复杂的数据如曲率,被指定,或者节点和、或权重也是系统的未知量,那么得到的系统是非线性的。
- 局部拟合算法——在本质上更几何化一些:分段地创建曲线或曲面,每一步只用到局部的数据。
以下仅介绍曲线曲面的全局拟合
首先给出非有理B样条曲线曲面的公式
B样条基函数
令
U
=
U
0
,
u
1
,
.
.
.
u
m
U={U_0,u_1,...u_m}
U=U0,u1,...um是一个单调不减的实例序列,即
u
i
≤
u
i
+
1
,
i
=
0
,
1
,
.
.
.
,
m
−
1
u_i\leq u_{i+1},i=0,1,...,m-1
ui≤ui+1,i=0,1,...,m−1.其中,
u
i
u_i
ui称为节点,
U
U
U称为节点矢量,用
N
i
,
p
(
u
)
N_{i,p}(u)
Ni,p(u)表示第
i
i
i个
p
p
p次
(
p
+
1
(p+1
(p+1阶
)
B
)B
)B样条基函数,其定义为
N
i
,
0
(
u
)
=
{
1
,
u
i
≤
u
≤
u
i
+
1
0
,其他
N
i
,
p
(
u
)
=
u
−
u
i
u
i
+
p
−
u
i
N
i
,
p
−
1
(
u
)
+
u
i
+
p
+
1
−
u
u
i
+
p
+
1
−
u
i
+
1
N
i
+
1
,
p
−
1
(
u
)
\\ N_{i,0}(u)= \begin{cases} 1,u_i\leq u \leq u_{i+1} \\ 0, 其他 \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)
Ni,0(u)={1,ui≤u≤ui+10,其他Ni,p(u)=ui+p−uiu−uiNi,p−1(u)+ui+p+1−ui+1ui+p+1−uNi+1,p−1(u)
其中,
1. N i , 0 ( u ) N_{i,0}(u) Ni,0(u)是一个阶梯函数,他在半开区间 u ∈ [ u i , u i + 1 ) u\in [u_i,u_{i+1}) u∈[ui,ui+1)外都为零;
2.当 p > 0 p>0 p>0时, N i , p ( u ) N_{i,p(u)} Ni,p(u)是两个 p − 1 p-1 p−1次基函数的线性组合;
3.计算一组基函数时需要事先指定节点矢量 U U U和次数 p p p;
4.规定 0 / 0 = 0 0/0=0 0/0=0;
B样条基函数的性质
1.局部支撑性:若 u ∉ [ u i , u i + p + 1 ) u\notin [u_i,u_{i+p+1}) u∈/[ui,ui+p+1),则 N i , p ( u ) = 0 N_{i,p}(u)=0 Ni,p(u)=0;
2.在任意给定的节点区间 [ u j , u j + 1 ) [u_j,u_{j+1}) [uj,uj+1)内,最多 p + 1 p+1 p+1个 N i , p N_{i,p} Ni,p是非零的,他们是 N j − p , p , . . . , N j , p N_{j-p,p},...,N_{j,p} Nj−p,p,...,Nj,p;
3.非负性:对于所有的 i , p i,p i,p和 u u u,有 N i , p ( u ) ≥ 0 N_{i,p}(u)\geq 0 Ni,p(u)≥0
4.规范性:对于任意的节点区间 [ u i , u i + 1 ) [u_i,u_{i+1}) [ui,ui+1),当 u ∈ [ u i , u i + 1 ) u\in[u_i,ui+1) u∈[ui,ui+1)时, ∑ j = i − p i ( u ) = 1 \sum^i_{j=i-p}(u)=1 ∑j=i−pi(u)=1.
5.等其他性质
1 全局插值
1 给定数据点的全局曲线插值
输入:
- 一组数据点, Q k , k = 0 , 1 , . . . , n {Q_k},k=0,1,...,n Qk,k=0,1,...,n
- 每一个数据的参数值 u ‾ k \overline{u}_k uk
- 节点矢量 U = u 0 , u 1 , . . . , u m U={u_0,u_1,...,u_m} U=u0,u1,...,um
输出:
- n + 1 n+1 n+1个控制点 P i P_i Pi
求解过程:
由输入,可以建立一个系数矩阵为
(
n
+
1
)
×
(
n
+
1
)
(n+1)\times(n+1)
(n+1)×(n+1)线性方程组
Q
k
=
C
(
u
‾
k
)
=
∑
i
=
0
n
N
i
,
p
(
u
‾
k
)
P
i
Q_k=C(\overline{u}_k)=\sum^n_{i=0}N_{i,p}(\overline{u}_k)P_i
Qk=C(uk)=i=0∑nNi,p(uk)Pi
- 确定参数值 u ‾ k \overline{u}_k uk
- 确定节点矢量 U U U
确定参数值 u ‾ k \overline{u}_k uk,假定 u ∈ [ 0 , 1 ] u\in [0,1] u∈[0,1]
-
均匀参数化:不推荐
-
弦长参数化:
令总弦长 d = ∑ k = 1 n ∣ Q k − Q k − 1 ∣ d=\sum^n_{k=1}|Q_k-Q_{k-1}| d=∑k=1n∣Qk−Qk−1∣,则 u ‾ 0 = 0 , u ‾ n = 1 \overline{u}_0=0,\overline{u}_n=1 u0=0,un=1,
u ‾ k = u ‾ k − 1 + ∣ Q k − Q k − 1 ∣ d , k = 1 , 2 , . . . , n − 1 \overline{u}_k=\overline{u}_{k-1}+\frac{|Q_k-Q_{k-1}|}{d},k=1,2,...,n-1 uk=uk−1+d∣Qk−Qk−1∣,k=1,2,...,n−1 -
向心参数化:
令总弦长 d = ∑ k = 1 n ∣ Q k − Q k − 1 ∣ d=\sum^n_{k=1}|Q_k-Q_{k-1}| d=∑k=1n∣Qk−Qk−1∣,则 u ‾ 0 = 0 , u ‾ n = 1 \overline{u}_0=0,\overline{u}_n=1 u0=0,un=1,
u ‾ k = u ‾ k − 1 + ∣ Q k − Q k − 1 ∣ d , k = 1 , 2 , . . . , n − 1 \overline{u}_k=\overline{u}_{k-1}+\frac{\sqrt{|Q_k-Q_{k-1}|}}{d},k=1,2,...,n-1 uk=uk−1+d∣Qk−Qk−1∣,k=1,2,...,n−1
确定节点矢量U
-
等距分布(不推荐,弦长参数和向心参数联立可能会产生奇异方程组):
u 0 = . . . = u p = 0 , u m − p = . . . = u m = 1 , u j + p = j n − p + 1 , j = 1 , 2 , . . . , n − p u_0=...=u_p=0,u_{m-p}=...=u_m=1,\\ u_{j+p}=\frac{j}{n-p+1},j=1,2,...,n-p u0=...=up=0,um−p=...=um=1,uj+p=n−p+1j,j=1,2,...,n−p -
取平均值法:
u 0 = . . . = u p = 0 , u m − p = . . . = u m = 1 , u j + p = 1 p ∑ i = j j + p + 1 u ‾ i , j = 1 , 2 , . . . , n − p u_0=...=u_p=0,u_{m-p}=...=u_m=1,\\ u_{j+p}=\frac{1}{p}\sum^{j+p+1}_{i=j}\overline{u}_i,j=1,2,...,n-p u0=...=up=0,um−p=...=um=1,uj+p=p1i=j∑j+p+1ui,j=1,2,...,n−p
2 端点导矢指定的全局曲线插值
输入:
- 一组数据点, Q k , k = 0 , 1 , . . . , n {Q_k},k=0,1,...,n Qk,k=0,1,...,n
- 每一个数据的参数值 u ‾ k \overline{u}_k uk
- 节点矢量 U = u 0 , u 1 , . . . , u m U={u_0,u_1,...,u_m} U=u0,u1,...,um
- 首末端点的一阶导矢 D 0 D_0 D0和 D n D_n Dn
输出:
- n + 1 n+1 n+1个控制点 P i P_i Pi
求解过程:
由于每个导矢都会增加一个节点和一个控制点,即增加一个线性方程。可以建立一个系数矩阵为
(
n
+
3
)
×
(
n
+
3
)
(n+3)\times(n+3)
(n+3)×(n+3)线性方程组
C
(
u
‾
k
)
=
∑
i
=
0
n
+
2
N
i
,
p
(
u
‾
k
)
P
i
C(\overline{u}_k)=\sum^{n+2}_{i=0}N_{i,p}(\overline{u}_k)P_i
C(uk)=i=0∑n+2Ni,p(uk)Pi
- 确定参数值 u ‾ k \overline{u}_k uk,这里只是多了两个节点。
- 确定节点矢量 U U U
3 三次样条曲线插值
上述方法使用于次数
p
>
1
p>1
p>1的情况,当p=3时有一个效率更高的方法生成常用的C^2连续的三次样条曲线。
u
‾
k
\overline{u}_k
uk的计算和前面一样。节点矢量如下确定
u
0
=
.
.
.
=
u
3
=
0
,
u
n
+
3
=
.
.
.
=
u
n
+
6
=
1
,
u
j
+
3
=
u
‾
j
,
j
=
1
,
2
,
.
.
.
,
n
−
1
u_0=...=u_3=0,u_{n+3}=...=u_{n+6}=1,\\ u_{j+3}=\overline{u}_j,j=1,2,...,n-1
u0=...=u3=0,un+3=...=un+6=1,uj+3=uj,j=1,2,...,n−1
曲线刚好在节点处插值于
Q
k
Q_k
Qk,开始的两个方程和最后两个方程分别为
P
0
=
Q
0
−
P
0
+
P
1
=
u
4
3
D
0
−
P
n
+
1
+
P
n
+
2
=
1
−
u
n
+
2
3
D
n
,
P
n
+
2
=
Q
n
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
P0=Q0−P0+P1=3u4D0−Pn+1+Pn+2=31−un+2Dn,Pn+2=Qn
在每一个内节点处,只有三个非零的三次基函数,所以剩下的 n − 1 n-1 n−1个方程可以写为如下形式
Q k = C ( u ‾ k ) = N k , 3 ( u ‾ k ) P k + N k + 1 , 3 ( u ‾ k ) P k + 1 + N k + 2 , 3 ( u ‾ k ) P k + 2 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} Qk=C(uk)=Nk,3(uk)Pk+Nk+1,3(uk)Pk+1+Nk+2,3(uk)Pk+2
其中 k = 1 , 2... , n − 1 k=1,2...,n-1 k=1,2...,n−1,令
a k = N k , 3 ( u ‾ k ) , b k = N k + 1 , 3 ( u ‾ k ) , c k = N k + 2 , 3 ( u ‾ k ) 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) ak=Nk,3(uk),bk=Nk+1,3(uk),ck=Nk+2,3(uk)
则得到三对角方程组系统
[ R 2 R 3 . . . R n − 1 R n ] = [ Q 1 − a 1 P 1 Q 2 . . . Q n − 2 Q n − 1 − c n − 1 P n = 1 ] = [ 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 ] [ P 2 P 3 . . . P n − 1 P n ] \left[ \begin{matrix} R_2\\ R_3\\ ...\\ R_{n-1}\\ R_n \end{matrix} \right]= \left[ \begin{matrix} Q_1-a_1P_1\\ Q_2\\ ...\\ Q_{n-2}\\ Q_{n-1}-c_{n-1}P_{n=1} \end{matrix} \right]= \left[ \begin{matrix} 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{matrix} \right] \left[ \begin{matrix} P_2 \\ P_3 \\ ...\\ P_{n-1} \\ P_n \end{matrix} \right] R2R3...Rn−1Rn = Q1−a1P1Q2...Qn−2Qn−1−cn−1Pn=1 = b1a2...00c1b2000c200............00an−2000bn−2an−100cn−2bn−1 P2P3...Pn−1Pn
4 全局曲面插值
给定
(
n
+
1
)
×
(
m
+
1
)
(n+1)\times(m+1)
(n+1)×(m+1)个数据点
Q
k
,
l
,
k
=
0
,
1
,
.
.
.
,
n
{Q_{k,l}},k=0,1,...,n
Qk,l,k=0,1,...,n和
l
=
0
,
1
,
.
.
.
,
m
l=0,1,...,m
l=0,1,...,m,创建一个非有理的
(
p
,
q
)
(p,q)
(p,q)次
B
B
B样条曲面使其插值与这些点,即
Q
k
,
l
=
S
(
u
‾
k
,
v
‾
l
)
=
∑
i
=
0
n
∑
i
=
0
m
N
i
,
p
(
u
‾
k
)
N
j
,
q
(
v
‾
l
)
P
i
,
j
Q_{k,l} =S(\overline{u}_k,\overline{v}_l) =\sum^n_{i=0}\sum^m_{i=0}N_{i,p}(\overline{u}_k)N_{j,q}(\overline{v}_l)P_{i,j}
Qk,l=S(uk,vl)=i=0∑ni=0∑mNi,p(uk)Nj,q(vl)Pi,j
计算合理的
(
u
‾
k
,
v
‾
l
)
(\overline{u}_k,\overline{v}_l)
(uk,vl)值以及节点矢量
U
U
U和
V
V
V,一个常用的方法是对每个
l
l
l,按照曲线插值的方法计算参数
u
‾
0
l
,
u
‾
1
l
.
.
.
,
u
‾
n
l
\overline{u}_0^l,\overline{u}_1^l...,\overline{u}_n^l
u0l,u1l...,unl,然后通过对所有的
u
‾
k
l
(
l
=
0
,
1
,
.
.
.
,
m
)
\overline{u}_k^l(l=0,1,...,m)
ukl(l=0,1,...,m)取平均值得到
u
‾
k
l
\overline{u}_k^l
ukl,即
u
‾
k
=
1
m
+
1
∑
l
=
0
m
u
‾
k
l
,
k
=
0
,
1
,
.
.
.
,
n
\overline{u}_k=\frac{1}{m+1}\sum^m_{l=0}\overline{u}_k^l,k=0,1,...,n
uk=m+11l=0∑mukl,k=0,1,...,n
考虑控制点的计算,
S
(
u
,
v
)
S(u,v)
S(u,v)是张量积曲面,
P
i
,
j
P_{i,j}
Pi,j可以通过一系列曲线插值来更简单和更高效地得到。对于固定的
l
l
l
Q
k
,
l
=
∑
i
=
0
n
N
i
,
p
(
∑
j
=
0
m
N
j
,
q
P
i
,
j
)
=
∑
i
=
0
n
N
i
,
p
(
u
‾
k
)
R
i
,
l
Q_{k,l} =\sum^n_{i=0}N_{i,p}(\sum^m_{j=0}N_{j,q}P_{i,j}) =\sum^n_{i=0}N_{i,p}(\overline{u}_k)R_{i,l}
Qk,l=i=0∑nNi,p(j=0∑mNj,qPi,j)=i=0∑nNi,p(uk)Ri,l
其中
R
i
,
l
=
∑
j
=
0
m
N
j
,
q
P
i
,
j
R_{i,l} =\sum^m_{j=0}N_{j,q}P_{i,j}
Ri,l=j=0∑mNj,qPi,j
Q
k
,
l
(
k
=
0
,
.
.
.
,
n
)
Q_{k,l}(k=0,...,n)
Qk,l(k=0,...,n)恰好是
R
i
,
l
R_{i,l}
Ri,l的数据点,
R
i
,
l
R_{i,l}
Ri,l恰好是
P
i
,
j
P_{i,j}
Pi,j的数据点。
算法步骤:
1.用节点矢量 U U U和参数 u ‾ k \overline{u}_k uk做 m + 1 m+1 m+1次曲线插值:对于 l = 0 , . . . , m l=0,...,m l=0,...,m,分别构造插值于点 Q 0 , l , . . . , Q n , l Q_{0,l},...,Q_{n,l} Q0,l,...,Qn,l的曲线,由此得到 R i , l R_{i,l} Ri,l;
2.用节点矢量 V V V和参数 v ‾ l \overline{v}_l vl做 n + 1 n+1 n+1次曲线插值:对于 i = 0 , . . . , n i=0,...,n i=0,...,n,分别构造插值于 R i , 0 , . . . , R i , m R_{i,0},...,R_{i,m} Ri,0,...,Ri,m的曲线,得到曲面的控制点 P i , j P_{i,j} Pi,j;
显然,此算法关于参数 u u u和 v v v是对称的,因而同样的曲面还可以按照如下方式得到:
1.做 n + 1 n+1 n+1次曲线插值(对于 k = 0 , . . . , n k=0,...,n k=0,...,n,分别构造插值于 Q k , 0 , . . . , Q k , m Q_{k,0},...,Q_{k,m} Qk,0,...,Qk,m的曲线)得到点 R k , j R_{k,j} Rk,j(等参曲线 S ( u ‾ k , v ) S(\overline{u}_k,v) S(uk,v)的控制点);
2.然后做 m + 1 m+1 m+1次曲线插值(对于 j = 0 , . . . , m j=0,...,m j=0,...,m,分别构造插值于 R 0 , j , . . . , R n , j R_{0,j},...,R_{n,j} R0,j,...,Rn,j的曲线)得到曲面的控制点 P i , j P_{i,j} Pi,j;
3 全局逼近
全局逼近的两种方式:
1 从最小的开始
- 从最少的或较少数量的控制点开始
- 使用一种全局拟合方法,构造一条逼近于所给数据点的曲线或曲面;
- 检查曲线或曲面与给定点数据的偏差;
- 如果每点处的偏差都小于E,返回;否则,增加控制点的数目然后转到步骤1b;
2 从最多的开始
1 最小二乘曲线逼近
输入:
- Q 0 , . . . Q m ( m > n , n ≥ p , p ≥ 1 ) Q_0,...Q_m(m>n,n \geq p,p\geq 1) Q0,...Qm(m>n,n≥p,p≥1)
试图找到一条
p
p
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_0=C(0),Q_m=C(1) Q0=C(0),Qm=C(1)
-
其余的数据点 Q k Q_k Qk在最小二乘的意义下被逼近,即
∑ k = 1 m − 1 ∣ Q k − C ( u ‾ k ) ∣ 2 \sum^{m-1}_{k=1}|Q_k-C(\overline{u}_k)|^2 k=1∑m−1∣Qk−C(uk)∣2
关于n+1个变量 P i P_i Pi达到最小; u ‾ k {\overline{u}_k} uk是预先计算好的参数值。注意生成的曲线一般不精确地通过数据点 Q k Q_k Qk,并且,一般 C ( u ‾ k ) C(\overline{u}_k) C(uk)不是曲面上与 Q k Q_k Qk最接近的点。令
R k = Q k − N 0 , p ( u ‾ k ) Q 0 − N n , p ( u ‾ k ) Q m , k = 1 , . . . , m − 1 R_k=Q_k-N_{0,p}(\overline{u}_k)Q_0-N_{n,p}(\overline{u}_k)Q_m,k=1,...,m-1 Rk=Qk−N0,p(uk)Q0−Nn,p(uk)Qm,k=1,...,m−1
然后令
f = ∑ k = 1 m − 1 ∣ Q k − C ( u ‾ k ) ∣ 2 = ∑ k = 1 m − 1 ∣ R k − ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ∣ 2 = ∑ k = 1 m − 1 [ R k ⋅ R k − 2 ∑ i = 1 n − 1 N i , p ( u ‾ k ) ( R k ⋅ P i ) + ( ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) ( ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) ] f=\sum^{m-1}_{k=1}|Q_k-C(\overline{u}_k)|^2\\ =\sum^{m-1}_{k=1}|R_k-\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i|^2 \\ =\sum^{m-1}_{k=1}[R_k\cdot R_k-2\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)(R_k\cdot P_i)+(\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i)(\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i)] f=k=1∑m−1∣Qk−C(uk)∣2=k=1∑m−1∣Rk−i=1∑n−1Ni,p(uk)Pi∣2=k=1∑m−1[Rk⋅Rk−2i=1∑n−1Ni,p(uk)(Rk⋅Pi)+(i=1∑n−1Ni,p(uk)Pi)(i=1∑n−1Ni,p(uk)Pi)]
f f f是关于 n − 1 n-1 n−1个变量 P 1 , . . . , P n − 1 P_1,...,P_{n-1} P1,...,Pn−1的标量值函数,应用标准的线性最小二乘拟合,为使目标函数 f f f最小,令 f f f关于 n − 1 n-1 n−1个未知控制点 P l P_l Pl的偏导都等于零。它的第 l l l个偏导为
∂ f ∂ P l = ∑ k = 1 m − 1 ( − 2 N l , p ( u ‾ k ) R k + 2 N l , p ( u ‾ k ) ∑ i = 1 n − 1 N i , p ( u ‾ k ) P i ) \frac{\partial{f}}{\partial{P_l}}=\sum^{m-1}_{k=1}(-2N_{l,p}(\overline{u}_k)R_k+2N_{l,p}(\overline{u}_k)\sum^{n-1}_{i=1}N_{i,p}(\overline{u}_k)P_i) ∂Pl∂f=k=1∑m−1(−2Nl,p(uk)Rk+2Nl,p(uk)i=1∑n−1Ni,p(uk)Pi)∑ i = 1 n − 1 ( ∑ k = 1 m − 1 N l , p ( u ‾ k ) N i , p ( u ‾ k ) ) P i = ∑ k = 1 m − 1 N l , p ( u ‾ k ) R k \sum^{n-1}_{i=1}(\sum^{m-1}_{k=1}N_{l,p}(\overline{u}_k)N_{i,p}(\overline{u}_k))P_i=\sum^{m-1}_{k=1}N_{l,p}(\overline{u}_k)R_k i=1∑n−1(k=1∑m−1Nl,p(uk)Ni,p(uk))Pi=k=1∑m−1Nl,p(uk)Rk
式中线性方程组以控制点 P 1 , . . . , P n − 1 P_1,...,P_{n-1} P1,...,Pn−1为未知量, l = 1 , 2 , . . . , n − 1 l=1,2,...,n-1 l=1,2,...,n−1,则得到一个含 n − 1 n-1 n−1个未知量和 n − 1 n-1 n−1个方程的线性方程组
( N T N ) P = R (N^TN)P=R (NTN)P=R
这里, N N N是由标量组成的 ( m − 1 ) × ( n − 1 ) (m-1)\times(n-1) (m−1)×(n−1)的矩阵
N = [ N 1 , p ( u ‾ 1 ) . . . N n − 1 , p ( u ‾ 1 ) . . . N 1 , p ( u ‾ m − 1 ) . . . N n − 1 , p ( u ‾ m − 1 ) ] N=\left[ \begin{matrix} 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{matrix} \right] N= N1,p(u1)N1,p(um−1).........Nn−1,p(u1)Nn−1,p(um−1)
R R R是由 n − 1 n-1 n−1个点组成的列向量:
R = [ N 1 , p ( u ‾ 1 ) R 1 + . . . + N 1 , p ( u ‾ 1 ) R m − 1 . . . N n − 1 , p ( u ‾ m − 1 ) R 1 + . . . + N n − 1 , p ( u ‾ m − 1 ) R m − 1 ] R= \left[ \begin{matrix} 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{matrix} \right] R= N1,p(u1)R1+...+N1,p(u1)Rm−1Nn−1,p(um−1)R1+...+Nn−1,p(um−1)Rm−1...