CN104680006A - 基于曲面精确表示的碰撞检测方法 - Google Patents

Abstract

一种基于曲面精确表示的碰撞检测方法，其步骤如下：一、包围盒粗判:若包围盒无交,则返回无碰撞；二、利用几何体的特殊几何特性,进行粗判排除无交情况,例如利用圆柱体的上下底面排除在上底面以上及在下底面以下的几何体；三、利用几何体的对称性及相对特殊位置关系，对某些特殊情况(比如可能转化为二维的情况)直接计算出碰撞点，判断出有无碰撞,例如两个轴线相交的圆环体，可以在两轴线所在平面内判断碰撞；四、排除几何体1和几何体2相互包含的情况；五、检测几何体的边界碰撞；六、几何体的表面内部碰撞检测；七、返回无碰撞。本发明实现了快速准确地对两个几何体进行干涉检查的方法，速度快，精度高。

Description

基于曲面精确表示的碰撞检测方法

技术领域

本发明提供一种基于曲面精确表示的碰撞检测方法，它涉及一种基于曲面精确表示的干涉检查方法，由于支持多种曲面类型，该方法由若干子算法组成，属于工程应用领域。

背景技术

具有强鲁棒性的干涉检查方法在机械、机器人仿真、虚拟现实、电脑游戏和虚拟外科手术等领域具有很大的应用价值，同时建立快速精确的干涉检查方法也是这些领域发展的一个技术瓶颈。在数控加工、机器人路径规划等领域中对碰撞检测的实时性和精确性要求越来越苛刻，因此，建立一种曲面之间的快速准确的干涉检查方法很有意义。

现有技术中，可以用于圆环体之间的干涉检查方法有多种。

现有技术一：基于包围盒的干涉检测法。该类方法通过构造形状简单且与被检测物体相差不大的包围盒，并对其进行碰撞检测，如果包围盒不相交，则被检测物体必不相交，如果包围盒相交，再进行被检物体进行精确检测。这类方法有轴向包围盒检测法，方向包围盒检测法等。

现有技术二：基于空间分解的干涉检查法。该方法通过将空间分解为体积相等的小单元格并且只对占据同一单元格或相临单元格的被检物体进行相交测试，这类方法有八叉树划分空间法，二叉空间剖分法等。

现有技术一不易建立起既简单又对被检物体包裹紧密的包围盒。

现有技术二存储空间占用大并且灵活性差。

现有技术在干涉检查过程中不可避免地要进行几何形体的三角片离散化操作或者曲面求交操作。本发明在干涉检查过程中不使用离散信息，而是基于曲面的原始几何信息进行快速精确干涉检查方法。

发明内容

1、目的：本发明的目的是提供一种基于曲面精确表示的碰撞检测方法，它是建立一套基于曲面原始几何信息的干涉检查方法，该方法能快速精确地判断空间中两曲面之间是否发生干涉。该干涉检查方法可以应用到工程实际中，如刀具的运动干涉检查及机器人的运动干涉检查等。

2、技术方案：本技术方案以共线法矢引理为理论依据，共线法矢引理描述为：如果两个几何体发生干涉，则必然存在一个点和一条直线，满足该点是两个几何体和直线三者的交集，且满足该直线同时垂直于这两个曲面，即：几何体T₁/几何体T₂有交点(碰撞)<＝＝>T₁和T₂在其极值位置处有一条T₁和T₂的共线法矢L，满足 其中T₁和T₂至少有一个是以下曲面的一种：平面(无穷大)、圆柱体、圆锥体、球体、圆环体、椭球体、多面体(相关算法不太完整)、裁剪曲面(不限于B样条曲面，目前实现了单张曲面情况的相关算法)，事实上，凡是等半径的球扫掠体(含等半径的管道体)都可以使用该判断准则。以上几何体是指由基本体素(方体、球体、圆柱体等)通过有限次并、交、差集合运算构成的点集。

下面给出该引理的证明：

假设Ω_i是一个球扫成体(sphere-swept body)，即一个等半径的球，半径为r_i＞0，其球心在曲面S_i(u,v)上运动形成的包络体，i＝1,2。在Ω₁和Ω₂有交(碰撞)，且曲面S₁(u,v)和S₂(u,v)无交的情况下，令λ是个变量，且0≤λ≤1，令Ω_i(λ)是球缩略体(sphere-swept body)，其半径变化规律为λ·r_i，i＝1,2，其球心在曲面S_i(u,v)上运动，令集合可以推出Φ₁至少含有1个元素，因而不为空，取Φ₁的下确界λ₁＝infΦ₁。Ω₁(λ₁)和Ω₂(λ₁)必然相切，(见附图1)，于是通过切点P就有共线法矢L。

(所述“下确界”的概念是数学分析中最基本的概念。考虑一个实数集合M.如果有一个实数S，使得M中任何数都大于或等于S,那么就称S是M的一个下界。

在所有那些下界中如果有一个最大的下界，就称为M的下确界。)

注意：

(1)Ω_i(λ)＝∪B_p,λ，其中B_p,λ表示以p为球心、λ为半径的球，p∈S_i(u,v)。若ξ＜η,则所以

(2)Ω_i(λ)是有界闭集，所有B_p,λ的内点构成Ω_i(λ)的内点，Ω_i(λ)的边界点必是某个B_p,λ的边界点。

