局部加密算法之NURBS曲面和隐式曲面求交

基于等值线法求取NURBS曲面与隐式曲面交线的原理，我们提出了一种局部加密的改进算法。通过局部加密算法减少正则网格单元顶点处^值计算数目。采用拟牛顿迭代法求交点、B样条曲线拟合参数域上的交线等改进算法，提高了NURBS曲面与隐式曲面求交算法的效率和精度，并通过MATLAB编程进行了验证。

一、全网格遍历法求解NURBS曲面与隐式曲面交线

1、基本原理

设NURBS曲面为x=g1(u，v），y=g2(u，v），z=g3(u，v），控制顶点Pij（i=0，1，…，n1；j=0，1，…，n1）呈拓扑矩形阵列，形成一个控制网格。ωij是与顶点Pij联系的权因子。Bik1(u)和色，Bjk2(v)分别为u向k1次和v向k2次规范B样条基。

设隐式曲面为f(x，y，Z）=0，其中厂的定义域为Ωf，则NURBS曲面与隐式曲面的交线方程为：f(g1，g2，g3）=h(u，v）=O的求解过程实际上就是曲面求交的过程。

2、全网格遍历求解多项式方程

（1）正则化网格的生成

曲面参数化后，将参数域Ωg沿u、v方向分别分为nu,nv份。每一个矩形单元△i,j的四个顶点分别为(ui,vj)， (ui,vj+1），(ui+1，vj+1), (ui+1，vj), (i=0,1，…，nu-1；j=0，1，…，nv-1)。对应的函数值为hij ，hi,j+1，hi+1J+1，hi+1J(i=0,1，…，nu-l;j=0,1，…，nv-1)。

(2)网格单元与等值线的交点计算

网格单元与等值线交点计算，主要是计算单元各边与等值线的交点，可以采用顶点判断和边上插值的方法计算交点，具体步骤如下：

1)若hiJ≤0，则顶点(ui,uj）记为“-”，否则记为“+”；

2)若单元的4个顶点都为“+"，或全为“-”，则网格单元与等值线h(u,v）=O无交点，否则转到3）；

3)对于两个顶点异号的单元边，采用拟牛顿迭代法计算等值线与这条边的交点，如图2中红色圆点，设hi+1,J+1，为“+”，hi+1,j为“-”，则交点(ut，vt)为：

vt=Vi+l（i=2，3，…，n）。该单元其它边上的交点可用相同方法求得。

3、等值线的生成

逐个计算出每一个网格单元与等值线的交点，用三次B样条曲线按顺序拟合这些交点，最终得到等值线h(u，v）=0。

需要注意的是，采用B样条曲线拟合单元内的交点，必须首先确定这些交点的先后次序。对于形状比较简单的等值线，可用其关键信息判断交点的先后次序。例如与u向传递的波浪面相交时，交点可按照其u坐标值的大小排序；与圆柱面相交时交点可按照交点和圆心的连线与水平线的夹角大小排序。对于形状比较复杂的等值线，目前还未找到完美的解决方法，可采用蚁群算法保证各点之间的总距离最小，以此为原则对数据点进行排序。

综上所述，全网格遍历法求解多项式方程h（u，v）=0的流程。

数据取自采用全网格遍历法编制的matlab程序求解NURBS曲面与圆柱曲面交线过程。本文所用NURBS曲面由201×401个型值点拟合而成，程序运行平台机型配置为AMD Athlon四核2.8GHz，8G内存。

当曲面曲率变化剧烈或者对交线精度要求较高时，往往需要较小的网格步长。全网格遍历要求对每个单元都要求解，由表1可知，当网格密度增加时，该算法计算量会成倍增加，导致其耗时也会成倍增加。在网格步长较小时，全网格遍历求解多项式h（u，v）=O效率很低。

从表1中不难看出，随着网格步长的减小，需要计算顶点处危值的数目也越多，而单独计算各顶点处h值所耗时间占整个求交程序耗时的绝大部分。由此可知，想要提高求交算法的效率，就需要在相同精度要求下，尽量减少尼值的计算数目。因此在求交过程中引入局部加密算法，期望通过减少危值计算数目，提高求交算法的效率。

二局部加密算法在NURBS曲面与隐式曲面求交中的应用

1、选定正则化基网格

局部加密算法首先需要选定一个网格密度较疏的正则化基网格，基网格的网格步长取△u=△v=△=O.1。求出每个网格单元四个顶点处的函数值h（u，v），运用全网格遍历求解算法求出网格单元与等值线的交点，如图4中红色星号所示。将顶点处h值储存在矩阵危。h0（i，j），(i=1，2，…，1/△+1)中。那么h0（i，j）=h((i- 1)△，(j-1）△)。

2、确定加密范围

以网格步长等于△为例，确定加密范围的方法如下：

1)将各单元与等值线的交点排序，排序方法与全网格遍历法中拟合等值线时交点排序的方法相同，排好序的交点分别记为al，a2，…，ak，ak+l，…，an(k=1，…，n）。对于单元i，各顶点编号如图5所示。以Ai为单元i的参考基点，则Ai的坐标为：uAi=fix [min(uk,uk+l)/△]△；：uAi=fix [min(vk,vk+l)/△]△，其中fix为取整函数。单元中Bi，Ci，Di的坐标可通过参考基点Ai的坐标求得，例如uBi=uAi+△， VBi =VAi；

2)如图6所示，按照交点顺序依次找到单元参考基点Al，A2，…，Ai，Ai+i，…，Am(=l，…，m)，由单元参考基点依次确定包含交点的每一个单元。

所求等值线必定在这个加密范围内。

3、网格加密

确定需要加密的网格范围后，将加密范围内的每个单元均分为四个小的单元。

4、数据处理

每一次加密网格后，需要更新顶点处h值，为避免重复计算，需要对数据做一些特殊处理：

1）确定加密区域所在的最小矩形域。

矩形域D1DZD3D4即为加密区域所在的最小矩形域。

2)生成判断矩阵。加密最小矩形域，将矩形域内每一个单元均分为四个小单元。

为便于后续计算，在加密后引入判断矩阵判断最小矩形域内的顶点是否在加密区域内。

3)更新网格顶点处h值，并将其储存在矩阵hK中。当K=1时，h1（i,j)=ho（i+VD1/△,j+VD1/△）。当K=2，3时，计算hK。

5、网格加密后单元与等值线交点的计算