基于曲面精確表示的距離極值點的計算及在刀具干涉檢測中的應用

田素凱， 寧 濤， 陳志同

(北京航空航天大學機械工程及自動化學院，北京 100191)

基于曲面精確表示的距離極值點的計算及在刀具干涉檢測中的應用

田素凱， 寧 濤， 陳志同

(北京航空航天大學機械工程及自動化學院，北京 100191)

針對基于曲面精確表示的剛體碰撞檢測中裁剪曲面距離極值點的求解問題，提出了利用平面向量場估計初始曲面距離極值點的方法，避免了曲面過度細分，討論了距離極值點滿足的微分幾何條件，給出了解析曲面/參數曲面、參數曲面/參數曲面、點/參數曲面和曲線/參數曲面的距離極值點迭代算法。實例驗證分析了該算法的高效性和可靠性。

碰撞檢測；平面向量場；距離極值點；迭代

碰撞檢測(collision detection)是指判斷多個剛體(零件、機械設備等)所占據的空間是否有重疊區域，一直是計算機仿真、干涉檢查、數控加工、機器人控制和虛擬現實等學科的關鍵技術?，F有文獻中的碰撞檢測算法分為基于圖形[1]和基于圖像[2]兩大類檢測算法?；趫D形的碰撞檢測包括模型分解(modal decomposition)和空間分解(space decomposition)兩種，前者常用方法有 AABB(axis aligned bounding box)層次樹[3]、包圍球(bounding ball)層次樹、OBB(oriented bounding box)層次樹和K-DOPs法(discrete orientation polytopes)[4]等；后者利用GPU硬件的加速功能對多面體模型進行凸分解，進而提高碰撞檢測的效率。這些算法都是基于模型離散或空間分解的方法，在精確性和實時性方面有一定的欠缺，無法解決最大過切、欠切位置計算等問題。近年來，針對曲面精確表示

的剛體碰撞檢測算法的研究有一定的進展，Choi[5]給出了空間橢球體的碰撞條件；Odehnal[6]研究了兩個圓環面的共線法矢；Lopes等[7]提出了二次曲面和超二次曲面的碰撞檢測數學框架；Zeng等[8]給出了支持多種幾何體的碰撞檢測算法庫。由于CAD系統支持的曲面有十多種之多，而以上方法僅適用于幾種簡單解析曲面，因此基于曲面精確表示的剛體碰撞檢測算法尚待研究。

在曲面精確表示的剛體碰撞檢測算法中，曲面之間的距離極值點對應最大干涉(過切、欠切)位置，相關計算是碰撞檢測算法的關鍵?，F有曲面間最小距離的求解主要有離散法和解析法，算法的實現效率及其精確度是兩個重要的性能指標。Gilbert和Foo[9]的算法首先將曲面離散為三角片，然后再計算兩三角片的最近距離，進而獲取兩曲面的近似距離。當精度要求比較高時，該算法的效率較低。Johnson和Cohen[10]提出了一種計算曲面間最小距離的框架，該框架以包圍盒細分為基礎，算法效率較低。Sohn等[11]提出了一種線性幾何(line geometry)的方法求解曲面間距離的計算，該算法減少方程組和未知量的數量，提高了算法效率，尤其適用于求解二次曲面距離極值點。劉浩等[12]提出了一種計算雙二次 NURBS曲面的最近距離點。文獻[11-12]都沒有給出一般曲面間距離求解的具體算法。陳麗萍等[13]從極值點的幾何條件出發，利用賦范空間投影法迭代，求取自由曲線到曲面的極值距離，但該算法并沒有給出初始點的計算方法。

本文結合平面向量場提出了一種求解裁剪曲面的初始距離極值點的算法[14-15]，該算法利用裁剪曲面的三角剖分網格，依據平面向量場的Poincare公式，估算出單張三角片內部是否有距離極值點。分別針對解析曲面/參數曲面、參數曲面/參數曲面等情況，給出了距離極值點迭代求精算法。同時，還討論了點到曲面距離、曲線到曲面距離的求解算法。

1 初始距離極值位置點的估計

在數控加工幾何仿真等領域，需要計算平面、二次曲面及圓環面等代數曲面與參數曲面距離極值位置。如何利用裁剪曲面的三角網格、基于平面向量場技術估計初始距離極值點，需引入幾個基本概念：