下面证明Ω₁(λ₁)和Ω₂(λ₁)必然相切。令集合Φ₀至少含有0，因而不为空，取Φ₀的上确界λ₀＝supΦ₀。容易看出λ₀≤λ₁。假设λ₀＜λ₁，取η＝(λ₀+λ₁)/2，Ω₁(η)∩Ω₂(η)等于空或不等于空都导致矛盾，所以λ₀＝λ₁，记η＝λ₀。

(1)如果由于Ω₁(η)、Ω₂(η)都是有界闭集，所以存在ε>0使||p-q||>ε,对任何的p∈Ω₁,q∈Ω₂都成立。总可以找到k>1，也就是把Ω₁(η)和Ω₂(η)都放大一点，使(见附图2)，这就与η＝supΦ₀矛盾，故

(2)任取Ω₁(η)∩Ω₂(η)＝中的一点p',于是存在p和q两点，使且如果p'是B_q,ξ或B_q,η的内点，则存在0<k<1使(见附图3),于是这与η＝supΦ₁矛盾，故p'只能是B_p,ξ和B_q,η边界上的点，p'是B_p,ξ和B_q,η的切点，即Ω₁(λ₁)和Ω₂(λ₁)相切。

根据共线法矢理论，我们可以得出几何体碰撞检测的通用算法

输入：几何体1和几何体2，容差ε>0

(该几何体1和几何体2是为了引用方便，将进行碰撞检测的两个几何体一个称为几何体1，另一个称为几何体2)

输出：是否碰撞及一个碰撞点

本发明一种基于曲面精确表示的碰撞检测方法，其步骤如下：

步骤一：包围盒粗判:若包围盒无交,则返回无碰撞；

步骤二：利用几何体的特殊几何特性,进行粗判排除无交情况,例如利用圆柱体的上下底面排除在上底面以上及在下底面以下的几何体；

步骤三：利用几何体的对称性及相对特殊位置关系，对某些特殊情况(比如可能转化为二维的情况)直接计算出碰撞点，判断出有无碰撞,例如两个轴线相交的圆环体，可以在两轴线所在平面内判断碰撞；

步骤四：排除几何体1和几何体2相互包含的情况：

(a)在几何体1内取一点，检测该点是否在几何体2内，若在几何体2内则返回有碰撞；

(b)在几何体2内取一点，检测该点是否在几何体1内，若在几何体1内则返回有碰撞；

步骤五：几何体的边界碰撞检测：

(a)几何体1的所有边界曲线与几何体2进行碰撞检测，若边界有碰撞即返回；

(b)几何体2的所有边界曲线与几何体1进行碰撞检测，若边界有碰撞即返回；

步骤六：几何体的表面内部碰撞检测：

遍历几何体1的所有表面，

遍历几何体2的所有表面，

计算面/面的所有极值位置，

遍历极值位置处的共线法矢L,

L在几何体1内裁断得L1，

若L1不为空，

则L1在几何体2内裁断得L2，

若L2不为空，则返回有碰撞，

结束遍历极值位置

结束遍历几何体2

结束遍历几何体1

对上述排版的解释：对于每一个小步骤，相对于其上一步骤向右缩进2个字符，若有两个小步骤的第一个字在同一竖直线上，则为一个遍历循环。

也可以这样表述：

对于几何体1和几何体2所有曲面的两两组合都做如下操作：计算两曲面的所有共线法矢，对于每一条共线法矢L，先在几何体1内裁断得L1，若

L1不为空，则L1在几何体2内裁断得L2，若L2不为空，则返回有碰撞，否则进行下一曲面对的判断，若所有曲面对都没有碰撞，则返回无碰撞。

步骤七：返回无碰撞

其中，在步骤六中所述的“遍历极值位置处的共线法矢L”，其含义是:先计算两几

何体1、2的共线法失，然后遍历，共线法失的计算是本算法的关键技术，一般来说共线法失的求解是根据几何体本身的特性而列出方程求解而得的，比如共线法矢过圆锥体/圆柱体轴线且与圆锥体/圆柱体的轴线成固定角度，共线法矢过圆环体的中心圆线且与圆环中心圆切线垂直等。

3、优点与功效

本发明实现了快速准确地对两个几何体进行干涉检查的方法。基于曲面的原始几何信息，既不使用曲面的离散信息，也不用计算曲面交线，有速度快和精度高的特点

附图说明

图1-a在Ω₁和Ω₂有交

图1-bΩ₁(λ₁)和Ω₂(λ₁)相切

图2存在k＞1使

图3p'是B_p,ξ或B_q,η的内点情况

图4两个圆锥共线法矢与轴线相交

图5-a圆锥的参数定义

图5-b圆环的参数定义

图6圆锥/圆环的共线法矢

图7α、π/2-α和π/2+α角度的几何意义

图8过中心圆1轴线的平面交中心圆2的轴线于Q₂点

图9当点P₁和P₂在极值位置时直线P₁P₂与中心圆及其轴线都相交

图10平面π与E的交线为

图11椭球体上一点的法矢与xy平面的交点位于一椭圆内

图12用于估计初始点的三角形

图13有向距离函数的定义

图14向量场中向量夹角的计算

图15平面/参数曲面迭代

图16球面/参数曲面迭代

图17圆柱面/参数曲面迭代

图18圆环面/参数曲面迭代

图19本发明所述方法流程图

图中序号、符号、代号说明如下

Ω_i为几何体i，i＝1,2

E为椭球体，平面π与E的交线

B_p,λ表示以p为球心、λ为半径的球，

