基于轮廓约束点的B样条曲面拟合算法

江本赤　韩　江　田晓青　夏　链

合肥工业大学，合肥，230009

基于轮廓约束点的B样条曲面拟合算法

江本赤韩江田晓青夏链

合肥工业大学，合肥，230009

提出了一种面向截面测量数据的B样条曲面拟合算法。首先对原始数据点列进行降噪处理，然后遴选出曲率优势点，并将其作为初始的轮廓约束点，得到插值于约束点的初始曲线。再在需改善拟合精度的区域增加约束点，直至获得满足精度要求的B样条曲线。最后以约束点数目最多的曲线为准，在其余的曲线上增加差额数目的约束点，并进行平均弦长参数化，构造出B样条曲线簇，最终获得B样条拟合曲面。仿真实验结果表明，该方法可显著压缩曲面模型的控制顶点数目，具有较高的曲面重构效率。

逆向工程；轮廓约束点；B样条；曲面拟合

0　引言

由截面测量数据进行曲面重构是科学计算可视化(visualization in scientific computing,VISC)的重要内容之一，在医学可视化、无损探伤和逆向工程等领域应用广泛[1]。B样条具有凸包性、连续性与光滑性等诸多优点，故在曲面重构中被广泛采用[2-3]。对于数目庞大的测量数据点，若采用直接插值的方法获取B样条曲面，将导致节点和控制顶点数量过大，给后续的曲面光顺与修改、数据交换等造成负担[4-5]。因此，在满足数据轮廓精度要求的前提下，如何使曲面的数学模型包含尽量小的数据量，一直是该领域研究热点问题之一。

文献[6-11]对B样条曲面拟合中的相关技术，包括散乱点数据参数化和数据压缩等进行了研究。早期有代表性的插值方法[6]是先插值各截面数据点列构造B样条曲线，再通过节点插入获得可相容的曲线簇，最后对曲线簇的控制顶点放样得到插值曲面。这种方法较直观，但由于引入了节点插入操作，势必会引起控制顶点数量的膨胀[2]。为此，Piegl等[7]构造了具有非退化特性的参数区间，以增加节点选择的灵活性，使所得的控制顶点数目明显减少；在此基础上，Park[8]采用能量极小化方法进一步压缩了曲面模型的数据量。王文珂等[10]借鉴传统蒙皮算法的思想，提出了一种由无序B样条曲线来拟合B样条曲面的算法，可得到满足误差要求的光滑拟合曲面，但当曲线数目较多时该算法会遇到数值稳定性问题。综合上述情况，对B样条曲面拟合相关技术展开进一步研究仍然非常必要。

本文提出一种面向截面测量数据的拟合算法。首先采用弦高法剔除原始数据点中的坏点，以得到较光滑的线形数据点列；然后求出离散曲率值，将超过平均曲率值的局部曲率极大值点定义为曲率优势点，并作为曲线轮廓的初始约束点，得到初始的插值B样条曲线；再计算出偏差分布情况，并将超过允许误差的局部偏差极大值点增补为轮廓约束点，直到满足拟合精度要求；然后以约束点数最多的曲线为准，将其余曲线的约束点补齐以得到拓扑矩形的数据点阵；最后计算出公用的节点矢量，获得最终的B样条单向蒙皮曲面。对一组截面测量数据的曲面拟合实验验证了算法的可行性。

1　截面B样条曲线的获取

对于各截面测量数据，先用弦高法对截面数据进行消噪处理，然后根据点列的曲率分布情况，确定曲率优势点并将优势点作为型值点，计算出各截面的初始B样条曲线。

1.1数据预处理

扫描数据点中一般会包含失真点，它们不能反映实际轮廓且对曲线光顺性存在较大影响，称作噪声点或坏点。在截面扫描数据点中，如果某个数据点偏离其相邻的两个点，且偏离程度超过预设值，则可视为坏点。关于剔除坏点的研究报道较多，针对截面测量数据的线形分布特点，本文采用文献[12]中的弦高差方法来去除坏点，对于考察点Pi，连接两个相邻点得到弦线Pi-1Pi+1，如果点Pi到该弦的距离大于允许误差，则将该点剔除。图1所示是某截面数据预处理效果。

图1　预处理后的某截面数据

下面以图1中的数据点为例，阐述获取截面B样条曲线的具体步骤。

1.2轮廓约束点的确定

1.2.1离散曲率计算

本文采用文献[13]给出的近似法，将数据点Pi处的曲率Ki近似为3个相邻点Pi-1、Pi和Pi+1所在圆弧的曲率。图1中数据点的曲率分布情况如图2所示。

