基于离散化的复杂曲面建筑网格划分方法

李铁瑞 吴慧 王奇胜 高博青

摘 要：为实现存在裁剪、孔洞的复杂自由曲面建筑网格划分，提出了一种基于离散的、以均匀性为目标的劃分方法.将复杂曲面离散并缝合，形成由大量面片组成的离散曲面，作为多个参数曲面的一体化表示.采用改进的误差扩散算法，在离散曲面上按一定的密度进行初始布点.采用基于空间距离的粒子动力松弛算法对点云进行初步均匀化，并应用基于曲面距离的k均值算法进行再次均匀化.对均匀的点云求曲面距离的Voronoi图，并获得相应网格.对网格进行拓扑优化和光顺优化.算例表明，本文算法可有效处理存在裁剪、孔洞的复杂自由曲面，并得到均匀光顺的三角网格.

关键词：复杂曲面；离散化；网格划分；均匀化；松弛

中图分类号：TU311.41 文献标志码：A

文章编号：1674—2974（2018）07—0048—06

Abstract： A grid generation method for complicated multiple surfaces with trimmings and holes is presented. This method is based on the discretization and concentrates on the aim of homogeneity. The multiple surfaces are discretized separately and seamed together to achieve a discrete surface. The points are distributed on the discrete surface according to the density applying improved error-diffusion method. The points are homogenized by particle dynamics method with Euclid distance and then homogenized once more by k-means algorithm with surface distance. The Voronoi diagram with surface distance is delivered on the discrete surface to obtain the grids. The topological and smooth relaxations are applied on the grids. Eventually， the case study indicates that this method can solve the problem of grid generation for complicated multiple surfaces effectively and achieve the homogeneous and smooth grids.

Key words： complex surfaces；discretization；grid generation；homogenization；relaxation

随着计算机辅助设计技术，特别是建模技术的发展，建筑师的创造力可以得到更好的发挥，自由曲面的建筑形式因其强力的视觉效果得到了人们的青睐.但是，对复杂的曲面进行网格划分十分困难，尚没有行之有效的方法，是目前空间结构领域研究的一个难点与热点.

对自由曲面，采用显式或隐式方程都无法表达，常采用NURBS技术[1]，或结合形态优化方法[2]建模.而对于特别复杂的曲面，用单个NURBS曲面表达也甚为困难，乃至无法表达，此时需采用多个曲面来建模.T-Splines[3]技术产生于2003年，并由T-Splines 公司.商业化，近年来也有大量的研究，该技术对多曲面建模提供了良好的支撑.由于NURBS曲面是由高阶非线性参数方程组表示的，一些经典的算法难以应用于NURBS曲面，而对于裁剪曲面特别是多曲面，则更为困难.

Owen[4]对经典网格划分算法进行了总结.早期的学者提出如波前推进法[5]、Delaunay法[6]、映射法[7]等.但是这些方法在自由曲面应用有其局限性.近年来，不少学者对该问题进行了深入研究，并提出了一些针对NURBS曲面的解决方案.现有的设计多采用参数域设计并映射到曲面，通过一些技术手段改善映射变形，或采用空间距离.熊英等[8]使用空椭圆准则代替传统的空圆准则解决曲面Delaunay问题；江存等[9]应用自定义单元法改善映射变形；苏亮等[10]采用基于主应力线的波前推进法，实现网格划分；危大结等[11]、潘炜等[12]采用等面积曲面展开方法，改善从平面到曲面的映射关系.对于多重曲面，现有的方法多是在各个曲面上分别划分网格，再进行调整.然而，此法曲面交界处可能会存在明显的弯折现象，且由于曲线分属多个曲面，难以共同优化调整，难以达到设计要求.

面对复杂曲面，先将其离散为多面体的离散化表示，是一种解决方案.得到离散化表达的曲面后，裁剪边界、洞口、多曲面，以及映射关系的负面影响将不再存在.一些基于曲面距离的算法也更容易应用.另一种情形是，建筑师首先制作实物模型，再采用3-D扫描技术数字化，得到离散化表示的曲面数据，就可以不必进行复杂的曲面拟合重构步骤[13]，而直接进行网格划分.

本文基于曲面离散，提出针对复杂自由曲面的网格划分方法.经过多个复杂自由曲面的验证，表明网格划分算法的有效性、可靠性.

1 算法概览

本算法的基本流程如下：

步骤1：将曲面离散为小面片组成的离散表示；

步骤2：设置密度控制函数与曲率控制函数；

步骤3：在离散曲面上按密度要求生成初始点；