n(u，v)为曲面法矢

具体实施方式：

如图1-a－图19，本发明一种基于曲面精确表示的碰撞检测方法(见附图19)，其步骤如下：

步骤一：包围盒粗判:若包围盒无交,则返回无碰撞；

步骤二：利用几何体的特殊几何特性,进行粗判排除无交情况；

步骤三：利用几何体的对称性及相对特殊位置关系，对某些特殊情况(比如可能转化为二维即2D的情况)直接计算出碰撞点，判断出有无碰撞,例如两个轴线相交的圆环体，可以在两轴线所在平面内判断碰撞；

步骤四：排除几何体1和几何体2相互包含的情况：

(a)在几何体1内取一点，检测该点是否在几何体2内，若在几何体2内则返回有碰撞；

(b)在几何体2内取一点，检测该点是否在几何体1内，若在几何体1内则返回有碰撞；

步骤五：几何体的边界碰撞检测：

(a)几何体1的所有边界曲线与几何体2进行碰撞检测，若边界有碰撞即返回；

(b)几何体2的所有边界曲线与几何体1进行碰撞检测，若边界有碰撞即返回；

步骤六：几何体的表面内部碰撞检测：

遍历几何体1的所有表面，

遍历几何体2的所有表面，

计算面/面的所有极值位置，

遍历极值位置处的共线法矢L,

L在几何体1内裁断得L1，

若L1不为空，

则L1在几何体2内裁断得L2，

若L2不为空，则返回有碰撞，

结束遍历极值位置

结束遍历几何体2

结束遍历几何体1

对上述排版的解释：对于每一个小步骤，相对于其上一步骤向右缩进2个字符，若有两个小步骤的第一个字在同一竖直线上，则为一个遍历循环。

也可以这样表述：

对于几何体1和几何体2所有曲面的两两组合都做如下操作：计算两曲面的所有共线法矢，对于每一条共线法矢L，先在几何体1内裁断得L1，若

L1不为空，则L1在几何体2内裁断得L2，若L2不为空，则返回有碰撞，否则进行下一曲面对的判断。

步骤七：返回无碰撞

其中，在步骤六中所述的“遍历极值位置处的共线法矢L”，其含义是:先计算两几何体的共线法失，然后遍历，共线法失的计算是本算法的关键技术，一般来说共线法失的求解是根据几何体本身的特性而列出方程求解而得的，比如共线法矢过圆锥体/圆柱体轴线且与圆锥体/圆柱体的轴线成固定角度，共线法矢过圆环体的中心圆线且与圆环中心圆切线垂直等。

对步骤六中排版的解释，对于每一个小步骤，相对于其上一步骤向右缩进2个字符，若有两个小步骤的第一个字在同一竖直线上，则为一个遍历循环。

下面详细介绍一下共线法矢的计算。

共线法矢的计算分为以下几种情况

1.两圆锥体共线法矢的计算

设两个圆锥的顶点、轴线及圆锥半角分别为O₁,O₂,Z₁,Z₂,θ₁,θ₂(见附图4)。对于两轴线Z₁,Z₂重合、平行、相交的情况，通过特殊方法可以得到共线法矢,对于轴线异面的情况，可以证明共线法矢与两轴相交，令共线法矢在轴线Z₁,Z₂上的交点分别为P₁,P₂,其长度为L。由于共线法矢与圆锥面垂直，所以共线法矢与圆锥轴线成固定的角度，所以满足一下方程：

