B-Spline Curve Library in Open Cascade

Size: px

Start display at page:

Download "B-Spline Curve Library in Open Cascade"

最充
5 years ago
Views:

1 B-Spline Curve Library in Open Cascade Open Cascade 中的 B 样条曲线库 eryar@163.com 摘要 Abstract: 简要介绍 Open Cascade 中的 B 样条曲线库 BSplCLib 的使用方法, 并且结合源程序来对 Open Cascade 中的 B 样条曲线的组成部分如节点矢量重复度等概念进行介绍, 以及通过对计算 B 样条基函数的算法进行分析, 加深对 B 样条曲线概念的理解关键字 Key Word:B Spline Curve Open Cascade Knot Vector Multiplicity 一概述 Overview 1946 年由 Schoenberg 提出了 B 样条理论, 给出了 B 样条的差分表达式 ;1972 年 de Boor 和 Cox 分别独立给出了关于 B 样条的标准算法 Gordon 和 Riesenfeld 又把 B 样条理论用于形状描述, 最终提出了 B 样条方法用 B 样条基替代了 Bernstein 基, 构造出 B 样条曲线, 这种方法继承了 Bezier 方法的一切优点, 克服了 Bezier 方法存在的缺点, 较成功地解决了局部控制问题, 又轻而易举地在参数连续性基础上解决了连接问题, 从而使自由曲线曲面形状的描述问题得到较好解决 p 次 B 样条曲线的定义为 : C( u ) = n i= 0 i,p ( u ) P 其中 : Pi 是控制顶点 (control point); i,p(u) 是定义在非周期节点矢量上的 p 次 B 样条基函数 ; i 有很多方法可以用来定义 B 样条基函数以及证明它的一些重要性质例如, 可以采用截尾幂函数的差商定义, 开花定义, 以及由 de Boor 和 Cox 等人提出的递推公式等来定义我们这里采用的是递推定义方法, 因为这种方法在计算机实现中是最有效的令 U=u0,u1,,um 是一个单调不减的实数序列, 即 ui<=u1,i=0,1,,m-1 其中,ui 称为节点,U 称为节点矢量, 用 i,p(u) 表示第 i 个 p 次 B 样条基函数, 其定义为 : i,0 i, p ( u ) = ( u ) = 1 ui u < u 0 otherwise u ui u u p B 样条基有如下性质 : i i,p 1 1 ( u ) + u u p+ 1 p+ 1 u u 1 1,p 1 ( u )

2 a) 递推性 ; b) 局部支承性 ; c) 规范性 ; d) 可微性 ; 根据 B 样条曲线定义可知, 给定控制顶点 Pi(control points), 曲线次数 p(degree) 及节点矢量 U(knot vectors),b 样曲线也就确定对于有理 B 样条曲线, 还需要参数权重 (weights)

3 二 OCC 中的 B 样条曲线库 BSplCLib in OCC 在 Open Cascade 中的工具箱 (Toolkit)TKMath 中的包 (package)bsplclib 是 B 样条曲线库, 为 B 样条曲线曲面的计算提供了支持它提供了三方面的功能 : 对节点矢量 (knot vectors) 及重复度 (multiplicities) 的管理 ; 对多维样条的支持, 即 B 样条方法中控制顶点的维数可以是任意维数 (dimension); 二维和三维样条曲线的方法 ; Open Cascade 中的 B 样条曲线由下列数据项定义 : 定义变量类型变量名称控制顶点 control points TColgp_Array1OfPnt Poles 权重 weights TColStd_Array1OfReal Weights 节点 knots TColStd_Array1OfReal Knots 重数 multiplicities TColStd_Array1OfInteger Mults 次数 degree Standard_Integer Degree 周期性 periodicity Standard_Boolean Periodic B 样条曲线库 BSplCLib 提供了一些基本几何算法 : B 样条基函数及其导数的计算 BSplCLib::EvalBsplineBasis(); 节点插入 BSplCLib::InsertKnot(); 节点去除 BSplCLib::RemoveKnot(); 升阶 BSplCLib::IncreaseDegree(); 降阶 ; 结合 The URBS Book 和 Open Cascade 中的 BSplCLib 的源程序, 可以高效的学习 URBS The URBS Book 中有详细的理论推导及算法描述, 而 Open Cascade 中有可以用来实际使用的程序理论联系实际, 有助于快速理解 URBS 的有关概念及其应用

