DOI: 10.11817/j.issn.1672-7207.2018.07.019

基于曲面拟合的复杂自由曲面网格划分

陈礼杰¹，吴慧²，李铁瑞¹，高博青¹

(1. 浙江大学 建筑工程学院，浙江 杭州，310058；

2. 浙江财经大学 公共管理学院，浙江 杭州，310018)

摘要：为了实现复杂自由曲面的建筑网格划分，将复杂自由曲面拟合成1个曲面，并利用映射法进行网格划分。首先提出基于非均匀有理B样条(NURBS)的线面求交和基于NURBS的面面求交2种方法，获得复杂自由曲面的N×M点云；将N×M点云作为控制点构造NURBS曲面，然后把修剪后的NURBS曲面边界按给定数目或给定长度划分，连接参数域上的边界节点后得到参数域网格，将其映射回空间曲面。对边界处网格进行投影处理以减少曲面拟合出现的边界褶皱及内缩问题。采用加密点阵技术减小曲面拟合的误差。研究结果表明：与多曲面网格划分方法相比，采用曲面拟合及映射法能够得到线条流畅、大小基本一致的网格，具有良好的建筑美学效果。

关键词：非均匀有理B样条(NURBS)；复杂自由曲面；曲面拟合；映射法；网格划分

中图分类号：TU393.3        文献标志码：A         文章编号：1672-7207(2018)07-1718-08

Grid generation on complex free-form surface based on surface fitting

CHEN Lijie¹, WU Hui², LI Tierui¹, GAO Boqing¹

(1. College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China;

2. Public Administration College, Zhejiang University of Finance & Economics, Hangzhou 310018, China)

Abstract: In order to obtain architectural mesh of complex free-form surface, the complex free-form surface was fitted into one surface and then grids were generated using mapping method. First, two methods called NURBS based on ray-surface intersection and NURBS based on plane-surface intersection were proposed to obtain the N×M point cloud of complex free-form surface. Secondly, a NURBS surface was constructed using the N×M point cloud as control points, then the boundaries of trimmed NURBS surface were divided by a given number or a given length, and the boundary nodes of the parametric domain were connected to get parametric domain grids. After that, the grids were mapped onto the free-form surface. Lastly, the grids near the boundaries were projected onto original surfaces to reduce boundary fold and shrinkage caused by surface fitting. To reduce the error of surface fitting, the number of point cloud was increased. The results show that compared with the method which respectively generates grids on multiple trimmed surfaces, the proposed methods produce fluent and uniform grids, with good architectural aesthetics.

Key words: non-uniform rational B-splines (NURBS); complex free-form surface; surface fitting; mapping method; grid generation

随着计算机科学技术的进步、建筑施工工艺的发展和人们建筑审美水平的提高，自由曲面结构以其优美的造型越来越多地出现在人们的视野当中。“自由曲面”指无法用单个或几个解析函数表达的曲面，也可理解为那些明显区别于传统建筑造型的曲面^[1]。由于解析曲面具有一定的周期性或规律性，其网格划分比较方便。而自由曲面的不规则性使传统的网格划分方法难以适用于自由曲面结构。因此，如何在自由曲面上生成线条流畅、大小基本一致的网格，是空间结构领域的热点和难题^[2]。网格划分起源于有限元领域^[3]，常用的方法有Delaunay法^[4]、波前法^[5]、映射法^[6]等。但有限元领域的网格划分方法难以满足建筑美学线条流畅、大小基本一致的要求。因此，近年来有研究者针对自由曲面的建筑网格划分进行了相关研究，如：WINSLOW等^[7]以杆件走向为优化变量，以结构性能为目标函数，用遗传算法实现自由曲面的多目标优化；SU等^[8]绘出自由曲面的主应力迹线，采用改进的波前法自动生成网格；潘炜等^[9]按照展开前、后面积变化最小的原则展开曲面，并利用线推进法生成平面网格，最后映射回空间曲面；丁慧等^[10]提取等参线并以杆件长度相等为条件分割等参线，连接等参线上的节点生成网格。复杂自由曲面由多个非均匀有理B样条(NURBS)裁剪曲面组成，且存在内部空洞等复杂情况。针对该类自由曲面，现有的网格划分方法是在每个裁剪曲面上分别划分网格，然后对交界附近网格进行调整。该方法的缺点是网格划分不易进行，交界处网格不好调整，难以达到建筑美学线条流畅的要求。本文作者提出基于NURBS的线面求交和面面求交2种方法，获得复杂自由曲面的N×M点云。第1种方法是构造矩形N×M平面点阵，过每个点作平行于Z坐标轴的射线与NURBS曲面求交，得到N×M点阵P_i,j；第2种方法是根据NURBS曲面的特征，构造N个平面并与NURBS曲面求交，对交线均匀采样M个点，修补边界后得到点阵P_i,j。利用向心参数化方法得到点阵P_i,j的U和V向节点矢量，构造以P_i,j为控制点阵的NURBS曲面。将修剪后的NURBS曲面边界按给定数目或给定长度划分，连接参数域上的边界节点后形成参数域网格，将其映射回空间曲面，对边界处网格作投影处理，即可得到线条流畅、大小基本一致的四边形网格。