图2　离散点的曲率求解

1.2.2轮廓约束点的遴选

经预处理之后的密集数据点，其曲率分布能大体反映出轮廓走势。图1中的点列共包含202个数据点，为了压缩拟合曲线的数据量，现通过曲率值来确定初始轮廓约束点，以下简称约束点。

首先考虑曲率极大值法。即数据点Pi处的曲率值Ki，需同时满足条件Ki>Ki-1和Ki>Ki+1。找出曲率极大值所对应的数据点，并将其作为约束点。则图1中的数据点包含67个约束点，其拟合情形如图3所示。由此得到数据点与拟合曲线间的最大偏差值为0.144 mm，平均偏差值为0.051 mm，均方差为0.028 mm，偏差分布情况见图4。

图3　曲率极大值法拟合效果

图4　图3对应的拟合偏差分布

从图3所示的约束点分布情况来看，该方法较有效地压缩了数据量，但即使在相对平直的局部区域，也分布了比较密集的约束点，显然没有必要。

再来考察曲率临界值法。事先设定曲率临界值Kcrt，将曲率值满足Ki>Kcrt的数据点作为约束点。不失一般性，此处将离散曲率的平均值Kavg设为临界值，即Kcrt=Kavg。对于图1中的数据点，用此方法确定的约束点为64个，所得拟合曲线如图5所示。最大拟合偏差值为0.969 mm，平均偏差值为0.100 mm，均方差为0.229 mm，偏差分布情况见图6。

图5　曲率临界值法拟合效果

图6　图5对应的拟合偏差分布

由图5可见，约束点可提取轮廓弯曲的局部区域，但由于该区域的离散曲率值一般都大于平均值，故约束点非常密集，而相对平直的区域约束点过少，这种分布显然也是不可取的。

上述两种确定约束点的方法，各有优劣。为利用二者的优点，本文选取超过平均曲率值的曲率极大值点作为约束点。即同时满足以下两个条件：①Ki>Ki-1且Ki>Ki+1；②Ki>Kavg。此处将符合上述条件的数据点称为曲率优势点，并将其作为轮廓约束点，对于开曲线，则需将首尾两个端点直接当作约束点。

采用曲率优势点法，则图1中的数据点只包含31个约束点，其拟合情形如图7所示。由此得到最大拟合偏差值为0.918 mm，平均偏差值为0.164 mm，均方差为0.261 mm，偏差分布情况见图8。

图7　曲率优势点法拟合效果图

图8　图7对应的拟合偏差分布

与图3和图5相比，图7中的约束点数目得到了进一步的缩减，同时能在一定程度上保证数据点的轮廓精度。下面将插值于曲率优势点的曲线作为初始曲线。

诚然，由曲率优势点所确定的曲线，尽管大大压缩了数据量，但并不能保证拟合精度。若初始曲线存在拟合偏差超过预设值的局部区域，则需进行局部优化。

1.3拟合精度的优化

由图8可以看出，只有少量数据点处的拟合偏差值为零，因此偏差曲线一般会出现多个峰值，记允许误差为[δ]，若初始曲线的某些峰值大于[δ]，则直接将偏差峰值对应的数据点增设为约束点。也就是说，新的约束点将由两部分构成，即曲率优势点和少数偏差峰值数据点。

以图7中的初始拟合曲线为例，设[δ]=0.5 mm，对照图8可知，需要增加2个约束点，新曲线的拟合效果如图9所示。相应的拟合偏差最大值为0.235 mm，平均偏差值为0.060 mm，均方差为0.057 mm，分布情况见图10。

图9　初始曲线的局部拟合精度优化

图10　图9对应的拟合偏差分布

可见，增加约束点之后，曲线局部区域的拟合精度得到了明显改善，而此时的轮廓约束点数为33个。该方法与曲率极值法和曲率临界值法相比，显然更加合理。当然，若偏差峰值仍大于[δ]，则需采用相同方法增加约束点，直到获得满足要求的B样条曲线。

2　B样条曲面的构造

利用上述方法，可以获得满足精度要求的各截面B样条曲线，但是这些曲线是孤立的，缺乏相容性。各截面数据轮廓曲率变化各异，因而各节点矢量也不尽相同。由于B样条曲面具有张量积特性，需要解决曲线之间的相容性问题。为此，下面设法获取各截面曲线公用的节点矢量。

2.1公用节点矢量的获取

2.1.1约束点的参数化

为计算公用节点矢量，首先必须保证各曲线具有相同数量的约束点。前面得到的孤立B样条曲线簇，尽管已满足拟合精度要求，其拟合偏差曲线仍会存在多个峰值。现将各峰值点对应的数据点，按其偏差值由大到小的顺序存入各自数组，称作备用约束点数组。