$\left\{\begin{array}{c}{P}_1=O_1+Z_1{t}_1...\left(1\right)\\ P_2=O_2+Z_2{t}_2...\\ (P_1-P_2,P_1-P_2)=L^{2...}\left(2\right)\\ (\frac{P_1-P_2}{L},Z_1)=\&PlusMinus;c_1...\left(4\right)\\ (\frac{P_1-P_2}{L},Z_2)=\&PlusMinus;c_2...\left(5\right)\end{array}\right.$

其中c₁＝sinθ₁，c₂＝sinθ₂(,)表示矢量的点积。由方程(4)和(5)可以得到，

$\begin{array}{c}{t}_1-(Z_1,Z_2)t_2=-(v,Z_1)\&PlusMinus;L{c}_1\\ (Z_1,Z_2)t_1-t_2=-(v,Z_2)\&PlusMinus;L{c}_2\\ v=O_1-O_2\end{array}$

注意上式中由于假设圆锥的轴线不平行，所以有(Z1,Z2)≠±1，另Δ＝(Z1,Z2)，则解上面两个二元一次线性方程组得，

$t_1=\frac{(v,Z_2\mathrm{\&Delta;}-Z_1)\&PlusMinus;L(c_1-c_2\mathrm{\&Delta;})}{1-{\mathrm{\&Delta;}}^2}$

$t_2=\frac{(v,Z_2-Z_1\mathrm{\&Delta;})\&PlusMinus;L(c_1\mathrm{\&Delta;}-c_2)}{1-{\mathrm{\&Delta;}}^2}$

将上面关于t₁,t₂的表达式带入方程(3)中，就得到关于L的一元二次方程，直接求解即得L，进而得到t₁,t₂和P₁,P₂。

当圆锥2退化为圆柱时，则取c₂＝0，上面的计算过程仍然有效，即上面的算法对于圆锥体/圆柱体也是成立的。

2.圆锥体与圆环体共线法矢的计算

假设圆锥的轴线上一点为P，轴线单位方向为D，其方程为L(t)＝P+Dt(见附图5-a),圆环的中心圆的圆心为O，半径为r，其方程为

C(θ)＝O+rXcosθ+rYsinθ(见附图5-b),易知中心圆的单位切矢为

-Xsinθ+Ycosθ。由于共线法矢与圆环的中心圆的切线垂直，所以有，

(P+Dt-O-rXcosθ-rYsinθ,-Xsinθ+Ycosθ)＝0

令P-O＝v，则根据上式可以推导出，

$t=\frac{(v,X)\sin \mathrm{\&theta;}-(v,Y)\cos \mathrm{\&theta;}}{(D,Y)\cos \mathrm{\&theta;}-(D,X)\sin \mathrm{\&theta;}}...\left(1\right)$

令共线法矢长度为L，则

(P+Dt-O-rXcosθ-rYsinθ，Dt-O-rXcosθ-rYsinθ)＝L²

上式展开可得，

(v,v)+2(v,D)t-2r(v,X)cosθ-2r(v,Y)sinθ-2rt(D,X)cosθ-2rt(D,Y)sinθ+t²+r²＝L²..........(2)

另一方面由于共线法矢与圆锥轴线成固定的角度，(见附图6)所以，

$(\frac{P+\mathrm{Dt}-O-\mathrm{rX}\cos \mathrm{\&theta;}-\mathrm{rY}\sin \mathrm{\&theta;}}{L},D)=\&PlusMinus;c$

其中c为常数，即c＝sinα。上式平方得，

[(v,D)+t-r(D,X)cosθ-r(D,Y)sinθ]²＝L²c²...............................(3)

利用(2)式和(3)式可以消去L，得到关于t和θ的方程，再根据(1)式消去t即可得到只含一个变量θ的三角函数方程，利用三角函数中的万能公式：

$\cos \mathrm{\&theta;}=\frac{1-{\mathrm{\&tau;}}^2}{1+{\mathrm{\&tau;}}^2},\sin \mathrm{\&theta;}=\frac{2\mathrm{\&tau;}}{1+{\mathrm{\&tau;}}^2},\mathrm{\&tau;}=\tan \frac{\mathrm{\&theta;}}{2},$ 此方程可转化为一元8次代数方程。

考虑到一元8次代数方程的直接求解方法数值稳定性差的问题，须另构造数值方法求解共线法矢。根据方程(1)可推导出，

$\frac{\mathrm{dt}}{\mathrm{d\&theta;}}\text{=}\frac{(v,X)(D,Y)-(v,Y)(D,X)}{{[(D,Y)\cos \mathrm{\&theta;}-(D,X)\sin \mathrm{\&theta;}]}^2}...\left(4\right)$

上式中的参数θ是圆环中心圆的参数，过圆环中心圆上对应参数θ的点且与该点处的切矢垂直的平面交圆锥轴线于一点，此点对应圆锥轴线的参数为t。根据上式可推出，在θ的有效参数域内t的变化是单调的，此特性决定了数值求解方法的稳定性。定义一个函数f(θ)，θ是圆环中心圆的参数，其参数域为[0,2π]，函数f(θ)返回值为共线法矢与圆锥轴线的夹角(假设p为共线法矢与圆锥轴线的交点，q为共线法矢与圆环中心圆上的交点，总是取共线法矢的方向为q-p)。数值求解的目标就是求解f(θ)＝π/2-α或π/2+α，其中α为圆锥的顶角。

数值算法的基本步骤是：

步骤一：将[0,2π]按函数f(θ)的奇异性、单调性分成4个子区间，

步骤二：在每个子区间内用区间对分法迭代计算满足f(θ)＝π/2-α或π/2+α的θ值(见附图7)。

上面算法的关键是如何将参数域[0,2π]分割成单调的子区间。下面给出将参数域[0,2π]近似分割成单调子区间的方法：

步骤一：

由于当(D,Y)cosθ-(D,X)sinθ＝0时t趋于正或负无穷大，求解此三角方程得θ₁和θ₂，

步骤二：由于直接对f(θ)求导数结果很复杂，考虑共线法矢与圆锥轴线单位矢量D的点积(P+Dt-O-rXcosθ-rYsinθ,D)，为求此点积的极值，对该点积求导并另其为0，则有

$(D\frac{\mathrm{dt}}{\mathrm{d\&theta;}}+\mathrm{rX}\sin \mathrm{\&theta;}-\mathrm{rY}\cos \mathrm{\&theta;},D)=0$ 即 $\frac{\mathrm{dt}}{\mathrm{d\&theta;}}=r[(D,Y)\cos \mathrm{\&theta;}-(D,X)\sin \mathrm{\&theta;}],$

再结合关系式(4)有， $(D,Y)\cos \mathrm{\&theta;}-(D,X)\sin \mathrm{\&theta;}=-\sqrt[3]{\frac{(v,X)(D,Y)-(V,Y)(D,Y)}{r}},$ 解此三角方程得θ₃和θ₄，

1)根据θ₁、θ₂、θ₃和θ₄再结合θ的周期性将[0,2π]分割成4个子区间。

3.两圆环体共线法矢的计算

圆环体/圆环体极值位置的计算等价于圆/圆极值位置的计算。

假设圆环体的中心圆为C_i(θ_i)＝O_i+r_iX_i cosθ_i+r_iY_i sinθ_i，i＝1,2，Z₁和Z₂是两圆所在平面的单位法矢，令P_i是C₁和C₂上的极值点对，根据极值理论有：

