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+tD t∈(-∞,+∞)
(1)
(2)NURBS曲线。一条p(p≥0)次NURBS曲线定义为
(2)
式中，Pi为控制点；wi为权因子；Ni,p(t)为定义在非均匀节点矢量T上的p次样条基函数。
基于De Boor-Cox递归公式，第i个p次样条基函数Ni,p(t)定义如下：
(3)
由NURBS曲线的规范性可知，对所有的
(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 区间分析
区间分析的主要思想是在整个数值计算过程以区间为运算对象，用区间代替取值不精确的浮点数进行计算，从而自动地记录计算机浮点运算中所产生的截尾误差和舍入误差，以区间形式输出计算结果，并且结果区间能够保证包含数据运算的真实结果，从而保证计算的可靠性。
区间与区间运算定义如下：实数集R上的一个连续子集
称为有界闭区间。所有R上有界闭区间的集合记作I(R)。区间X的绝对值、宽度、中点分别定义为
1.2.1 Newton算子
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.2 Krawczyk算子
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)在不同迭代区间存在唯一解或多个解的情况示意图。
(1)如图1a所示，f(x)在该区间满足单调性，由牛顿迭代序列
及F′(X)=[f′(a),f′(b)],F′(X)为正区间，
及由区间除法定义N(Xk)=[A,B], 可得新的迭代区间为Xk+1=N(Xk)∩Xk=[A,B],区间宽度W(Xk+1)=B-A<W(Xk)=b-a,故由区间迭代解的存在性定理也可判断出f(x)在该区间单调，存在唯一解，再利用牛顿迭代法或拟牛顿迭代法即可快速迭代出真解
(2)如图1b所示，f(x)在该区间不满足单调性，由牛顿迭代序列
⊃
及区间除法定义N(Xk)⊂(-∞,A]∪[B,+∞),可得新的迭代区间为Xk+1=N(Xk)∩Xk⊂[a,A]∪[B,b],区间宽度W(Xk+1)=b-a=W(Xk),故由区间迭代解的存在性定理也可判断出f(x)在该区间非单调，可能存在多个解也可能无解，下次迭代需要进一步分别对两个区间进行判断，直到找到在某区间内f(x)满足单调性，f(x)在N(X)中存在唯一解或无解，则停止迭代。
(a)迭代区间有唯一解(b)迭代区间有多个解
图1 区间迭代的几何意义
Fig.1 Geometrical interpretation of the univariate interval Newton method
1.3 仿射算术
区间算术的缺点在于其保守性，经区间算术运算得到的结果区间往往被很大程度地放大，放大倍数取决于区间运算计算链的长度，理想的实际区间为结果区间的子区间，而且结果区间往往比理想的实际区间大得多，这一问题在长计算链的区间运算中尤为突出，因此会产生“误差爆炸”现象，而这样的长计算链在实际计算中经常出现[12]。仿射算术为以解决区间运算中相关性问题为目的而建立的有效的数值计算模型，该模型可以解决区间算术过于保守的问题。与区间算术一样，仿射算术能够自动记录浮点数的截尾和舍入误差，此外，仿射算术还能保持输入和计算中各个不确定量之间的依赖关系，并自动应用于原始运算中。基于此优点，仿射算术能得到比区间算术区间范围更小更合理的区间，特别是在实际计算的长计算链中优势更加明显。
在仿射算术中，一个不确定量x的仿射形式为
其中,εi为噪声元，假定εi落在[-1,1]范围内,对应的系数xi为实数，决定了噪声元εi的大小和符号，x0为仿射形式
的中心值。εi表示一个对x的总的不确定性起一定作用的独立误差源, 如输入数据的不确定性、公式的截断误差、运算中的舍入误差以及其他不确定性因素。不同的仿射形式如果有相同的噪声元εi，则说明这些不确定仿射形式对应的不确定量之间具有一定的相关性。
区间和仿射形式之间可以相互转换：对于一给定区间
其对应的仿射形式为
反过来，给定一仿射形式
其对应的区间为
给定两个仿射形式进行运算：
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)的区间运算中需要多次调用递归运算，属于长计算链的区间运算，在计算过程中如果迭代区间多次累计放大，会导致迭代次数呈非线性增长，区间迭代效率降低或出现不收敛的情况；此外，对不包含真实解的无效区间进行判断是影响区间迭代效率的另一主要原因，如何快速准确地定位包含真实解的有效区间直接影响算法的效率。
针对长计算链的区间运算导致结果区间过度放大的问题，笔者将仿射算术运用到样条基函数Ni,p(t)的区间运算中，在每次迭代时将当前步要迭代的区间分量变为仿射形式
代入样条基函数Ni,p(t)中进行计算，在一定程度上可以降低结果区间的放大倍数，减少迭代的次数。针对不包含真实解的无效区间的加速判断，本文采用曲线离散和直线与包围盒求交的方法(图2)，快速准确定位包含真实解的有效区间，很大程度上减少了对不包含真实解的无效区间的判断，提高了区间迭代算法的效率。
(a)NURBS曲线 (b)细分步骤(1)
(c)细分步骤(2) (d)构建包围盒
(e)直线与参数空间NURBS曲线求交
图2 曲线离散和直线与包围盒求交
Fig.2 Curve discretization and intersection of line and boxes
算法具体步骤描述如下：
(1) 对NURBS曲线进行边离散。采用基于曲率的二分法对目标曲线进行离散，为简化计算采用弦长比作为衡量曲线曲率的标准，弦长比
如图2a所示。若R不满足预指定的边离散曲率参数，将该弧
平分为弧
和弧
如图2b、图2c所示，再分别对弧和弧
计算弦长比R，不满足要求继续细分，直到满足预定参数终止。
(2)对离散的曲线段创建包围盒(图2d、图2e)。如图2d、图2e所示，判断直线与包围盒是否相交，若不相交则去除这些包围盒对应曲线段的区间，避免对这些不包含交点的区间进行无效迭代，同时筛选出与直线相交的包围盒作为下一步进行区间迭代的预迭代区间。
(3)对筛选的区间进行区间Newton迭代。若存在多个解，迭代过程中迭代区间会进行自适应划分，进一步对细分区间进行判断；若无解，则停止迭代,终止计算；若存在唯一解，则在这一区间内给定一个初值，通过牛顿迭代法或拟牛顿迭代法求出真解。
3 直线与NURBS曲面快速求交算法
3.1 构建直线与NURBS曲面求交目标函数
为了表述清晰，本文在局部坐标系下进行直线与NURBS曲面求交[13]，将直线L和NURBS曲面所在的原坐标系转换为以L为Z轴的新坐标系。在新坐标系下，直线L与Z轴重合，直线L的参数方程L(t)=O+tDZ,t∈[-∞,+∞],DZ为Z轴的单位向量,则NURBS曲面的控制点Pi,j变换为
此时，NURBS曲面方程为
(14)
0≤u,v≤1
图3 绝对坐标系下的直线与NURBS曲面
Fig.3 Line and NURBS surface in the absolute coordinate system
如图3、图4所示，经过空间坐标变换后，直线L和NURBS曲面S(u,v)求交转化为Z轴与NURBS曲面S′(u,v)求交。在新坐标系下交点必然既位于NURBS曲面S′(u,v)上，又同时位于Z轴上，由此可知，交点的坐标x、y均为零，坐标z未受约束为自由项。又由NURBS曲面的规范性可知，对所有的(u,v)∈[0,1]×[0,1],有
故在新坐标系下求解直线与NURBS曲面交点的目标函数为
图4 局部坐标系下的直线与NURBS曲面
Fig.4 Line and NURBS surface in the local coordinate system
(15)
式中，
分别为NURBS曲面控制点的坐标。
3.2 基于仿射算术的Krawczyk算子运算
直线与NURBS曲面求交中同样存在直线与NURBS曲线求交遇到的问题，另外，在Krawczyk算子计算中也需要考虑相关性问题。
对于在NURBS曲面基函数的计算中需多次调用样条基函数存在长计算链的区间运算问题，处理方式与前述内容相同，此处不作赘述。对于如何快速准确定位包含真实解的有效区间，本文提出了面曲率的二叉树细分对曲面进行自适应区域划分方法，并基于直线与体包围盒求交的思想，快速准确定位包含真实解的有效区间，从而在很大程度上缩小了求交搜索范围，减少了不必要的迭代次数，极大提高了区间迭代算法的效率。另外，在Krawczyk算子计算中需要考虑相关性问题。根据前文对仿射算术限制的讨论，Y(k)为实系数矩阵，F′(X(k))为区间变量矩阵。基于Krawczyk算子公式，在每次迭代时将当前步要迭代的区间分量变为仿射形式
代入区间变量矩阵F′(X(k))得到矩阵的仿射形式，再与实数矩阵Y(k)进行仿射形式的乘法运算得到仿射矩阵，再将仿射矩阵转化为新的区间矩阵，代入Krawczyk算子公式计算求得K(X(k))。
算法具体步骤描述如下：
(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)
其中，Y(k)=[m(F′(X(k)))]-1,M=(m1,m2,…,mn)T为X的中点，R=(r1,r2,…,rn)T为X的半径的向量形式，
(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
表1 不同直线与参数空间NURBS曲线求交时误差与计算时间比较(迭代容差为10-12 mm)
Tab.1 Comparison of error and calculation time between intersection of straight line and parameter space NURBS curve(iteration tolerance is 10-12 mm)
直线误差(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
表2 不同直线与参数区间NURBS曲面求交时误差与计算时间比较(迭代容差为10-12 mm)
Tab.2 Comparison of error and calculation time between intersection of straight line and NURBS surface(iteration tolerance is 10-12 mm)
直线误差(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
表3 直线与复杂NURBS曲面求交时误差与计算时间比较(迭代容差为10-12 mm)
Tab.3 Comparison of error and calculation time between intersection of straight line and complex NURBS surface(iteration tolerance is 10-12 mm)
误差(mm) 计算时间(ms) M1M2M37.114×10-1394.361.7平均时间65～70
表4 不同直线与弯管曲面求交时的误差与计算时间比较(迭代容差为10-9 mm)
Tab.4 Comparison of error and calculation time between intersection of straight line and the NURBS surface of canal(iteration tolerance is 10-9 mm)
直线误差(mm)计算时间(ms)M1M2M414.379×10-1056.942.128.692×10-1087.361.7平均时间50～60
5 结论
(1)本文将仿射算术与区间运算相结合，一定程度上解决了传统区间运算的“保守性”，减少了迭代过程中不必要的求交判断，并放宽了对初始区间的要求。
(2)将基于边曲率或面曲率的子域分解方法应用到求交算法中，快速定位预迭代区间，在相同的容差要求下，比传统区间迭代求交算法更快，提高了迭代求交算法的效率。
(3)通过计算区间算子判断给定直线与参数空间NURBS曲线或直线与NURBS曲面有无交点和存在交点时的交点数目，可保证求解交点精度，且不会遗漏交点。
[1] MOORE R E. Interval Analysis[M]. Englewood Clifs: Prentice Hall, 1966.
[2] PIEGL L, TILLER W. The NURBS Book[M]. Berlin: Springer Science & Business Media, 2012.
[3] NELLES O. Nonlinear System Identification: from Classical Approaches to Neural Networks and Fuzzy Models[M]. Berlin: Springer Science & Business Media, 2013.
[4] OWEN S J. An Introduction to Automatic Mesh Generation Algorithms—Part Ⅱ[R].Albuquerque: Sandia National Lab, 2016.
[5] CHEN J , XIAO Z , ZHENG Y, et al. Automatic Sizing Functions for Unstructured Surface Mesh Generation[J]. International Journal for Numerical Methods in Engineering, 2017, 109(4):577-608.
[6] DEVENDRA K E A . Guaranteed Ray Intersections with Implicit Surfaces[J]. ACM Siggraph Computer Graphics, 1989, 23(23):297-306.
[7] GOLDBERG P W, PAPADIMITRIOU C H, SAVANI R. The Complexity of the Homotopy Method, Equilibrium Selection, and Lemke-Howson Solutions[J]. ACM Transactions on Economics and Computation, 2013, 1(2):9.
[8] PIRNIA M, CANIZARES C A, BHATTACHARYA K, et al. A Novel Affine Arithmetic Method to Solve Optimal Power Flow Problems with Uncertainties[J]. IEEE Transactions on Power Systems, 2014, 29(6):2775-2783.
[9] 李庆扬, 莫孜中, 祁力群. 非线性方程组的数值解法[M] .北京:科学出版社, 1987.
LI Qingyang, MO Zizhong, QI Liqun. Numerical Methods for No-linear Systems of Equations[M]. Beijing: Science Press, 1987.
[10] KRAWCZYK R. A Class of Interval-newton-operators[J]. Computing, 1986, 37(2):179-183.
[11] NEUMAIER A, ZUHE S. The Krawczyk Operator and Kantorovich’s Theorem[J]. Journal of Mathematical Analysis and Applications, 1990, 149(2):437-443.
[12] 寿华好, 王国瑾, 沈杰. 区间算术和仿射算术的研究与应用[J].中国图象图形学报, 2006, 11(10): 1351-1358.
SHOU Huahao, WANG Guojin, SHEN Jie. A Survey on Research and Applications of Interval Arithmetic and Affine Arithmetic[J]. Journal of Image and Graphics Image Graph, 2006, 11(10): 1351-1358.
[13] 王保庆, 张 俐, 李东升. 逆向工程中 NURBS 曲面与直线交点快速计算[J]. 工程图学学报, 2010, 31(2): 149-152.
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): 149-152.
[14] 朱建宁, 王敏杰, 魏兆成,等. 直线与高精度细分曲面交点快速计算方法[J]. 大连理工大学学报, 2013(3):376-381.
ZHU Jianning, WANG Minjie, WEI Zhaocheng. Rapid Calculation of Intersection of Line and Catmull-Clark Surface[J]. Journal of Dalian University of Technology, 2013(3):376-381.
[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.
