NURBS曲面拟合

NURBS曲面拟合（共4篇）

NURBS曲面拟合 篇1

0引言

在逆向工程中,自由曲面重构是零件三维CAD模型重建中的核心技术环节。目前相对具有代表性的自由曲面重构技术,主要是以三角网格面为基础的近似曲面重构和以样条曲面为基础的自由曲面重构两种方案。三角网格面的研究相对较早,其技术和方法已经相对比较成熟,并取得了一些比较成熟的研究成果。但是三角网格面在光顺性以及与CAD/CAM系统在模型表示方法的统一上存在缺陷。

由于非均匀有理B样条(NURBS)曲面本身所具有的良好光顺性、局部可编辑和标准的数学表达形式等优点,引起了在逆向工程领域研究中学者们的重视[1]。尤其是近几年来,研究人员在采用NURBS进行曲面重构方面进行了许多研究和实验,也取得了大量的NURBS曲面拟合方法和结果。而在对各种方法的优劣和拟合结果的评价过程中,如果仅靠视觉直观等简单的手段显然是不科学的,必须给出一种定量的计算方法对NURBS曲面拟合的精度作出客观和准确的评测。

本文在采用NURBS对散乱点云进行自由曲面重构的基础上,设计了一种区域划分搜索迭代的拟合精度评测算法,通过此算法可以较快速地求解出每一个原始数据点到NURBS曲面的误差,并以最终的统计结果反映拟合的精度。该方法有效地避免了求解高次方程计算费时,而且容易导致计算发散等问题。

1NURBS曲面拟合过程

1.1B样条及NURBS的基本概念[1,2]

B样条方法是继Bezier方法之后,由美国通用汽车公司的Gordon等人出于外形设计的需求,于1974年提出。该方法克服了Bezier曲线曲面不能作局部修改、曲线阶次受控制顶点个数决定等缺点,同时继承了Bezier方法的一切优点。B样条曲面的定义如式(1)所示:

undefined (1)

其中,顶点Pi,j为控制顶点,k、l为B样条曲线的次数。Bi,k(u)和Bj,l(v)分别为k和l次B样条基函数。NURBS在B样条基础上加入了权因子,使其不仅能表示自由曲线曲面形状,而且能精确表示圆锥曲线、二次曲面等规则曲线曲面。同时NURBS设置节点矢量两端重复度为曲线次数加1,解决了型值点分布不均匀时,用均匀B样条曲线插值求得的曲线两端产生拐点的问题。

1.2基于NURBS的散乱点云曲面重构

本文针对散乱点云数据,在上述基本定义的基础上,采用NURBS曲线插值的方式进行自由曲面重构。该过程主要有两个阶段:第一阶段,通过对型值点的一次插值计算得到单条NURBS曲线,由交互方式经过多次选取数据点和插值计算,得到多条NURBS曲线。第二阶段,由选定的多条NURBS曲线进一步拟合生成NURBS曲面[3,4],处理流程如图1所示。

本文方法采用了修正弦长参数化的方法计算并确定节点矢量,同时根据已有型值点和计算得到的对应节点矢量,代入NURBS曲线的定义公式,通过求解方程组得到控制顶点。

在统一节点矢量的过程中,本文方法采用基于二分插值的节点统一方案,尽可能地保留了曲线的细节部分[5]。在样条次数为3次的条件下,本文通过实验得到了如图2和图3的拟合结果,其中(a)、(b)分别为拟合生成的NURBS曲面网格和打开光照的效果。

2拟合曲面的精度评测

2.1模型精度评价

逆向工程中,精度反映着逆向模型和实物以及模型和产品差距的大小,其评价指标分为量化指标和非量化指标[6]。非量化指标主要用于曲面模型表面光顺性等结果的评价,通过光照效果、法矢和主曲率,并结合人的感官进行评价。量化指标则以定量和准确的数值计算结果来反映模型与实物的实际偏差程度,具体可通过实物上的原始点到模型表面的偏差距离来表示。曲面拟合的精度对模型重构过程中的误差控制和曲面重构方法的选择与改进有着重要指导意义,因此有必要研究出一种实用、有效的精度评测方法。

2.2传统计算方法

对于量化指标的计算,传统的方法是通过求解非线性方程组获得实物点到拟合结果曲面间的距离。该方法的基本原理是曲面外一点Q和曲面上最近的点P(u,v)的距离矢量undefined与曲面在P(u,v)的法矢方向相同,如图4所示。这样将点到曲面的最短距离的计算转换为计算点在参数曲面r(u,v)上的投影,投影方向为曲面的法矢。用方程表示如(2)式。

undefined (2)

其中R为点Q到曲面r(u,v)的最短距离,undefined为参数曲面的偏导数。对于复杂曲面,式(2)所示的方程是一个高次非线性方程,通常采用Newton-Rephson迭代法进行求解。该方法的缺陷是对初值要求严格,选取不当会导致计算发散,同时Newton-Rephson法比较费时,在实践中不是理想的方法。

在实际处理过程中,为了避免求解非线性方程组,常采用几何分割法将曲面离散成一系列平面片组成的多面体,然后计算点到各个平面片的距离,取其中最小值为要求的最短距离。该方法的计算精度与曲面片的离散程度成正比,但如果要得到较高的精度,就必须提高曲面片的离散精度,但计算速度也随之下降。虽然到目前为止,在计算点到Bezier曲面的最短距离方面有一些研究成果,但针对NURBS曲面的拟合精度方面仍然缺乏十分行之有效的计算方法[7]。

2.3基于区域划分的搜索迭代算法