(P₂-P₁,-X₁ cosθ₁+Y₁ sinθ₁)＝0…………………………(1)

(P₂-P₁,-X₂cosθ₂+Y₂sinθ₂)＝0…………………………(2)

记v＝O₂-O₁，根据方程(1)和方程(2)可以分别推出，

r₂[(Y₁,X₂)cosθ₁-(X₁,X₂)sinθ₁]cosθ₂

+r₂[(Y₁,Y₂)cosθ₁-(X₁,Y₂)sinθ₁]sinθ₂＝………………………(3)

(v,X₁)sinθ₁-(v,Y₁)cosθ₁

[r₁(X₁,Y₂)cosθ₁+r₁(Y₁,Y₂)sinθ₁-(v,Y₂)]cosθ₂

...............(4)

-[r₁(X₁,X₂)cosθ₁+r₁(Y₁,X₂)sinθ₁-(v,X₂)]sinθ₂＝0

将公式(3)、(4)写成如下形式，

$\left\{\begin{array}{c}C\cos {\mathrm{\&theta;}}_2-D\sin {\mathrm{\&theta;}}_2=0\\ E\cos {\mathrm{\&theta;}}_2+F\sin {\mathrm{\&theta;}}_2=G\end{array}\right.$

其中

E＝r₂[(Y₁,X₂)cosθ₁-(X₁,X₂)sinθ₁]

F＝r₂[(Y₁,Y₂)cosθ₁-(X₁,Y₂)sinθ₁]

G＝(v,X₁)sinθ₁-(v,Y₁)cosθ₁

C＝r₁(X₁,Y₂)cosθ₁+r₁(Y₁,Y₂)sinθ₁-(v,Y₂)

D＝r₁(X₁,X₂)cosθ₁+r₁(Y₁,X₂)sinθ₁-(v,X₂)。

解此关于θ₂的方程得，

$\cos {\mathrm{\&theta;}}_2=\frac{\mathrm{DG}}{\mathrm{FC}+\mathrm{DE}}$

$\sin {\mathrm{\&theta;}}_2=\frac{\mathrm{CG}}{\mathrm{FC}+\mathrm{DE}}$

所以有

${\left(\frac{\mathrm{DG}}{\mathrm{FC}+\mathrm{DE}}\right)}^2+{\left(\frac{\mathrm{CG}}{\mathrm{FC}+\mathrm{DE}}\right)}^2=1$

即

(D²+C²)G²＝(FC+DE)²……………………………(5)

方程(5)仅含有cosθ₁,sinθ₁，且次数最高为4，所以方程(5)等价一个8次多项式方程。本文采用Newton迭代的方法求解方程(5)，其中关键是确定初始解。下面针对两中心圆轴线异面的情况，讨论如何划分区间[0,2π]以便于初始解的确定。下面针对两中心圆轴线异面的情况，讨论如何划分区间[0,2π]以便于初始解的确定。

当点P₁沿中心圆1运动，对应参数为θ₁，令平面π(过中心圆1的轴线)与中心圆2的轴线交于Q₂(见附图8)。假设中心圆2的轴线方程为P(t)＝O₂+Z₂·t，则Q₂点对应的参数t满足，

$t=\frac{(v,X_1)\sin {\mathrm{\&theta;}}_1-(v,Y_1)\cos {\mathrm{\&theta;}}_1}{(Z_2,Y_1)\cos {\mathrm{\&theta;}}_1-(Z_2,X_1)\sin {\mathrm{\&theta;}}_1}...\left(6\right)$

解方程(Z₂,Y₁)cosθ₁-(Z₂,X₁)sinθ₁＝0得ξ₁和ξ₂作为分段区间[0,2π]的两个角度(当θ₁取ξ₁和ξ₂时，平面π平行于中心圆2的轴线)。

由圆的几何性质易知，当直线P₁Q₂与两中心圆及其轴线都相交时，点P₁和P₂就是两中心圆的极值位置(见附图9)。直线P₁Q₂与中心圆2所在平面平行是一种特殊情况，此时P₁Q₂与中心圆2及其轴线可能无交。P₁Q₂与中心圆2所在平面平行等价于(P₁Q₂,Z2)＝0，即

(O₂+tZ₂-O₁-r₁X₁cosθ₁-r₁Y₁sinθ₁,Z₂)＝0

将方程(6)带入上式得

$(v,Z_2)+\frac{(v,X_1)\sin {\mathrm{\&theta;}}_1-(v,Y_1)\cos {\mathrm{\&theta;}}_1}{(Z_2,Y_1)\cos {\mathrm{\&theta;}}_1-(Z_2,X_1)\sin {\mathrm{\&theta;}}_1}=r_1(X_1,Z_2)\cos {\mathrm{\&theta;}}_1+r_1(Y_1,Z_2)\sin {\mathrm{\&theta;}}_1...\left(7\right)$

解方程(7)等价于4次代数方程，解之得ξ₃、ξ₄、ξ₅和ξ₆。定义关于变量θ函数f(θ)＝(D²+C²)G²-(FC+DE)²，根据ξ_i，i＝1,2,3,4,5,6，将区间[0,2π]分成子区间，采用区间对分、f(θ)和f′(θ)极值估计方法，结合Newton迭代求解方程f(θ)＝0，得到两中心圆的极值位置，这些极值位置也是对应圆环体的极值位置。

4.圆/椭球体碰撞检测算法

由于圆环体可以看做是球沿中心圆运动产生的包络面，所以只要将检测容差改为圆环体的半径，此算法亦适用于圆环体/椭球体碰撞检测。

令空间圆C的方程为：P(θ)＝O+rXcosθ+rYsinθ，O＝(x₀,y₀,z₀)^T，X＝(A₁,B₁,C₁)，Y＝(A₂,B₂,C₂)，Z＝X×Y为圆C所在平面的单位法矢，记该平面为π；椭球体E的方程为：a,b,c皆为正实数。圆/椭球体碰撞检测算法如下：

输入：圆C，椭球体E及检测容差ε。

步骤一:若圆C和椭球体E的包围盒无交则转步骤六。

步骤二:平面π与椭球体E求交。平面与椭球体如果相交，交线一定是椭圆。如果相交转步骤三，相切转步骤四，无交则转步骤五。

步骤三：(相交)令交线为椭圆E～(见附图10)。计算圆C/E～的交线，如果有交，则圆C/椭球体E有碰撞；否则找到圆C与E～的最近点，通过迭代计算圆C与E的最近点，根据最近距离判断是否碰撞。转步骤六。

步骤四：(相切)判断切点是否在圆C上，若在则碰撞，否则以切点为初始点通过迭代计算圆C与E的最近点，根据最近距离判断是否碰撞。转步骤六。

步骤五：(无交)求平面π与椭球体E的最近点(此点唯一存在)及最近距离d。若d大于ε则C/E无碰撞，转步骤六；否则，令与平面π平行且距离为ε的平面(与E在π同侧)为π'，计算π'与E的交线，记该椭圆为E'。易知E'的点到E的距离小于ε。计算E'/C的交线：

(1)若E'/C相交或相切，则圆C/椭球体E有碰撞，转步骤六。

(2)若E'/C无交，可分3中情况处理：

(a)E'/C相离，则以E'/C的最近点为初始点，迭代计算圆C/E的最近点，进而判断出是否碰撞，转步骤六；

(b)E’包含C，则C/E碰撞，转步骤六；

(c)C包含E'，处理方法同(a)，转步骤六。

步骤六：算法结束。

以上算法中用到了平面/椭球体求交方法，为提高收敛速度，本文根据初始点，利用Newton迭代求解非线性特征方程。对于椭球体退化为球或旋转体的情况，可用特殊方法处理，下面仅讨论对a>b>c的情况。初始点的估计是基于下面的引理：

引理：令椭球体E的方程为满足a＞b＞c＞0，则E表面上所有z坐标大于0的点，其法矢与xy平面的交点皆位于xy平面上的一个椭圆内，该椭圆方程为 $\frac{x^2}{{(a-\frac{c^2}{a})}^2}+\frac{y^2}{{(b-\frac{c^2}{b})}^2}=1$ (见附图11)。

证明：由于E上一点(x₀,y₀,z₀)^T(z₀>0)的法矢方向为该点处对应的法线方程为令法线方程中z坐标为0，可解出t＝-c²，得到法线方程与xy平面交点为该点满足 $\frac{x^2}{e^2}+\frac{y^2}{f^2}<1,$ 其中 $e=a-\frac{c^2}{a},f=b-\frac{c^2}{b}.$

上面引理显示出椭球体的法矢分布规律，根据这个规律可以估计最近距离点的初始值。假设R(x₀,y₀,z₀)是空间一点，下面仅考虑x₀＞0,y₀＞0,z₀＞0情况，取点P₁,P₂,P₃和Q₁,Q₂,Q₃，构成两个三角形(见附图12)，其坐标依次为(a,0,0)、(0,b,0)、(0,0,c)、(e,0,0)和(f,0,0)。下面算法目的是在椭球体上找到一点，距离R越近越好。基本思路是实现一种映射，满足P_i→Q_i，为此引入重心坐标，需要找到这样的系数w_i，使∑w_i＝1，w_i≥0，且P＝∑P_iw_i，Q＝∑Q_iw_i，R＝(1-t)P+tQ。这里P在ΔP₁P₂P₃内，Q在ΔQ₁Q₂Q₃内，P,Q,R共线。于是t和w_i满足如下方程：

$\left[\begin{array}{ccc}(1-t)a+\mathrm{te}& & \\ & (1-t)b+\mathrm{tf}& \\ & & (1-t)c\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]=\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]$