步骤4：对初始点进行密度与曲率控制的均匀化；

步骤5：形成网格；

步骤6：对网格进行拓扑与光顺优化.

1.1 离散方法

对参数曲面，如NURBS曲面，若曲面完整，则采用参数域均匀布点，形成网格，映射至曲面.对裁剪曲面，先形成均匀点，再将边界外的点舍弃，接着在参数平面上求剖分，并映射至曲面.对于多重曲面或环面，需要在曲面交界处进行缝合.

1.2 改进误差扩散算法

误差扩散算法（Error-Diffusion，简称E-D算法）由Floyd和Steinberg[14]于1976年发明.最初被应用于计算机图像处理领域的灰度图打印.一张图片上每一个像素点都有一个灰度信息，而打印时只存在有墨点和无墨点两种状态，如何分布这些墨点使得图片打印后看起来与原图一致即是误差扩散算法所要解决的问题.

离散曲面布点和灰度图生成情形类似，每一小三角形面片都具有一个点密度值，需要解决的问题是怎样分布这些点，以满足密度分布要求.

本文采用密度修正函数fρ与曲率修正函数fk对离散曲面上点密度的分布进行调整.该函数根据用户的美学需求，可任意定义.

根据前人的实践经验，误差扩散算法容易出现重复模式（artifacts）缺陷.现有的研究采用以蛇形[15]或空间填充曲线[16]替代逐行扫描，以及精心设计误差分配的范围与系数[17]，甚至采用变系数法[18]等方法规避这一缺陷.然而这些研究的对象都为像素矩阵，而本文的应用情形不存在行列概念，故难以采用现成的改进方法.

本文参考前人的思路对E-D算法进行改进.传统算法中，当某点处密度累积到达阈值t后，即布点.作者在该步骤中加入了概率因子p，布点前需要进行随机判断，p概率布点，解决了重复模式问题.概率因子p建议值受布点数k与离散曲面顶点数 n的比值k/n影响，该比值越大p的取值应越大，反之亦然.

传统算法中阈值t一般取常量0.5，即四舍五入.而根据作者的实验结果，t取常量0.5将会有起始点附近出现空洞的现象.本文针对此现象进行改进，采用变化的阈值，t将随着访问点数的增加而逐渐增大至0.5，即初始時更容易布点，以解决起始点附近的空洞现象.t的增长率取值根据比值k/n确定，该比值越大，则增长率取值应越高；反之亦然.根据作者实验，修正后的算法，起始点选择对最终结果影响不大，一般取0号点即可.

图1为该布点算法的平面算例，该算法可以有效地满足密度分布，且并未出现重复模式现象，为下一步均匀性优化提供了一个可行的初值.

1.3 基于空间距离的均匀化算法

前文生成的点阵虽然在曲面上基本满足密度分布要求，但是均匀性不良，需要应用松弛算法进行优化.本文采用粒子动力松弛算法[19].

本算法基本思想为，将各点视为具有质量mi、电量qi的粒子，质量与电量一般取1.0.任意两个粒子Vi，Vj间存在斥力：

其中，kq为用于控制斥力强度的系数，该系数对布点结果影响较大，需要试算并调整；ri，j为节点Vi与Vj的距离，此处采用欧几里得距离.本文中不考虑引力，距离大于临界距离rc后，即认为不存在相互作用力，临界距离rc一般可取为1.5 ～ 2.5倍目标杆长，亦可根据计算结果调整.其中，kq与rc可受用户自定义函数fρ与fk的调整，以控制不同区域的点密度.点将在斥力作用下运动并稳定在平衡位置附近.

图2为该松弛算法的平面算例.实验表明，该算法可以有效地均匀化点阵， 且对初始布点不敏感.

但应用于曲面时，由于动态更新两点间的曲面距离复杂度较高，本算法采用欧几里得距离，存在不足.因此仅作为初步优化，为后文步骤提供初值.

1.4 基于曲面距离的均匀化算法

Lloyd算法（又称Voronoi迭代、Voronoi松弛）由Lloyd[20]于1957年发明，1982年发表.Macqueen[21]改进后用于离散对象，也称k均值算法（k-means）.用于机器学习领域的聚类分析.

对于空间中需要聚类（cluster）的n个d维向量（x1，x2，…，xn），k均值算法的目标是将其分为k个簇 S = {S1，S2，…，Sk}，且簇内各点到其中心的距离之和最小.即：

直接求解优化函数（2）较困难.k均值算法采用迭代思想，交替进行以下两个步骤直至收敛[22]：分配步骤：根据距离函数确定xi属于哪个簇Si.更新步骤：用各个簇的中心更新其均值.