上述通过求解方程组的传统方法中的关键处理环节是求得点在曲面上的空间距离最近的点。从整张曲面上搜索最近的点可以避免直接求解方程组,但又会带来计算时间过长的弊端。本文所研究方法是根据NURBS曲面网格和曲面节点矢量的对应关系,将空间点对应到曲面上的每一个曲面网格,并在曲面网格内按设定的步长迭代搜索离空间点最近的点,这种方法避免了对最近点的基于全局策略的盲目搜索。最后统计出所有空间点到NURBS曲面的最长、最短和平均距离,完成对NURBS曲面拟合精度的量化评测。

在上述算法基本思想下,基于区域划分的搜索迭代算法步骤如下:

Step1 从原始点云数据中读入所有点到链表Points中,并给每一个节点添加三个属性,横向节点矢量编号uNum,纵向节点矢量编号vNum,到曲面的最短距离minDis。

Step2 假设生成的NURBS曲面为S(u,v),横向和纵向节点向量数组分别为U[uControlNum+k+1],V[vControlNum+l+1],其中uControlNum,vControlNum分别为S在横向和纵向的控制顶点个数,k、l为NURBS曲线的次数。实验中采用3次NURBS曲线,故这里k=l=3。

Step3 从Points中依次取出一个点P,初始化其P.uNum=uControlNum,P.vNum=vControlNum。设置循环变量i=k～uControlNum,j=1～vControlNum。不从0开始的目的是为去除节点向量数组中的端点重复值。

Step4 循环比较P的横向坐标值(若横向为x方向,则为P.x;若为y方向,则为P.y)和S(U[i],0)的大小,若小于S(U[i],0),则使P.uNum=i,结束本次循环。

Step5 循环比较P的纵向坐标值和S(U[P.uNum],V[j])的大小,若大于S(U[P.uNum],V[j],则使P.vNum=j,结束本次循环。跳转到Step3继续执行,直至Points中所有的点都完成区域划分。

Step6 从完成区域划分的Points中依次取出一个点P,在u=U[P.uNum-1]～U[P.uNum],v=V[P.vNum-1]～V[P.vNum]区域内迭代搜索曲面S离P最近的点,搜索的步长取du=(U[P.uNum]-U[P.uNum-1])/M,dv=(V[P.vNum]-V[P.vNum-1])/N。显然,搜索的次数为M×N。

Step7 完成点P的搜索后,记下P到曲面S最近点的距离作为P到S的最短距离,保存到P的属性minDis。跳转到Step6继续执行,直至所有的点P都完成迭代搜索。

Step8 最后计算得到所有点中到曲面距离的最大值、最小值和平均值。

该算法处理步骤中,在对点云数据中每个点的区域划分过程中,时间复杂度为O(uControlNum×vControlNum),且容易分析出是由NURBS曲面的控制顶点个数决定。在迭代搜索每个点到曲面上最近点的过程中,时间复杂度是O(M×N)。搜索的精确程度与M和N的取值大小有关,M、N取值越大,搜索越精细,搜索得到与最近点距离越能代表空点和NURBS曲面的最短距离。但随着M、N取值的增大,算法的计算时间也会随抛物线级增长,影响算法执行的效率。故需要在尽可能精确的条件下给M、N设定一个合适的值。

3实验结果

作者以本文第2节中所述及的“管状”零件模型和人脸模型为例,采用本文算法对其经重构生成的NURBS曲面进行拟合精度计算。实验的硬件环境为P4 3.0GHZ CPU,512MB 内存,编程环境为Visual C++ 6.0,图形开发环境为OpenGL。

经过多次测试和试验,在实验中取M和N的值为10,表1给出了采用本文算法计算得到的NURBS曲面拟合结果的量化精度计算结果。为了使数值结果具有可参照性,表1中同时给出了模型平均宽度和高度。

为了验证表1中评测结果是否真实反映NURBS曲面拟合精度,作者采用本文NURBS方法对标准球体模型进行曲面拟合,对拟合得到的曲面采样获取一定数量曲面上的点,计算这些点到球心的距离与半径之间的偏差,计算结果如表2所示。表2中各偏差相对于半径的百分比与表1中各偏差分量相对于模型宽度和高度的百分比接近,从而间接证明表1中评测结果是NURBS曲面拟合偏差的真实体现,也证明了本文算法的正确性。

4结语

本文针对散乱点云数据采用NURBS进行了自由曲面重构,得到了良好的拟合效果。在此基础上设计了一种基于区域划分的搜索迭代算法,该算法为NURBS曲面拟合精度分析过程提供了一种量化的计算方法。实验结果验证了该算法在时间能够被接受的前提下得到了当前曲面与原始模型的偏差,同时也证明了该算法的可行性和正确性。

参考文献

[1]张曦煌,杜俊利.计算机图形学[M].北京:邮电大学出版社,2006.

[2]Donald Hearn,M Pauline Baker.计算机图形学[M].蔡士杰,等译.北京:电子工业出版社,2005.

[3]慈瑞梅,李东波.逆向工程中NURBS曲面重构技术研究[J].南京理工大学学报,2004,28(4):391-394.

[4]王田苗,曹宇男,陈友东,等.基于de Boor算法的NURBS曲线插补和自适应速度控制研究[J].中国机械工程,2007,18(21):2608-2613.

[5]林致浩,姜晓峰.逆向工程中曲面重构技术的研究与实现[D].苏州大学,2006.

[6]金涛,童水光,等.逆向工程技术[M].北京:机械工业出版社,2003.

[7]Hu SM,Wallner J.A second order algorithm for orthogonal projectiononto curves and surfaces[J].Computer Aided Geometric Design,2005,22(3):251-260.

NURBS曲面拟合 篇2

1 曲线节点插入技术在曲面中的推广

1.1 NURBS曲面

假设一张次NURBS曲面如下形式

$\overline{r}(u,v)=\frac{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\omega }}{}_{ij}d{}_{ij}{\rm N}{}_{ik}(\mathit{u}){\rm N}{}_{jh}(\mathit{v})}{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\omega }}{}_{ij}{\rm N}{}_{ik}(\mathit{u}){\rm N}{}_{jh}(\mathit{v})}(1)$

