基于四边形网格参数细分的平面与自由曲面求交算法*

李慧莹，陈良骥

（郑州航空工业管理学院 机电工程学院，河南 郑州 450015）

0 引 言

一直以来，平面与曲面间求解交线的问题都是计算机辅助设计与制造（computer-aided design/manufac⁃turing，CAD/CAM）编程领域中最普遍的工程问题。在曲面造型与裁剪、加工刀具轨迹计算、加工几何图形验证等实际应用中，常常需要对平面与曲面进行求交运算［1-3］。国内外在求交计算方面做了大量的研究，大致包括牛顿迭代法［4］、曲面离散法［5-6］、区间算术法［7］、光线跟踪法［8-9］以及近年来出现的等值线法［10］、三角网格方法［11-12］等。牛顿迭代法通常是在给定初值后将求交问题转化为求线性方程组或常微分方程解的问题，但收敛与否以及收敛的速度与迭代初值的选取有很大关系。曲面离散算法是用小平面片近似逼近曲面的一种几何求交算法，有较高的可靠性，但缺点是运算量较大、效率低、精度低。区间算术法的提出则只是为了解决直线与隐曲面求交的问题。

鉴于以上各种求交算法各自的缺点，本研究将首先提出一种四边形参数曲面片模型，在此基础上可将平面与自由曲面的求交问题简化为直线段与平面的位置关系问题。与现行常见的求交方法相比，本研究算法的特点是计算稳定可靠、精度高而且具有较强的普适性。

1 自由曲面的NURBS表示

一张笛卡尔空间自由曲面可表示为如下非均匀有理B样条形式［13-14］:

式中:Pi,j—三维控制点；wi,j—Pi,j对应权重；（n+1）,（m+1）—u向和v向控制点的数目；Ni,p(u)，Nj,q(v)—沿u向的p次和v向的q次B样条基函数。

由自由曲面的NURBS定义可知，一张自由曲面对应于参数平面内一正方形区域｛（u,v）|0≤u≤1，0≤v≤1｝，笛卡尔空间与参数空间的映射关系如图1所示。

图1 笛卡尔空间与参数空间的映射关系

2 空间点、线与平面的关系

2.1 空间点V与平面的关系

已知平面单位法矢为｛A，B，C｝、平面常数为D，平面与曲面相交如图2所示。平面为所有点V=（x，y，z）的集合｛V|F（V）=F（x，y，z）=Ax+By+Cz+D=0｝。该平面将空间分为3个部分，即｛V|F（V）＜0｝、｛V|F（V）=0｝和｛V|F（V）＞0｝。因此，空间任一点V与平面有如下关系:若点V在平面上，则F（V）=0；反之，若V使得F（V）≠0，则V定不在平面上。

图2 平面与曲面相交

2.2 空间线段V1V2与平面的关系

空间线段V1V2与平面的关系可由线段两个端点V1和V2分别与平面的关系来确定。具体位置关系可以有以下4种情形:①若V1和V2分布在平面相同一侧，则线段V1V2与平面不相交，此时F(V1)、F(V2)皆非零且具有相同的符号；②若V1和V2分布在平面相异一侧，则线段V1V2与平面交于一点，此时F(V1)、F(V2)皆非零且具有相异的符号；③若V1（或V2）在平面上而V2（或V1）不在平面上，则线段与平面交于V1（或V2）；④若V1和V2都在平面上，则线段位于平面上。

将以上线段与平面之间的关系归纳如下:

3 四边形参数面片

如图1所示，位于笛卡尔空间的曲面的每个空间四边形都对应于参数平面内的一个小长方形。四边形参数面片模型如图3所示，本研究将每个长方形分别按照一定顺序标记出4个角点（1、2、3、4）和4条边（Ⅰ、Ⅱ、Ⅲ、Ⅳ），并定义这样的长方形为四边形参数面片。参数平面内，平面与曲面交线的参数变化曲线v=f(u)，如图4所示。本研究用前一小节中介绍的两种关系可以判定出各个空间四边形每个边与平面的关系，进而可以判断出四边形参数面片与v=f(u)的位置关系。曲面所有四边形参数面片与v=f(u)间的位置关系分别会有如图3所示的a、b、c、d、e、f、g等7种状态。

图3 四边形参数面片模型

本研究对这7种状态分别做如下定义:①空间四边形的每个边均在平面的相同一侧为a态；②空间四边形的仅有一个顶点在平面上为b态；③空间四边形的两个相邻边与平面相交为c态；④空间四边形的一个顶点在平面上、一条边与平面相交为d态；⑤空间四边形的两个对边与平面相交为e态；⑥空间四边形的两个对角顶点在平面上为f态；⑦空间四边形的某边位于平面上时为g态。

4 交线链表生成

对图1中的每个小四边形，顺次计算各角点坐标，应用前面所介绍的两种关系和几种状态的定义方法，可以得到如图4所示的一张状态表。

