鄭志浩，汪國昭

（浙江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，浙江 杭州310027）

將式（9）及ω（t）的表達式代入式（10）和（11），整理后得到以下方程組：

在計算機輔助幾何設(shè)計中，圓錐曲線可用帶權(quán)因子的有理多項式表示，權(quán)因子的易控特點使得圓錐曲線具有很強的幾何靈活度，因而被廣泛用于航空和汽車工業(yè)設(shè)計、機械加工、路徑規(guī)劃及藝術(shù)字體的造型中.利用圓錐曲線進行幾何造型會頻繁涉及放樣與拼接問題，其中蘊含了大量的曲線求交工作，因此，須進行有理多項式與多項式的系統(tǒng)轉(zhuǎn)化，其主要方法就是用多項式對圓錐曲線進行高精度逼近［1-6］.同時，幾何造型會頻繁處理等距線的構(gòu)造及弧長計算問題.由于目前圓錐曲線（除圓弧外）的弧長或等距線往往不能用多項式或有理多項式表示，不兼容于現(xiàn)有的NURBS系統(tǒng)，因此，對圓錐曲線進行弧長表示的有理轉(zhuǎn)化是一種非常有效的方法.PH曲線是Farouki［7］引入的一種新穎的多項式曲線，即平面多項式參數(shù)曲線P（t）＝（x（t），y（t））是PH 曲線的充要條件如下：存在多項式σ（t）使得關(guān)系式x′2（t）＋y′2（t）＝σ2（t）成立.PH 曲線主要特征是其弧長和等距線可用多項式及有理多項式表示.自引入PH 曲線后，許多學(xué)者對其進行了廣泛研究，并用來對曲線作插值逼近.主要成果包括：由曲線端點條件插值生成G1的三次PH 曲線或PH 樣條［8-10］；Pelosi等［11］利用解非線性復(fù)方程的方法構(gòu)造C2五次PH 樣 條 曲 線；Moon 等［12］用 復(fù) 分 析 中 的Minkowski幾何代數(shù)法研究了五次PH 插值曲線的形狀與端點插值條件的關(guān)系；Sir等［13］由端點C1的Hermite插值條構(gòu)造五次PH 樣條曲線，逼近誤差階數(shù)為4；Juettler［14］由基曲線端點，基曲線端點一階導(dǎo)數(shù)及端點曲率為插值條件，插值方程通過適當?shù)淖鴺俗儞Q得到相關(guān)代數(shù)參量的四次方程，生成G2連續(xù)的七次PH 逼近曲線；Sir等［15］構(gòu)造了次數(shù)更高的九次PH 樣條曲線，逼近誤差階數(shù)為6；Wang等［16］對四次PH 曲線的控制多邊形的幾何結(jié)構(gòu)作了深入研究，給出了四次Bézier曲線成為PH 曲線的邊角分離的幾何條件.張偉紅等［17］利用五次PH曲線對圓弧作了等弧逼近.方林聰?shù)龋?8］在C1Hermite插值條件下構(gòu)造了2類六次PH 曲線，并能自由選擇奇異點的位置.楊平等［19-20］討論了七次Bézier曲線為PH 曲線的幾何條件.此外，PH 曲線還應(yīng)用于過渡曲線的構(gòu)造和逆向工程［21-25］.

在現(xiàn)有PH 插值及逼近結(jié)果基礎(chǔ)上，本文進一步探討較為新穎的四次PH 曲線及一類特殊的五次PH 曲線（具有保弧長的特點）的插值生成算法.分別從幾何結(jié)構(gòu)及復(fù)分析的角度具體構(gòu)造G1插值的PH 曲線來對圓錐曲線作逼近.由于大量涉及圓錐曲線逼近的研究所采用誤差度量均是Hausdorff距離［1-6］，本文同樣給出這些逼近的Hausdorff距離誤差估計.

1 圓錐曲線及其逼近誤差定義

平面內(nèi)的圓錐曲線b（t）＝（bx（t），by（t））可用二次有理Bézier曲線表示為

則X＝τ0b0＋τ1b1＋τ2b2，τ0＋τ1＋τ2＝1.