4 三 OCC 中 B 样条曲线库的节点和重数 Knots and Multiplicity in BSplCLib 由 B 样条曲线的可微性可知, 节点的重数与 B 样条曲线的连续性相关在节点区间内部,i,p(u) 是无限次可微的, 因在每个节点区间内部, 它是一个多项式在节点处 i,p(u) 是 p-k 次连续的, 其中 k 是节点的重复度 (multiplicity, 有时也称为重数 ) 因此, 增加次数 p 将提高曲线的连续性, 而增加节点的重复度则使连续性降低重复度 (multiplicity, 有时也称为重数 ) 有两种不同的理解方式 : 节点在节点矢量中的重复度 ; 节点相对于一个特定的基函数的重复度 ; 在 Open Cascade 中对重复度的理解是前者, 即节点在节点矢量中的重复度下面结合源程序来进行说明函数 BSplCLib::Knots() 用来将给定的节点矢量 ( 节点序列 knot sequence) 转换为节点的重复度不大于 1 的 Knots 数组和每个节点对应的重复度 Mults 数组, 且数据 Knots 和 Mults 的长度必由函数 BSplCLib::KnotsLength() 得到 Knots() 函数的源程序如下所示 : ======================================================================= function : Knots purpose : Computes the sequence of knots Knots without repetition of the knots of multiplicity greater than 1. Length of <Knots> and <Mults> must be KnotsLength(KnotSequence,Periodic). ======================================================================= void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots, TColStd_Array1OfReal &knots, TColStd_Array1OfInteger &mult, const Standard_Boolean Periodic) const Standard_Boolean ) Standard_Real val = SeqKnots(1); Standard_Integer kk=1; knots(kk) = val; mult(kk) = 1; for (Standard_Integer jj=2;jj<=seqknots.length();jj++) test on strict equality on nodes if (SeqKnots(jj)!=val) val = SeqKnots(jj); kk++; knots(kk) = val; mult(kk) = 1; else

5 mult(kk)++; 从上述代码可知, 直接使用了不等于来判断两个节点的值是否相同, 而没有采用误差处理, 即严格的相等比较程序将节点重复度不大于 1 的节点及其相应的重复度分别保存到 knots 和 mult 中

6 四 B 样条曲线的分类 B Spline Curve Type B 样条曲线一般按定义基函数的节点序列是否等距 ( 均匀 ) 分为均匀 B 样条曲线 (Uniform B-Spline Curve) 和非均匀 B 样条曲线 (on Uniform B-Spline Curve) B 样条曲线按节点序列中节点分布情况不同, 又分为四种类型 : 均匀 B 样条曲线准均匀 B 样条曲线分段 Bezier 曲线一般非均匀 B 样条曲线设给定特征多边形顶点 Vi, i=0,1,,n, 曲线次数 k, 则有 : 均匀 B 样条曲线 (uniform B-Spline curve): 节点序列中节点沿参数轴均匀或等距分布, 即所有节点区间长度为大于零的常数 (constant): i = t 1 ti = constant > 0, i = 0,1,, n + k 准均匀 B 样曲线 (quasi-uniform B-Spline curve): 其节点序列中两端节点具有重复度 k+1, 而所有内节点均匀分布, 具有重复度 1 分段 Bezier 曲线 (piecewise Bezier curve): 其节点序列中两端节点重复度与准均匀 B 样条曲线的相同, 所不同的是所有内节点重复度为 k 非均匀 B 样条曲线 (general non-uniform B-Spline curve): 这是对任意分布的节点序列, 只要在数学上成立, 即节点序列非递减, 都可取在基础类模块 (Module FoundationClasses) 的工具箱 (Toolkit TKMath) 中的包 (GeomAbs) 中有对 B 样样条曲线类型的定义, 源程序如下所示 :! This enumeration is used to note specific curve form. <br> enum GeomAbs_BSplKnotDistribution GeomAbs_onUniform, GeomAbs_Uniform, GeomAbs_QuasiUniform, GeomAbs_PiecewiseBezier ; 而类 BSplCLib 主要是用来管理节点和重复度的, 所有将节点和重复度也进行了分类根据节点矢量是否均匀分布, 将节点分配方式 ( Knot Distribution ) 分为 : 均匀 (BSplCLib_Uniform) 和非均匀 (BSplCLib_onUniform) 源程序如下所示: enum BSplCLib_KnotDistribution BSplCLib_onUniform, BSplCLib_Uniform ; 根据重复度数组将重复度的分配方式分为如下三种类型 : BSplCLib_Constant: 重复度都相同 ; BSplCLib_QuasiConstant: 首尾节点的重复度与内部节点的重复度不同 ; BSplClib_onConstant: 其它情况 ; 源程序如下所示 : enum BSplCLib_MultDistribution BSplCLib_onConstant, BSplCLib_Constant, BSplCLib_QuasiConstant ; 判断节点矢量和重复度矢量类型分别由下列函数实现 :