(1) 平面向量場[15]v是連續函數 v:D →R2，即對于平面區域D中的任何點p，v定義平面向量v(p)與之對應。向量場在自然界廣泛存在，例如引力場、重力場、電磁場和流體的速度場等。在微分動力系統中，微分方程的穩定性用向量場刻劃[16]。對于平面區域D中一點 p，若v(p)=0，則稱p為平面向量場v的臨界點(critical point)。假設p為平面向量場v的孤立臨界點，N為p的領域，且在N及其邊界?N上沒有其他臨界點。對于邊界?N 上的每個點q，平面向量場v都定義非0向量v(q)與之對應。假設一個手持指南針的人，沿逆時針方向在邊界上行走，指南針總保持與向量場同向。當回到出發點時，指南針亦回到初始位置。所以在此過程中，指南針相對于底盤必轉動整數圈。逆時針轉一圈記為1，順時針轉一圈記為–1。逆時針和順時針所轉圈數的代數和稱為旋轉數，如圖 1所示。向量場的旋轉數可以用 Poincare公式[16]計算：，其中C是平面區域中的封閉曲線，在 C上向量場▽φ = (φu,φv)非退化。有定理：設v是連續的平面向量場，C為 v上封閉簡單光滑曲線，若旋轉數W(v，C)不為0，則閉曲線C內至少有一個臨界點。臨界點和旋轉數反映向量場v的拓撲結構。當向量場v發生微小變動時，其拓撲結構一般保持不變，所以向量場的拓撲結構被認為是一種穩定的結構，能反映向量場的內在性質。

圖1 平面向量場的旋轉數

Kriezis等[14]用有向距離函數定義平面向量場。如圖2所示，假設點p是曲面 S1(u,v)上的任意一點，則q為曲面 S2(s,t)上距離點p最近的一點。當曲面 S1(u,v)上的點 p變動時，得到映射M(p )=q。對于曲面 S1(u,v)上任意一點p，S1(u,v)

有向距離函數定義為：φ(u,v) =(n2(M (p)), p -M(p) )，其中( , )表示內積n2，n2表示曲面S2(s,t)的單位法向量。有向距離函數 φ( u,v)的梯度定義了u、v參數域上的向量場，在 φ( u,v)的有效定義域內，有 ▽φ (u,v) =((n2,S1u),(n2,S1v))。所以對于 φ( u,v)的任何臨界點(u, v)，有 (n2,S1u) =0, (n2 ,S 1v)=0。

圖2 有向距離函數的定義

(2) 討論距離極值點的估算方法。在以上有向距離函數的定義中，假設 S1(u,v)、 S2(s,t)是裁剪參數曲面， ?1=(V1,T1)、?2=(V2,T2)分別是兩裁剪曲面的三角片集合，其中 Vi表示所有頂點集合，Ti表示所有三角片集合，i=1,2。利用 S2(s,t)的三角片集合，使用點/三角片最近距離算法計算出頂點集合 V1中每個頂點對應的 φ(u,v)。對于三角片集合T1中的任何三角形，假設v1,v2,v3是其3個頂點對應的平面向量場φ(u,v)的3個矢量，如果對某個i=1,2,3成立，則該三角形內可能存在距離極值點，否則計算v1和v2，v2和v3，及v3和v1的夾角之和α，若，則該三角形內可能存在距離極值點。這里“v1和 v2夾角”是指向量v1按逆時針轉到另向量v2的夾角，如圖3所示。

如果 S2(s,t)退化為一個三維點，則無須利用?2=(V2,T2)，可直接計算出有向距離φ( u,v)在三角片網格 ?1=(V1,T1)上的取值，進而得到對應的平面向量場；如果 S2(s,t)是二次曲面或圓環面，則由于 S1(u,v)上的任意點p到 S2(s,t)上最近點q容易計算，所以可以用解析方法得到▽φ(u,v) =((n2,S1u),(n2,S1v))。后續對距離極值點的估計同上給出的參數曲面到參數曲面距離極值點的算法。

圖3 向量場中向量夾角的計算

2 距離極值點的迭代算法

以上提出了一種針對代數曲面和裁剪曲面的初始極值位置點的估計方法。表 1針對具體的代數曲面，構造對應的迭代公式，可實現初始極值位置點的迭代求精。

說明：

(1) 圖示中S(u,v)代表參數曲面S(u,v)，n(u,v)代表曲面上某點處的法向量，G(t)為參數曲線。

(2) 對于球面與參數曲面的情況，給定球心為O(x0,y0,z0)，及參數曲面S(u,v)，曲面上各點處單位法矢為 n( u,v)。曲面 S(u,v)的法矢方程為L(s) = S(u,v) +s ·n( u,v)，要找到這樣的u, v參數，使L(s)通過球心O。在L(s)上找一點距離O最近，由于