1  曲面拟合

曲面拟合技术在逆向工程和计算机辅助设计等领域中得到了广泛应用^[11]。利用计算机视觉技术(如激光扫描)对自由曲面进行数字化，可得到大量散乱数据点，根据数据点集可得到实物的CAD模型^[12]。

在建筑工程领域，逆向工程用于古建测绘和古建保护、快速建立城市模型、空间网格结构的逆向建模等方面。李娜^[13]借鉴逆向工程的操作原理，利用激光三维扫描系统进行实物样件的点云数据采集，采用蒙皮法进行曲面造型，然后进行网格划分。该技术基于对实物进行激光扫描得到大规模散乱数据点。在复杂自由曲面网格划分时，一般建筑师已经提供符合其设计意图的曲面CAD模型；为了获得线条流畅、大小基本一致的网格，本文作者提出基于NURBS的线面求交和面面求交2种方法，获得复杂自由曲面的N×M点云，构造以N×M点云为控制点阵的NURBS曲面，为映射法划分网格提供基础。

图1所示为某工程CAD模型。下面以图1为例，介绍基于NURBS的线面求交法和基于NURBS的面面求交法。

图1  某工程CAD模型

Fig. 1  CAD model of a project

1.1  基于NURBS的线面求交法

根据NURBS曲面的强凸包性^[14]和控制点阵为矩形拓扑^[12]的特点，构造1个N×M矩形平面点阵作射线与NURBS曲面求交，提取原曲面信息。算法流程如图2所示。

1) 延伸曲面上边界。由于N×M空间矩形点阵在包围盒底面，使得以P_i,j为控制点阵的拟合曲面在上边界处产生内缩。因此，适当延伸曲面上边界，得到拟合曲面后对其进行裁剪，便可避免内缩。延伸曲面上边界如图3所示。

2) 构造曲面包围盒。CAD模型有K个NURBS曲面(其中，k=1, 2, …, K；为第k个曲面的控制点；为控制点对应的有理基函数)。第k个曲面有(nk+1)×(mk+1)个点的控制点阵。找出所有控制点阵中控制点最大X坐标值，记为X_max；同理，找出曲面模型的X_min，Y_max，Y_min，Z_max和Z_min。利用这6个值即可在模型空间内绘出曲面的包围盒，如图4所示。根据NURBS曲面的强凸包性，第k个NURBS曲面一定位于其控制点阵形成的凸包内，因此，K个NURBS曲面一定包含在曲面的包围盒内。

图2  基于NURBS的线面求交法流程图

Fig. 2  Flow chart of NURBS based ray-surface intersection method

图3  延伸曲面上边界

Fig. 3  Extend upper boundaries of surfaces

3) 以包围盒底面为基础，构造N×M矩形点阵。将包围盒底面在X和Y轴的正、负2个方向上适当延伸，向Z轴负向平移一定距离，根据CAD模型的曲面个数及复杂程度确定矩形点阵X和Y方向的点数N和M，在延伸和平移后的底面上构造1个N×M矩形点阵(见图4)。

图4  N×M矩形点阵及曲面包围盒

Fig. 4  N×M rectangular lattice and surface bounding box

4) 遍历所有平面点作平行于Z轴的射线与曲面求交。

