基于仿射算术和区间运算的直线与NURBS曲线/曲面求交

池宝涛 张见明 鞠传明

1.湖南大学机械与运载工程学院，长沙，4100822.湖南大学汽车车身先进设计制造国家重点实验室，长沙，410082

0 引言

直线与参数空间NURBS(non-uniform rational B-splines)曲线、直线与NURBS曲面求交问题是计算机辅助几何设计(computer aided geometric design, CAGD)的一个基本问题，也是其中非常重要和复杂的问题之一。曲线、曲面求交算法是CAD/CAE系统中的一个核心算法，广泛应用于实体造型、曲面裁剪等问题中。区间分析[1]是数值计算领域一个新的分支，它是以“区间数”为研究对象的一种新的数学分析方法。直线与参数空间NURBS曲线、直线与NURBS曲面[2]求交的本质问题就是对非线性方程或非线性方程组进行求解，求交算法的稳定性、精度和效率直接影响算法的实用性[3]。

当前CAD几何造型系统中，直线与参数空间曲线或曲面的求交方法主要归纳为两大类：一类为离散分割近似方法，其特点是易于实现，计算精度较低，存储量较大，效率不高；另一类为数值计算方法，其优点是计算精度高，几何数据存储量较少，缺点是一般只适用于给定曲面方程的曲面，对于较为复杂的方程组几乎无法求解。按照方法来分，离散分割的近似方法主要包括离散细分法[4]、网格划分法[5]等。数值计算方法主要包括代数解析法、点迭代法、跟踪求交法[6]、同伦迭代法[7]等。

本文提出了改进的基于仿射算术[8]和区间运算的直线与参数空间NURBS曲线、直线与NURBS曲面求交的算法。与传统的离散分割近似方法相比，该方法由于采用区间运算，故细分的阈值可以任意给定，且无需考虑细分的曲线段或曲面片是否近似线性，从而细分层数得到很大程度的减少；与传统的点迭代法相比，该方法由于采用了区间运算，迭代过程不需给定合适的迭代初始值，具有更好的灵活性；与传统的区间迭代法相比，该方法放宽了对初始区间的要求，采用基于线曲率和面曲率的子域分解方法，可以快速筛选预迭代区间，提高迭代效率，另外，通过运用仿射算术，考虑计算过程中数据的相关性，可以大大减少迭代次数，有效弥补了区间算法的局限性，在相同的容差要求下，可提高迭代求交的效率。

1 算法基本理论

1.1 直线、NURBS曲线与NURBS曲面的定义

(1)直线。设直线的起始点为O，直线的单位方向矢量为D，则直线L的参数方程可表示为

L(t)=O+tDt∈(-∞,+∞)

(1)

(2)NURBS曲线。一条p(p≥0)次NURBS曲线定义为

(2)

式中，Pi为控制点；wi为权因子；Ni,p(t)为定义在非均匀节点矢量T上的p次样条基函数。

基于De Boor-Cox递归公式，第i个p次样条基函数Ni,p(t)定义如下：

(3)

(3)NURBS曲面。一张在u方向p次和v方向q次的NURBS曲面或T-Spline曲面定义为

(4)

0≤u,v≤1

式中，Pi,j为控制点；wi,j为权因子；Ni,p(u)、Nj,q(v)分别为定义在非均匀节点矢量U和V上的非有理样条基函数。

1.2 区间分析

区间分析的主要思想是在整个数值计算过程以区间为运算对象，用区间代替取值不精确的浮点数进行计算，从而自动地记录计算机浮点运算中所产生的截尾误差和舍入误差，以区间形式输出计算结果，并且结果区间能够保证包含数据运算的真实结果，从而保证计算的可靠性。

1.2.1Newton算子

Newton算子定义如下：

(5)

其中，X为区间向量。

假定f:Rn→Rn，F为f的区间扩展[9]，F′为F的一阶导数，则新区间X(k+1)及Newton算子区间迭代序列[10]为

(6)

k=0,1,2,…

由区间除法的定义可知，区间Newton法的区间划分判断是自适应的，这是区间Newton法的最大优势。

Newton算子具有如下性质。

假定映像f:X⊂Rn→Rn,求解非线性方程组f(x)=0,x∈X⊂Rn。

(1)若x*为f(x)=0的解，且x*∈X，则x*∈N(X)。

(2)若X∩N(X)=∅，则f(x)在X中无解。

(3)若N(X)⊂X，则f(x)在X中必有解存在，当且仅当W(N(X))⊂W(X)时，f(x)在X中存在唯一解。

1.2.2Krawczyk算子

Krawczyk算子定义如下：

K(X)=y-Yf(y)+[I-YF′(X)](X-y)

(7)

Krawczyk算子[11]区间迭代规则如下：

(8)

k=0,1,2,…

其中，y∈X，I为n阶单位矩阵，Y为n阶非奇异矩阵。特别地，若在式(8)中可取y=m(X),Y=[m(F′(X))]-1。