定義1 同一參數(shù)區(qū)間［0，1.0］內(nèi)的2條平面曲線p（t）及q（t）的Hausdorff距離定義為

Hausdorff距離已成為計算曲線間誤差的最常用的測度方式.特別針對圓錐曲線的逼近，F(xiàn)loater［4］還給出以下結(jié)論.

定理1 包含于Δb0b1b2的連續(xù)曲線r（t），t∈［0，1］，滿足r（0）＝b0，r（1）＝b2，｛τj｝2j＝0是r（t）關(guān)于Δb0b1b2的三重心坐標，則r（t）與圓錐曲線b（t）的Hausdorff誤差距離上界有如下估計式：

式中：

2 用四次PH 逼近圓錐曲線

對式（1）所表示的圓錐曲線，本章將根據(jù)其端點及端點單位切向量來構(gòu)造G1四次PH 插值曲線，把該PH 曲線作為圓錐曲線的插值逼近曲線，并利用定理1進行誤差估計.

2.1 G1 四次PH 曲線的插值構(gòu)造

在復(fù)平面中，四次Bézier曲線的參數(shù)多項式P（t）＝（x（t），y（t））可 表 示 成 復(fù) 數(shù) 形 式P（t）＝x（t）＋i y（t），P（t）是四次PH 曲線的充分必要條件是：存在一次實多項式w（t）＝α（1-t）＋t及一次復(fù)多項式G（t）＝β（1-t）＋γt（α為實系數(shù)，β與γ 為復(fù)系數(shù)），使得P′（t）＝w（t）G2（t）.基于復(fù)分析的方法，Wang等［16］給出四次PH 曲線的控制多邊形幾何特征，如定理2所述.

定理2 四次Bézier曲線

圖1 四次Bézier曲線的控制多邊形Fig.1 Control polygon of quartic Bézier curve

圖2 圓錐曲線及四次PH 曲線的控制多邊形Fig.2 Conic section and control polygon of PH quartic

式中：

因此，P（t）的控制頂點可用上述有關(guān)邊長及角度幾何量表示：

式中：λ可自由選取.當λ為某一特定值時，只要求出L＊0、L＊1、L＊2就能同時算出L0、L3、R1、R3，由式（4）產(chǎn)生四次PH 曲線的控制頂點坐標.為此，利用插值條件P4＝b2，對應(yīng)以下方程組：

因此，在設(shè)定λ值后，能直接解出L＊0、L＊1、L＊2，得到PH 曲線的控制頂點.

由上述邊長L＊0、L＊1、L＊2的關(guān)系可知，k＝1 即λ＝3，就有

若取P＊0＝b0、P＊1＝E、P＊2＝F、P＊3＝b2，所構(gòu)造的三次Bézier曲線

四次PH 曲線可退化成三次PH 曲線.

2.2 插值逼近的誤差

當以圓錐曲線的端點及端點單位切向為插值條件構(gòu)造了四次PH 逼近曲線后，進一步來估計逼近誤差.不妨設(shè)

則四次PH 曲線的控制頂點｛Pj與點b0、b1、b2就有以下線性組合關(guān)系：

式中：n0、n1、n2是 點P2關(guān) 于Δb0b1b2的 重 心 坐 標，通過具體計算得到

將式（7）代入

整理后化成τ0b0＋τ1b1＋τ2b2的線性組合形式，系數(shù)｛τj就是P（t）關(guān)于Δb0b1b2的重心坐標，計算結(jié)果如下：

進一步計算可得

這里Bernstein基函數(shù)的系數(shù)分別為

總結(jié)以上的分析及結(jié)論，得到如下誤差定理和推論.

定理3 圓錐曲線b（t）與由其端點G1Hermite插值產(chǎn)生的四次PH 曲線P（t）的Hausdorff距離誤差有以下估計：

對于以上估計式，為獲得最大值，可能會涉及解六次多項式方程，但精度較高.如為了避免解方程，由定理3的結(jié)論進一步可得以下推論.

