徐水龍，徐周波，古天龍

（1.桂林電子科技大學(xué) 機(jī)電工程學(xué)院，廣西 桂林 541004；2.桂林電子科技大學(xué) 計算機(jī)科學(xué)與工程學(xué)院，廣西 桂林 541004）

隨著計算機(jī)輔助圖形設(shè)計理論的發(fā)展，非均勻有理B樣 條（non－uniform rational B－spline，簡 稱NURBS）技術(shù)以其獨(dú)特的優(yōu)點(diǎn)，從眾多形狀描述曲線曲面中脫穎而出，成為定義工業(yè)產(chǎn)品幾何形狀的唯一數(shù)學(xué)方法，并被廣泛應(yīng)用于CAD／CAM中［1］。但由于控制算法的不足及CPU處理能力的限制，NURBS技術(shù)在CAM中的應(yīng)用不盡人意，其中NURBS曲線曲面的表示問題顯得尤為突出。

在NURBS曲線曲面插補(bǔ)中，為完成一維參數(shù)空間到三維坐標(biāo)空間的映射，需計算大量插補(bǔ)點(diǎn)的坐標(biāo)；同時為了控制刀具的速度，需對曲線上各點(diǎn)的曲率進(jìn)行計算。傳統(tǒng)的Deboor算法［2］在求解NURBS曲線上點(diǎn)的坐標(biāo)和導(dǎo)數(shù)的遞推過程中存在大量的重復(fù)計算，從而增加了許多額外的時間開銷。近年來，矩陣?yán)碚摰陌l(fā)展和計算機(jī)對矩陣運(yùn)算處理水平的提高，使得NURBS的矩陣表示算法越來越受到國內(nèi)外學(xué)者的重視。Choi等［3］和王學(xué)福等［4］均提出了用于計算NURBS曲線的符號矩陣表示算法，但都未給出任意階NURBS曲線顯式矩陣表達(dá)式的遞推過程。Grabowski等［5］給出了B樣條基的系數(shù)矩陣元素間的遞推公式，但并沒有給出系數(shù)矩陣的矩陣遞推公式。秦開懷［6］采用Toeplitz矩陣，導(dǎo)出了任意階NURBS曲線的顯式矩陣遞推公式，但該公式中所含的矩陣相乘運(yùn)算將增加乘法的計算量，且算法不易于用代碼實(shí)現(xiàn)。潘日晶等［7］引入差商的概念，提出了NURBS曲線的顯式矩陣表示算法，但該算法的推導(dǎo)過程較為繁瑣，時間復(fù)雜性隨NURBS曲線階數(shù)的增加而大大增加，不利于保證直接插補(bǔ)的實(shí)時性。王國勛等［8］通過對B樣條的導(dǎo)數(shù)公式進(jìn)行逐階積分，給出了系數(shù)矩陣的元素遞推關(guān)系，但該遞推關(guān)系分類復(fù)雜，同樣不利于用代碼實(shí)現(xiàn)。

鑒于此，基于文獻(xiàn)［5］提出一種新型的系數(shù)矩陣遞推公式。該算法具有計算量小、計算穩(wěn)定、易于用代碼實(shí)現(xiàn)等優(yōu)點(diǎn)，可用于任意階NURBS曲線和B樣條曲線的矩陣表示及其升階、降階。在NURBS曲線直接插補(bǔ)可通過該算法進(jìn)行預(yù)處理，求出曲線的系數(shù)矩陣并存儲，便于后續(xù)軌跡規(guī)劃和速度控制時直接調(diào)用該系數(shù)矩陣進(jìn)行曲線求值求導(dǎo)，從而保證了插補(bǔ)的實(shí)時性。

1 NURBS曲線的矩陣表示形式

1.1 B樣條曲線的表示式

k次B樣條曲線的表示式為：