7 BSplCLib::KnotForm(); BSplCLib::MultForm(); 具体的判断方法可以查看源程序将节点分布方式与重复度的分布方式进行组合, 可以得出 B 样条曲线的那几种类型

8 五 B 样条基函数的计算 Evaluate the B-Spline Basis B 样条基函数的计算主要使用了 B 样条基了函数的递推公式 (Cox-deBoor 公式 ) 的局部支撑性质, 如下所示 : i,0 i, p ( u ) = ( u ) = 1 ui u < u 0 otherwise u ui u u p i i,p 1 1 ( u ) + u u p+ 1 p+ 1 u u 1 1,p 1 ( u ) 直接由定义可知 : i,0(u) 是一个阶梯函数, 它在半开区间 u [ui, u1) 外都为零 ; 当次数 p>0 时,i,p(u) 是两个 p-1 次基函数的线性组合 ; 计算一组基函数需要事先指定节点矢量 U 和次数 p; 半开区间 [ui,u1) 称为第 i 个节点区间 (knot span), 它的长度可以为零, 因为相邻节点可以是相同的 ; 计算 p 次基函数的过程可以生成一个如下形式的三角形陈列 : 0,0 1,0 2,0 3,0 4,0 0,1 1,1 2,1 3,1 0,2 1,2 2,2 0,3 1,3 B 样条有局部支撑性, 即若 u 不在区间 [ui, up+1), 则 i,p(u)=0 可从下面的三角形中看出 1,3 是 1,0 2,0 3,0 和 4,0 的线性组合, 而 1,0 在区间 [u1, u2) 上非零,2,0 在区间 [u2,u3) 上非零,3,0 在区间 [u3,u4) 上非零,4,0 在区间 [u4,u5) 上非零, 所以 1,3 仅在区间 [u1,u5) 上非零在任意给定的节点区间 [uj,uj+1) 内, 最多有 p+1 个是非零的, 它们是 j-p,p j-p+1 j,p 例如, 在 [u3,u4) 上, 零次基函数中只有 3,0 是非零的, 一次基函数只有 2,1 和 3,1 是非零的, 非零的三次基函数只有 0,3 1,3 2,3 3,3 这个性质如下图所示:

9 上面两幅图中右边的图中所示的推算过程表明, 给定节点序列 U 及 B 样条曲线的次数 p, 给出任意一个 u 值, 找出其所在的节点区间 [ui,u1) 上, 最多有 i-p,p,i-p+1,p,,i,p 个非零的基函数例如我们根据递推公式写出二次基函数的一般形式, 如下所示 : 当给定的 u 值在区间 [u3,u4) 上即 (i=3) 时, 根据上面的三角形, 得出下列重要结论 : i 2,1 3 2,1 3+ 1,1 ( u) ( u) ( u) = = = 1,1 1,1 4,1 ( u) ( u) ( u) = 0 = 0 = 0 即这两项不需要计算另外一个重要结论就是图中用相同颜色框中的部分是相同的, 也就是下面程序中的变量 temp 表示的内容我们引入下面符号 :