特別地，當λ＝3時，插值構(gòu)造的四次PH 曲線退化為由式（6）所表示的三次PH 曲線.誤差計算可通過次數(shù)較低的三次PH 曲線P＊（t）計算.設(shè)

那么，P＊（t）關(guān)于Δb0b1b2的重心坐標為τ0＝3（1-u）（1-t）2t＋（1-t）3，τ1＝3u（1-t）2t＋3v（1-t）t2，τ2＝t3＋3（1-v）（1-t）t2.

考慮函數(shù)

式中：

總結(jié)以上過程又有以下推論.

推論2 當λ＝3時，四次PH 曲線P（t）對b（t）的逼近退化為三次PH 曲線（6）對b（t）逼近.

兩者間的Hausdorff距離誤差估計如下：

3 用五次PH 曲線逼近圓錐曲線

復(fù)數(shù)形式的五次Bézier曲線為

若點P0確定，則其余控制頂點為

3.1 逼近的誤差估計

在△b0b1b2內(nèi)構(gòu)造插值圓錐曲線端點及其單位切向的G1五次PH 曲線.根據(jù)要求可知，點P1、P4應(yīng) 分 別 在b0-b1和b1-b2的 邊 上，點P2、P3為△b0b1b2的內(nèi)點，因此，PH 曲線的控制頂點與點b0、b1、b2存在以下線性關(guān)系：

式中：

總結(jié)以上結(jié)果得到定理4.

定理4 圓錐曲線b（t）與插值逼近的五次PH

曲線P（t）的Hausdorff距離估計式為

定理4的最大值可通過八次多項式方程求解得到，精度較高.同樣分析定理4 的結(jié)論進一步得到推論3.

3.2 插值曲線的類型

要構(gòu)造對圓錐曲線插值逼近的G1五次PH 曲線，即只要確定復(fù)系數(shù)｛ωj｝2j＝0，用定理4或推論3可估計逼近誤差.關(guān)于復(fù)系數(shù)的確定可采用以下2種方式.

1）基于C1插值逼近.

為使插值產(chǎn)生的五次PH 曲線在圓錐曲線點達到C1，參照Moon［12］的結(jié)論，取系數(shù)

2）G1等弧方式逼近圓錐曲線.

為了得到更高的逼近精度，主要采用等弧逼近的方式.如果逼近曲線與基曲線保持弧長的一致性，那么曲線的內(nèi)在幾何特征就能充分得以體現(xiàn).比如：當在路徑規(guī)劃設(shè)計中使用逼近算法時，基曲線與逼近曲線弧長的對應(yīng)就很重要，逼近曲線應(yīng)?；￠L.

在此要解決的問題是：由式（1）所表示的圓錐曲線的弧長為S，現(xiàn)根據(jù)其端點及其端點的單位切向量（一階幾何Hermite插值條件），構(gòu)造與之等弧長的五次PH 逼近曲線P（t）.參照圖2 的坐標系，以點b0為原點且b0的切向為x 軸正向，末端點b2的切向與x 軸的夾角為2θ.由于P（t）插值端點b0、b2及其單位切向T0、T1，故

從而ω0為實數(shù)，記為ω0，為減少計算參數(shù)，設(shè)定五次PH 曲線的兩端點切向量模長相等，即取ω2＝ω0exp（iθ），另外ω1＝ω1x＋iω1y.由P（t）插值兩端點b0、b2得到

及PH 曲線的弧長同為S，則

將式（9）及ω（t）的表達式代入式（10）和（11），整理后得到以下方程組：

由式（12）、（14）得到

取

由式（13）解得

再將ω1x、ω1y的表達式代入式（14），經(jīng)整理后產(chǎn)生了關(guān)于ω0的四次多項式方程：

解方程式（17），特別取

4 實 例

為了更好地說明本文中各類PH 插值曲線對圓錐曲線的逼近應(yīng)用，列舉如下2個實例.

圖3 橢圓及1／4段的控制頂點Fig.3 Ellipse and control polygon of quarter section

圖4 例1的誤差主項函數(shù)Fig.4 Key function for error in example 1

表1 用PH曲線逼近1／4橢圓段的Hausdorff距離誤差估計Tab.1 Hausdorff distance evaluation between quarter section of ellipse and PH quartics