其中：u∈［0，1］；di＝（xi，yi，zi），i＝0，1，…，n為曲線的控制頂點(diǎn)；Ni，k（u）為由節(jié)點(diǎn)矢量U根據(jù)Deboor－Cox遞推公式（式（2））決定的k次規(guī)范B樣條基函數(shù)，簡稱B樣條基。節(jié)點(diǎn)矢量U＝［u0，u1，…，un＋k＋1］，0≤ui≤1且ui≤ui＋1是經(jīng)過規(guī)范化處理［7］的節(jié)點(diǎn)矢量。

根據(jù)B樣條基的局部支承性，即當(dāng)u∈［ui，ui＋1）時，有

于是可通過判斷參數(shù)u所處區(qū)間，對曲線進(jìn)行分段處理。故式（1）可寫成如下分段形式：

其中，由于任意l次B樣條基u∈［ui，ui＋1）可展開成如下冪級數(shù)形式：

其中：（i）（與參數(shù)u無關(guān)，但與i有關(guān)）為l次B樣條基Ni－h(huán)，l中uj項的系數(shù)；l＝1，2，…，k；h＝l，l－1，…，0；j＝0，1，…l。

記

則式（3）可寫成

其中，（k＋1）×（k＋1）階矩陣Mk（i）稱為k次B樣條基的系數(shù)矩陣。

1.2 NURBS曲線的表示式

由于專門用于自由型曲線曲面設(shè)計的B樣條方法無法滿足構(gòu)造初等曲線曲面的要求，于是學(xué)者們對B樣條方法進(jìn)行改進(jìn)，提出了非均勻有理B樣條（NURBS）方法。k次NURBS曲線C（u）的表達(dá)式

其中u∈［0，1］，wi（i＝0，1，…，n）稱為權(quán)因子，其他各量的含義與式（1）相同。由式（6）可知，當(dāng)wi＝1時，NURBS曲線便成為一條B樣條曲線?；谑剑?）的推導(dǎo)過程，式（6）可寫成如下矩陣形式：

其 中：Q（i）＝（wi－kdi－k，wi－k＋1di－k＋1，…，widi）T；W（i）＝（wi－k，wi－k＋1，…，wi）T；u∈［ui，ui＋1）。

2 系數(shù)矩陣的遞推算法

由式（7）可見，k次NURBS曲線的矩陣表示的關(guān)鍵是求其系數(shù)矩陣Mk（i）。文獻(xiàn)［5］基于式（2）給出了系數(shù)矩陣元素的遞推關(guān)系式，從關(guān)系式中可看出，每個高階元素由4項低階元素的加權(quán)項相加而得。但在遞推過程中，這些加權(quán)項的權(quán)系數(shù)存在重復(fù)計算，且在計算系數(shù)矩陣中位于矩陣邊緣的元素時，需要判斷這4個加權(quán)項中是否含有為0的項。為彌補(bǔ)上述不足，通過對該遞推關(guān)系式中的權(quán)系數(shù)進(jìn)行變量代換，進(jìn)一步對得到的關(guān)系式進(jìn)行形式變換，找出某一低階系數(shù)矩陣元素對哪些高階系數(shù)矩陣元素中的哪些項有影響，提取受影響的項并組合疊加，最終推導(dǎo)出系數(shù)矩陣顯式遞推公式。

2.1 系數(shù)矩陣顯式遞推過程的推導(dǎo)

當(dāng)u∈［ui，ui＋1）時，令 式（4）中＝ui－h(huán)＋l－為變量代換后的權(quán)系數(shù)。由此可得：＝1－；則根據(jù)式（2）有：

結(jié)合式（4）、（8）并展開可得

其中：l＝1，2，…，k；h＝l，l－1，…，0；j＝0，1，…，l。

將式（9）兩邊進(jìn)行多項式展開，由等式兩邊含uj項的系數(shù)對應(yīng)相等，可得遞推關(guān)系：

其中

1）u∈［ui，ui＋1），＝1；

2）當(dāng)j＝0時，

3）當(dāng)j＝l時，