解此方程得

$w_1=\frac{x_0}{(1-t)a+\mathrm{te}}$

$w_2=\frac{y_0}{(1-t)b+\mathrm{tf}}$

$w_3=\frac{z_0}{(1-t)c}$

将上面3个w_i的表达式代入∑w_i＝1，即得到关于t的三次方程，解此方程得t，进而确定了w_i，按w_i的比例估计出椭球体上一点作为R的初始值，该点坐标为：用随机布点的方式进行测试，本文算法平均迭代7次。

5.代数曲面/参数曲面极值位置的计算

代数曲面/参数曲面极值位置的计算基于平面向量场技术估计初始极值位置点，根据不同代数曲面的具体几何特性构造迭代公式，求解极值位置点。下面先介绍如何利用平面向量场估计初始极值位置点。

平面向量场v是连续函数v:D→R²，即对于平面区域D中的任何点p，v定义平面向量v(p)与之对应。对于平面区域D中一点p，若v(p)＝0，则称p为平面向量场v的临界点(critical point)。假设p为平面向量场v的孤立临界点，N为p的领域，且在N及其边界上没有其它临界点。对于边界上的每个点q，平面向量场v都定义非0向量v(q)与之对应。

假设对于曲面S¹(u,v)上的任意点p，在曲面S²(s,t)上都唯一存在点q与p最近。当p在曲面S¹(u,v)上变动时，得到映射M(p)＝q。对于曲面S¹(u,v)上任意一点p，S¹(u,v)＝p，有向距离函数定义为:,φ(u,v)＝(n₂(M(p)),p-M(p))，其中(,)表示内积,n₂表示曲面S²(s,t)的单位法向量。有向距离函数φ(u,v)的梯度定义了u、v参数域上的向量场，在φ(u,v)的有效定义域内，有所以对于φ(u,v)的任何临界点(u,v)，有