10 由二次基函数推出的三个公式可写为 : 上述推导过程为 The URBS Book 中的算法, 算法代码如下所示 : 理解了变量 temp 的意义之后, 整个程序就很好理解了将 Open Cascade 中计算基函数的算法是不同的, 将其源程序摘抄如下所示 : ======================================================================= function : Build BSpline Matrix

11 purpose : Builds the Bspline Matrix ======================================================================= Standard_Integer BSplCLib::EvalBsplineBasis (const Standard_Integer Side, = 1 rigth side, -1 left side (const Standard_Integer, = 1 rigth side, -1 left side const Standard_Integer DerivativeRequest, const Standard_Integer Order, const TColStd_Array1OfReal& FlatKnots, const Standard_Real Parameter, Standard_Integer& FirstonZeroBsplineIndex, math_matrix& BsplineBasis) the matrix must have at least DerivativeRequest + 1 row and Order columns the result are stored in the following way in the Bspline matrix Let i be the FirstonZeroBsplineIndex and t be the parameter value, k the order of the knot vector, r the DerivativeRequest : B (t) B (t) B (t) i 1 k-1 (1) (1) (1) B (t) B (t) B (t) i 1 k-1 (r) (r) (r) B (t) B (t) B (t) i 1 k-1 Standard_Integer ReturnCode, ii, pp, qq, ss, umpoles, LocalRequest ;,Index ;

12 Standard_Real ewparameter, Inverse, Factor, LocalInverse, Saved ;, *FlatKnotsArray ; ReturnCode = 0 ; FirstonZeroBsplineIndex = 0 ; LocalRequest = DerivativeRequest ; if (DerivativeRequest >= Order) LocalRequest = Order - 1 ; if (BsplineBasis.LowerCol()!= 1 BsplineBasis.UpperCol() < Order BsplineBasis.LowerRow()!= 1 BsplineBasis.UpperRow() <= LocalRequest) ReturnCode = 1; goto FIISH ; umpoles = FlatKnots.Upper() - FlatKnots.Lower() Order ; BSplCLib::LocateParameter(Order - 1, FlatKnots, Parameter, Standard_False, Order, umpoles+1, ii, ewparameter) ; FirstonZeroBsplineIndex = ii - Order + 1 ; BsplineBasis(1,1) = 1.0e0 ; LocalRequest = DerivativeRequest ; if (DerivativeRequest >= Order) LocalRequest = Order - 1 ; for (qq = 2 ; qq <= Order - LocalRequest ; qq++) BsplineBasis(1,qq) = 0.0e0 ; for (pp = 1 ; pp <= qq - 1 ; pp++)

13 this should be always invertible if ii is correctly computed Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ; Saved = Factor * BsplineBasis(1,pp) ; BsplineBasis(1,pp) *= (1.0e0 - Factor) ; BsplineBasis(1,pp) += BsplineBasis(1,qq) ; BsplineBasis(1,qq) = Saved ; for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) for (pp = 1 ; pp <= qq - 1 ; pp++) BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ; BsplineBasis(1,qq) = 0.0e0 ; for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) BsplineBasis(Order - ss + 2,qq) = 0.0e0 ; for (pp = 1 ; pp <= qq - 1 ; pp++) Inverse = 1.0e0 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ; Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ; Saved = Factor * BsplineBasis(1,pp) ; BsplineBasis(1,pp) *= (1.0e0 - Factor) ; BsplineBasis(1,pp) += BsplineBasis(1,qq) ; BsplineBasis(1,qq) = Saved ; LocalInverse = (Standard_Real) (qq - 1) * Inverse ; for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ; BsplineBasis(Order - ss + 2, pp) *= - LocalInverse ; BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2,qq) ; BsplineBasis(Order - ss + 2,qq) = Saved ; FIISH : return (ReturnCode) ; 函数的作用是用来计算所有的基函数及其导数, 并将结果以矩阵 ( 数组 ) 的形式保存结合