记第i条曲线的现有约束点数目为Ci，Ci的最大值记为Cmax。为确保各曲线具有相同的节点数目，现从各备用数组中依次取出必要数目的约束点，对于第i条曲线，需取出Cmax-Ci个数据点作为约束点。至此，所有曲线的约束点数目都为Cmax。

接下来实现各约束点的参数化。为尽量考虑型值点的分布对曲线形状的影响，此处采用平均弦长参数化方法。设共有k条截面曲线，对于每一条曲线，依次连接各约束点得到Cmax-1段弦线，记第i条曲线的第j个弦长值为Li,j，则曲线簇第j段的弦长平均值Lj为

(1)

j=1,2,…,Cmax-1

因此，用来计算节点矢量的总弦长d为

(2)

而其余参数

(3)

t=1,2,…,Cmax-2

2.1.2节点矢量的计算

记曲线次数为p，m=Cmax+p。下面计算曲线簇公用的节点矢量U。其中，u0=u1=…=up=0，um-p=um-p+1=…=um=1，中间的节点采用平均值法，用下式求出：

(4)

j=1,2,…,Cmax-p-1

这种方法得到的控制顶点与约束点的数目相等，避免了非奇异矩阵的产生，便于后续控制顶点的计算。

2.2控制网络的反算

(5)

图11所示的是B样条曲线簇的矩形拓扑控制网络示意图。

图11　控制顶点网络

至此，曲面的节点矢量和控制顶点均已求出，即B样条曲面已被定义。

3　仿真实验与分析

为验证本文算法的可行性与有效性，下面以一组截面点云数据为例，用一张双三次B样条曲面对其进行拟合。图12所示的是经预处理之后的截面数据点，共21条线性点列，包含2557个数据点，分布在40mm×50mm×12mm的空间里，平均每条点列包含约120个数据点。

图12　截面数据举例

按照上述步骤，通过遴选曲率优势点获得了初始拟合曲线，并以[δ]=0.5 mm为逼近偏差约束条件，对局部区域进行了精度优化，最终得到了V向B样条曲线簇(即各曲线的参数u为常数)，图13所示的是其中的一条曲线。单条曲线包含的控制顶点数目最多为22个，最少为13个。

图13　三次B样条曲线

为满足曲线簇的相容性要求，将所有曲线的轮廓约束点都增加至22个。然后运用平均弦长参数法，得到公用节点矢量为

U=(0，0，0，0，0.047，0.095，0.142，0.190，0.238，0.285，0.333，0.381，0.428，0.476，0.523，0.571，0.619，0.667，0.714，0.761，0.809，0.857，0.904，0.952，1，1，1，1)

V=(0，0，0，0，0.067，0.153，0.179，0.186，0.229，0.243，0.309，0.381，0.398，0.434，0.515，0.551，0.588，0.629，0.667，0.695，0.715，0.750，0.811，0.923，1，1，1，1)

所得的控制网络为21行×22列，即共有462个控制顶点，如图14所示。最终重构的B样条蒙皮曲面见图15。与原数据点数目相比，减少了2095个，即减少81.97%。

图14　重构曲面的控制网络

图15　B样条曲面重构

4　结论

(1)通过曲率优势点法所确定的轮廓优势点，可在提取截面数据基本形状的前提下，大幅缩减曲面模型的数据量，为后续的数据运算奠定基础。

(2)将初始拟合的偏差峰值对应点增设为B样条曲线型值点，可快速有效地提高局部逼近精度。

(3)与曲率极大值法和曲率临界值法相比，利用曲率优势点法所得型值点的分布更加合理，可避免曲率较大区域的型值点过于密集，同时对于较平直的区域亦具有较高的轮廓保真度。

[1]王瑜,郑津津,周洪军,等.断层轮廓的双三次非均匀B样条曲面重构[J].中国科学技术大学学报, 2011, 41(5):384-391.

Wang Yu, Zheng Jinjin, Zhou Hongjun, et al. Bi-cubic Non-uniform B-spline Surface Reconstruction for Slice Contours[J].Journal of University of Science and Technology of China,2011,41(5):384-391.

[2]张丽艳,周来水,蔡炜斌,等.基于截面测量数据的B样条曲面重建[J].应用科学学报,2002, 20(2): 173-177.

Zhang Liyan, Zhou Laishui, Cai Weibin, et al. B-spline Surface Reconstruction Based on Digitized Section Data[J]. Journal of Applied Sciences, 2002, 20 (2): 173-177.