4）當(dāng)h＝0時，

5）當(dāng)h＝l時

為進(jìn)一步得出系數(shù)矩陣的直接顯式遞推公式，對式（10）作如下形式變換：

式（11）～（14）中，p、q＝0，1，…，l－1；l＝1，2，…，k。

從上可看出，l－1次B樣條基系數(shù)矩陣中的元素，僅對l次B樣條基系數(shù)矩陣中的、產(chǎn)生影響，于是可提取式（11）～（14）中包含的項。最后將系數(shù)矩陣Ml－1（i）中的各個元素對系數(shù)矩陣Ml（i）中元素的影響進(jìn)行組合疊加，可得u∈［ui，ui＋1）時Ml（i）的顯式遞推公式：

其中：

1）右邊的矩陣中其余元素均為0；

2）＝1；

6）p、q＝0，1，…，l－1；l＝1，2，…，k。

從式（15）可看出，等式右邊由l2個l＋1階方陣相加，因而存在大量的加法運(yùn)算，尤其是零加零的運(yùn)算。為省去這些零加零運(yùn)算，降低算法的時間復(fù)雜度，可作如下優(yōu)化：先計算并存儲對Ml（i）中各個元素的影響，再將對Ml（i）中各個元素的影響累加到先前存儲的元素中，依此類推，便可遞推出系數(shù)矩陣Ml（i）。

2.2 遞推算法的代碼實(shí)現(xiàn)

為實(shí)現(xiàn)上述算法，首先需要定義并初始化曲線的相關(guān)參數(shù)（曲線次數(shù)k、控制點(diǎn)個數(shù)n＋1、節(jié)點(diǎn)矢量U［n＋k＋2］）。將上述系數(shù)矩陣遞推算法寫成預(yù)處理函數(shù)，其C程序如下：

double M［k＋1］［k＋1］［k＋1］［n－2］＝0；

／／定義系數(shù)矩陣并初始化為0

void pretreat（）

｛

int p，q，L，i；

for（i＝0；i＜＝n－k；i＋＋）

｛／／對曲線分段處理

M［0］［0］［0］［i］＝1；

／／M［p］［q］［L］［i］即為文中的（i）

for（L＝1；L＜＝k；L＋＋）

｛for（q＝0；q＜＝L－1；q＋＋）

｛double beta，C，D；

beta＝U［i＋k＋q＋1］－U［i＋k＋q－L＋1］；

／／由于u0～uk均為0，故第i段曲線對應(yīng)于

／／u∈［ui＋k，ui＋k＋1）

if（beta＝＝0）

｛C＝0； D＝0；｝

else

｛D＝－1／beta； C＝U［i＋k＋q＋1］／beta；｝

for（p＝0；p＜＝L－1；p＋＋）

｛double var1，var2；

／／定義2個中間變量，以減少重復(fù)計算

var1＝C＊M［p］［q］［L－1］［i］；

var2＝D＊M［p］［q］［L－1］［i］；

／／將式（15）中的各項用累加的方法表示，

／／以省去零加零運(yùn)算

M［p］［q］［L］［i］＝M［p］［q］［L］［i］＋var1；

M［p］［q＋1］［L］［i］＝M［p］［q＋1］［L］［i］

＋M［p］［q］［L－1］［i］－var1；

M［p＋1］［q］［L］［i］＝M［p＋1］［q］［L］［i］

＋var2；

M［p＋1］［q＋1］［L］［i］＝M［p＋1］［q＋1］

［L］［i］－var2；

｝

｝

｝

｝

｝

由上可看出本算法的代碼實(shí)現(xiàn)簡單，過程清晰。而文獻(xiàn)［6］的算法涉及矩陣的相關(guān)運(yùn)算（如矩陣擴(kuò)行、雙對角矩陣的構(gòu)造、矩陣相乘等），故代碼實(shí)現(xiàn)將降低算法的執(zhí)行效率。利用文獻(xiàn)［8］的算法求系數(shù)矩陣中的元素，其遞推公式分為5類，且遞推公式需要調(diào)用同階系數(shù)矩陣的元素，不利于編程實(shí)現(xiàn)和提高算法效率。綜上所述，從代碼實(shí)現(xiàn)角度分析，本算法具有明顯優(yōu)勢。