式中:dij是控制网格(m+1列n+1行)中第i列第j行的点。ωij是与顶点对应的权因子,Nik与Njh分别是u向k次和v向次的规范B样条基,分别由节点矢量u,v按Cox-DeBoor递推公式决定。

1.2 曲面上的节点插入技术

假设NURBS曲面的向节点矢量,u规范化后为$\mathit{u}=[\underset{k+1}{{\displaystyle 0,0,\cdots ,0}},u{}_{k+1},\cdots ,\underset{r}{{\displaystyle u{}_{i-r+1}=\cdots =u{}_i}},\cdots ,u{}_m,\underset{k+1}{{\displaystyle 1,1,\cdots ,1}}]$,其中节点ui(k<i<m+1)的重复度为r,将ui重复插入l次后,生成一组新的节点矢量,重新编号为$\mathit{u}{}^{\mathit{l}}=[\underset{k+1}{{\displaystyle 0,0,\cdots ,0}},u{}_{k+1}^l,\cdots ,\underset{r}{{\displaystyle u{}_{i-r+1}^l=\cdots =u{}_i^l}},\cdots ,u{}_m^l,\underset{k+1}{{\displaystyle 1,1,\cdots ,1}}]$。

假定u为网格横向,v为纵向,则这个新的节点矢量ul决定了一组新的B样条基N${}_{i^{\prime },k}^l$(u)(i′=0,1,…,m+l),则原来的NURBS曲面可以表示为

$\overline{\mathit{r}}(\mathit{u},\mathit{v})$$=\frac{{\displaystyle \sum _{i^{\prime }=1}^{m+l}{\displaystyle \sum _{j=1}^n\omega }}{}_{i^{\prime },j}^{l,0}d{}_{i^{\prime },j}^{l,0}{\rm N}{}_{i^{\prime },k}^l(\mathit{u}){\rm N}{}_{j,h}(\mathit{v})}{{\displaystyle \sum _{i^{\prime }=1}^{m+l}{\displaystyle \sum _{j=1}^n\omega }}{}_{i^{\prime },j}^{l,0}{\rm N}{}_{i^{\prime },k}^l(\mathit{u}){\rm N}{}_{j,h}(\mathit{v})}(2)$

以控制网格中的第j行为例,新控制顶点d${}_{i^{\prime },j}^{l,0}$和权值ω${}_{i^{\prime }j}^{l,0}$可由老控制顶点和权值按下面的公式计算