NURBS曲面与射线求交是一种基础的几何算法，常用于雷达散射截面计算(RCS)的射线追踪计算^[15]、逆向工程中曲面的制造误差计算^[16]等方面。

过平面点作平行于Z轴的射线与NURBS曲面求交点即求解二元非线性方程组：

式中：和分别为平面点在X和Y轴上的坐标值；和分别为曲面控制点在X和Y轴上的坐标值；为控制点对应的有理基函数。

首先利用Krawczyk算子判断式在某一区间是否存在解，然后对存在唯一解的区间给定1个初值，利用牛顿迭代法求解。详细流程参见文献[16]。

复杂自由曲面中可能存在2种曲面：一种是完整NURBS曲面，另一种是裁剪NURBS曲面。对于完整的NURBS曲面，上述方程求得的交点即为有效交点。而裁剪NURBS曲面只是将NURBS曲面的删除部分不予显示，因此，射线与曲面的交点可以分为有效交点和无效交点，如图5所示(其中，有效交点位于保留部分，无效交点位于删除部分)。

作出平面点的射线L后，分别与CAD模型中的K个NURBS曲面求交。假定射线与整个模型只有1个有效交点，当射线L与第k个曲面求得有效交点P后，令P_i,j=P，进行下一个平面点的求交。若射线L与K个NURBS曲面都没有有效交点，则令。由此，得到点阵P_i,j。

图5  裁剪NURBS曲面的有效交点与无效交点

Fig. 5  Valid intersections and invalid intersections of trimmed NURBS surface

1.2  基于NURBS的面面求交法

前面提出一种提取原曲面信息的方法，现通过构造N个平面与NURBS曲面求交提取原曲面信息。算法流程如图6所示。

1) 补齐曲面空洞。当复杂自由曲面内部存在空洞时，为了保证平面与曲面的交线在平面上能够拟合为一条曲线，用平面补齐空洞。补齐空洞及平面切割方式如图7所示。

图6  基于NURBS的面面求交法流程图

Fig. 6  Flow chart of NURBS based plane-surface intersection method

图7  补齐空洞及平面切割方式

Fig. 7  Patterns of filling holes and plane cutting

2) 选择平面切割方式。复杂自由曲面由多个NURBS裁剪曲面组成，应根据曲面各部分的走向和特征，选择合适的切割方式，如直线式和圆弧式切割(见图7)。原曲面由11个NURBS曲面组成(不包括补齐空洞的平面)，分为5部分切割。第1部分由XY平面上点1和点2及平面数N₁确定。首先过点1和点2作直线L，然后以L为正法线，等距离地作N₁个平面与曲面求交。第2部分由XY平面上点3、点4、点5及平面数N₂确定。首先过点4以直线L为正法线作平面P，然后以点4为旋转中心，Z轴为旋转轴，将平面P旋转N₂-1次，每次旋转角度为∠345/(N₂-1)，即可得到N₂个平面与曲面相交。其余3个部分依此类推。

3) NURBS曲面与平面求交。用平面与CAD模型的K个曲面求交，如图8所示。平面与4个曲面相交，共有5条交线，其中与曲面2有2条交线即交线2和4。

4) 处理同一平面上的交线(见图8)。首先处理裁剪曲面与平面求交的情况。曲面1是裁剪曲面，其删除部分不显示，因此，把位于曲面1删除部分且与平面相交的交线删除。对曲面4进行同样处理。平面与NURBS曲面共有5条交线，每条交线均匀采样一定数目点，然后构造1条插值于所有采样点的NURBS曲线。每个平面都得到1条与原曲面相交的曲线。

图8  同一平面上交线处理

Fig. 8  Treatment of intersecting lines on the same plane

5) 拟合曲线采样。曲面拟合中NURBS曲面的型值点或控制点必须为矩形拓扑，故得到平面和原曲面拟合后的交线后，对该交线均匀采样M个点，由此保证NURBS曲面与N个平面求交后得到的点阵P_i,j为矩形拓扑。

6) 补齐边界。对于与NURBS曲面相交的第1个和最后1个平面(下面称边界平面)，总存在一部分曲面在边界平面之外，其大小不能忽略。取这部分曲面的边界线，并均匀采样M个点，可以保证完整地提取原曲面的信息。补齐边界方式如图9所示。

图9  补齐边界