在以上有向距离函数的定义中，假设S¹(u,v)是一般参数曲面，S²(s,t)是代数曲面,(见附图13)。由于S¹(u,v)上的任意点p到S²(s,t)上最近点q容易计算，所以可以用解析方法得到假设S¹(u,v)是裁剪曲面，Ω＝(V,T)是该裁剪曲面的三角片集合，其中V表示所有顶点集合，T表示所有三角片集合。利用S²(s,t)的几何特性，计算出顶点集合V中每个顶点对应的φ(u,v)。对于三角片集合T中的任何三角形，假设v₁,v₂,v₃是其三个顶点对应的φ(u,v)，如果||v_i||＜ε对某个i＝1,2,3成立，则该三角形内可能存在临界点，否则计算v₁和v₂,v₂和v_3，及v₃和v₁的夹角之和α，若则该三角形内可能存在临界点。这里“v₁和v₂夹角”是指向量v₁按逆时针转到另向量v₂的夹角(见附图14)：

以上提出了一种针对代数曲面和裁剪曲面的初始极值位置点的估计方法。下面针对具体的代数曲面，构造对应的迭代公式，利用这些迭代公式可实现初始极值位置点的迭代求精。

(1)平面与参数曲面极值位置的计算

假设给定参数曲面S(u,v)和平面π，曲面上各点处单位法矢为n(u,v)，平面的单位法矢为N(见附图15)，定义目标函数f(u,v)＝(S_u,N)，g(u,v)＝(S_v,N)，对u,v求导得，

f_u＝(S_uu,N)

f_v＝(S_vv,N)

g_u＝(S_uv,N)

g_v＝(S_vv,N)

利用导数f_u,f_v,g_u和g_v构造迭代Newton公式如下,