[3]李涛,刘浩.B样条曲面的双向插值造型算法[J].计算机工程与应用, 2012,48(35):178-181.

Li Tao, Liu Hao. Bidirectional Interpolation Modeling Algorithm for B-spline Surface[J]. Computer Engineering and Applications, 2012, 48(35):178-181.

[4]杜兴吉,周晓军,吴瑞明,等.超声检测中复杂曲面数字化方法研究[J].中国机械工程,2005, 16(8):686-689.

Du Xingji,Zhou Xiaojun,Wu Ruiming,et al. Research on Sculptured Surface Digitizing on Ultrasonic Inspection[J]. China Mechanical Engineering,2005,16(8):686-689.

[5]Woodward C. Skinning Techniques for Interactive B-spline Surface Interpolation[J]. Computer-Aided Design，1988，20(8):441-451.

[6]李运平,周来水,王志国,等.基于刚架模型的多张B样条曲面变形[J].中国机械工程, 2013，24(24):3355-3365.

Li Yunping, Zhou Laishui, Wang Zhiguo,et al. Shape Modification of Multi B-spline Surfaces Via Rigid Frame[J]. China Mechanical Engineering, 2013, 24(24):3355-3365.

[7]Piegl L，Tiller W. Reducing Control Points in Surface Interpolation[J].IEEE Computer Graphics and Applications, 2000, 20(5):70-74.

[8]Park H. Lofted B-spline Surface Interpolation by Linearly Constrained Energy Minimization[J].Computer-Aided Design, 2003, 35(14): 1261-1268.

[9]Weiss V, Andor L, Renner G, et al. AdvancedSurface Fitting Techniques[J].Computer Aided Geometric Design, 2002, 19(1):19-42.

[10]王文珂,李思昆.无序B样条曲线的曲面拟合算法[J].计算机辅助设计与图形学学报,2013, 25(5):679-685.

Wang Wenke, Li Sikun. B-spline Surface Fitting to Unorganized Curves[J].Journal of Computer-Aided Design & Computer Graphic, 2013, 25(5):679-685.

[11]刘俊,王启富,陈立平,等.B样条曲线曲面的二阶预测压缩算法[J].中国机械工程,2008,19 (3):304-307.

Liu Jun,Wang Qifu,Chen Liping,et al. Compression of B-spline Based on Local Coordinate Second-order Prediction[J]. China Mechanical Engineering, 2008, 19(3):304-307.

[12]何丽辉.反求工程中截面点云数据的预处理[J].鞍山师范学院学报, 2012,14(6):53-56.

He Lihui.Pre-process of Reversing Section Point Cloud Data in Engineering[J]. Journal of Anshan Normal University, 2012, 14(6):53-56.

[13]Zhou K, Wang G J, Jin H Z, et al. NURBS Interpolation Based on Exponential Smoothing Forecasting[J]. Int. J. Adv. Manuf. Technol., 2008, 39:1190-1196.

[14]Les P, Wayne T. 非均匀有理B样条[M].2版.赵罡,穆国旺,王拉柱,译.北京:清华大学出版社, 2010.

(编辑王艳丽)

B-spline Surface Fitting Algorithm Based on Contour’s Constraint Points

Jiang BenchiHan JiangTian XiaoqingXia Lian

Hefei University of Technology,Hefei,230009

Aiming at cross-sectional measured data, a B-spline surface fitting method was proposed. The original data points were denoised firstly, then vantage points of curvature were picked out and taken as the initial contour constraint points, and an initial curve interpolated with the constraint points was gained. After that, constraint points were added to the region whose fitting accuracy needed to be improved, until the B-spline curve could reach preset accuracy. Finally, subjecting to curve which contained the most constraint points, constraint points of the balance number were added to the rest of the curve, and average chord-length parameterization was conducted to construct a family of B-spline curves, and a B-spline fitting surface was obtained at last. Simulation results show that, this method can compress the number of control vertices of the surface model significantly with a high efficiency of surface reconstruction.

reverse engineering;contour’s constraint point; B-spline; curved surface fitting

2014-09-28

国家自然科学基金资助项目(51275147)

TP391.72< class="emphasis_italic">DOI

：10.3969/j.issn.1004-132X.2015.15.006

江本赤，男，1979年生。合肥工业大学机械与汽车工程学院博士研究生。主要研究方向为数控技术与装备。发表论文10余篇。韩江，男，1963年生。合肥工业大学机械与汽车工程学院教授、博士研究生导师。田晓青，女，1987年生。合肥工业大学机械与汽车工程学院讲师、博士。夏链(通信作者)，女，1964年生。合肥工业大学机械与汽车工程学院教授、博士。