所以L(s)與球心O的最近距離為 (V,V )-(n,V )2，其中 V =S-O?；谝陨戏治隽畹繕撕瘮禐閒(u,v) = (V,V ) -(n,V )2，求導得

表1 迭代公式

(3) 對于圓柱面與參數曲面的情況，令給定的圓柱面軸線為 P(t)=P+tD，其中P是圓柱軸線上一點 (x0,y0,z0)，D是圓柱軸線單位方向，給定參數曲面 S(u,v)，曲面上各點處單位法矢為 n(u,v)。曲面 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)上兩點距離之平方，令

其中， V =S-P ，對g(s,t)求偏導得

所以兩直線的最近距離點對滿足

解此方程得

其中， Δ=1 -(D,n)2，將上式中的s, t代入g(s, t)得直線L(s)與P(t)的最近距離為

所以迭代目標函數定義為

f(u,v)為L(s)與P(t)最近距離的平方加上L(s) 與P(t)夾角余弦的平方。實際上，當(n,D )→0時，Δ→ 1，于是迭代目標函數 f(u,v)可簡化為

對于 f(u,v)求偏導得

利用上面的導數構造迭代公式：

(4) 對于圓環面與參數曲面情況，令給定的圓環面中心圓方程為： Q (θ) =O+ rX cosθ + rY sinθ，給定參數曲面 S(u,v)，曲面上各點處單位法矢為n( u,v)。由于圓環面與參數曲面的極值位置點對連線垂直于曲面一階偏導數Su和Sv，且與中心圓的切矢垂直，所以定義目標函數為

3 距離極值點算法應用

在數控加工幾何仿真領域，為計算過切量和確定過切位置，需要計算平面、二次曲面及圓環面等代數曲面與參數曲面極值位置。本文實現了平面、圓柱面、球面和圓環面對參數曲面的極值位置計算算法，并結合平面向量場技術和裁剪曲面剖分，實現了三坐標加工過切仿真功能。測試實例為：B樣條曲面u, v參數的次數都是5次，控制頂點是 30×30個，x, y, z最小最大范圍分別為[–134,134]×[–133,136]×[–172,77]，單位為mm，如圖4所示。以φ10球頭刀為刀具采用變步長方法計算刀位軌跡，其逼近 B樣條曲面等距面的精度為0.001，刀位點總計64 976個。裁剪曲面的三角片總計8 818個，三角片結點共4 458個。

圖4 B樣條曲面及其裁剪部分

在仿真算法中，兩個相鄰刀位點之間的刀具包絡體分為球體和圓柱體兩部分處理，如圖 5所示。由于此算法中三角片和刀位點數據量很大，為提高搜索與兩個相鄰刀位點之間的刀具包絡體鄰近的三角片，在加工坐標系xoy坐標面上劃分正方形網格單元(邊長等于刀具半徑)，構建三角片相鄰關系索引表。該方法大幅減少了三角片的查詢計算量和時間。

圖5 單條插補段的刀具運動包絡體

仿真環境是Windows 7，處理器是Intel Core Duo P8600，主頻是2.4 GHz，內存2 G。仿真結果如圖6~7所示。最大過切量是0.000 86 mm，耗時

46.4 s，見表2，過切基本出現在參數曲面曲率大的位置，且為刀具包絡的圓柱部分與零件曲面干涉，刀具包絡的球面部分與零件無干涉。仿真結果表明本文實現的代數曲面/參數曲面極值位置計算算法穩定可靠，能夠實現加工過切仿真等CAD/CAM功能。

圖6 裁剪曲面及其刀位軌跡

圖7 放大10 000倍的過切情況

表2 測試結果

4 結 論

現有成熟的碰撞檢測算法都是基于模型離散或空間分解的算法，難以解決給出最大干涉量和干涉位置等問題，而基于曲面精確表示的碰撞檢測則可以解決這樣的問題，但相關算法更為復雜。本文研究了點、線、面到裁剪曲面的距離極值點求解問題，給出了一種通用算法思路。該算法基于裁剪曲面的三角片網格、利用平面向量場估算初始距離極值點，再通過Newton迭代精確的距離極值點。對于二次曲面及圓環面/參數曲面的情況，由于二次曲面及圓環面的幾何特性，相關平面向量場以及Newton迭代的計算大為簡化。后續研究內容主要包括二次曲面及圓環面之間的距離極值點的求解問題，以及基于曲面精確表示的碰撞檢測系統的框架設計。