14 二次基函数的推导方法, 将述代码写成公式的形式函数的参数及其描述如下表所示 : 变量描述 DerivativeRequest 导数的次数 Order B 样条基函数的阶数 ( 次数 +1) FlatKnots 节点矢量 Parameter 参数 FirstonZeroBspline 第一个非零基函数的索引值 BsplineBasis 基函数值矩阵当导数次数 DerivativeRequest 大于 B 样条基的阶数 Order 时, 将计算导数的次数设置为 B 样条基的次数 (Order-1) 程序代码如下所示: LocalRequest = DerivativeRequest ; if (DerivativeRequest >= Order) LocalRequest = Order - 1 ; 对 B 样条基数计算结果矩阵 BsplineBasis 存储空间进行检查若存储空间不足, 则会退出, 程序代码如下所示 : if (BsplineBasis.LowerCol()!= 1 BsplineBasis.UpperCol() < Order BsplineBasis.LowerRow()!= 1 BsplineBasis.UpperRow() <= LocalRequest) ReturnCode = 1; goto FIISH ; 确定参数 Parameter 所在的节点区间的下标 ( 索引值 ), 程序代码如下所示 : umpoles = FlatKnots.Upper() - FlatKnots.Lower() Order ; BSplCLib::LocateParameter(Order - 1, FlatKnots, Parameter, Standard_False, Order, umpoles+1, ii, ewparameter) ; 确定参数 Parameter 所在区间的算法是用二分法搜索得到程序代码如下所示 : ======================================================================= function : Hunt purpose : ======================================================================= void BSplCLib::Hunt (const Array1OfReal& XX, const Standard_Real X, Standard_Integer& Ilc) replaced by simple dichotomy (RLE) Ilc = XX.Lower();

15 const Standard_Real *px = &XX(Ilc); px -= Ilc; if (X < px[ilc]) Ilc--; return; Standard_Integer Ihi = XX.Upper(); if (X > px[ihi]) Ilc = Ihi + 1; return; Standard_Integer Im; while (Ihi - Ilc!= 1) Im = (Ihi + Ilc) >> 1; if (X > px[im]) Ilc = Im; else Ihi = Im; 确定参数所在区间 [ui,u1) 后, 可得到第一个非零基函数的索引值为 i-p; FirstonZeroBsplineIndex = ii - Order + 1 ; 基函数计算的主要算法代码如下所示 : BsplineBasis(1,1) = 1.0e0 ; for (qq = 2 ; qq <= Order - LocalRequest ; qq++) BsplineBasis(1,qq) = 0.0e0 ; for (pp = 1 ; pp <= qq - 1 ; pp++) this should be always invertible if ii is correctly computed Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ; Saved = Factor * BsplineBasis(1,pp) ; BsplineBasis(1,pp) *= (1.0e0 - Factor) ; BsplineBasis(1,pp) += BsplineBasis(1,qq) ; BsplineBasis(1,qq) = Saved ; 其中 : Factor = u ui u u p 为节点区间 [ui,u1) 上的基函数左边部分的系数 ; i

16 (1 Factor) = u u p p u u 为节点区间 [ui-1,ui) 上的基函数右边部分的系数 ; i 六程序示例 Sample Codes 将上述内容以一个简单示例程序来验证, 程序代码如下所示 : /* * Copyright (c) 2013 eryar All Rights Reserved. * * File : Main.cpp * Author : eryar@163.com * Date : * Version : * * Description : Learn the B-Spline Curve library in the Open Cascade. * */ #include <BSplCLib.hxx> #include <math_matrix.hxx> #include <TColStd_Array1OfReal.hxx> #include <TColStd_Array1OfInteger.hxx> #include <Geom2d_BSplineCurve.hxx> #pragma comment(lib, "TKernel.lib") #pragma comment(lib, "TKMath.lib") #pragma comment(lib, "TKG2d.lib") int main(int argc, char* argv[]) Knot vector: [0,0,0,1,2,3,4,4,5,5,5] TColStd_Array1OfReal knotseq(1, 11); knotseq.init(0); knotseq.setvalue(1, 0); knotseq.setvalue(2, 0); knotseq.setvalue(3, 0); knotseq.setvalue(4, 1); knotseq.setvalue(5, 2); knotseq.setvalue(6, 3); knotseq.setvalue(7, 4); knotseq.setvalue(8, 4);