3 算法效率的比較

設(shè)Tadd、Tsub、Tmul、Tdiv分別代表所需要的加、減、乘、除的運(yùn)算次數(shù)，并假定它們滿足如下關(guān)系：Tadd≈Tsub，Tmul≈Tdiv，Tmul≈1.12Tadd

［6］。將系數(shù)矩陣遞推算法效率與相關(guān)文獻(xiàn)的算法效率進(jìn)行比較，其各系數(shù)矩陣遞推算法時間復(fù)雜度見表1所示。

k次B樣條基與文獻(xiàn)［6］中的k＋1階B樣條基對應(yīng)。從表1可看出，本算法的時間復(fù)雜度優(yōu)于另外幾種算法。

為驗證本算法的高效性，將Deboor算法［2］、遞推矩陣算法［6］與本算法均用C語言代碼實(shí)現(xiàn)，然后在主頻為2.99GHz的計算機(jī)運(yùn)行比較。以計算5000次NURBS曲線上的點(diǎn)所需時間為基準(zhǔn)，得出3種算法的平均時間開銷分別為3.160 222、1.224 315、0.917 387ms。由此可見，本算法在效率上具有一定的優(yōu)勢，能為NURBS曲線的高速高精度插補(bǔ)提供有力保證。

表1 各系數(shù)矩陣遞推算法時間復(fù)雜度Tab.1 Time complexity of different coefficient matrix recursion algorithms

4 NURBS矩陣表示在直接插補(bǔ)中的應(yīng)用

近年來，隨著NURBS曲線曲面被廣泛應(yīng)用于產(chǎn)品設(shè)計中，NURBS曲線的直接插補(bǔ)技術(shù)也越來越受到國內(nèi)外學(xué)者的重視。2007年，王田苗等［9］將Deboor算法用于NURBS曲線的直接插補(bǔ)中，但由于Deboor算法在求值求導(dǎo)過程中存在大量重復(fù)遞推計算，與NURBS矩陣表示算法相比，不能通過預(yù)處理減小插補(bǔ)周期內(nèi)的計算量，難以保證插補(bǔ)的實(shí)時性。

4.1 NURBS曲線上點(diǎn)的計算

在NURBS曲線的直接插補(bǔ)中需要完成一維參數(shù)空間到三維坐標(biāo)空間的映射，計算大量插補(bǔ)點(diǎn)的坐標(biāo)，以得到下一個插補(bǔ)周期中刀具的進(jìn)給增量。可先調(diào)用上面的算法進(jìn)行預(yù)處理，求得k次B樣條基的系數(shù)矩陣Mk（i），進(jìn)而得到不同分段區(qū)間內(nèi)NURBS曲線的矩陣表示式。對給定參數(shù)u，要計算其在三維坐標(biāo)空間的映射點(diǎn)，只需判斷u所處的分段區(qū)間并調(diào)用相應(yīng)的系數(shù)矩陣，再將u的值代入式（7）即可求得NURBS曲線的坐標(biāo)（Cx，Cy，Cz）。

4.2 NURBS曲線上各點(diǎn)處的導(dǎo)數(shù)及曲率計算

在NURBS曲線的直接插補(bǔ)中，為滿足一定的加工精度和機(jī)床動力性能，常常需要根據(jù)曲線的曲率對刀具的速度、加速度等狀態(tài)進(jìn)行實(shí)時控制，這就需要計算曲線在某點(diǎn)處的導(dǎo)數(shù)。首先通過對式（7）進(jìn)行變換［8］可得到：）