[1] Lin M C, Gottschalk S. Collision detection between geometric models: a survey [C]//Proceedings of IMA Conference on Mathematics of Surfaces, Cambridge: Andalucia En La Historia, 1998: 37-56.

[2] 范昭煒, 萬華根, 高曙明. 基于圖像的快速碰撞檢測算法[J]. 計算機輔助設計與圖形學學報, 2002, 14(9): 805-809.

[3] 陳 華. 確定任意形狀物體最小包圍盒的一種方法[J].圖學學報, 2010, 31(2): 49-53.

[4] Klosowski J T, Held M, Mitchell J S B, et al. Efficient collision detection using bounding volume hierarchies of k-DOPs [J]. IEEE Transactions on Visualization & Computer Graphics, 1998, 4(1): 21-36.

[5] Choi Y K. Collision detection for ellipsoids and other quadrics [D]. Hong Kong: The University of Hong Kong, 2008.

[6] Odehnal B. Common normals of two tori [J]. Journal for Geometry & Graphics, 2005, (9): 51-65.

[7] Lopes D S, Silva M T, Ambrósio J A, et al. A mathematical framework for rigid contact detection between quadric and superquadric surfaces [J]. Multibody System Dynamics, 2010, 24(3): 255-280.

[8] Zeng L, Ning T, Xi P. Research and application of general collision detection simulation platform [M]//The 19th International Conference on Industrial Engineering and Engineering Management. Berlin: Springer Berlin Heidelberg, 2013: 1315-1323.

[9] Gilbert E G, Foo C P. Computing the distance between general convex objects in three-dimensional space [J]. IEEE Transactions on Robotics & Automation, 1990, 6(1): 53-61.

[10] Johnson D E, Cohen E. A framework for efficient minimum distance computations [C]//IEEE International Conference on Robotics & Automation. New York: IEEE Press, 1998: 3678-3684.

[11] Sohn K A, Juttler B, Kim M S, et al. Computing distances between surfaces using line geometry [C]// Computer Graphics and Applications, 2002. Proceedings. Pacific Conference on. New York: IEEE Press, 2002: 236-245.

[12] 劉 浩, 唐月紅, 廖文和. 雙二次NURBS曲面間的最短距離[J]. 計算機輔助設計與圖形學學報, 2003, 15(10): 1298-1302.

[13] 陳麗萍, 陳 燕, 胡德金. 一種快速完備的自由曲線和曲面間最短距離求取算法[J]. 上海交通大學學報, 2003, 37(S1): 41-44.

[14] Kriezis G A, Patrikalakis N M, Wolter F E. Topological and differential-equation methods for surface intersection [J]. Computer-Aided Design, 1992, 24(1): 41-55.

[15] 寧 濤, 馬德昌, 王亞平, 等. 平面向量場與曲率分析在曲面求交中的應用[J]. 計算機學報, 1997, 20(12): 1074-1080.

[16] 張錦炎, 錢 敏. 微分動力系統導引[M]. 北京: 北京大學出版社, 1991: 41-47.

On Computation of Distance Extremum Points Based on Exact Surface Representation and Its Application in Tool Interferance Detection

Tian Sukai, Ning Tao, Chen Zhitong

(School of Mechanical Engineering and Automation, Beihang University, Beijing 100191, China)

Aiming at the problem of solving distance extreme points on trimmed surfaces in collision detection of rigid bodies which are represented by exact surface, a method is presented to obtain the initial distance extreme points using the plane vector field, which can avoid the excessive subdivision of the surface. The differential geometry conditions of the distance extreme points are discussed, and iterative formulas for obtaining distance extreme points of the analytic surface/parametric surface and parametric surface/parametric surface are given in this paper. The efficiency and reliability of the algorithm are verified by an example.

collision detection; plane vector field; distance extreme points; iteration

V 260.5；TP 391

10.11996/JG.j.2095-302X.2016050620

A

2095-302X(2016)05-0620-06

2016-03-24；定稿日期：2016-05-13

國家科技重大專項–高檔數控機床與基礎制造裝備科技重大專項課題(2013ZX04011031)

田素凱(1989–)，男，山東棗莊人，碩士研究生。主要研究方向為計算機輔助設計、CAD技術。E-mail：tskdqq@163.com

寧 濤(1966–)，男，遼寧丹東人，副教授，博士。主要研究方向為計算機輔助設計、CAD/CAM技術。E-mail：ningtao@buaa.edu.cn