Fig. 9  Filling vacancies of boundaries

经过以上步骤，可以得到矩形拓扑点阵P_i,j。

1.3  拟合NURBS曲面

利用曲面拟合中常用的向心参数化^[14]方法构造矩形拓扑点阵P_i,j的节点矢量U和V。获得点阵P_i,j和节点矢量U和V之后，进行曲面拟合。

对于矩形拓扑点阵P_i,j，传统的曲面拟合方法是将其作为型值点阵Q_i,j，创建1个非有理的(p，q)次B样条曲面S(u,v)，使其插值于这些点。

式中：Q_i,j为型值点；u_i和v_j分别为其对应的U和V向参数值；R_i,j(u_i, v_j)为其对应的有理基函数；P_k,l为待求的控制点。

另一种方法用1个非有理的(p，q)次B样条曲面S(u，v)逼近这些点，以控制点为优化变量，满足如下条件：

为了准确地表示复杂自由曲面，其矩形拓扑点阵P_i,j的点数N×M不能太少。若使用插值或逼近方法，则至少需要解1个(N×M)×(N×M)的线性方程组，这将使计算时间大大增加。

针对复杂自由曲面拟合问题，采用控制点逼近方法。控制多面体是对B样条曲面的平面片逼近，次数越低，逼近效果越好；当节点矢量无限加细时，控制多面体将收敛于曲面^[17]。因此，以P_i,j为控制点阵，在节点矢量U和V上构造非有理的双二次B样条曲面S(u, v)，根据CAD模型的原有边界情况对S(u, v)进行裁剪。基于NURBS的线面求交法和面面求交法的拟合曲面与裁剪曲面分别如图10和图11所示，2种求交方法的拟合曲面虽不尽相同，但裁剪曲面都能反应复杂自由曲面的特征，可为映射法网格划分提供基础。

图10  基于NURBS的线面求交法的拟合曲面与裁剪曲面

Fig. 10  Fitting surface and trimmed surface of NURBS based ray-surface intersection method

图11  基于NURBS的面面求交法的拟合曲面与裁剪曲面

Fig. 11  Fitting surface and trimmed surface of NURBS based plane-surface intersection method

2  映射法划分网格

图12所示为映射法划分网格示意图。建筑师提供的CAD模型由3个曲面构成(见图12(a))，利用基于NURBS的线面求交法，将3个NURBS裁剪曲面拟合为1个曲面(见图12(b))，利用映射法划分网格。

图12  映射法划分网格示意图

Fig. 12  Diagram of grid generation through mapping method

1) 将空间曲面边界分段。以二分法为基础，将空间曲面的边界按弦长分段。分段方式有给定长度、给定段数、等比例或指定长度数组等(见图12(b))。

2) 将分段点投影回参数域，连接参数域上的分段点后形成参数域网格。分段后的边界点反向映射回参数域中，在参数域中将边界点按照建筑师期望的走向依次连直线。

3) 将参数域网格映射回空间曲面，得到空间网格(见图12(c))。

4) 处理边界处网格。由图12(c)可知：用映射法在拟合曲面上画出的网格在原有曲面交界处线条流畅，但在拟合曲面的边界处质量不佳，为此需处理边界处网格。下面介绍该具体方法。

3  算例与比较

3.1  2种求交方法的比较

前面以图1为例详细介绍了基于NURBS的线面求交法(以下简称方法一)及基于NURBS的面面求交法(以下简称方法二)的算法流程。表1所示为2种方法的比较结果(其中，t为将复杂自由曲面拟合为1个曲面所需的时间，△_U为拟合曲面上边界最大误差，△_L为拟合曲面下边界最大误差)。

方法一耗时最久的步骤是曲面与直线求交，需进行49万(700×700)次求交运算；方法二耗时最久的步骤是曲面与平面求交，只需进行900次求交运算。故方法二的速度远快于方法一的速度。

表1  2种方法的比较

Table 1  Comparison of two methods

图13所示为方法一边界处褶皱。方法一中，当曲面边界走向不严格为平面点阵走向时，边界附近的控制点不能很好地还原曲面的边界信息，出现“褶皱”问题，导致上、下边界最大误差较大。此时，可以通过加密平面点阵的方法来减小褶皱，但会使时间大大增加。