其中，Mk（u）、Q（u）、W（u）與參數(shù)u無關(guān)，簡寫為Mk、Q、W。再對式（16）兩邊求r階導(dǎo)數(shù)：

由式（17）可遞推出曲線C（u）在某點(diǎn)處的r階導(dǎo)數(shù)C（r）（u）。一般在NURBS曲線插補(bǔ)中只涉及曲線的一階和二階導(dǎo)數(shù)，其計算公式為：

在插補(bǔ)過程中，由于進(jìn)給步長很小，可用圓弧近似各步長間的曲線，將該圓弧半徑作為對應(yīng)插補(bǔ)點(diǎn)的曲率半徑。利用式（18）、（19）可計算曲率半徑ρ，

此外，實(shí)時插補(bǔ)過程中還要不斷通過計算進(jìn)給步長確定下一時刻插補(bǔ)點(diǎn)的位置，通常采用的方法是一階或二階泰勒展開公式。式（20）同樣需要反復(fù)采用曲線在插補(bǔ)點(diǎn)處的導(dǎo)數(shù)信息，可見，NURBS曲線求值求導(dǎo)算法的效率對保證其插補(bǔ)實(shí)時性具有非常重要的作用。

5 結(jié)束語

基于文獻(xiàn)［5］的思想，將Deboor－Cox遞推公式中的B樣條基用冪函數(shù)表示，利用公式兩邊次數(shù)相同項的系數(shù)相等的關(guān)系，得出系數(shù)矩陣元素之間的遞推關(guān)系。本研究對該關(guān)系式進(jìn)行變量代換，避免了中間量的重復(fù)計算。再通過對該關(guān)系式進(jìn)行形式變換、元素提取、組合疊加，最終得出系數(shù)矩陣的顯式遞推矩陣公式，進(jìn)一步減小了計算量。與傳統(tǒng)的Deboor算法相比，本算法采用矩陣表示NURBS曲線，便于對曲線進(jìn)行預(yù)處理，避免了重復(fù)遞推，更有利于保證插補(bǔ)的實(shí)時性。與其他矩陣遞推算法相比，本算法的時間復(fù)雜度更優(yōu)，且易于用代碼實(shí)現(xiàn)。本算法不僅可用于任意階NURBS曲線和B樣條曲線的矩陣表示及其升階、降階，同時，還可為進(jìn)一步研究NURBS曲面表示和直接插補(bǔ)奠定基礎(chǔ)。

［1］施法中.計算機(jī)輔助幾何設(shè)計與非均勻有理B樣條［M］.北京：高等教育出版社，2001：211－310.

［2］Deboor C.A practical guide to splines［M］.New York：Springer－Verlag，1978：115－117.

［3］Choi B K，Yoo W S，Lee C S.Matrix representation for NURBS curves and surfaces［J］.Computer Aided Design，1990，22（4）：235－240.

［4］王學(xué)福，孫家廣，秦開懷.NURBS的符號矩陣表示及其應(yīng)用［J］.計算機(jī)學(xué)報，1993，16（1）：28－34.

［5］Grabowski H，Li X.Coefficient formula and matrix of nonuniform B－spline functions［J］.Computer Aided Design，1992，24（12）：637－642.

［6］秦開懷.NURBS曲線和曲面的遞推矩陣及其應(yīng)用［J］.計算機(jī)學(xué)報，1996，19（12）：941－947.

［7］潘日晶.NURBS曲線曲面的顯式矩陣表示及其算法［J］.計算機(jī)學(xué)報，2001，24（4）：358－366.

［8］王國勛，舒啟林，王軍，等.NURBS直接插補(bǔ)技術(shù)中快速求值求導(dǎo)算法［J］.東北大學(xué)學(xué)報：自然科學(xué)版，2012，33（7）：1021－1024.

［9］王田苗，曹宇男，陳友東，等.基 于de Boor算 法 的NURBS曲線插補(bǔ)和自適應(yīng)速度控制研究［J］.中國機(jī)械工程，2007，18（21）：2608－2613.