$\omega {}_{i^{\prime },j}^{s,0}=\{\begin{array}{l}\omega {}^{s-1,0}0\le i^{\prime }\le i-k+s-1\\ \alpha {}_{i^{\prime },j}^{s,0}\omega {}_{i^{\prime },j}^{s-1,0}+(1+\alpha {}_{i^{\prime },j}^{s,0})\omega {}_{i^{\prime }-1,j}^{s-1,0}i-k+s\le i^{\prime }\le i-r\\ \omega {}_{i^{\prime },j}^{s-1,0}i-r+1\le i^{\prime }\le n+s\end{array}(3)$

$d{}_{i^{\prime },j}^{s,0}=\{\begin{array}{l}d{}_{i^{\prime },j}^{s-1,0}0\le i^{\prime }\le i-k+s-1\\ \frac{\alpha {}_{i^{\prime },j}^{s,0}\omega {}_{i^{\prime },j}^{s-1,0}d{}_{i^{\prime },j}^{s-1,0}+(1-\alpha {}_{i^{\prime },j}^{s,0})\omega {}_{i^{\prime }-1,j}^{s-1,0}d{}_{i^{\prime }-1,j}^{s-1,0}}{\omega {}_{i^{\prime },j}^{s,0}}i-k+s\le i^{\prime }\le i-r\\ d{}_{i^{\prime },j}^{s-1,0}i-r+1\le i^{\prime }\le n+s\end{array}(4)$

其中s=1,2,…,l,ω${}_{i^{\prime },j}^{0,0}$=ωi′,j,d${}_{i^{\prime },j}^{0,0}$=di′,j

以上对控制顶网格中的第j行沿u向l级递推,由图1所示,将虚线框中k-1个老顶点用实框中新顶点所替换,其余的老控制顶点不变。权值的处理方法同理可得。

当l+r=k时,实框中最后一行的控制顶点就合成一个点d${}_{i-r,j}^{l,0}$。

2 曲面的NURBS向分段有理Bezier表示形式的转化

2.1 B样条曲面转化为Bezier曲面的经典算法

对于k次B样条曲线的节点矢量满足端点重复度为k+1,内部节点的重复度为k时,则B样条曲线就是一分段有理Bezier曲线,重复节点之间的一段曲线就是一个有理Bezier曲线。基于此对于曲面的情况Boehm提出过这样的思路:对控制网格(k×k)逐行利用曲面节点插入技术直到每一行对应的u向内部节点(uk+1…,…um)的重复度全部达到k,则将生成一张新的控制网格定义其为中间网格。再对这中间网格逐列执行插入算法直到每一列对应的v向内节点全达到h重,则整个中间网格最终可转换成Bezier网格。经过研究,文中提出了一种改进的优化算法,可简化整个转换过程。

2.2 电磁散射计算中曲面的NURBS到Bezier转换的优化算法

在电大尺寸的RCS计算中,要求在NURBS曲面转化过程简单易于实现计算量尽量的小,并且转化后的Bezier曲面片最好都能够拥有单独的控制网格信息和规范的数据存储方式,以便于后续的计算。上述“逐行”,再“逐列”的经典算法直接从控制网格入手,不但过程繁琐,计算量大,节点区间存在着重复处理而且对结果结果还需要进一步的处理来寻找相应的控制网格信息。

由式(3),式(4)知NURBS控制网格的u向插值运算中的递推过程在不同的控制行里的步骤相同,即递推涉及的控制顶点及权因子在每一行中对应相同的序号,且由式(5)知α${}_{i^{\prime },0}^{s,0}$=α${}_{i^{\prime },1}^{s,0}$=…=α${}_{i^{\prime },j}^{s,0}$=…=α${}_{i^{\prime },n}^{s,0}$。基于这一特点处理NURBS曲面时如果从节点矢量入手,先根据内部非空节点区间的情况将控制网格剖分,这样可在转化过程中实现分片处理的方法,同时也可使结果直接反映出曲面片的个数,有利于后面的电磁计算的展开。这样的方法也可以形象的比喻为切蛋糕的过程。

先对u向寻找非空节点区间段,最多可找出m-k个区间,根据找到区间将控制网格沿u向划分成对应的网格片。由于在“片”内每一行的曲面节点插入过程均一样,因此分片处理可以简并相同递推量的计算,将整“片”对应的节点矢量区间两端节点插值k为次重复度后即完成了对一个“片”的转化。同理依次从左到右直到所有的“片”都被处理完,如图2所示,则生成过渡网格。再对整个过渡网格进行v向分“片”,同理处理所有的片,最后将生成一系列Bezier控制网格。这样一种处理方式下使插入算法原来要处理的(m-k)(n+1)+(n-h)(m+1)个节点区间的情况缩减为(m-k)+(n-h)个,减少了节点插入的次数,避免重复的递推过程,另外在对行执行节点插入计算时,进一步再将B样条曲线递推算法中不同级之间递推系数的关系推广到曲面上后,有α${}_{b-k+i,j}^{1,0}$=α${}_{b-k+i+1,j}^{2,0}$=α${}_{b-r,j}^{k-r-i,0}$,同理对v向,这样一来可使需要计算的$(n+1)\times {\displaystyle \sum _{i=k=1}^m\frac{(k-r{}_i)\times (k-r{}_i+1)}{2}}$个横向递推系数α${}_{i^{\prime },j}^{s,0}$和$(m+1)\times {\displaystyle \sum _{j=h=1}^n\frac{(h-r{}_j)\times (h-r{}_j+1)}{2}}$个纵向递推系数α${}_{i,j^{\prime }}^{0,s}$的计算减小到:${\displaystyle \sum _{i=k=1}^m(}k-r{}_i)$和${\displaystyle \sum _{j=h=1}^n(}h-r{}_j)$次。其中ri,rj分别为u向的第i个节点和v向的第j个节点的重复度。

3 算例

以文献[6]中的图形数据为例,利用Matlab,按其NURBS表示生成的图形见图3(a),应用文中的算法将其转化为分段有理Bezier形式后生成的图形,如图3(b),分片的图3(c)所示,与原NURBS曲面相同。但曲面拼接处的光滑性会降低,对于RCS计算而言,对目标的光滑性要求不高,不会影响计算结果。

4 结束语

文中提出的算法在曲面次数高、形状复杂、控制顶点数大、面片数量多的情况下优势尤为明显。由于Bezier在数值计算方面的稳定性和简洁性,因此采用NURBS技术来进行几何建模,再用文中的算法将其转化为分段Bezier表示形式来进行数值计算。在不损失模型精度的同时,提高了整个工程计算的速度和可靠性。另外,文中对控制顶点“分块”处理的方法可应用于NURBS曲面的生成算法中,提高效率。

参考文献

[1]施法中.计算机辅助几何设计与非均匀有理B样条[M].北京:高等教育出版社,2001.

[2]Perez J,Catedra M F.Application of Physics Optics to the RCS Computation of Bodies Modeled with NURBS Surfaces[J].IEEE Trans on Antennas Propagation,1994,42(10):1404-1411.

[3]Boehm W.Inserting New Knots Into B-spline Curves[J].Computer Aided Design,1980,12(4):199-201.

[4]Boehm W.Generating the Bezier Points of B-spline Curves and Surfaces[J].Computer Aided Design,1981,13(6):365-366.

[5]王鸿.NURBS曲线转化为Bezier表示形式的快速算法[J].西安公路交通大学学报,2000,20(3):120-122.

NURBS曲面拟合 篇3

1 点云数据的采集与预处理

通常把三维空间内点的集合称为点云, 点云的数据个数从几百、几千到几万个不等, 排列方式也有散乱点云、扫描线点云、网格化点云等。不同数据采集装备所产生的点云是不同的, 所采用的处理方法也是不同的。

1.1 数据采集与精简

反求工程采用的测量方法主要有两种:一是传统的接触测量法, 二是无接触测量法[2]。本文的数据采集装置是激光扫描仪, 通过规划路径, 可以获得一系列点的Z坐标值, 进而进行后续的点云数据处理与对比分析。采用均匀网格法可以减少数据量, 其原理是:首先把所得的数据点进行均匀网格划分 (通常由用户指定) , 然后从每个网络中提取点根据Z值排序, 如果某个点位于各个点中间, 那么这个点被选中保留, 其他点滤除[3]。

1.2 数据平滑

经过数据精简的点云数据若直接进行曲面重建, 可能出现局部曲率过大的情况, 容易增加后期曲面重建时运算负担, 所以要进行平滑处理。高斯滤波法以高斯滤波器在指定域内的权重为高斯分布 (正态分布) , 其平均效果较小, 在数据平滑的同时, 能较好地保持原数据形貌, 因而常被使用。高斯滤波算法按式 (1) 计算。

式中pi-3~pi+3--连续相邻7个点的坐标值 (mm×mm×mm) 。

1.3 边缘识别与分割

采集装置采集得到的点云数据不仅包括被测零件点云数据, 还包括零件周围台面及台面上其他物体的点云数据。必须进行边缘识别与分割, 将被测零件从台面中区分出来。被测零件与台面间的边界主要属于褶皱边界 (切矢不连续) , 可以通过检查曲率半径的方式来判定是否为边界点。

2基于NURBS技术的曲面重建

NURBS曲面可采用有理分式方法表示为:

这里控制定点di, j (i=0, 1, ……, m;j=0, 1, ……, n) 呈拓扑矩形阵列, 形成一个控制网格。ωi, j是与定点di, j联系的权因子, 规定四角顶点处用正权因子, 即ω0, 0、ωm, 0、ω0, n、ωm, n>0, 其余ωi, j≥0且顺序k×l个权因子不同时为0。Ni, k (u) (i=0, 1, ……, m) 和Nj, l (v) (j=0, 1, ……, n) 分别为u向k次和v向l次的规范B样条基。它们分别由u向和v向德节点矢量U=[u0, u1, ……, um+j+1]与V=[u0, u1, ……un+l+1]按德布尔递推公式决定。控制网格的生成是重建NURBS曲面十分重要的一步。在生成m×n的控制网格时, 由于点云数据在沿x方向上的密度远大于沿y方向上的密度, 因而令n等于扫描路径的数目。沿x方向上控制点数m由用户设定, 按式 (3) 计算跨距S:

在x方向上每隔跨距S取一点的坐标值作为控制点的坐标, 从而生成m×n的均匀矩形控制网格。要确定一个NURBS曲面的形状, 不光要知道控制点信息, 还要知道节点向量和权值。为使NURBS曲面边界与点云数据边界具有更高的重合度, k阶NURBS曲面节点向量的前k+1个参数为0, 最后k+1个参数为1。由于点云数据基本均匀分布, 权值可以都设为1。在确定控制点网格、节点向量和权值后, 就可以通过插值生成NURBS曲面。

3 算法实现

采用上述方法, 对如图1所示的点云进行了曲面重建实验, 边缘分割结果和曲面重建结果分别如图2和图3所示。在数据平滑过程中, 如果一次平滑的效果不理想, 可以进行多次数据平滑。对于一个m×n的控制网格, m由扫描路径条数决定, 在参数输入区输入n的值, 就能根据这个m×n的控制网格分段NURBS曲面, n值的大小决定了重建曲面的精度。从实验结果可以看出, 本文算法成功地实现了数据的预处理、分割以及NURBS曲面重建。

4 结语

在点云数据曲面重建过程中, 原始点云数据不能直接用来重建曲面, 需要通过数据精简、数据平滑、边缘识别与分割才能生成适于重建曲面的点云数据。从该点云数据中选取适当点作为NURBS曲面的控制点, 并赋予适当的权值和节点向量, 就能生成NURBS曲面。

摘要：在反求工程中, 点云数据的曲面重建质量直接关系到结果的精确性和实用性。本文基于N U R B S技术对点云数据的曲面重建过程进行了研究, 建立了点云数据预处理和曲面重建的算法, 通过实验验证了文中算法的正确性。

关键词：曲面重建,NURBS,点云,反求工程

参考文献

[1]WEIYIN MA, PEIREN HE.B-spline surface local updating with unorganized points[J].Computer-Aided Design, 1998, (11) :853-862.

[2]苏海.反求工程在工业设计中的应用研究[D].昆明:昆明理工大学, 2000.

NURBS曲面拟合 篇4

虽然自由曲线曲面直接插补方式已成为数控加工的主流,但以微小直线段离散化曲线加工路径的微段加工方式仍被广泛地应用。该方式借助CAM软件,通过后置处理将连续加工路径离散化为大量的微小直线段,并生成数控程序文件后,数控系统根据由微小直线段组成的折线进行插补和加工。但是,微段加工方式下加工程序量过大和需要频繁加减速来满足加工精度要求饱受诟病。

为了实现微小直线段的高速平滑加工,提高加工效率,国内外学者做了大量的研究工作。徐志明等[1]用递归算法求得微段转接点处的最大允许进给速度;王宇晗等[2]利用微小线段速度衔接数学模型求得进给速度限制的近似最优解;许海峰等[3]通过改进和完善上述模型,给出了离散化下的速度限制值计算方法;彭芳瑜等[4]进一步考虑了机床的动力学约束;黄昕等[5]提出了一种以最大进给速度为目标,通过双向扫描算法,获得拐点处的最优进给速度的算法;叶伟等[6]和李小清等[7]也根据各自的微小段加工速度、加速度衔接模型改进了微段加工方法;冷洪滨等[8]进一步对单轴的运动速度进行限制来改善加工;杨开明等[9]则通过在相邻程序段间加入过渡段,在牺牲局部精度的前提下实现速度的平滑。

但是,以上算法无法满足高速、高精加工的需求。为此,将微小直线段用多项式曲线、三次样条和B样条等进行拟合,然后进行直接插补成为更好的改进方向。Erkorkmaz等[10]和Yau等[11]以弦长为参数,分别采用五次样条曲线和Bézier曲线来拟合微小直线段从而生成新的光滑加工路径来解决微段加工出现的问题;任锟等[12]对离散加工路径采用三次样条曲线拟合,分析三次样条发现高曲率点,预估高曲率点处的最优速度,并通过S形加减速控制实现微段优化加工;张园等[13]采用五次样条曲线拟合微小直线段生成新的光滑加工路径,并附加考虑弓高误差的限定,对整个样条拟合曲线进行速度规划,实现了进给速度的自适应控制。

本文采用基于NURBS曲线拟合的微段高速平滑加工算法(简称微段平滑算法)来实现微段小直线的平滑加工。

1 微段平滑算法流程和自适应分区

微段平滑算法的基本思路是将解析得到的刀位点进行自适应分区,并分区拟合为NURBS曲线,再进行NURBS曲线直接插补。即先进行一定段数的预读,然后根据给定的限制条件将微段自适应分为各个区域,最后对符合拟合条件的区域逐区进行带权因子和一阶导数约束的最小二乘NURBS曲线逼近,并进行NURBS曲线实时插补,与此对应,将不符合拟合条件的区域划分为直线段区域,保持线性插补方式,如图1所示。

1.1 自适应分区的原则

通过前瞻,可以得到关于微段的两个基本信息:相邻微段间的夹角θi(i=0,1,2,…)和微段的长度Li。拟合的时间与被拟合的点数相关,因此,自适应分区受拟合最大刀位点数、微段间夹角和微段长度的限制,三者为自适应分区的限制条件。将临界处的刀位点作为断点,以断点为界将微段划分为不同的区域,即断点前为上一区域(断点为终点),断点后为下一区域(断点为下一区域的起点)。

1.2 拟合最大刀位点数限制

针对微段长度小和数量多的特点,首先对最大刀位点数nlim进行限制,将其作为自适应分区的首要约束条件。nlim的大小与数控系统控制器的运算能力和对实时性的要求(插补周期T)有关。不失一般性,不妨设上一断点为Pr(r=0,1,2,…),当前刀位点为Pr+n。当n≥nlim时,将刀位点Pr+n作为新的断点,断点前的从r到r+n的n+1个刀位点划为一个适合拟合区域。

1.3 微段夹角限制

在确定点数的情况下,微段的几何特性决定了拟合的计算量和拟合曲线的形状,而微段间的夹角关系是一个主要表征。因此,先对微段夹角的大小θi进行限制,并对θi≤θlim(θlim为最大夹角)的微段进一步进行夹角差Δθi和夹角累加值${\displaystyle \sum _{j=0}^i\theta }{}_j$的判定和约束。

根据给定进给速度vf、最大加速度amax和插补周期T估算得到最大夹角θlim,即

θlim=arccos(amaxT/vf) (1)

当θI>θlim(I=2,3,…)时,PI为一尖点,将PI-2到PI+2的4个微小直线段划分为直线段区,PI-2为断点,PI+2为下一区域的起点。当θI≤θlim时,若ΔθI>Δθlim或${\displaystyle \sum _{j=0}^n\theta }{}_j>\theta {}_{\text{Σ}}^{\max }(\text{Δ}\theta {}_{\lim }$为最大夹角差,θmaxΣ为最大夹角累加值),则当前刀位点PI为断点,否则,继续读入下一个刀位点,并重复微段夹角限制流程。

1.4 微段长度限制

微段是由后置处理程序将刀轨按照给定误差离散化得到的,反映的是刀轨的几何特性。当刀轨曲率变化较小时,微段较长;相反,当刀轨曲率变化较大时,微段较短。换言之,若微段较短(接近1~2个插补长度),刀轨曲率变化就很大,加入拟合环节后得到的加工精度较直接线性插补得到的加工精度变化不大,而只会增加计算的时间和复杂度;与此类似,若微段较长(超过100个插补长度)时,刀轨更趋近于直线,对这些微段进行拟合,反而会使拟合后得到的加工路径偏离原微段所描述的刀轨,使加工精度变差。因此,必须对微段长度加以约束,即规定待拟合微段长度必须介于最小微段长度Lmin和最大微段长度Lmax之间,一方面,使拟合得到的加工路径可以准确反映刀轨的几何特性,另一方面,使拟合计算尽可能简化。

不失一般性,不妨设上一断点为Ps(s=0,1,2,…),并进一步假设相邻的两个刀位点为Pi和Pi+1(i≥s),Li为$\overline{{\rm P}{}_i{\rm P}{}_{i+1}}$的长度。若Lmin≤Li≤Lmax,则微段$\overline{{\rm P}{}_i{\rm P}{}_{i+1}}$为适合拟合微段,并读入下一个刀位点Pi+2,继续判定微段长度;否则Pi为断点,并将从Ps到Pi的刀位点划为适合拟合区域,将微段$\overline{{\rm P}{}_i{\rm P}{}_{i+1}}$划为直线段区域,进行直接线性插补。默认取Lmin=vfT,Lmax=100vfT(vf为进给速度,T为插补周期)。

2 带权因子和一阶导数约束的NURBS曲线最小二乘逼近算法

根据后置处理的原理,微小直线段是刀轨按给定误差离散化得到的。因此,微段平滑算法的基础是利用最小二次逼近法,按照给定误差对适合拟合区域的微段进行拟合,近似得到连续的加工路径来逼近原有刀轨。在拟合误差允许范围内,为了进一步提高拟合曲线的平滑度,本文提出的算法在NURBS曲线最小二乘逼近算法的基础上,对各数据点增加了权因子约束和进给速度方向限制(一阶导数约束),使拟合得到曲线更逼近原始加工路径的同时,增加曲线的平滑度。

2.1 算法基本流程

如图2所示,带约束NURBS最小二乘逼近算法基本流程为:首先按照修正弦长法对读入的数据点进行参数化,然后根据各点的权因子和一阶导数约束进行最小二乘拟合,得到最大逼近误差范围内的控制点,生成新的NURBS曲线加工路径。

2.2 算法实现

假设适合拟合的r+1个数据点为qj(j=1,2,…,r),qj受权因子hj(hj≥1)或一阶导数fj(进给速度方向)的约束,拟合得到的NURBS曲线为C(u)(u∈[ui,ui+1]⊂[u3,un+1])[14]。

默认权因子hj=1(若hj通过手工编程指定,则hj为编程值),当数据点qj最小二乘数逼近误差εj超出给定最大逼近误差εmax时,增大该点处的权因子(hj←chhj,ch为权因子修正系数,默认ch=5),然后重新迭代计算,直到满足εj≤εmax。

按数据点是否受fj约束,将数据点分为导数约束点和导数非约束点,分别简称为约束点和非约束点,并假设约束点的个数为s+1(1≤s≤r),非约束点的个数为t+1(-1≤t≤r-s-1),拟合后的NURBS曲线包含m+1(m≥s+1)个控制点。根据实际情况,一般非约束点远多于约束点,现规定:不计首末端点,非约束点比约束点至少多两个。将约束点标识为qcj(k)(k=0,1,…,s),对应的一阶导数标识为fcj(k),对应的权因子为ωcj(k);非约束点标识为qucj(k),对应的一阶导数标识为fucj(k),对应的权因子为ωucj(k)。

进一步令S为一个含有2(s+1)个约束数据项qcj(k)和fcj(k)的列向量,S=(qcj(0),qcj(1),…,qcj(s),fcj(0),fcj(1),…,fcj(s));T为一个含有2(t+1)个非约束数据项qucj(k)和fucj(k)的列向量,T=(qucj(0),qucj(1),…,qcj(t),fucj(0),fucj(1),…,fucj(t));D为一个包含m+1个控制点(dl=(xl,yl,zl),l=0,1,…,m)的矩阵,D=[d0d1 … dm]T;Wc为一个对角线元素为约束点权因子ωcj(k)的r+1维方阵,Wc=diag(ωcj(0),ωcj(1),…,ωcj(s));Wuc为一个对角线元素为非约束点权因子ωucj(k)的r+1维方阵,Wuc=diag(ωucj(0),ωucj(1),…,ωucj(s));Mc为一个包含约束点p次基函数和对应p次fcl(k),p的(r+1)×(m+1)维矩阵;Muc为一个包含非约束点p次基函数和对应p次fucl(k),p的(r+1)×(m+1)维矩阵:

$\begin{array}{l}\mathit{{\rm M}}{}^{\text{c}}=\\ \left[\begin{array}{cccccccc}{\rm N}{}_{0(0),p}& {\rm N}{}_{0(1),p}& \cdots & {\rm N}{}_{0(r),p}& \mathit{f}{}_{0(0),p}^{\text{c}}& \mathit{f}{}_{0(1),p}^{\text{c}}& \cdots & \mathit{f}{}_{0(r),p}^{\text{c}}\\ {\rm N}{}_{1(0),p}& {\rm N}{}_{1(1),p}& \cdots & {\rm N}{}_{1(r),p}& \mathit{f}{}_{1(0),p}^{\text{c}}& \mathit{f}{}_{1(1),p}^{\text{c}}& \cdots & \mathit{f}{}_{1(r),p}^{\text{c}}\\ ⋮& ⋮& & ⋮& ⋮& ⋮& & ⋮\\ {\rm N}{}_{m(0),p}& {\rm N}{}_{m(1),p}& \cdots & {\rm N}{}_{m(r),p}& \mathit{f}{}_{m(0),p}^{\text{c}}& \mathit{f}{}_{m(1),p}^{\text{c}}& \cdots & \mathit{f}{}_{m(r),p}^{\text{c}}\end{array}\right]\\ \mathit{{\rm M}}{}^{\text{u}\text{c}}=\\ \left[\begin{array}{cccccccc}{\rm N}{}_{0(0),p}& {\rm N}{}_{0(1),p}& \cdots & {\rm N}{}_{0(r),p}& \mathit{f}{}_{0(0),p}^{\text{u}\text{c}}& \mathit{f}{}_{0(1),p}^{\text{u}\text{c}}& \cdots & \mathit{f}{}_{0(r),p}^{\text{u}\text{c}}\\ {\rm N}{}_{1(0),p}& {\rm N}{}_{1(1),p}& \cdots & {\rm N}{}_{1(r),p}& \mathit{f}{}_{1(0),p}^{\text{u}\text{c}}& \mathit{f}{}_{1(1),p}^{\text{u}\text{c}}& \cdots & \mathit{f}{}_{1(r),p}^{\text{u}\text{c}}\\ ⋮& ⋮& & ⋮& ⋮& ⋮& & ⋮\\ {\rm N}{}_{m(0),p}& {\rm N}{}_{m(1),p}& \cdots & {\rm N}{}_{m(r),p}& \mathit{f}{}_{m(0),p}^{\text{u}\text{c}}& \mathit{f}{}_{m(1),p}^{\text{u}\text{c}}& \cdots & \mathit{f}{}_{m(r),p}^{\text{u}\text{c}}\end{array}\right]\end{array}$

式中,Nm(k),p为p次基函数。

p次fcl(k)和fucl(k)(统一简化表示为fl(k))为

$\begin{array}{l}\mathit{f}{}_{l,0}(u)=\{\begin{array}{l}\mathit{f}{}_{i(0)}(u)u\in [u{}_l,u{}_{l+p})\\ 0\text{其}\text{他}\end{array}\\ \mathit{f}{}_{l,p}(u)=\frac{u-u{}_l}{u{}_{l+p}-u{}_l}\mathit{f}{}_{l,p-1}(u)+\\ \frac{u{}_{l+p+1}-u}{u{}_{l+p+1}-u{}_{l+1}}\mathit{f}{}_{l+1,p-1}(u)u\in [u{}_l,u{}_{l+p})\end{array}\}(2)$

首末端点的一阶导数f0和fr可以根据给定的首末端点边界条件计算得到,例如q0处的边界条件为切矢条件时,f0为上一段加工程序在q0处的切矢。中间点的一阶导数按照下列规则得到:若未给定,则该数据点为非约束点;若手工编程指定,则fj为编程值;若已给定法向矢量,则fj通过法向矢量近似计算,近似计算方法如下:

假设法向矢量为nj,qjqj+1为点qj到qj+1的矢量,fj可近似计算得

$\begin{array}{l}\mathit{f}{}_j=\text{Δ}\mathit{q}{}_j=\left[\begin{array}{ccc}x{}_{j+1}-x{}_j& y{}_{j+1}-y{}_j& z{}_{j+1}-z{}_j\end{array}\right]\times \\ \left[\begin{array}{ccc}\sin \alpha & -\cos \alpha & 0\\ \cos \alpha & \sin \alpha & 0\\ 0& 0& 1\end{array}\right](3)\end{array}$

$\alpha =\arccos \frac{\mathit{q}{}_j\mathit{q}{}_{j+1}\cdot \mathit{n}{}_j}{|\mathit{q}{}_j\mathit{q}{}_{j+1}||\mathit{n}{}_j|}$

式中,(xj,yj,zj)、(xj+1,yj+1,zj+1)分别为点qj和qj+1的坐标;α为qjqj+1与nj的夹角。

通过引入拉格朗日乘数法,即引入2(s+1)个拉格朗日乘数因子,带约束的NURBS曲线最小二乘逼近的算法实现问题衍化为解一个系数矩阵为2s+m+3维方阵的线性方程组。设拉格朗日乘数因子为λa(a=0,1,…,2s+m+3),令A=[λ0λ1 … λ2s+m+3],可得非约束点满足的方程:

MucD=S (4)

约束点满足的方程:

McD=T (5)

非约束点的误差可以表示为S-MucD,按照最小二乘法逼近原理,在约束条件McD=T下,使误差S-MucD的加权平方达到最小。

因此,根据拉格朗日乘数法,应使下式达到最小值:

Y=(S-MucD)TWuc(S-MucD)+

AT(McD-T)TWc(McD-T) (6)

对Y求A和D求的偏导数,并令其等于0,得

$\begin{array}{l}-2(\mathit{S}{}^{\text{Τ}}-\mathit{{\rm M}}{}^{\text{u}\text{c}}\mathit{D})\mathit{W}{}^{\text{u}\text{c}}\mathit{{\rm M}}{}^{\text{u}\text{c}}+\mathit{A}{}^{\text{Τ}}\mathit{{\rm M}}{}^{\text{c}}\mathit{W}{}^{\text{c}}=0\\ \mathit{{\rm M}}{}^{\text{c}}\mathit{W}{}^{\text{c}}\mathit{D}=\mathit{{\rm T}}\end{array}\}(7)$

解式(7),并用A/2代替A,得

T=McWc[(Muc)TWucMuc]-1(MucWucS-McWcA) (8)

进一步可推得

A=(McWc)-1×

[MucWucS-(Muc)TWucMuc(McWc)-1T] (9)

通过式(9)可以得到拉格朗日乘数因子矩阵A,将A代入式(8),即可解得所需的控制点矩阵D,得到所需的p次NURBS曲线。

3 仿真与实验

为了验证算法的正确性和有效性,对图3所示的原始图形(三瓣花)的连续加工路径进行后置处理,借助CAM软件生成微小直线段,微段平均长度为1.0923mm。按照上述算法分别进行仿真与实验。

拟合限制参数如下:最大数据点数nlim=30,微段最小长度Lmin=2μm,微段最大长度Lmax=20mm,最大夹角θlim=π/9,最大夹角差Δθlim=π/18,最大夹角累加值θmaxΣ=π/2,给定逼近内外误差e∈[-2μm,2μm]。加工参数如下:进给速度vf=300mm/s,最大加速度amax=100m/s2,插补周期T=2ms,允许最大弓高误差E=2μm。分别利用微段加工方式和本文算法对三瓣花微段加工路径进行测试和实验,并借助MATLAB对加工过程进行仿真,可得三瓣花微段加工路径NURBS曲线逼近的仿真结果(图4a),以及两种加工方式下进给速度对比图(图4b)、加速度对比图(图4c)和加工误差对比图(图4d)。

分析图4a可知,三瓣花逼近的实际最大误差为3.24μm,位于给定的逼近内外公差范围内,验证了本拟合算法的可行性。由图4b、图4c可知,本文算法的加速度没有像微段加工方式那样发生突变,比较平稳,进给速度较微段加工方式也更平滑,且具有更高的加工效率。由图4d可知,本文算法的加工误差也更小,且更均匀。

(d)加工误差对比图

4 结语

为了消除微段加工方式下加工程序量过大和需要频繁加减速的诟病,实现微小直线段的高速平滑加工,提高加工效率,本文提出了基于带权因子和一阶导数约束的NURBS曲线最小二乘逼近算法,并在此基础上进一步给出了微段平滑加工算法。将离散的微段数据点拟合成一条连续的NURBS曲线并将其作为新的加工路径,然后根据该加工路径进行NURBS实时插补,从而实现微小直线段的平滑加工。验证结果表明,该算法有助于实现微小直线段的平滑加工。