NURBS曲線矩陣表示及其在插補中的應(yīng)用①

徐水龍，徐周波，古天龍

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

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

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

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

1 NURBS曲線的矩陣表示形式

1.1 B樣條曲線的表示式

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

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

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

于是可通過判斷參數(shù)u所處區(qū)間，對曲線進行分段處理。故式（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樣條方法進行改進，提出了非均勻有理B樣條（NURBS）方法。k次NURBS曲線C（u）的表達式

其中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）。文獻［5］基于式（2）給出了系數(shù)矩陣元素的遞推關(guān)系式，從關(guān)系式中可看出，每個高階元素由4項低階元素的加權(quán)項相加而得。但在遞推過程中，這些加權(quán)項的權(quán)系數(shù)存在重復(fù)計算，且在計算系數(shù)矩陣中位于矩陣邊緣的元素時，需要判斷這4個加權(quán)項中是否含有為0的項。為彌補上述不足，通過對該遞推關(guān)系式中的權(quán)系數(shù)進行變量代換，進一步對得到的關(guā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）兩邊進行多項式展開，由等式兩邊含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時

為進一步得出系數(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）中元素的影響進行組合疊加，可得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階方陣相加，因而存在大量的加法運算，尤其是零加零的運算。為省去這些零加零運算，降低算法的時間復(fù)雜度，可作如下優(yōu)化：先計算并存儲對Ml（i）中各個元素的影響，再將對Ml（i）中各個元素的影響累加到先前存儲的元素中，依此類推，便可遞推出系數(shù)矩陣Ml（i）。

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

為實現(xiàn)上述算法，首先需要定義并初始化曲線的相關(guān)參數(shù)（曲線次數(shù)k、控制點個數(shù)n＋1、節(jié)點矢量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）中的各項用累加的方法表示，

／／以省去零加零運算

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；

｝

｝

｝

｝

｝

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

3 算法效率的比較

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

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

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

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

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

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

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

4.1 NURBS曲線上點的計算

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

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

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

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

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

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

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

5 結(jié)束語

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

［1］施法中.計算機輔助幾何設(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］.計算機學(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］.計算機學(xué)報，1996，19（12）：941－947.

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

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

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