曲线散乱数据的有序化处理
大多数情况下,三维测量得到的数据点是散乱点,在三维空间上不具备邻接关系,这样未经处理的散乱数据很难直接用于数据参数化和曲线拟合,因此必须对他们进行有序化的预处理,使其按照既定的空间位置顺序排列.
假设散乱曲线点所在曲线至少是 G 1 G^1 G1连续的,则可给出如下排序策略:
- 选择任一点 p 0 p_0 p0作为开始点,将其连接到它的最近点 p s p_s ps,将 p 0 p s p_0p_s p0ps作为搜索方向;
- 更新开始点为 p s p_s ps, 继续向前搜索最近点 p w p_w pw,并计算 p s p w p_sp_w pspw与 p 0 p s p_0p_s p0ps之间的夹角
θ s = a r c c o s ( p 0 p s ∣ ∣ p 0 p s ∣ ∣ ⋅ p s p w ∣ ∣ p s p w ∣ ∣ ) \theta_s=arccos(\frac{p_0p_s}{||p_0p_s||}\cdot\frac{p_sp_w}{||p_sp_w||}) θs=arccos(∣∣p0ps∣∣p0ps⋅∣∣pspw∣∣pspw)
如果 θ s < θ t h , θ \theta_s<\theta_{th},\theta θs<θth,θ为给定阈值,一般可选择为 60 ° − 90 ° 60°-90° 60°−90°,那么 p s p_s ps将被作为新的开始点,同时将搜索矢量更新为 p s p w p_sp_w pspw,继续上述过程,直至到达一点 p e 1 p_e1 pe1,其最近点不再满足 θ s < θ t h \theta_s<\theta_{th} θs<θth,则认为已搜索到曲线端点;
- 将开始点重新设置为 p 0 p_0 p0,搜索方向设为 p 0 p s p_0p_s p0ps的反方向 p s p 0 p_sp_0 psp0,继续上述搜索过程,直至另一个端点也被找到.
将数据点按新的顺序存储,得到序化数据.
曲线有序数据点参数化
给定一组有序数据点 p i ( i = 0 , 1 , . . . , m ) p_i(i=0,1,...,m) pi(i=0,1,...,m), 要确定一条插值或拟合这些点的B样条曲线,首先要为这些数据点确定一个参数值 u i ‾ \overline{u_i} ui.通常有三种参数化方法.
- 均匀参数化
u ‾ 0 = 0 , u ‾ m = 1 , u ‾ i = u ‾ i − 1 + 1 m , i = 1 , . . . , m − 1 \overline{u}_0=0,\overline{u}_m=1,\overline{u}_i=\overline{u}_{i-1}+\frac{1}{m},i=1,...,m-1 u0=0,um=1,ui=ui−1+m1,i=1,...,m−1
仅适用于数据点分布较均匀的情况,一旦分布不均,生成的曲线会产生扭曲变形甚至出现尖点或自交现象.
- 累加弦长参数化
u ‾ 0 = 0 , u ‾ m = 1 , u ‾ i = u ‾ i − 1 + ∣ ∣ p i − p i − 1 ∣ ∣ / ∑ j = 1 m ∣ ∣ p j − p j − 1 ∣ ∣ \overline{u}_0=0,\overline{u}_m=1,\overline{u}_i=\overline{u}_{i-1}+||p_i-p_{i-1}||/\sum_{j=1}^{m}||p_j-p_{j-1}|| u0=0,um=1,ui=ui−1+∣∣pi−pi−1∣∣/j=1∑m∣∣pj−pj−1∣∣
目前最常用的参数化方法,反映了数据点按弦长分布的情况.
- 向心参数化
u ‾ 0 = 0 , u ‾ m = 1 , u ‾ i = u ‾ i − 1 + ∣ ∣ p i − p i − 1 ∣ ∣ / ∑ j = 1 m ∣ ∣ p j − p j − 1 ∣ ∣ i = 1 , . . . , m − 1 \overline{u}_0=0,\overline{u}_m=1,\overline{u}_i=\overline{u}_{i-1}+\sqrt{||p_i-p_{i-1}||}/\sqrt{\sum_{j=1}^{m}||p_j-p_{j-1}||}\\ i=1,...,m-1 u0=0,um=1,ui=ui−1+∣∣pi−pi−1∣∣/j=1∑m∣∣pj−pj−1∣∣i=1,...,m−1
由波音公司Lee提出的方法,当数据点急剧转弯变化时,该方法能得到比弦长参数化更好的结果[1].
function T = CumuPara(P)
%累加弦长参数化
%P输入数据点,T累加弦长参数化得到的参数
m=size(P,1);%数据点个数
T=zeros(1,m);
sum_chord=0;
for j=1:m-1
sum_chord=sum_chord + norm(P(j+1,:)-P(j,:),2);
end
chord=0;
for i=2:m
chord=chord+norm(P(i,:)-P(i-1,:),2);
T(i)=chord/sum_chord;
end
% T(m)=1;
end
曲面阵列数据点参数化
每一行或每一列含有相同数目的数据点称为规则的曲面数据点
P
i
,
j
(
i
=
0
,
.
.
.
,
m
;
j
=
0
,
1
,
.
.
.
,
n
)
{P_{i,j}}(i=0,...,m;j=0,1,...,n)
Pi,j(i=0,...,m;j=0,1,...,n). 对于沿
u
u
u向的第
j
j
j行数据点
P
i
,
j
(
i
=
0
,
.
.
.
,
m
)
{P_{i,j}}(i=0,...,m)
Pi,j(i=0,...,m)可利用上面曲线参数化方法进行参数化,设对应参数为
u
i
,
j
(
i
=
0
,
.
.
.
,
m
)
u_{i,j}(i=0,...,m)
ui,j(i=0,...,m),则曲面数据点
u
u
u向的参数化可取所有行数据点的参数值的算术平均值
u
‾
i
=
1
n
+
1
∑
j
=
0
n
u
i
,
j
,
i
=
0
,
.
.
.
,
m
\overline{u}_i=\frac{1}{n+1}\sum_{j=0}^n u_{i,j},i=0,...,m
ui=n+11j=0∑nui,j,i=0,...,m
类似地,
v
v
v向数据点参数化:
v
‾
j
=
1
m
+
1
∑
i
=
0
m
v
i
,
j
,
j
=
0
,
1
,
.
.
.
,
n
\overline{v}_j=\frac{1}{m+1}\sum_{i=0}^{m} v_{i,j},j=0,1,...,n
vj=m+11i=0∑mvi,j,j=0,1,...,n
v
i
,
j
v_{i,j}
vi,j为
v
v
v向第
i
(
i
=
0
,
.
.
.
,
m
)
i(i=0,...,m)
i(i=0,...,m)列数据点
P
i
,
j
P_{i,j}
Pi,j经曲线点参数化方法得到的参数值.
曲面散乱数据点的参数化
参数化基面的构造
散乱数据并不能直接进行参数化,而是通过事先构造基面的方法. 基面参数化是依据数据点和基面上点的对应关系确定每个数据点参数值的过程[2~6]. 最理想的基面是与实际曲面最相近的曲面,但实际处理中,要根据数据点分布情况确定切合实际的基面.
- Shepard插值曲面
- 双线性Coons曲面
由于圆柱刀机械加工过程中关注的是刀具轴线形成的轴迹面对设计曲面等距面的逼近程度,因此这里重点考虑刀具轴迹面作为基面的构造方法,利用刀轴确定的一族直线段插值形成一张B样条曲面表达的直纹面(详见蒙面法).
散乱数据点的直接基面参数化
基面构造完成后,就可以把数据点投影到基面上,以对应投影点的 ( u , v ) (u,v) (u,v)参数作为该点参数值,完成散乱数据点的参数化. 最常用的投影方式是计算散乱数据点到基面的最近点.
数据点
P
i
{P_i}
Pi与基面
r
(
u
,
v
)
r(u,v)
r(u,v)间的向量可表示为
d
(
u
,
v
)
=
P
i
−
r
(
u
,
v
)
d(u,v)=P_i-r(u,v)
d(u,v)=Pi−r(u,v)
如果直线向量
d
(
u
,
v
)
d(u,v)
d(u,v)与基面上一点
q
i
=
r
(
u
i
,
v
i
)
q_i=r(u_i,v_i)
qi=r(ui,vi)处切平面垂直,则其必与曲面在该点处偏导矢
r
u
,
r
v
r_u,r_v
ru,rv垂直,即满足
{ f ( u , v ) = r u ⋅ [ P i − r ( u , v ) ] = 0 g ( u , v ) = r v ⋅ [ P i − r ( u , v ) ] = 0 \begin{cases} f(u,v)=r_u\cdot[P_i-r(u,v)]=0 \\ g(u,v)=r_v\cdot[P_i-r(u,v)]=0 \end{cases} {f(u,v)=ru⋅[Pi−r(u,v)]=0g(u,v)=rv⋅[Pi−r(u,v)]=0
为找到点 P i P_i Pi在基面 r ( u , v ) r(u,v) r(u,v)上最近点 q i q_i qi, 可采用 N e w t o n Newton Newton迭代法进行求解.
- 估算曲面上 n × n n\times n n×n个采样点,从中找到距离 P i P_i Pi最近的采样点,将其 ( u , v ) (u,v) (u,v)参数作为投影点迭代计算的初始值.
- 在第 i i i步 N e w t o n Newton Newton迭代中,需求解如下线性方程系统
[ f u ( u i . v i ) f v ( u i , v i ) g u ( u i , v i ) g v ( u i , v i ) ] = [ σ u σ v ] [ f ( u i , v i ) g ( u i , v i ) ] \begin{bmatrix} &f_u(u_i.v_i) & f_v(u_i,v_i)\\ &g_u(u_i,v_i) &g_v(u_i,v_i) \end{bmatrix}= \begin{bmatrix} \sigma u \\ \sigma v \end{bmatrix} \begin{bmatrix} f(u_i,v_i)\\ g(u_i,v_i) \end{bmatrix} [fu(ui.vi)gu(ui,vi)fv(ui,vi)gv(ui,vi)]=[σuσv][f(ui,vi)g(ui,vi)]
则第
i
+
1
i+1
i+1次迭代的起始点
(
u
i
+
1
,
v
i
+
1
)
(u_{i+1},v_{i+1})
(ui+1,vi+1)为
{
u
i
+
1
=
u
i
+
σ
u
,
v
i
+
1
=
v
i
+
σ
v
.
\begin{cases} u_{i+1}=u_i+\sigma u,\\ v_{i+1}=v_i+\sigma v. \end{cases}
{ui+1=ui+σu,vi+1=vi+σv.
Piegl等利用两个容差判断迭代的收敛性[7]:
ε
1
\varepsilon_1
ε1度量欧几里得距离是否为0,
ε
2
\varepsilon_2
ε2度量余弦是否为0,收敛准则如下:
- 点是否重合:
∣ ∣ P i − r ( u i , v i ) ≤ ε 1 ||P_i-r(u_i,v_i)\leq \varepsilon_1 ∣∣Pi−r(ui,vi)≤ε1
- 余弦是否重合:
∣ ∣ r u ⋅ [ P i − r ( u i , v i ] ∣ ∣ ∣ ∣ r u ∣ ∣ ∣ ∣ P i − r ( u i , v i ) ∣ ∣ ≤ ε 2 \frac{||r_u\cdot[P_i-r(u_i,v_i]||}{||r_u|| ||P_i-r(u_i,v_i)||}\leq \varepsilon_2 ∣∣ru∣∣∣∣Pi−r(ui,vi)∣∣∣∣ru⋅[Pi−r(ui,vi]∣∣≤ε2
- 参数是否在定义域内
a ≤ u i ≤ b , c ≤ v i ≤ d a\leq u_i \leq b,\quad c\leq v_i \leq d a≤ui≤b,c≤vi≤d
- 参数是否不再显著改变
∣ ∣ ( u i + 1 − u i ) r u ( u i , v i ) + ( v i + 1 − v i ) r v ( u i , v i ) ∣ ∣ ≤ ε 1 ||(u_{i+1}-u_i)r_u(u_i,v_i)+(v_{i+1}-v_i)r_v(u_i,v_i)||\leq \varepsilon_1 ∣∣(ui+1−ui)ru(ui,vi)+(vi+1−vi)rv(ui,vi)∣∣≤ε1
如果条件1,2或4满足,则迭代停止,将投影点参数值 ( u i , v i ) (u_i,v_i) (ui,vi)作为数据点的参数值,遍历所有散乱数据点后,完成参数化过程. 也可以采用其他关于曲面上最近点的计算方法.
Reference
孙玉文,徐金亭,任斐,郭强著. 复杂曲面高性能加工技术与方法[M]. 北京:科学出版社, 2014:45-51.