17 knotseq.setvalue(9, 5); knotseq.setvalue(10, 5); knotseq.setvalue(11, 5); cout<<"knot Sequence: [ "; for (Standard_Integer i = 1; i <= knotseq.length(); +) cout<<knotseq.value(i)<<" "; cout<<"]"<<endl; Standard_Integer knotslen = BSplCLib::KnotsLength(knotSeq); TColStd_Array1OfReal knots(1, knotslen); TColStd_Array1OfInteger mults(1, knotslen); Test Knots, Mults and Knot sequence of BSplCLib. BSplCLib::Knots(knotSeq, knots, mults); cout<<"knots: [ "; for (Standard_Integer i = 1; i <= knots.length(); +) cout<<knots.value(i)<<" "; cout<<"]"<<endl; cout<<"multiplicity: [ "; for (Standard_Integer i = 1; i <= mults.length(); +) cout<<mults.value(i)<<" "; cout<<"]"<<endl; if (BSplCLib::KnotForm(knots, 1, knotslen) == BSplCLib_Uniform) cout<<"knots is uniform."<<endl; else cout<<"knots is non-uniform."<<endl; Standard_Real rvalue = 2.5; Standard_Integer iorder = 2+1; Standard_Integer ifirstonzeroindex = 0;

18 math_matrix bsplinebasis(1, 1, 1, iorder, 0); BSplCLib::EvalBsplineBasis(1, 0, iorder, knotseq, rvalue, ifirstonzeroindex, bsplinebasis); cout<<"first on-zero Basis index: "<<ifirstonzeroindex<<endl; cout<<bsplinebasis<<endl; return 0; 上述代码对节点矢量重复度的概念的验证, 并以一个实例计算所有非零基函数的值程序输出为 : Knot Sequence: [ ] Knots: [ ] Multiplicity: [ ] Knots is uniform. First on-zero Basis index: 3 math_matrix of Rowumber = 1 and Columber = 3 math_matrix ( 1, 1 ) = math_matrix ( 1, 2 ) = 0.75 math_matrix ( 1, 3 ) = Press any key to continue...

19 七结论 Conclusion 通过学习 The URBS Book 并给合 Open Cascade 的源程序, 理论联系实际, 使对 URBS 的学习更轻松根据 B 样条基的递推公式,B 样条曲线的局部性是通过节点来具体实现的与 Bezier 曲线不同的就是增加了节点这个参数根据 Cox-deBoor 递推公式亲自推导出一次二次三次 B 样条基函数, 可以加深对 B 样条曲线的理解计算给定节点矢量次数及参数, 计算参数所在区间上所有非零基函数算法的步骤为 : 通过二分法查找出参数所在的节点区间 ; 根据 B 样条基的局部支撑性, 计算出所在节点区间上所有非零基函数 ; 八致谢 Acknowledgments 感谢晓天的支持与鼓励九参考文献 Bibliography 1. 赵罡, 穆国旺, 王拉柱译 Les Piegl,Wayne Tiller The URBS Book(Second Edition) 2010 清华大学出版社 2. 莫容, 常智勇计算机辅助几何造型技术 2009 科学出版社 3. 王仁宏, 李崇君, 朱春钢计算几何教程 2008 科学出版社

绘制OpenCascade中的曲线

绘制OpenCascade中的曲线在 OpenSceneGraph 中绘制 OpenCascade 的曲线 Draw OpenCascade Geometry Curves in OpenSceneGraph eryar@163.com 摘要 Abstract: 本文简要说明 OpenCascade 中几何曲线的数据, 并将这些几何曲线在 OpenSceneGraph 中绘制出来关键字 KeyWords:OpenCascade Geometry