Krawczyk算子具有如下性质。

假定映像f:X⊂Rn→Rn,求解非线性方程组f(x)=0,x∈X⊂Rn。

(1)若x*为f(x)=0的解，且x*∈X，则x*∈K(X)。

(2)若X∩K(X)=∅，则f(x)在X中无解。

(3)若K(X)⊂X，则f(x)在X中必有解存在，当且仅当W(K(X))⊂W(X)，f(x)在X中存在唯一解。

1.2.3区间迭代几何意义

对一任意函数f(x)，设初始迭代区间为[a,b]，图1给出了f(x)在不同迭代区间存在唯一解或多个解的情况示意图。

(a)迭代区间有唯一解(b)迭代区间有多个解图1 区间迭代的几何意义Fig.1 Geometrical interpretation of the univariate interval Newton method

1.3 仿射算术

区间算术的缺点在于其保守性，经区间算术运算得到的结果区间往往被很大程度地放大，放大倍数取决于区间运算计算链的长度，理想的实际区间为结果区间的子区间，而且结果区间往往比理想的实际区间大得多，这一问题在长计算链的区间运算中尤为突出，因此会产生“误差爆炸”现象，而这样的长计算链在实际计算中经常出现[12]。仿射算术为以解决区间运算中相关性问题为目的而建立的有效的数值计算模型，该模型可以解决区间算术过于保守的问题。与区间算术一样，仿射算术能够自动记录浮点数的截尾和舍入误差，此外，仿射算术还能保持输入和计算中各个不确定量之间的依赖关系，并自动应用于原始运算中。基于此优点，仿射算术能得到比区间算术区间范围更小更合理的区间，特别是在实际计算的长计算链中优势更加明显。

在仿射算术中，一个不确定量x的仿射形式为

给定两个仿射形式进行运算：

2 直线与NURBS曲线快速求交算法

2.1 构建参数空间直线与NURBS曲线求交目标函数

在二维参数空间中，直线与NURBS曲线交点(u,v)应满足

(9)

和

(10)

其中，Pui、Pvi分别为NURBS曲线控制点的坐标，(u0,v0)及(u1,v1)为直线上两个给定点的参数坐标。在二维参数空间求解直线与NURBS曲线交点的目标函数如下。

(1) 当u1≠u0且v1≠v0时

(11)

(2)当u1≠u0且v1=v0时

(12)

(3)当u1=u0且v1≠v0时

(13)

2.2 基于仿射算术的Newton算子运算

对于Newton区间迭代法，利用区间除法的性质，对迭代区间进行自适应划分，这是Newton区间迭代法的优势，而在样条基函数Ni,p(t)的区间运算中需要多次调用递归运算，属于长计算链的区间运算，在计算过程中如果迭代区间多次累计放大，会导致迭代次数呈非线性增长，区间迭代效率降低或出现不收敛的情况；此外，对不包含真实解的无效区间进行判断是影响区间迭代效率的另一主要原因，如何快速准确地定位包含真实解的有效区间直接影响算法的效率。

(a)NURBS曲线 (b)细分步骤(1)

(c)细分步骤(2) (d)构建包围盒

(e)直线与参数空间NURBS曲线求交图2 曲线离散和直线与包围盒求交Fig.2 Curve discretization and intersection of line and boxes

算法具体步骤描述如下：

(2)对离散的曲线段创建包围盒(图2d、图2e)。如图2d、图2e所示，判断直线与包围盒是否相交，若不相交则去除这些包围盒对应曲线段的区间，避免对这些不包含交点的区间进行无效迭代，同时筛选出与直线相交的包围盒作为下一步进行区间迭代的预迭代区间。

(3)对筛选的区间进行区间Newton迭代。若存在多个解，迭代过程中迭代区间会进行自适应划分，进一步对细分区间进行判断；若无解，则停止迭代,终止计算；若存在唯一解，则在这一区间内给定一个初值，通过牛顿迭代法或拟牛顿迭代法求出真解。

3 直线与NURBS曲面快速求交算法

3.1 构建直线与NURBS曲面求交目标函数

(14)

0≤u,v≤1

图3 绝对坐标系下的直线与NURBS曲面Fig.3 Line and NURBS surface in the absolute coordinate system

图4 局部坐标系下的直线与NURBS曲面Fig.4 Line and NURBS surface in the local coordinate system

(15)

3.2 基于仿射算术的Krawczyk算子运算

直线与NURBS曲面求交中同样存在直线与NURBS曲线求交遇到的问题，另外，在Krawczyk算子计算中也需要考虑相关性问题。

算法具体步骤描述如下：

(1)对NURBS曲面进行区域划分。在二维参数空间考虑黎曼度量和面曲率生成二叉树网格，对于裁剪曲面无需将二叉树网格与边界进行求交处理。为简化计算，采用弦长比作为衡量曲面曲率的标准。如图5所示，曲线ab为细分单元内一条等参线所对应曲线的剖视图，根据面曲率对目标面进行细分，然后判断各子域是否满足给定细分阈值，若满足则将此子域细分，否则停止细分，直到单元内的面曲率满足预给定阈值，然后将细分结果映射到三维空间。