圖5 離散后的兩等弧逼近曲線及控制多邊形Fig.5 Two approximation curves preserving arc-lengths and their control polygons after subdivision

若要獲得更高的逼近精度，可將表示橢圓段的二次有理Bézier曲線在肩點（t＝1／2）處離散成2段子曲線b1（t）與b2（t）（見圖5）：

b1（t）與b2（t）對應(yīng)的控制頂點及權(quán)因子為

如果對這兩離散段分別構(gòu)造五次PH 曲線作等弧逼近，則第一段誤差估計值為σ＝0.127 9×10-3，而第二段誤差估計值為σ＝0.001 2.由此可見，適當?shù)貙A錐曲線進行離散可以使逼近精度迅速提高.當離散段數(shù)較多時，誤差估計式采用推論1及推論3可避免多次解方程.

圖6 拋物線及控制頂點Fig.6 Parabola and control polygon

圖7 例2的誤差主項函數(shù)Fig.7 Key function for error in example 2

表2 例2的誤差主項函數(shù)Hausdorff距離誤差估計Tab.2 Hausdorff distance evaluation between parabola and PH quartics in example 2

5 結(jié) 語

采用四次PH 曲線及具有保弧長特征的五次PH 曲線可對圓錐曲線進行插值逼近.基于對控制多邊形邊角分離的幾何結(jié)構(gòu)分析，得到四次PH 曲線退化為三次PH 曲線的條件.利用不同類型的PH 曲線逼近圓錐曲線的Hausdorff距離誤差的估計結(jié)果，可選取適當類型的PH 曲線，在用戶給定的誤差范圍內(nèi)，將圓錐曲線轉(zhuǎn)化成PH 曲線.若要提高轉(zhuǎn)化精度，可利用二次有理Bézier曲線的中點離散公式將圓錐曲線離散成多段曲線，再分別對各離散段構(gòu)造插值逼近的PH 曲線，而對應(yīng)的PH 曲線段將組成C1、G1的樣條曲線.采用PH 曲線轉(zhuǎn)化的方法，不僅實現(xiàn)了圓錐曲線的多項式的表示，還能得到弧長的多項式表示形式及等距線的有理形式，對減少計算量、壓縮數(shù)據(jù)量都有一定的貢獻.為獲得更高連續(xù)階的逼近曲線，在今后工作中將研究用六次及七次PH 曲線構(gòu)造G2插值逼近曲線，并繼續(xù)將本文的逼近思想應(yīng)用到對其他類型曲線的逼近.

（

）：

［1］AHN Y J.Approximation of conic section by curvature continuous quartic Bézier curves［J］.Computers and Mathematics with Applications，2010，60（7）：1986-1993.

［2］FANG L.G3approximation of conic sections by quintic polynomial curves［J］.Computer Aided Geometric Design，1999，16（8）：755-766.

［3］FLOATER M.High-order approximation of conic sections by quadratic splines［J］.Computer Aided Geometric Design，1995，12（6）：617-637.

［4］FLOATER M.An O（h2n）Hermite approximation for conic sections［J］.Computer Aided Geometric Design，1997，14（2）：135-151.

［5］KIM S H，AHN Y J.An approximation of circle arcs by quartic Bézier curves［J］.Computer-Aided Design，2007，39（6）：490-493.

［6］AHN Y J.Helix approximations with conic and quadratic Bézier curves［J］.Computer Aided Geometric Design，2005，22（6）：551-565.

［7］FAROUKI R T.The conformal map z→z2of the hodograph plane ［J］.Computer Aided Geometric Design，1994，11（4）：363-390.

［8］LI Y J，DENG C Y.2012.C-shaped C2Hermite interpolation with circular precision based on cubic PH curve interpolation ［J］.Computer-Aided Design，2012：44（11），1056-1061.

［9］MEEK D S，WALTON D J.Geometric Hermite interpolation with Tschirnhausen cubics［J］.Journal of Computational and Applied Mathematics，1997，81（2）：299-309.

