OpenCASCADE Conic to BSpline Curves-Parabola
eryar@163.com
Abstract. Rational Bezier Curve can represent conic curves such as circle, ellipse, hyperbola, .etc. But how to convert a conic curve to BSpline curve is still question, i.e. Represent a conic curve in BSpline form. Parabola curve is the most simple conic curve, that the parabola does not require rational functions. Let’s begin from the simplest one...
Key Words. OpenCASCADE, Convert, Parabola, BSplineCurve, Conic Curve
1. Introduction
圓錐截線（Conic或稱為二次曲線）和圓在CAD/CAM中有著廣泛應(yīng)用。毫無疑問NURBS的一個(gè)最大優(yōu)點(diǎn)就是既能精確表示圓錐截線和圓，也能精確表示自由曲線曲面。這個(gè)優(yōu)點(diǎn)的意義是方便編程，使所有的曲線可以采用統(tǒng)一的數(shù)據(jù)結(jié)構(gòu)來表示。通過有理的方式可以精確來表示這些二次曲線，那么給定一個(gè)二次曲線的相關(guān)參數(shù)（如圓的圓心和半徑等），如何構(gòu)造出對(duì)應(yīng)的NURBS曲線呢？
在圓錐截線中，拋物線（Parabola）是不需要用有理函數(shù)來表示的，所以是形式最簡單的二次曲線。先從簡單的著手，來學(xué)習(xí)如何將一個(gè)拋物線從隱式方程的形式轉(zhuǎn)換成NURBS曲線形式。
先簡要回顧一下高中數(shù)學(xué)中關(guān)于拋物線的相關(guān)知識(shí)點(diǎn)，如下圖所示為我上高中時(shí)數(shù)學(xué)課本中的關(guān)于拋物線方程及焦點(diǎn)坐標(biāo)（focus）和準(zhǔn)線（directrix）方程一個(gè)圖表：
Figure 1.1 Parabola Fuction
OpenCASCADE中對(duì)應(yīng)拋物線的隱式方程表示的類是gp_Parab/gp_Parab2d。本文主要介紹OpenCASCADE中如何將gp_Parab2d轉(zhuǎn)換為NURBS曲線。
2. Parametric Representations
在CAD/CAM的應(yīng)用中，圓錐截線有兩種重要的參數(shù)表示形式：有理形式和最大內(nèi)接面積形式（Rational and maximum inscribed area forms）。表示拋物線的這兩種形式是相同的，如下所示：
圓錐截線的有些有理參數(shù)表示形式可能是有相當(dāng)差的參數(shù)化，即均勻分布的參數(shù)值對(duì)應(yīng)于曲線上分布很不均勻的點(diǎn)。利用線性有理函數(shù)對(duì)有理曲線進(jìn)行重新參數(shù)化可以改變（因而可能改善）其參數(shù)化。
假設(shè)C(u)=(x(u), y(u))是一條在標(biāo)準(zhǔn)位置的圓錐截線的參數(shù)表示。現(xiàn)在我們對(duì)拋物線給出的參數(shù)方程也是上式，它是一個(gè)好的參數(shù)化：對(duì)于任意給定的整數(shù)n和參數(shù)邊界a與b，取n個(gè)等間隔分布的參數(shù)：
點(diǎn)C(u1)，C(u2), ..., C(un)形成曲線上n-1邊多邊形，它的閉合多邊形具有最大的內(nèi)接面積。
根據(jù)最大內(nèi)接面積表示法可以求出節(jié)點(diǎn)矢量。
3. Conversion Algorithm
將隱式表示的拋物線方程轉(zhuǎn)換為NURBS（有理Bezier是NURBS的特例）曲線需要確定NURBS的以下信息：節(jié)點(diǎn)矢量，權(quán)因子，次數(shù)，控制頂點(diǎn)。
因?yàn)閽佄锞€是二次曲線，所以對(duì)應(yīng)的NURBS曲線的次數(shù)也為2。因?yàn)槭怯糜欣淼腂ezier曲線來表示的，所以需要的控制頂點(diǎn)數(shù)為3。
其中節(jié)點(diǎn)矢量可由最大內(nèi)接面積表示法來確定。下面來確定剩余的所需數(shù)據(jù)。由有理Bezier曲線的公式得二次有理Bezier曲線弧的表示形式為：
稱k為形狀不變因子，公式如下所示：
可以證明同一組控制頂點(diǎn)選取不同 的權(quán)因子，只要形狀因子k相等，則由它們決定的二次有理Bezier曲線是同一條曲線段，不同的權(quán)因?qū)?yīng)不同的參數(shù)化，而且可以根據(jù)形狀不變因子對(duì)二次曲線進(jìn)行分類：
v K=0； 表示退化的二次曲線：一對(duì)直線段P0P1和P1P2；
v K∈[0，1]； 表示雙曲線；
v K=1； 表示拋物線；
v K∈[1, +∞]； 表示橢圓；
v K=+∞； 表示連接P0和P2的直線段；
習(xí)慣上我們選擇ω0=ω2=1稱為標(biāo)準(zhǔn)參數(shù)化。此時(shí)由形狀因子k公式得拋物線的ω1=1。所以拋物線的權(quán)因子也因此而確定，即ω0=ω1=ω2=1。
Figure 3.1 不同的權(quán)因子ω1 定義的圓錐截線
由二次有理Bezier曲線公式可知，當(dāng)u=0和u=1時(shí)，C(0)=P0, C(1)=P2，即曲線通過特征多邊形的首末頂點(diǎn)。由此可確定拋物線的兩個(gè)控制頂點(diǎn)P0和P2，現(xiàn)在只剩下最后一個(gè)P1頂點(diǎn)未確定。
P1點(diǎn)可以由二次有理Bezier曲線的公式列方程計(jì)算得出，但這并不是一個(gè)方便的方法。一種更方便的方法是：在這條曲線上指定第三個(gè)點(diǎn)，該點(diǎn)對(duì)應(yīng)于某個(gè)特定的參數(shù)，例如u=1/2。點(diǎn)S＝C(1/2)稱為圓錐截線的肩點(diǎn)(shoulder point)，如圖3.1所示。將u=1/2代入二次有理Bezier公式得：
其中M是弦P0P2的中點(diǎn)。肩點(diǎn)S＝C（1/2）。所以由上式可確定P1點(diǎn)。下面將各個(gè)控制頂點(diǎn)的計(jì)算列出如下所示：
設(shè)拋物線的起止區(qū)間為[UF，UL]，則因?yàn)榍€通過特征多邊形的頂點(diǎn)，所以通過首點(diǎn)P0＝（X0，Y0）＝（X0，UF）。又由拋物線的參數(shù)方程可知：
同理可計(jì)算出通過末點(diǎn)P2的坐標(biāo)，將他們分別列出如下：
由計(jì)算肩點(diǎn)S的公式，將Sy＝Y(jié)1＝C（u1）＝（UF+UL）/2代入可得：
所以P1點(diǎn)的坐標(biāo)為：
至此，拋物線的三個(gè)控制頂點(diǎn)P0，P1，P2都已計(jì)算出來了。即拋物線的NURBS表示所需的數(shù)據(jù)都已經(jīng)得到了。下面看看OpenCASCADE中的實(shí)現(xiàn)代碼。
4. Code Analysis
OpenCASCADE的Math工具集中有個(gè)包Covert用來將圓錐曲線曲面轉(zhuǎn)換為NURBS曲線曲面。其中轉(zhuǎn)換拋物線的類為：Convert_ParabolaToBSplineCurve，實(shí)現(xiàn)代碼如下所示：
//=======================================================================
//function : Convert_ParabolaToBSplineCurve
//purpose :
//=======================================================================
Convert_ParabolaToBSplineCurve::Convert_ParabolaToBSplineCurve
(const gp_Parab2d& Prb,
const Standard_Real U1 ,
const Standard_Real U2 )
: Convert_ConicToBSplineCurve (MaxNbPoles, MaxNbKnots, TheDegree)
{
Standard_DomainError_Raise_if( Abs(U2 - U1) < Epsilon(0.),
"Convert_ParabolaToBSplineCurve");
Standard_Real UF = Min (U1, U2);
Standard_Real UL = Max( U1, U2);
Standard_Real p = Prb.Parameter();
nbPoles = 3;
nbKnots = 2;
isperiodic = Standard_False;
knots->ChangeArray1()(1) = UF; mults->ChangeArray1()(1) = 3;
knots->ChangeArray1()(2) = UL; mults->ChangeArray1()(2) = 3;
weights->ChangeArray1()(1) = 1.;
weights->ChangeArray1()(2) = 1.;
weights->ChangeArray1()(3) = 1.;
gp_Dir2d Ox = Prb.Axis().XDirection();
gp_Dir2d Oy = Prb.Axis().YDirection();
Standard_Real S = ( Ox.X() * Oy.Y() - Ox.Y() * Oy.X() > 0.) ? 1 : -1;
// poles expressed in the reference mark
poles->ChangeArray1()(1) =
gp_Pnt2d( ( UF * UF) / ( 2. * p), S * UF );
poles->ChangeArray1()(2) =
gp_Pnt2d( ( UF * UL) / ( 2. * p), S * ( UF + UL) / 2. );
poles->ChangeArray1()(3) =
gp_Pnt2d( ( UL * UL) / ( 2. * p), S * UL );
// replace the bspline in the mark of the parabola
gp_Trsf2d Trsf;
Trsf.SetTransformation( Prb.Axis().XAxis(), gp::OX2d());
poles->ChangeArray1()(1).Transform( Trsf);
poles->ChangeArray1()(2).Transform( Trsf);
poles->ChangeArray1()(3).Transform( Trsf);
}
由上面的代碼可知，先設(shè)置曲線次數(shù)為2，再設(shè)置節(jié)點(diǎn)矢量為[UF，UF，UF，UL，UL，UL]，即首參數(shù)UF和末參數(shù)UL的重?cái)?shù)皆為3，由節(jié)點(diǎn)矢量可知轉(zhuǎn)換后的NURBS曲線為Bezier曲線。
三個(gè)控制頂點(diǎn)對(duì)應(yīng)的權(quán)因也都設(shè)置為1。三個(gè)控制頂點(diǎn)計(jì)算方法按上一節(jié)中所述。關(guān)鍵是P1點(diǎn)的計(jì)算。上面的計(jì)算方法僅為個(gè)人觀點(diǎn)，歡迎討論交流。
最后根據(jù)有理Bezier曲線的仿射不變性：對(duì)有理Bezier曲線進(jìn)行旋轉(zhuǎn)、平移和縮放變換，其表達(dá)式不變，只是控制點(diǎn)發(fā)生了改變。新的控制點(diǎn)可以通過對(duì)原控制點(diǎn)作變換得到。即要對(duì)有理Bezier曲線進(jìn)行仿射變換，只需對(duì)其控制點(diǎn)作變換即可。
圓錐截線的轉(zhuǎn)換類的使用是很簡單的，且計(jì)算都是在構(gòu)造函數(shù)中完成。下面給出一個(gè)將拋物線轉(zhuǎn)換為NURBS曲線的具體示例來說明其用法。
/*
* Copyright (c) 2014 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-10-02 20:46
* Version : 1.0v
*
* Description : OpenCASCADE conic to BSpline curve-Parabola.
*
* Key words : OpenCascade, Parabola, BSpline Curve, Convert
*/
#define HAVE_CONFIG_H
#include <gp_Parab2d.hxx>
#include <Convert_ParabolaToBSplineCurve.hxx>
void DumpConvertorInfo(const Convert_ConicToBSplineCurve &theConvertor)
{
Standard_Integer aCounter = 0;
std::cout << "Convert Result" << std::endl;
std::cout << "Degree: " << theConvertor.Degree() << std::endl;
std::cout << "Periodic: " << (theConvertor.IsPeriodic() ? "yes" : "no") << std::endl;
std::cout << "Knots: " << std::endl;
for (Standard_Integer i = 1;i <= theConvertor.NbKnots(); ++i)
{
for (Standard_Integer j = 1; j <= theConvertor.Multiplicity(i); ++j)
{
std::cout << ++aCounter << ": " << theConvertor.Knot(i) << std::endl;
}
}
std::cout << "Poles(Weight): " << std::endl;
for (Standard_Integer i = 1; i <= theConvertor.NbPoles(); ++i)
{
gp_Pnt2d aPole = theConvertor.Pole(i);
std::cout << i << ": " << aPole.X() << ", " << aPole.Y()
<< " W(" << theConvertor.Weight(i) << ")" << std::endl;
}
}
void TestParabolaConversion(void)
{
gp_Parab2d aParabola(gp::OX2d(), 1.0);
Convert_ParabolaToBSplineCurve aConvertor(aParabola, 1.0, M_PI);
DumpConvertorInfo(aConvertor);
}
int main(int argc, char **argv)
{
TestParabolaConversion();
return 0;
}
程序開始部分需要定義HAVE_CONFIG_H，這與在Windows中編程定義WNT有所不同。程序輸出結(jié)果如下圖所示：
Figure 4.1 Convert Parabola to BSpline Curve Result
5. Conclusion
NURBS的一個(gè)優(yōu)勢(shì)就是統(tǒng)一了曲線曲面的表示方法，即不僅可以表示自由曲線曲面，還可精確表示圓錐曲線曲面。其中拋物線是最簡單的圓錐截線，從簡單的拋物線轉(zhuǎn)換成NURBS曲線著手，學(xué)習(xí)NURBS是如何來表示圓錐截線的。
第一次使用Debian系統(tǒng)來編程，選用了一個(gè)輕量級(jí)的IDE開發(fā)工具Codelite，用法與Visual Studio類似。期間也遇到包含引用目錄及程序運(yùn)行時(shí)找不到動(dòng)態(tài)庫的問題。其中在Codelite中添加引用庫路徑的方法如下圖所示：
Figure 5.1 Set Library Path and Add Libraries in Codelite
添加上述引用庫后，解決了鏈接錯(cuò)誤，但是直接運(yùn)行程序，發(fā)現(xiàn)也有與在Windows中類似的問題，即缺少依賴的動(dòng)態(tài)庫，如下圖所示為使用lld命令檢查程序的依賴項(xiàng)。最后通過向ld.so.conf中添加動(dòng)態(tài)庫所在的目錄解決了問題。相關(guān)命令如下所示：
# cat /etc/ld.so.conf
include ld.so.conf.d/*.conf
# echo "/home/eryar/opencascade-6.7.1/lib" >> /etc/ld.so.conf
# ldconfig
注：以上命令需要有root權(quán)限。
Figure 5.2 Use lld command to check depends
通過解決以上的問題，來適應(yīng)在Debian中使用Codelite開發(fā)程序。
6. References
1. 人民教育出版社中學(xué)數(shù)學(xué)室. 數(shù)學(xué)第二冊(cè)（上）. 人民教育出版社. 2000
2. 趙罡，穆國旺，王拉柱譯. 非均勻有理B樣條. 清華大學(xué)出版社. 2010
3. 王仁宏，李崇君，朱春鋼. 計(jì)算幾何教程. 科學(xué)出版社. 2008