k均值算法目标与本文均匀分布离散曲面点的目标一致.Du Qiang等[23]采用该算法进行二维点的变密度分布，周炎涛等[24]将其应用于复杂三维点的聚类，取得良好的效果.

本文将该算法应用于离散化曲面.且可受函数fρ与fk的调整.实验表明该算法可以有效地实现空间曲面上点的变密度分布.

1.5 离散曲面Voronoi图计算方法

Voronoi图的基本思想是存在k个站点（site），平面分成k个区域，使得每个区域内的点到它所属区域站点的距离最近.

本文将曲面离散化后，可以应用Voronoi图的定义求解.将Voronoi图定义推广到离散曲面.离散曲面可视作一个以面片为顶点，相邻面片之间存在一条边的图.当小面片的相对尺度足够小时，即可认为该图两点间的最短路径即是对应的曲面距离.

求得Voronoi图后，即可形成网格.

1.6 拓扑优化

然而，三角形质量最优的网格并非建筑上要求的网格.在1.5节步骤之后，需要对网格进行拓扑优化.Frey等[25]认为，对于三角网格，内部节点度数为6是较优的.基本思想为遍历每一个四边形，若交换该四边形的对角线能使内部顶点度数更接近6，则进行交换.

如图3所示，左图圈中四边形箭头所指对角线对应着质量更高的两个三角形，但四个顶点皆为奇异点，其度数分别为5，7，5，7.进行边交换后，得到右图，对应的四个顶点度数都为6，规整性提高.

1.7 光顺优化

本文采用弹簧质点法[26].其基本思想是将杆件视为具有刚度ks的弹簧，节点视为具有质量的质点.其中， m一般取常量1.0即可；ks的大小影响计算速度与收敛性，需要根据经验调整；计算过程中通过调整各杆件的原长控制弹簧力的大小与方向，并影响最终结果.原长可受密度修正函数fρ与曲率修正函数fk的修正，以控制不同区域的密度.节点将在弹簧力作用下运动，并稳定在平衡位置.实验表明，最终可得到较为光顺的网格.

2 算例分析

2.1 完整单曲面算例

作者创建了一个月牙形曲面模型，如图4（a）所示.长约100 m，跨度约40 m，矢高约20 m.该模型建模时采用一个完整曲面，无裁剪及孔洞，作为基准算例.将其离散为320 000个小三角形面片.点布置以及Voronoi图如图4（c）所示，网格效果如图4（d）所示.本算法生成的网格均匀性较好，但由于曲面变化较大，为保证均匀，网格流线存在发散和收缩现象，存在少量奇异节点，但过渡较为流畅.

作为对照的映射法网格如图4（b）所示.由于该模型为单曲面，且无裁剪及孔洞，适合应用映射法，由于曲面形态、建模方式等原因，映射变形较大，生成的网格流畅性较好，但是存在明显均匀性问题.

均匀性指标如表1所示，杆长相近的情况下，本算法网格均匀性显著优于映射法.

2.2 裁剪孔洞曲面算例

作者对2.1节曲面进行了裁剪，如图5（a）所示.模型尺寸同2.1节，洞口直径为12 m，远大于网格尺寸，因此将会影响网格生成.

算例表明，本算法可有效处理裁剪和孔洞，优化后的点布置及Voronoi图如图5（c）所示，生成的网格如图5（d）所示.洞口附近通过少量奇异点过渡，无明显不流畅现象.洞口网格质量较好，整体流畅性尚可.均匀性指标如表2所示，杆长相近的情况下，本算法网格均匀性显著优于映射法.

作为对照的映射法网格如图5（b）所示.从图中，可以發现，映射法受孔洞影响较大，孔洞和裁剪边界周围网格难以处理，出现明显的不流畅、不均匀现象.且由于过于不均匀，1.7节的优化算法难以应用，优化过程中易出现网格折叠现象.

2.3 复杂曲面算例

复杂的建筑模型，采用多重曲面建模可以简化建模的难度，提高建模质量，提高建模效率[3].

然而对于复杂曲面，网格划分困难.现有的方法多是在各个曲面上分别划分网格，再对多个网格进行缝合，最后进行调整和优化.由于一条曲线分属两个曲面，难以共同调整优化.流畅性与均匀性都难以保证，曲面交界处往往存在明显弯折.