$\left\{\begin{array}{c}{f}_u\mathrm{\&Delta;u}+f_v\mathrm{\&Delta;v}=-f\\ g_u\mathrm{\&Delta;u}+g_u\mathrm{\&Delta;v}=-g\end{array}\right.$

(2)球面与参数曲面极值位置的计算

对于球面与参数曲面的情况，给定球心为O(x₀,y₀,z₀)，及参数曲面S(u,v)，曲面上各点处单位法矢为n(u,v)(见附图16)。曲面S(u,v)的法矢方程为L(s)＝S(u,v)+s·n(u,v)，要找到这样的u,v参数，使L(s)通过球心O。在L(s)上找一点距离O最近，由于

(S+s·n-O,S+s·n-O)

＝(V+s·n,V+s·n)

＝(V,V)+2s(n,V)+s²

＝(V,V)-(n,V)²+[(n,V)-s]²

所以L(s)与球心O的最近距离为(V,V)-(n,V)²，这里V＝S-O。基于以上分析令迭代目标函数为f(u,v)＝(V,V)-(n,V)²，求导得

f_u＝2(V_u,V)-2(n,V)[(n_u,V)+(n,V_u)]

f_v＝2(V_v,V)-2(n,V)[(n_v,V)+(n,V_v)]

这里V_u＝S_u，V_v＝S_v。利用根据裁剪曲面三角网格估计出的初始极值位置点，结合导数f_u和f_v即可实现球面/参数曲面极值位置点的Newton迭代。

(3)圆柱面与参数曲面极值位置的计算

对于圆柱面/参数曲面情况，令给定的圆柱面轴线为P(t)＝P+tD，其中P是圆柱轴线上一点(x₀,y₀,z₀)，D是圆柱轴线单位方向，给定参数曲面S(u,v)，曲面上各点处单位法矢为n(u,v)(见图17)。曲面S(u,v)的法矢方程为L(s)＝S(u,v)+s·n(u,v)，要使L(s)也是圆柱面的法矢，则L(s)必须通过圆柱面轴线P(t)，且满足(n(u,v),D)＝0。考虑到两直线两直线L(s)和P(t)上两点距离之平方，令

g(s,t)＝(S+s·n-P-tD,S+s·n-P-tD)

＝(V+s·n-tD,V+s·n-tD)

＝(V,V)+s²+t²+2s(n,V)-2t(D,V)-2st(D,n)

这里V＝S-P，对g(s,t)求导得，

$\frac{\&PartialD;g}{\&PartialD;s}=2(n,V)+2s-2t(D,n)$

$\frac{\&PartialD;g}{\&PartialD;t}=2t-2(D,V)-2s(D,n)$

所以两直线的最近距离点对满足，

$\left\{\begin{array}{c}s-(D,n)t=-(n,V)\\ (D,n)s-t=-(D,V)\end{array}\right.$

解此方程得，

$s=\frac{(D,n)(D,V)-(n,V)}{\mathrm{\&Delta;}}$

$t=\frac{(D,V)-(n,V)(D,n)}{\mathrm{\&Delta;}}$

其中Δ＝1-(D,n)²，将上式中的s,t代入g(s,t)得直线L(s)与P(t)的最近距离为，

$g(s,t)=(V,V)+\frac{2(D,n)(D,V)(n,V)-{(n,V)}^2-{(D,V)}^2}{\mathrm{\&Delta;}}$

所以迭代目标函数定义为，

$f(u,v)=(V,V)+\frac{2(D,n)(D,V)(n,V)-{(n,V)}^2-{(D,V)}^2}{\mathrm{\&Delta;}}+{(n,D)}^2$

f(u,v)为L(s)与P(t)最近距离的平方加上L(s)与P(t)夹角余弦的平方。实际上，当(n,D)→0时，Δ→1，于是f(u,v)可简化为，

f(u,v)＝(V,V)-(n,V)²-(D,V)²+(n,D)²

对于f(u,v)求导得，

f_u＝2(V_u,V)-2(n,V)[(n_u,V)+(n,V_u)]-2(D,V)(D,V_u)+2(n,D)(n_u,D)

f_v＝2(V_v,V)-2(n,V)[(n_v,V)+(n,V_v)]-2(D,V)(D,V_v)+2(n,D)(n_v,D)

利用上面的导数构造迭代公式：f_uΔu+f_vΔv＝-f，u_n+1＝u_n+Δu，v_n+1＝v_n+Δv， $\mathrm{\&Delta;u}=\frac{-f_{u}f}{f_u^2+f_v^2},\mathrm{\&Delta;v}=\frac{-f_{v}f}{f_u^2{f}_v^2}.$

(4)圆环面与参数曲面极值位置的计算

对于圆环面/参数曲面情况，令给定的圆环面中心圆方程为：Q(θ)＝O+rXcosθ+rYsinθ，给定参数曲面S(u,v)，曲面上各点处单位法矢为n(u,v)(见附图18).由于圆环面与参数曲面的极值位置点对连线垂直于曲面一阶偏导数S_u和S_v，且与中心圆的切矢垂直，所以定义目标函数为：f(u,v)＝(S-Q,S_u)，g(u,v)＝(S-Q,S_v)，h(u,v)＝(S-Q,-Xsinθ+Ycosθ)＝(S-O,-Xsinθ+Ycosθ)，

将此三个函数对u,v,θ求导得，

f_u＝(S_u,S_u)+(S-Q,S_uu)

f_v＝(S_v,S_u)+(S-Q,S_uv)

f_θ＝r(Xsinθ-Ycosθ,S_u)

g_u＝f_v

g_v＝(S_v,S_v)+(S-Q,S_vv)

g_θ＝r(Xsinθ-Ycosθ,S_v)

h_u＝(S_u,-Xsinθ+Ycosθ)

h_v＝(S_v,-Xsinθ+Ycosθ)

h_θ＝-(S-O,Xcosθ+Ysinθ)

求解如下线性方程组，

$\left\{\begin{array}{c}{f}_u\mathrm{\&Delta;u}+f_v\mathrm{\&Delta;v}+f_{\mathrm{\&theta;}}\mathrm{\&Delta;\&theta;}=-f\\ g_u\mathrm{\&Delta;u}+g_v\mathrm{\&Delta;v}+g_{\mathrm{\&theta;}}\mathrm{\&theta;}=-g\\ h_u\mathrm{\&Delta;u}+h_v\mathrm{\&Delta;v}+h_{\mathrm{\&theta;}}\mathrm{\&Delta;\&theta;}=-h\end{array}\right.$

得到Δu,Δv,Δθ，构造迭代公式u_n+1＝u_n+Δu v_n+1＝v_n+Δvθ_n+1＝θ_n+Δθ。

以上所述仅为本发明较佳的实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化和替换都应涵盖在本发明的保护范围之内，另外本发明提供的方法可以集成到任何包含干涉检查的软件中。

Claims (1)

1.一种基于曲面精确表示的碰撞检测方法，其特征在于：其步骤如下：

步骤一：包围盒粗判:若包围盒无交,则返回无碰撞；

步骤二：利用几何体的特殊几何特性,进行粗判排除无交情况；

步骤三：利用几何体的对称性及相对特殊位置关系，对某些特殊情况直接计算出碰撞点，判断出有无碰撞；

步骤四：排除几何体1和几何体2相互包含的情况：

(a)在几何体1内取一点，检测该点是否在几何体2内，若在几何体2内则返回有碰撞；

(b)在几何体2内取一点，检测该点是否在几何体1内，若在几何体1内则返回有碰撞；

步骤五：几何体的边界碰撞检测：

(a)几何体1的所有边界曲线与几何体2进行碰撞检测，若边界有碰撞即返回；

(b)几何体2的所有边界曲线与几何体1进行碰撞检测，若边界有碰撞即返回；

步骤六：几何体的表面内部碰撞检测：

对于几何体1和几何体2所有曲面的两两组合都做如下操作：计算两曲面的所有共线法矢，对于每一条共线法矢L，先在几何体1内裁断得L1，若L1不为空，则L1在几何体2内裁断得L2，若L2不为空，则返回有碰撞，否则进行下一曲面对的判断，若所有曲面对都没有碰撞，则返回无碰撞；

步骤七：返回无碰撞。