图14所示为方法二下边界误差。原曲面下边界出现直角转折，在下边界附近的3个采样点由此也产生较大转折(见方法二步骤5))。因为NURBS曲面的强凸包性，拟合曲面在此处出现最大误差为250.5 mm，可以通过加密采样点的方法减小误差。而上边界处转折较小，最大误差仅为2.1 mm。

图13  方法一边界处褶皱

Fig. 13  Boundary fold of method 1

图14  方法二下边界误差

Fig. 14  Lower boundary error of method 2

2种方法除边界处出现误差外，其余部分误差均小于5 mm。用映射法网格划分时需将边界处网格投影至原曲面以消除误差，方法一还需手动处理因“褶皱”而畸形的网格。

通过2种方法的比较可以得出以下结论：1) 方法二比方法一复杂；2) 方法二的速度远快于方法一的速度，最大误差也比方法一的小；3) 利用映射法划分网格时，2种方法均需将边界处网格投影回原曲面，方法一还需手动处理因“褶皱”而畸形的网格。

3.2  分曲面划分网格与曲面拟合划分网格比较

以某工程CAD模型(见图12(a))为例，用分曲面的方法划分网格，如图15所示。首先将曲面1的边界1等长度分段，根据欲划分的网格形式确定边界2的分段点，最后连接边界1分段点或同时连接边界1、边界2分段点(见图15(a))。划分好曲面1的网格后，凭经验确定曲面2上边界3的分段点，连接边界2分段点或边界2、边界3分段点(见图15(b))。同理画出曲面3的网格(见图15(c))。但上述过程依赖于工程师的经验，且网格划分容易失败。

在边界2附近，曲线1为边界1分段点相连所得曲线，曲线2为边界1和2分段点相连所得曲线，因此，导致网格大小不均匀(见图15(a))。由于分别在3个曲面上划分网格，导致边界处网格不均匀，线条不流畅(见图15(c))。对边界处质量不佳的曲线进行一定调整，调整后网格均匀性有一定提高，但交界处曲线流畅性仍然不佳(见图15(d))。

图15  分曲面划分网格

Fig. 15  Respective generation of grids on multiple trimmed surfaces

现用曲面拟合的方法对图12(a)进行网格划分。该工程平面投影为椭圆形，平面面积为146.0 m×116.0 m。用方法一获得600×600点阵，拟合为1个曲面，等分外边界，在参数域上将外边界分段点依次连为直线后将参数域网格映射回空间曲面(见图12(c))。将边界处网格投影回原曲面，并处理极少数因“褶皱”而畸形的网格，得到最终的网格划分(见图12(d))。

图15中在原有的3个曲面上分别划分网格，因此，在交界附近网格不均匀，线条过渡不流畅，且难以调整，网格划分过程依赖工程师个人经验，且容易失败；图12(d)中线条流畅，网格大小基本一致，且过渡均匀。虽然在边界处出现一定误差，但可以通过加密点阵减小误差，或将边界处网格投影回原曲面，以消除边界处误差。

4  结语

1) 提出2种点云生成方法：方法1中构造1个N×M矩形平面点阵并作射线与复杂自由曲面求交，从而获得矩形拓扑点阵；方法2需合理选择平面切割方式，将复杂自由曲面分为几个部分，构造N个平面与复杂自由曲面求交，并将多条交线拟合后均匀采样M个点，修补边界后获得矩形拓扑点阵。

2) 将二者进行比较，令点云为控制点阵，把复杂自由曲面拟合为1个NRUBS曲面，以映射法划分空间网格。

3) 与分曲面划分方法相比，曲面拟合和映射法能够得到线条流畅、大小基本一致的四边形网格。

参考文献：

[1] 沈世钊, 武岳. 结构形态学与现代空间结构[J]. 建筑结构学报, 2014, 35(4): 1-10.

SHEN Shizhao, WU Yue. Structural morphology and modern space structures [J]. Journal of Building Structures, 2014, 35(4): 1-10.

[2] 李承铭, 卢旦. 自由曲面单层网格的智能布局设计研究[J]. 土木工程学报, 2011, 44(3): 1-7.

LI Chengming, LU Dan. Study of intelligent layout design of single-layer lattice shell of free form surface[J]. China Civil Engineering Journal, 2011, 44(3): 1-7.