图6是一个船形曲面的局部.由图6（a）可见，由于模型形式复杂，建筑师建模时采用了两个曲面组合来表达.如图6（a）的曲面结构线所示，上部曲面是裁剪曲面，且一端存在急剧收缩的现象，这都给传统网格划分方法造成了很大困难.

对于该算例，采用1.1节中所述的方法，先分别将两个曲面离散为约150 000、210 000个小三角形面片，再对两个网格进行缝合，最终形成共360k个面片.

实验证明本算法可以处理裁剪曲面、多重曲面等形式复杂的模型.如图6（b）所示，生成的网格均匀性较好，杆长均值为2.313 m，杆长方差仅为0.044 m2.但由于模型变化较大，为保证上下网格密度相同，则网格流线不可避免存在发散和收缩现象，个别内部节点所连杆件数量不为6，优化后的网格整体流畅性尚可.

作者对上下两个曲面分别采用映射法划分网格，如图6（c）所示，下部曲面是完整曲面，映射法效果较好；而上部曲面由于是裁剪曲面，且一端急剧收缩，映射法仅保证了拓扑流畅，实际空间效果较差.由图6（c）可见，映射法存在严重的映射变形问题.对上下两曲面分别应用1.7节的均匀化后，如图6（d）所示，上部曲面的映射变形问题得到一定改善.但整体仍存在下部杆件过于密集，上部杆件过长的问题.而且网格流线在两曲面分界线处存在明显的弯折现象，而流线分属两个曲面，难以调整.而图6（b）所示的本算法网格并没有受到曲面分界线的影响.

作者根据标高设置了分段密度控制函数，上部疏下部密，并得到了图6（e）和（f）的最终结果，网格流畅，且满足密度分布，仅存在少量奇异点.实验表明本算法可以根据用户的美学需求进行变密度网格生成.

3 结 论

针对复杂曲面，本文提出了一种以密度控制为目标的网格划分方法.首先，对复杂曲面离散并缝合，一体化为由大量三角面片组成的离散曲面，可消除裁剪边界、孔洞、多曲面拼接、闭合等负面影响.其次，采用改进的误差扩散算法进行初始点的分布，该算法布点结果虽不够均匀，但满足密度分布，为之后的步骤提供了可靠的初值.然后，采用基于空间距离的粒子动力松弛算法，对点云进行初步均匀化，应用基于曲面距离的k均值算法对点云进行再次均匀化，得到曲面上满足密度要求的均匀点阵.接着求离散曲面Voronoi图，由此得到曲面Delaunay图.之后，进行拓扑优化，提高规整性.最后，采用弹簧质点松弛算法对该三角网格进行光顺优化，得到均匀流畅的最终网格.

算例表明，本文算法可以处理裁剪、开洞、变化剧烈、环面、多曲面拼接等复杂情况，适应性好.最终形成的网格能满足密度要求，但网格流线可能存在发散和收缩情况，存在内部节点连有5根或7根杆件的现象.

奇异节点不可避免，本文算法的缺陷在于，奇异节点的位置并不能控制，本文采用的拓扑优化算法并无全局优化能力，可能需要在形成网格后交互地修改拓扑，消除奇异点的负面影响.对此仍需要进行更多的研究探索，在形成网格前就交互确定奇异点，或在算法中全局地优化奇异点的数量与位置，以获得更规整的网格.

参考文献

[1] PIEGLL L，TILLER W. The NURBS Book[M] .2nd ed. Berlin Heidelberg： Springer-Verlag，1997.

[2] 王磊，杨彬，张其林. 非均匀有理B样条曲面形状优化方法[J]. 湖南大学学报（自然科学版）， 2012，39（7）：14—19.

WANG L， YANG B，ZHANG Q L. Shape optimization of non-uniform rational B-spline surface[J].Journal of Hunan University（Natural Sciences），2012，39（7）：14—19. （In Chinese）

[3] SEDERBERG T W， ZHENG J， BAKENOV A， et al. T-splines and T-NURCCs[J]. ACM Transactions on Graphics，2003，22（3）： 477 —484.

[4] OWEN S. A survey of unstructured mesh generation[C]// International Meshing Roundtable. Dearborn， USA： IMR， 1998： 239—267.

[5] LHNER R. Progress in grid generation via the advancing front technique[J]. Engineering with Computers，1996，12（3）：186—210.

[6] BOROUCHAKI H，HECHT F，SALTEL E，et al.Reasonably efficient Delaunay based mesh generator in 3 dimensions[C]// International Meshing Roundtable. South Lake Tahoe， USA： IMR， 1999： 3—14.