［10］BYRTUS M，BASTL B.G1Hermite interpolation by PH cubics revisited［J］.Computer Aided Geometric Design，2010，27（8）：622-630.

［11］PELOSI F，SAMPOLI M L，F(xiàn)AROUKI R T，et al.A control polygon scheme for design of planar C2PH quintic spline curves［J］.Computer Aided Geometric Design，2007，24（1）：28-52.

［12］MOON H P，F(xiàn)AROUKI R T.Construction and shape analysis of PH quintic Hermite interpolants［J］.Computer Aided Geometric Design，2001，18（2）：93-115.

［13］SIR Z，F(xiàn)EICHTINGER.F，JUETTLER B.Approximating curves and their offsets using baircs and Pythagorean hodograph quintic［J］.Computer-Aided Design，2006，38（6）：608-618.

［14］JUETTLER B.Hermite interpolation by Pythagorean hodograph curves of degree seven［J］.Mathematics of Computation，2000，70（235）：1089-1111.

［15］SIR Z，JUETTLER B.Constructing acceleration continuous tool paths using Pythagorean Hodograph curves［J］.Mechanism and Machine Theory，2005，40（11）：1258-1272.

［16］WANG G Z，F(xiàn)ANG L C.On control polygon of quartic Pythagorean-hodograph curves［J］.Computer Aided Geometric Design，2009，26（9）：1006-1015.

［17］張偉紅，蔡亦青，馮玉瑜.圓弧的五次PH 曲線等弧長逼近［J］.計算機輔助設(shè)計與圖形學(xué)報，2010，22（7）：1082-1086.ZHANG Wei-hong，CAI Yi-Qin，F(xiàn)ENG Yu-yu.Arclength-preserving approximation of circular arcs by quintic PH quintic PH curves［J］.Journal of Computer-Aided Design and Computer Graphics，2010，22（7）：1082-1086.

［18］方林聰，汪國昭.六次PH 曲線C1Hermite插值［J］.中國科學(xué)：數(shù)學(xué)，2014，44（7）：799-804.Fang Lin-cong，Wang Guo-zhao.C1Hermite interpolation using sextic PH curves［J］.SCIENCE CHINA Mathematics，2014，44（7）：799-804.

［19］楊平，汪國昭.7次PH 曲線的控制多邊形的幾何性質(zhì)［J］.計算機輔助設(shè)計與圖形學(xué)報，2014，26（3）：378-384.YANG Ping，WANG Guo-zhao.Geometric properties of control polygon of septic PH curve［J］.Journal of Computer-Aided Design and Computer Graphics，2014，26（3）：378-384.

［20］楊平，汪國昭.C3連續(xù)的七次PH 樣條閉曲線插值［J］.浙江大學(xué)學(xué)報：工學(xué)版，2014，48（5）：934-941.YANG Ping，WANG Guo-zhao.C3Spline interpolation by Pythagorean hodograph closed curve of degree seven［J］.Journal of Zhejiang University：Engineering Science，2014，48（5）：934-941.

［21］鄭志浩，汪國昭.三次PH 曲線的曲率單調(diào)性及過渡曲線構(gòu)造［J］.計算機輔助設(shè)計與圖形學(xué)報，2014，26（8）：1219-1224.ZHENG Zhi-hao，WANG Guo-zhao.On curvature monotony for a PH cubic curve and constructing transition curve［J］.Journal of Computer-Aided Design and Computer Graphics，2014，26（8）：1219-1224.

［22］WALTON D J，MEEK D S.G2curve design with a pair of Pythagorean hodograph quintic spiral segments［J］.Computer Aided Geometric Design，2007，24（5）：267-285.

［23］HABIB Z，SAKAI M.On PH quintic spirals joining two circles with one circle inside the other［J］.Computer-Aided Design，2007，39（2）：125-132.

［24］FAROUKI R T，GIANNELLI C，SESTINI A.Identification and“reverse engineering”of Pythagorean-hodograph curves［J］.Computer Aided Geometric Design，2015，34（1）：21-36.

［25］FAROUKI R T，SIR Z.Rational Pythagorean-hodograph space curves［J］.Computer Aided Geometric Design，2011，28（2）：75-88.