[3] LO S H. Finite element mesh generation and adaptive meshing[J]. Progress in Structural Engineering & Materials, 2002, 4(4): 381-399.

[4] SI Hang. TetGen, a delaunay-based quality tetrahedral mesh generator[J]. ACM Transactions on Mathematical Software, 2015, 41(2): 11.

[5] PIRZADEH S. Structured background grids for generation of unstructured grids by advancing-front method[J]. AIAA Journal, 2012, 31(2): 285-289.

[6] LISEIKIN V. Grid generation methods[M]. Amsterdam, Netherlands: Springer, 2010: 10-14.

[7] WINSLOW P, PELLEGRINO S, SHARMA S B. Multi-objective optimization of free-form grid structures[J]. Structural & Multidisciplinary Optimization, 2010, 40(1/2/3/4/5/6): 257-269.

[8] SU Liang, ZHU Shunlai, XIAO Nan, et al. An automatic grid generation approach over free-form surface for architectural design [J]. Journal of Central South University, 2014, 21(6): 2444-2453.

[9] 潘炜, 吴慧, 李铁瑞, 等. 基于曲面展开的自由曲面网格划分[J]. 浙江大学学报(工学版), 2016, 50(10): 1973-1979.

PAN Wei, WU Hui, LI Tierui, et al. Grid generation on free-form surface based on surface flattening[J]. Journal of Zhejiang University: Engineering Science, 2016, 50(10): 1973-1979.

[10] 丁慧, 罗尧治. 自由形态网壳结构网格生成的等参线分割法[J]. 浙江大学学报(工学版), 2014, 48(10): 1795-1801.

DING Hui, LUO Yaozhi. Isoparametric line dividing grid generation method for free-form lattice shells[J]. Journal of Zhejiang University: Engineering Science, 2014, 48(10): 1795-1801.

[11] BRUJIC D, AINSWORTH I. Fast and accurate NURBS fitting for reverse engineering[J]. The International Journal of Advanced Manufacturing Technology, 2011, 54(5): 691-700.

[12] 来新民, 黄田. 基于NURBS的散乱数据点自由曲面重构[J]. 计算机辅助设计与图形学学报, 1999, 11(5): 433-436.

LAI Xinmin, HUANG Tian. NURBS based approach of surface reconstruction from scattered data[J]. Journal of Computer Aided Design and Computer Graphics, 1999, 11(5): 433-436.

[13] 李娜. 空间网格结构几何形态研究与实现[D]. 浙江大学建筑工程学院, 2009: 77-92.

LI Na. Research and implementation of geometric configuration of space grid structures[D]. Hangzhou: Zhejiang University. College of Civil Engineering and Architecture, 2009: 77-92.

[14] PIEGL L, TILLER W. The NURBS book[M]. New York, USA: Springer, 1995: 94-95.

[15] 官火梁, 吴强, 席平. RCS计算中NURBS曲面和射线求交的快速计算[J]. 图学学报, 2006, 27(1): 87-91.

GUAN Huoliang, WU Qiang, XI Ping. A fast algorithm for intersection calculation of ray and NURBS surface in predicting radar cross section calculation[J]. Journal of Engineering Graphics, 2006, 27(1): 87-91.

[16] 王保庆, 张俐, 李东升. 逆向工程中NURBS曲面与直线交点快速计算[J]. 图学学报, 2010, 31(2): 155-158.

WANG Baoqing, ZHANG Li, LI Dongsheng. Rapid calculation of intersection points between NURBS surface and line in Reverse Engineering[J]. Journal of Engineering Graphics, 2010, 31(2): 155-158.

[17] BOOR C D. Cutting corners always works[J]. Computer Aided Geometric Design, 1987, 4(S1/2): 125-131.

(编辑  伍锦花)

收稿日期：2017-07-02；修回日期：2017-08-20

基金项目(Foundation item)：国家自然科学基金资助项目(51378457，51678521)；浙江省自然科学基金资助项目(LY15E080017) (Projects(51378457, 51678521) supported by the National Natural Science Foundation of China; Project(LY15E080017) supported by the Natural Science Foundation of Zhejiang Province)

通信作者：高博青，博士，教授，从事空间结构研究；E-mail: bqgao@zju.edu.cn