[7] COOK W A，OAKES W R. Mapping method for generating three-dimensional meshes： past and present[C]// International Computer Engineering Conference.San Diego，USA：ASME，1982.

[8] 熊英，胡于進，赵建军. 基于映射法和Delaunay方法的曲面三角网格划分算法[J]. 计算机辅助设计与图形学学报， 2002， 14（1）： 56—60.

XIONG Y， HU Y J， ZHAO J J. An algorithm of surface triangulation based on mapping and Delaunay method[J]. Journal of Computer-Aided Design & Computer Graphics，2002，14（1）： 56—60.（In Chinese）

[9] 江存，高博青. 基于自定义单元法的自由曲面网格划分研究[J]. 建筑结构， 2015， 45（5）： 44—48.

JIANG C，GAO B Q. Research on free-form surface meshing based on self-defined element method[J]. Building Structure， 2015， 45（5）： 44—48. （In Chinese）

[10] SU L， ZHU S， XIAO N， 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.

[11] 危大结，舒赣平. 自由曲面网格的划分与优化方法[J]. 建筑结构， 2013， 43（19）： 48—53.

WEI D J， SHU G P. Mesh generation and optimization method for free-form surface grid[J]. Building Structure， 2013， 43（19）： 48—53. （In Chinese）

[12] 潘炜，吴慧，李铁瑞，等.基于曲面展开的自由曲面网格划分[J]. 浙江大学学报（工学版）， 2016， 50（10）： 1973—1979.

PAN W， WU H， LI T R，et al. Grid generation on free-form surface based on surface flattening[J]. Journal of Zhejiang University（Engineering Science），2016，50（10）：1973—1979.（In Chinese）

[13] SCHALL O， SAMOZINO M， FALCIDIENO B， et al. Surface from scattered points： a brief survey of recent developments[C]// Proccedings of the first international workshop towards Semantic Virtual Environments. Villars， Switzerland： SVE， 2005： 138—147.

[14] FLOYD R W， STEINBERG L. Adaptive algorithm for spatial greyscale[C]// Proceeding SID. New York， USA： SID， 1976：75—77.

[15] ULICHNEY R. Digital halftoning[M].Cambrige，MA：MIT Press， 1987：12—13.

[16] VELHO L， GOMES J D M. Digital halftoning with space filling curves[J]. ACM SIGGRAPH Computer Graphics， 1991， 25（4）：81—90.

[17] FAN Z. Set of easily implementable coefficients in error diffusion with reduced worm artifacts[C]// Proceeding SPIE. San Jose， USA： SPIE，1996：222—225.

[18] OSTROMOUKHOV V. A simple and efficient error-diffusion algorithm[C]// Proceeding SIGGRAPH. Los Angeles， USA： ACM， 2001： 567—572.

[19] BRONSON J R， LEVINE J A， WHITAKER R T. Particle systems for adaptive， isotropic meshing of CAD models[J]. Engineering with Computers， 2012， 28（4）： 331—344.

[20] LLOYD S. Least squares quantization in PCM[J]. IEEE Transactions on Information Theory， 1982， 28（2）： 129—137.

[21] MACQUEEN J. Some METHODS for classification and analysis of multi variate observations[C]//Proceedings of Berkeley Symposium

on Mathematical Statistics and Probability. Oakland， USA： University of California Press， 1967：281—297.

[22] MACKAY D.Information theory，inference，and learning algorithms[M]. Cambridge， UK： Cambridge University Press，2003： 640.

[23] DU Q， FABER V， GUNZBURGER M. Centroidal Voronoi tessellations： applications and algorithms[J]. Siam Review， 1999， 41（4）： 637—676.

[24] 周炎濤，吴正国，易兴东. 基于网格带有参考参数的扩展聚类算法[J]. 湖南大学学报（自然科学版），2009，36（2）：48—52.

ZHOU Y T， WU Z G， YI X D. Extended grid-based clustering algorithm with referential parameters[J]. Journal of Hunan University（Natural Sciences），2009，36（2）：48—52.（In Chinese）

[25] FREY W H， FIELD D A. Mesh relaxation： a new technique for improving triangulations[J]. International Journal for Numerical Methods in Engineering， 1991， 31（6）： 1121—1133.

[26] 李基拓，陆国栋. 基于边折叠和质点-弹簧模型的网格简化优化算法[J].计算机辅助设计与图形学学报，2006，18（3）：426—432.

LI J T， LU G D. Mesh simplification and optimization with edge collapse and mass-spring model[J]. Journal of Computer-Aided Design & Computer Graphics，2006，18（3）：426—432.（In Chinese）