(2)对目标面的剖分区域创建体包围盒。将直线与体包围盒进行相交判断，筛选出与直线相交的包围盒作为下一步进行区间迭代的预迭代区间。

(3)对筛选的区间进行区间迭代。若在u,v的某一区域有多个解，对U、V采用四叉树方法进行细分，分成四个小区域，对各个子区域计算Krawczyk算子，并对这些子区域进行判断。若无解则终止计算；若存在唯一解，则在这一区域内给定一初值，通过牛顿迭代法或拟牛顿迭代法迭代求出真解。

(4)在计算过程中对Krawczyk算子进行改进：

X(0)=X=M+Rε= (m1+r1ε1,m2+r2ε2,…,mn+rnεn)′
X(k+1)=X(k)∩K(X(k))

将X(k+1)转化为仿射形式，应用仿射算术计算G(X(k))=Y(k)F′(X(k)),将仿射矩阵G(X(k))转化为新的区间矩阵H(X(k))。代入Krawczyk算子公式，应用区间算术计算

K(X(k))=m(X(k))-Y(k)f(m(X(k)))+
[I-H(X(k))](Rε)(k)

(a)NURBS曲线 (b)细分步骤(1)

(c)细分步骤(2) (d)细分步骤(3)

(e)NURBS曲面细分实例 (f)曲面细分俯视图图5 基于曲率的NURBS曲面细分Fig.5 NURBS surface subdivision based on curvature

4 算法实例验证

本文提出一种改进的基于仿射算术和区间运算的直线与参数空间NURBS曲线及直线与NURBS曲面快速求交算法，分别验证了三种情况下的直线与5次NURBS参数空间曲线求交，并将本文改进的方法与传统的区间迭代求交方法进行比较，如图6及表1所示。其中，M1为传统区间迭代求交方法，M2为本文改进的区间迭代求交方法。为验证本文算法的可行性与一般性，将本文算法与传统的区间迭代求交算法进行对比，验证了直线与控制点为7×3的双三次NURBS曲面、直线与控制点为12×3的双三次NURBS曲面，在相应容差要求下的计算结果和时间的对比情况，如图7及表2所示。将本文改进方法与直线与复杂曲面求交或实际工程算例进行对比，在相应容差要求下的计算结果和时间对比情况如图8、图9及表3、表4所示。其中，M3为文献[13]中直线与复杂曲面求交的方法，M4为文献[14]中直线与弯管曲面求交的方法。图8为直线与控制点为13×3的u方向四次、v方向三次的复杂NURBS曲面求交算例[14]。图9为直线与控制点为14×43的u方向四次、v方向三次的复杂弯管曲面求交的工程算例[15]。计算结果表明采用本文方法计算多个交点的求解时间与其他算例中的平均时间相当或优于其他算例求解时间，证明了本文算法的可行性。

(a)直线u=u0与参数空间NURBS曲线求交

(b)直线v=v0与参数空间NURBS曲线求交

(c)一般直线与参数空间NURBS曲线求交图6 直线与参数空间NURBS曲线求交Fig.6 Intersection between the line and the NURBS curve

直线误差(mm) 计算时间(ms)M1M211.208×10-1427.517.126.758×10-1448.226.437.933×10-1452.628.3

(a)直线与NURBS曲面求交存在一个交点情况

(b)直线与NURBS曲面求交存在两个交点情况

(c)直线与NURBS曲面求交存在多个交点情况图7 直线与NURBS曲面求交Fig.7 Intersection between the line and the NURBS surface

直线误差(mm) 计算时间(ms)M1M213.238×10-1457.837.126.712×10-1469.745.934.693×10-1371.449.6

图8 直线与复杂NURBS曲面求交Fig.8 Intersection between the line and the complex NURBS surface

(a)直线与弯管曲面求交存在两个交点情况

(b)直线与弯管曲面求交存在多个交点情况图9 直线与弯管曲面求交Fig.9 Intersection between the line and the NURBS surface of canal

误差(mm) 计算时间(ms) M1M2M37.114×10-1394.361.7平均时间65～70

表4 不同直线与弯管曲面求交时的误差与计算时间比较(迭代容差为10-9 mm)

5 结论

(1)本文将仿射算术与区间运算相结合，一定程度上解决了传统区间运算的“保守性”，减少了迭代过程中不必要的求交判断，并放宽了对初始区间的要求。

(2)将基于边曲率或面曲率的子域分解方法应用到求交算法中，快速定位预迭代区间，在相同的容差要求下，比传统区间迭代求交算法更快，提高了迭代求交算法的效率。

(3)通过计算区间算子判断给定直线与参数空间NURBS曲线或直线与NURBS曲面有无交点和存在交点时的交点数目，可保证求解交点精度，且不会遗漏交点。