图4 参数变化曲线与状态表

求交线的过程为:①建立一张空的四边形面片结构的链表并对照状态表将所有处于非a态和非b态的四边形按一定顺序（如从下到上顺次取出每一行、对取出的每行按从左至右顺次取出每个四边形）存入该链表中，得到一张初始链表；②对初始链表中所有处于非f态以及非g态的各结点（小四边形）按同样的方法再进行细分，并将细分后得到的所有非a态和非b态的四边形按与前面相同的顺序插入该链表中，得到一张新链表；③重复第②步的操作直到链表中所有的结点状态都为f态和g态为止。最后得到的四边形链表能很好地逼近交线且交线必然位于平面内，称该链表为交线链表。以图4表示的状态表来说明该计算过程，求交线链表的过程如图5所示。

图5 求交线链表的过程

将交线链表中各四边形位于平面内的角点作为型值点，用NURBS样条进行拟合，即可得到平面与自由曲面的交线。

5 实例计算

本节将以一张双三次NURBS曲面为例对所提出的方法进行验证。该曲面由16个控制点形成，各控制点权重取为1。为说明问题，先找出曲面上已知的3个点P、P1和P2，通过这三点确定一个平面，求这个平面与自由曲面的交线。

取P=S(0.5,0.5)=[30,67.5,45]T，P1=S(0.15,0.35)=[9,34.528 3,31.5]T，P2=S(0.85,0.65)=[51,34.321 7,58.5]T，经计算，由点P、P1和P2所确定的平面方程各系数分别为:A＝0.540 757 6，B＝0.0，C＝-0.841 178 5，D＝21.630 304。利用VC++和OpenGL在计算机上编程实现了本研究所提出的求交算法，从初始链表开始，对链表中每个四边形分别进行10×10细分，共经过6层细分可以得到图5中的交线链表，对其进行曲线拟合，求得的平面与自由曲面的交线如图6所示。

图6 平面与自由曲面的交线

以下通过求直线与交线的交点来进一步验证算法的精确性。理论上，P1、P2应该在交线上，但实际计算平面内直线P1P2与交线的交点时采用将曲线变折线段的方式进行，因而只能得到它们的近似值和，求得的和分别为:

与P1和P2的距离误差分别为:

6 结束语

本研究所提出的平面与自由曲面求交线方法的思路是:根据笛卡尔空间与参数平面之间的映射关系，将空间四边形与平面的位置关系转化为参数平面内四边形与交线参数变化曲线间的关系，进而求得交线。最后，本研究结合实例进行了计算，研究结果表明，求得的交线能够与理论交线精确符合，完全可以在几何造型、曲面裁剪以及五轴编程刀位计算等诸多方面得到实际应用。

（References）:

［1］张和明，柯映林，程耀东.参数曲面与平面求交的一种新方法［J］.工程图学学报，1995（2）:31-37.

［2］张和明，张玉云，熊光楞，等.参数曲面与平面的精确求交及其应用［J］.机械工程学报，1997，33（5）:31-36.

［3］吕晓倩，赵玉刚，周海安.空间曲面与平面交线的一种插补算法［J］.组合机床与自动化加工技术，2008（3）:13-15.

［4］ROY U，DASARI R.Implementation of polygonal algorithm for surface-surface intersections［J］.Computers Industry Engineering，1998，34（2）:399-412.

［5］马 翔，周儒荣.自由曲面与平面的一种分割、跟踪求交方法［J］.南京航空航天大学学报，1994，26（1）:75-79.

［6］郑立垠，张 丽，张云鹏.细分曲面求交交线计算方法的研究［J］.微计算机应用，2008，29（1）:78-81.

［7］余正生，李启炎，肖少拥，等.一种直线与隐式曲面求交的方法［J］.工程图学学报，2000，21（3）:20-23.

［8］ROTH S D.Ray casting for modeling solid［J］.Computer Graphics&Image Process，1982，18（2）:109-144.

［9］余正生，楚广琳.一种跟踪隐式曲面交线的算法［J］.计算机应用研究，2008，25（7）:2235-2237.

［10］宋宏勋，韩 毅，吴初娜.一种基于等值线法的NURBS曲面与平面的求交算法［J］.数字技术与应用，2011（7）:103-105.

［11］孙殿柱，孙永伟，田中朝，等.三角网格曲面模型快速求交算法［J］.北京工业大学学报，2012，38（8）:1121-1124.

［12］孙殿柱，康新才，李延瑞，等.三角Bézier曲面快速求交算法［J］.机械工程学报，2011，47（3）:89-94.

［13］PIEGL L.On NURBS:a survey［J］.IEEE Computer Graphics&Application，1991，11（1）:55-71.

［14］隆 强，谢延敏，杨 川.基于Foleg参数法反算三次NURBS曲线的算法研究［J］.机械，2012，39（7）:5-8，40.