// Copyright (c) 1995-1999 Matra Datavision // Copyright (c) 1999-2014 OPEN CASCADE SAS // // This file is part of Open CASCADE Technology software library. // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License version 2.1 as published // by the Free Software Foundation, with special exception defined in the file // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT // distribution for complete text of the license and disclaimer of any warranty. // // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. // xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem // Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix, // EvalPolynomial, RationalDerivatives. // xab : 22-Mar-94 : fixed rational problem in BuildCache // xab : 12-Mar-96 : added MovePointAndTangent #include #include #include #include #include #include #include //======================================================================= //struct : BSplCLib_DataContainer //purpose: Auxiliary structure providing buffers for poles and knots used in // evaluation of bspline (allocated in the stack) //======================================================================= struct BSplCLib_DataContainer { BSplCLib_DataContainer(Standard_Integer Degree) { (void)Degree; // avoid compiler warning Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25, "BSplCLib: bspline degree is greater than maximum supported"); } Standard_Real poles[(25+1)*(Dimension_gen+1)]; Standard_Real knots[2*25]; Standard_Real ders[Dimension_gen*4]; }; //======================================================================= //function : Reverse //purpose : //======================================================================= void BSplCLib::Reverse(Array1OfPoints& Poles, const Standard_Integer L) { Standard_Integer i, l = L; l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1); Array1OfPoints temp(0,Poles.Length()-1); for (i = Poles.Lower(); i <= l; i++) temp(l-i) = Poles(i); for (i = l+1; i <= Poles.Upper(); i++) temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i); for (i = Poles.Lower(); i <= Poles.Upper(); i++) Poles(i) = temp(i-Poles.Lower()); } // // CURVES MODIFICATIONS // //======================================================================= //function : RemoveKnot //purpose : //======================================================================= Standard_Boolean BSplCLib::RemoveKnot (const Standard_Integer Index, const Standard_Integer Mult, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults, const Standard_Real Tolerance) { Standard_Boolean rational = Weights != NULL; Standard_Integer dim; dim = Dimension_gen; if (rational) dim++; TColStd_Array1OfReal poles(1,dim*Poles.Length()); TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length()); if (rational) PLib::SetPoles(Poles,*Weights,poles); else PLib::SetPoles(Poles,poles); if (!RemoveKnot(Index,Mult,Degree,Periodic,dim, poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance)) return Standard_False; if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights); else PLib::GetPoles(newpoles,NewPoles); return Standard_True; } //======================================================================= //function : InsertKnots //purpose : insert an array of knots and multiplicities //======================================================================= void BSplCLib::InsertKnots (const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, const TColStd_Array1OfReal& AddKnots, const TColStd_Array1OfInteger* AddMults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults, const Standard_Real Epsilon, const Standard_Boolean Add) { Standard_Boolean rational = Weights != NULL; Standard_Integer dim; dim = Dimension_gen; if (rational) dim++; TColStd_Array1OfReal poles(1,dim*Poles.Length()); TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length()); if (rational) PLib::SetPoles(Poles,*Weights,poles); else PLib::SetPoles(Poles,poles); InsertKnots(Degree,Periodic,dim,poles,Knots,Mults, AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add); if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights); else PLib::GetPoles(newpoles,NewPoles); } //======================================================================= //function : InsertKnot //purpose : //======================================================================= void BSplCLib::InsertKnot(const Standard_Integer , const Standard_Real U, const Standard_Integer UMult, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights) { TColStd_Array1OfReal k(1,1); k(1) = U; TColStd_Array1OfInteger m(1,1); m(1) = UMult; TColStd_Array1OfReal nk(1,Knots.Length()+1); TColStd_Array1OfInteger nm(1,Knots.Length()+1); InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults, k,&m,NewPoles,NewWeights,nk,nm,Epsilon(U)); } //======================================================================= //function : RaiseMultiplicity //purpose : //======================================================================= void BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex, const Standard_Integer Mult, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights) { TColStd_Array1OfReal k(1,1); k(1) = Knots(KnotIndex); TColStd_Array1OfInteger m(1,1); m(1) = Mult - Mults(KnotIndex); TColStd_Array1OfReal nk(1,Knots.Length()); TColStd_Array1OfInteger nm(1,Knots.Length()); InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults, k,&m,NewPoles,NewWeights,nk,nm,Epsilon(k(1))); } //======================================================================= //function : IncreaseDegree //purpose : //======================================================================= void BSplCLib::IncreaseDegree (const Standard_Integer Degree, const Standard_Integer NewDegree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults) { Standard_Boolean rational = Weights != NULL; Standard_Integer dim; dim = Dimension_gen; if (rational) dim++; TColStd_Array1OfReal poles(1,dim*Poles.Length()); TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length()); if (rational) PLib::SetPoles(Poles,*Weights,poles); else PLib::SetPoles(Poles,poles); IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults, newpoles,NewKnots,NewMults); if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights); else PLib::GetPoles(newpoles,NewPoles); } //======================================================================= //function : Unperiodize //purpose : //======================================================================= void BSplCLib::Unperiodize (const Standard_Integer Degree, const TColStd_Array1OfInteger& Mults, const TColStd_Array1OfReal& Knots, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, TColStd_Array1OfInteger& NewMults, TColStd_Array1OfReal& NewKnots, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights) { Standard_Boolean rational = Weights != NULL; Standard_Integer dim; dim = Dimension_gen; if (rational) dim++; TColStd_Array1OfReal poles(1,dim*Poles.Length()); TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length()); if (rational) PLib::SetPoles(Poles,*Weights,poles); else PLib::SetPoles(Poles,poles); Unperiodize(Degree, dim, Mults, Knots, poles, NewMults,NewKnots, newpoles); if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights); else PLib::GetPoles(newpoles,NewPoles); } //======================================================================= //function : Trimming //purpose : //======================================================================= void BSplCLib::Trimming(const Standard_Integer Degree, const Standard_Boolean Periodic, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const Standard_Real U1, const Standard_Real U2, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults, Array1OfPoints& NewPoles, TColStd_Array1OfReal* NewWeights) { Standard_Boolean rational = Weights != NULL; Standard_Integer dim; dim = Dimension_gen; if (rational) dim++; TColStd_Array1OfReal poles(1,dim*Poles.Length()); TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length()); if (rational) PLib::SetPoles(Poles,*Weights,poles); else PLib::SetPoles(Poles,poles); Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2, NewKnots, NewMults, newpoles); if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights); else PLib::GetPoles(newpoles,NewPoles); } //-------------------------------------------------------------------------- //ELEMENTARY COMPUTATIONS //-------------------------------------------------------------------------- // // All the following methods are methods of low level used in BSplCLib // but also in BSplSLib (but not in Geom) // At the creation time of this package there is no possibility // to declare private methods of package and to declare friend // methods of package. It could interesting to declare that BSplSLib // is a package friend to BSplCLib because it uses all the basis methods // of BSplCLib. //-------------------------------------------------------------------------- //======================================================================= //function : BuildEval //purpose : builds the local array for evaluation //======================================================================= void BSplCLib::BuildEval(const Standard_Integer Degree, const Standard_Integer Index, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, Standard_Real& LP) { Standard_Real w, *pole = &LP; Standard_Integer PLower = Poles.Lower(); Standard_Integer PUpper = Poles.Upper(); Standard_Integer i; Standard_Integer ip = PLower + Index - 1; if (Weights == NULL) { for (i = 0; i <= Degree; i++) { ip++; if (ip > PUpper) ip = PLower; const Point& P = Poles(ip); PointToCoords (pole,P,+0); pole += Dimension_gen; } } else { for (i = 0; i <= Degree; i++) { ip++; if (ip > PUpper) ip = PLower; const Point& P = Poles(ip); pole[Dimension_gen] = w = (*Weights)(ip); PointToCoords (pole, P, * w); pole += Dimension_gen + 1; } } } //======================================================================= //function : PrepareEval //purpose : stores data for Eval in the local arrays // dc.poles and dc.knots //======================================================================= static void PrepareEval (Standard_Real& u, Standard_Integer& index, Standard_Integer& dim, Standard_Boolean& rational, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, BSplCLib_DataContainer& dc) { // Set the Index BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u); // make the knots BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots); if (Mults == NULL) index -= Knots.Lower() + Degree; else index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults); // check truly rational rational = (Weights != NULL); if (rational) { Standard_Integer WLower = Weights->Lower() + index; rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree); } // make the poles if (rational) { dim = Dimension_gen + 1; BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles); } else { dim = Dimension_gen; BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles); } } //======================================================================= //function : D0 //purpose : //======================================================================= void BSplCLib::D0 (const Standard_Real U, const Standard_Integer Index, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, Point& P) { // Standard_Integer k,dim,index = Index; Standard_Integer dim,index = Index; Standard_Real u = U; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc); BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles); if (rational) { Standard_Real w = dc.poles[Dimension_gen]; CoordsToPoint (P, dc.poles, / w); } else CoordsToPoint (P, dc.poles, ); } //======================================================================= //function : D1 //purpose : //======================================================================= void BSplCLib::D1 (const Standard_Real U, const Standard_Integer Index, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, Point& P, Vector& V) { Standard_Integer dim,index = Index; Standard_Real u = U; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc); BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles); Standard_Real *result = dc.poles; if (rational) { PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders); result = dc.ders; } CoordsToPoint (P, result, ); CoordsToPoint (V, result + Dimension_gen, ); } //======================================================================= //function : D2 //purpose : //======================================================================= void BSplCLib::D2 (const Standard_Real U, const Standard_Integer Index, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, Point& P, Vector& V1, Vector& V2) { Standard_Integer dim,index = Index; Standard_Real u = U; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc); BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles); Standard_Real *result = dc.poles; if (rational) { PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders); result = dc.ders; } CoordsToPoint (P, result, ); CoordsToPoint (V1, result + Dimension_gen, ); if (!rational && (Degree < 2)) NullifyPoint (V2); else CoordsToPoint (V2, result + 2 * Dimension_gen, ); } //======================================================================= //function : D3 //purpose : //======================================================================= void BSplCLib::D3 (const Standard_Real U, const Standard_Integer Index, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, Point& P, Vector& V1, Vector& V2, Vector& V3) { Standard_Integer dim,index = Index; Standard_Real u = U; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc); BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles); Standard_Real *result = dc.poles; if (rational) { PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders); result = dc.ders; } CoordsToPoint (P, result, ); CoordsToPoint (V1, result + Dimension_gen, ); if (!rational && (Degree < 2)) NullifyPoint (V2); else CoordsToPoint (V2, result + 2 * Dimension_gen, ); if (!rational && (Degree < 3)) NullifyPoint (V3); else CoordsToPoint (V3, result + 3 * Dimension_gen, ); } //======================================================================= //function : DN //purpose : //======================================================================= void BSplCLib::DN (const Standard_Real U, const Standard_Integer N, const Standard_Integer Index, const Standard_Integer Degree, const Standard_Boolean Periodic, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger* Mults, Vector& VN) { Standard_Integer dim,index = Index; Standard_Real u = U; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc); BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles); if (rational) { Standard_Real v[Dimension_gen]; PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False); CoordsToPoint (VN, v, ); } else { if (N > Degree) NullifyPoint (VN); else { Standard_Real *DN = dc.poles + N * Dimension_gen; CoordsToPoint (VN, DN, ); } } } //======================================================================= //function : Solves a LU factored Matrix //purpose : //======================================================================= Standard_Integer BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, const Standard_Integer UpperBandWidth, const Standard_Integer LowerBandWidth, Array1OfPoints& PolesArray) { Standard_Real *PArray ; PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ; return BSplCLib::SolveBandedSystem(Matrix, UpperBandWidth, LowerBandWidth, Dimension_gen, PArray[0]) ; } //======================================================================= //function : Solves a LU factored Matrix //purpose : if HomogeneousFlag is 1 then the input and the output // will be homogeneous that is no division or multiplication // by weigths will happen. On the contrary if HomogeneousFlag // is 0 then the poles will be multiplied first by the weights // and after interpolation they will be divided by the weights //======================================================================= Standard_Integer BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, const Standard_Integer UpperBandWidth, const Standard_Integer LowerBandWidth, const Standard_Boolean HomogeneousFlag, Array1OfPoints& PolesArray, TColStd_Array1OfReal& WeightsArray) { Standard_Real *PArray, * WArray ; PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ; WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ; return BSplCLib::SolveBandedSystem(Matrix, UpperBandWidth, LowerBandWidth, HomogeneousFlag, Dimension_gen, PArray[0], WArray[0]) ; } //======================================================================= //function : Evaluates a Bspline function : uses the ExtrapMode //purpose : the function is extrapolated using the Taylor expansion // of degree ExtrapMode[0] to the left and the Taylor // expansion of degree ExtrapMode[1] to the right // if the HomogeneousFlag == True than the Poles are supposed // to be stored homogeneously and the result will also be homogeneous // Valid only if Weights //======================================================================= void BSplCLib::Eval(const Standard_Real Parameter, const Standard_Boolean PeriodicFlag, const Standard_Boolean HomogeneousFlag, Standard_Integer& ExtrapMode, const Standard_Integer Degree, const TColStd_Array1OfReal& FlatKnots, const Array1OfPoints& PolesArray, const TColStd_Array1OfReal& WeightsArray, Point& aPoint, Standard_Real& aWeight) { Standard_Real Inverse, P[Dimension_gen], *PArray, *WArray ; Standard_Integer kk ; PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ; WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ; if (HomogeneousFlag) { BSplCLib::Eval(Parameter, PeriodicFlag, 0, ExtrapMode, Degree, FlatKnots, Dimension_gen, PArray[0], P[0]) ; BSplCLib::Eval(Parameter, PeriodicFlag, 0, ExtrapMode, Degree, FlatKnots, 1, WArray[0], aWeight) ; } else { BSplCLib::Eval(Parameter, PeriodicFlag, 0, ExtrapMode, Degree, FlatKnots, Dimension_gen, PArray[0], WArray[0], P[0], aWeight) ; Inverse = 1.0e0 / aWeight ; for (kk = 0 ; kk < Dimension_gen ; kk++) { P[kk] *= Inverse ; } } for (kk = 0 ; kk < Dimension_gen ; kk++) aPoint.SetCoord(kk + 1, P[kk]) ; } //======================================================================= //function : CacheD0 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplCLib::CacheD0(const Standard_Real Parameter, const Standard_Integer Degree, const Standard_Real CacheParameter, const Standard_Real SpanLenght, const Array1OfPoints& PolesArray, const TColStd_Array1OfReal* WeightsArray, Point& aPoint) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 Standard_Real NewParameter, Inverse; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ; Standard_Real * myPoint = (Standard_Real *) &aPoint ; NewParameter = (Parameter - CacheParameter) / SpanLenght ; PLib::NoDerivativeEvalPolynomial(NewParameter, Degree, Dimension_gen, Degree * Dimension_gen, PArray[0], myPoint[0]) ; if (WeightsArray != NULL) { const TColStd_Array1OfReal& refWeights = *WeightsArray; Standard_Real * WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ; PLib::NoDerivativeEvalPolynomial(NewParameter, Degree, 1, Degree, WArray[0], Inverse) ; Inverse = 1.0e0 / Inverse; ModifyCoords (myPoint, *= Inverse); } } //======================================================================= //function : CacheD1 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplCLib::CacheD1(const Standard_Real Parameter, const Standard_Integer Degree, const Standard_Real CacheParameter, const Standard_Real SpanLenght, const Array1OfPoints& PolesArray, const TColStd_Array1OfReal* WeightsArray, Point& aPoint, Vector& aVector) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 Standard_Real LocalPDerivatives[Dimension_gen << 1] ; // Standard_Real LocalWDerivatives[2], NewParameter, Inverse ; Standard_Real LocalWDerivatives[2], NewParameter ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ; Standard_Real * myPoint = (Standard_Real *) &aPoint ; Standard_Real * myVector = (Standard_Real *) &aVector ; NewParameter = (Parameter - CacheParameter) / SpanLenght ; PLib::EvalPolynomial(NewParameter, 1, Degree, Dimension_gen, PArray[0], LocalPDerivatives[0]) ; // // unormalize derivatives since those are computed normalized // ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght); if (WeightsArray != NULL) { const TColStd_Array1OfReal& refWeights = *WeightsArray; Standard_Real * WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ; PLib::EvalPolynomial(NewParameter, 1, Degree, 1, WArray[0], LocalWDerivatives[0]) ; // // unormalize the result since the polynomial stored in the cache // is normalized between 0 and 1 // LocalWDerivatives[1] /= SpanLenght ; PLib::RationalDerivatives(1, Dimension_gen, LocalPDerivatives[0], LocalWDerivatives[0], LocalPDerivatives[0]) ; } CopyCoords (myPoint, LocalPDerivatives); CopyCoords (myVector, LocalPDerivatives + Dimension_gen); } //======================================================================= //function : CacheD2 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplCLib::CacheD2(const Standard_Real Parameter, const Standard_Integer Degree, const Standard_Real CacheParameter, const Standard_Real SpanLenght, const Array1OfPoints& PolesArray, const TColStd_Array1OfReal* WeightsArray, Point& aPoint, Vector& aVector1, Vector& aVector2) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 Standard_Integer ii,Index,EndIndex; Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ; // Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ; Standard_Real LocalWDerivatives[3], NewParameter, Factor ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ; Standard_Real * myPoint = (Standard_Real *) &aPoint ; Standard_Real * myVector1 = (Standard_Real *) &aVector1 ; Standard_Real * myVector2 = (Standard_Real *) &aVector2 ; NewParameter = (Parameter - CacheParameter) / SpanLenght ; PLib::EvalPolynomial(NewParameter, 2, Degree, Dimension_gen, PArray[0], LocalPDerivatives[0]) ; // // unormalize derivatives since those are computed normalized // Factor = 1.0e0/ SpanLenght ; Index = Dimension_gen ; EndIndex = Min (2, Degree) ; for (ii = 1 ; ii <= EndIndex ; ii++) { ModifyCoords (LocalPDerivatives + Index, *= Factor); Factor /= SpanLenght ; Index += Dimension_gen; } Index = (Degree + 1) * Dimension_gen; for (ii = Degree ; ii < 2 ; ii++) { NullifyCoords (LocalPDerivatives + Index); Index += Dimension_gen; } if (WeightsArray != NULL) { const TColStd_Array1OfReal& refWeights = *WeightsArray; Standard_Real * WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ; PLib::EvalPolynomial(NewParameter, 2, Degree, 1, WArray[0], LocalWDerivatives[0]) ; for (ii = Degree + 1 ; ii <= 2 ; ii++) { LocalWDerivatives[ii] = 0.0e0 ; } // // unormalize the result since the polynomial stored in the cache // is normalized between 0 and 1 // Factor = 1.0e0 / SpanLenght ; for (ii = 1 ; ii <= EndIndex ; ii++) { LocalWDerivatives[ii] *= Factor ; Factor /= SpanLenght ; } PLib::RationalDerivatives(2, Dimension_gen, LocalPDerivatives[0], LocalWDerivatives[0], LocalPDerivatives[0]) ; } CopyCoords (myPoint, LocalPDerivatives); CopyCoords (myVector1, LocalPDerivatives + Dimension_gen); CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2); } //======================================================================= //function : CacheD3 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplCLib::CacheD3(const Standard_Real Parameter, const Standard_Integer Degree, const Standard_Real CacheParameter, const Standard_Real SpanLenght, const Array1OfPoints& PolesArray, const TColStd_Array1OfReal* WeightsArray, Point& aPoint, Vector& aVector1, Vector& aVector2, Vector& aVector3) { // // the CacheParameter is where the cache polynomial was evaluated in homogeneous // form // the SpanLenght is the normalizing factor so that the polynomial is between // 0 and 1 Standard_Integer ii, Index, EndIndex; Standard_Real LocalPDerivatives[Dimension_gen << 2] ; // Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ; Standard_Real LocalWDerivatives[4], Factor, NewParameter ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ; Standard_Real * myPoint = (Standard_Real *) &aPoint ; Standard_Real * myVector1 = (Standard_Real *) &aVector1 ; Standard_Real * myVector2 = (Standard_Real *) &aVector2 ; Standard_Real * myVector3 = (Standard_Real *) &aVector3 ; NewParameter = (Parameter - CacheParameter) / SpanLenght ; PLib::EvalPolynomial(NewParameter, 3, Degree, Dimension_gen, PArray[0], LocalPDerivatives[0]) ; Index = (Degree + 1) * Dimension_gen; for (ii = Degree ; ii < 3 ; ii++) { NullifyCoords (LocalPDerivatives + Index); Index += Dimension_gen; } // // unormalize derivatives since those are computed normalized // Factor = 1.0e0 / SpanLenght ; Index = Dimension_gen ; EndIndex = Min(3,Degree) ; for (ii = 1 ; ii <= EndIndex ; ii++) { ModifyCoords (LocalPDerivatives + Index, *= Factor); Factor /= SpanLenght; Index += Dimension_gen; } if (WeightsArray != NULL) { const TColStd_Array1OfReal& refWeights = *WeightsArray; Standard_Real * WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ; PLib::EvalPolynomial(NewParameter, 3, Degree, 1, WArray[0], LocalWDerivatives[0]) ; // // unormalize the result since the polynomial stored in the cache // is normalized between 0 and 1 // Factor = 1.0e0 / SpanLenght ; for (ii = 1 ; ii <= EndIndex ; ii++) { LocalWDerivatives[ii] *= Factor ; Factor /= SpanLenght ; } for (ii = (Degree + 1) ; ii <= 3 ; ii++) { LocalWDerivatives[ii] = 0.0e0 ; } PLib::RationalDerivatives(3, Dimension_gen, LocalPDerivatives[0], LocalWDerivatives[0], LocalPDerivatives[0]) ; } CopyCoords (myPoint, LocalPDerivatives); CopyCoords (myVector1, LocalPDerivatives + Dimension_gen); CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2); CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3); } //======================================================================= //function : BuildCache //purpose : Stores theTaylor expansion normalized between 0,1 in the // Cache : in case of a rational function the Poles are // stored in homogeneous form //======================================================================= void BSplCLib::BuildCache (const Standard_Real U, const Standard_Real SpanDomain, const Standard_Boolean Periodic, const Standard_Integer Degree, const TColStd_Array1OfReal& FlatKnots, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, Array1OfPoints& CachePoles, TColStd_Array1OfReal* CacheWeights) { Standard_Integer ii, Dimension, LocalIndex, index = 0 ; Standard_Real u = U, LocalValue ; Standard_Boolean rational; BSplCLib_DataContainer dc(Degree); PrepareEval(u, index, Dimension, rational, Degree, Periodic, Poles, Weights, FlatKnots, (BSplCLib::NoMults()), dc); // // Watch out : PrepareEval checks if locally the Bspline is polynomial // therefore rational flag could be set to False even if there are // Weights. The Dimension is set accordingly. // BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles); LocalValue = 1.0e0 ; LocalIndex = 0 ; if (rational) { for (ii = 1 ; ii <= Degree + 1 ; ii++) { CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue); LocalIndex += Dimension_gen + 1; LocalValue *= SpanDomain / (Standard_Real) ii ; } LocalIndex = Dimension_gen; LocalValue = 1.0e0 ; for (ii = 1 ; ii <= Degree + 1 ; ii++) { (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ; LocalIndex += Dimension_gen + 1; LocalValue *= SpanDomain / (Standard_Real) ii ; } } else { for (ii = 1 ; ii <= Degree + 1 ; ii++) { CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue); LocalIndex += Dimension_gen; LocalValue *= SpanDomain / (Standard_Real) ii ; } if (Weights != NULL) { for (ii = 1 ; ii <= Degree + 1 ; ii++) (*CacheWeights)(ii) = 0.0e0 ; (*CacheWeights)(1) = 1.0e0 ; } } } void BSplCLib::BuildCache(const Standard_Real theParameter, const Standard_Real theSpanDomain, const Standard_Boolean thePeriodicFlag, const Standard_Integer theDegree, const Standard_Integer theSpanIndex, const TColStd_Array1OfReal& theFlatKnots, const Array1OfPoints& thePoles, const TColStd_Array1OfReal* theWeights, TColStd_Array2OfReal& theCacheArray) { Standard_Real aParam = theParameter; Standard_Integer anIndex = theSpanIndex; Standard_Integer aDimension; Standard_Boolean isRational; BSplCLib_DataContainer dc(theDegree); PrepareEval(aParam, anIndex, aDimension, isRational, theDegree, thePeriodicFlag, thePoles, theWeights, theFlatKnots, (BSplCLib::NoMults()), dc); // // Watch out : PrepareEval checks if locally the Bspline is polynomial // therefore rational flag could be set to False even if there are // Weights. The Dimension is set accordingly. // Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension; BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles); Standard_Real aCoeff = 1.0; Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol())); Standard_Real* aPolyCoeffs = dc.poles; for (Standard_Integer i = 0; i <= theDegree; i++) { for (Standard_Integer j = 0; j< aDimension; j++) aCache[j] = aPolyCoeffs[j] * aCoeff; aCoeff *= theSpanDomain / (i+1); aPolyCoeffs += aDimension; aCache += aDimension; if (aCacheShift > aDimension) { aCache[0] = 0.0; aCache++; } } if (aCacheShift > aDimension) theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0); } //======================================================================= //function : Interpolate //purpose : //======================================================================= void BSplCLib::Interpolate(const Standard_Integer Degree, const TColStd_Array1OfReal& FlatKnots, const TColStd_Array1OfReal& Parameters, const TColStd_Array1OfInteger& ContactOrderArray, Array1OfPoints& Poles, Standard_Integer& InversionProblem) { Standard_Real *PArray ; PArray = (Standard_Real *) &Poles(Poles.Lower()) ; BSplCLib::Interpolate(Degree, FlatKnots, Parameters, ContactOrderArray, Dimension_gen, PArray[0], InversionProblem) ; } //======================================================================= //function : Interpolate //purpose : //======================================================================= void BSplCLib::Interpolate(const Standard_Integer Degree, const TColStd_Array1OfReal& FlatKnots, const TColStd_Array1OfReal& Parameters, const TColStd_Array1OfInteger& ContactOrderArray, Array1OfPoints& Poles, TColStd_Array1OfReal& Weights, Standard_Integer& InversionProblem) { Standard_Real *PArray, * WArray ; PArray = (Standard_Real *) &Poles(Poles.Lower()) ; WArray = (Standard_Real *) &Weights(Weights.Lower()) ; BSplCLib::Interpolate(Degree, FlatKnots, Parameters, ContactOrderArray, Dimension_gen, PArray[0], WArray[0], InversionProblem) ; } //======================================================================= //function : MovePoint //purpose : Find the new poles which allows an old point (with a // given u as parameter) to reach a new position //======================================================================= void BSplCLib::MovePoint (const Standard_Real U, const Vector& Displ, const Standard_Integer Index1, const Standard_Integer Index2, const Standard_Integer Degree, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& FlatKnots, Standard_Integer& FirstIndex, Standard_Integer& LastIndex, Array1OfPoints& NewPoles) { // calculate the BSplineBasis in the parameter U Standard_Integer FirstNonZeroBsplineIndex; math_Matrix BSplineBasis(1, 1, 1, Degree+1); Standard_Integer ErrorCode = BSplCLib::EvalBsplineBasis(0, Degree+1, FlatKnots, U, FirstNonZeroBsplineIndex, BSplineBasis); if (ErrorCode != 0) { FirstIndex = 0; LastIndex = 0; for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) { NewPoles(i) = Poles(i); } return; } // find the span which is predominant for parameter U FirstIndex = FirstNonZeroBsplineIndex; LastIndex = FirstNonZeroBsplineIndex + Degree ; if (FirstIndex < Index1) FirstIndex = Index1; if (LastIndex > Index2) LastIndex = Index2; Standard_Real maxValue = 0.0; Standard_Integer i, kk1=0, kk2, ii; for (i = FirstIndex-FirstNonZeroBsplineIndex+1; i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) { if (BSplineBasis(1,i) > maxValue) { kk1 = i + FirstNonZeroBsplineIndex - 1; maxValue = BSplineBasis(1,i); } } // find a kk2 if symetriy kk2 = kk1; i = kk1 - FirstNonZeroBsplineIndex + 2; if ((kk1+1) <= LastIndex) { if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) { kk2 = kk1+1; } } // compute the vector of displacement Standard_Real D1 = 0.0; Standard_Real D2 = 0.0; Standard_Real hN, Coef, Dval; for (i = 1; i <= Degree+1; i++) { ii = i + FirstNonZeroBsplineIndex - 1; if (Weights != NULL) { hN = Weights->Value(ii)*BSplineBasis(1, i); D2 += hN; } else { hN = BSplineBasis(1, i); } if (ii >= FirstIndex && ii <= LastIndex) { if (ii < kk1) { Dval = kk1-ii; } else if (ii > kk2) { Dval = ii - kk2; } else { Dval = 0.0; } D1 += 1./(Dval + 1.) * hN; } } if (Weights != NULL) { Coef = D2/D1; } else { Coef = 1./D1; } // compute the new poles for (i=Poles.Lower(); i<=Poles.Upper(); i++) { if (i >= FirstIndex && i <= LastIndex) { if (i < kk1) { Dval = kk1-i; } else if (i > kk2) { Dval = i - kk2; } else { Dval = 0.0; } NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ); } else { NewPoles(i) = Poles(i); } } } //======================================================================= //function : MovePoint //purpose : Find the new poles which allows an old point (with a // given u as parameter) to reach a new position //======================================================================= //======================================================================= //function : MovePointAndTangent //purpose : //======================================================================= void BSplCLib::MovePointAndTangent (const Standard_Real U, const Vector& Delta, const Vector& DeltaDerivatives, const Standard_Real Tolerance, const Standard_Integer Degree, const Standard_Integer StartingCondition, const Standard_Integer EndingCondition, const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const TColStd_Array1OfReal& FlatKnots, Array1OfPoints& NewPoles, Standard_Integer & ErrorStatus) { Standard_Real *delta_array, *delta_derivative_array, *poles_array, *new_poles_array ; Standard_Integer num_poles ; num_poles = Poles.Length() ; if (NewPoles.Length() != num_poles) { throw Standard_ConstructionError(); } delta_array = (Standard_Real *) &Delta ; delta_derivative_array = (Standard_Real *)&DeltaDerivatives ; poles_array = (Standard_Real *) &Poles(Poles.Lower()) ; new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ; MovePointAndTangent(U, Dimension_gen, delta_array[0], delta_derivative_array[0], Tolerance, Degree, StartingCondition, EndingCondition, poles_array[0], Weights, FlatKnots, new_poles_array[0], ErrorStatus) ; } //======================================================================= //function : Resolution //purpose : //======================================================================= void BSplCLib::Resolution(const Array1OfPoints& Poles, const TColStd_Array1OfReal* Weights, const Standard_Integer NumPoles, const TColStd_Array1OfReal& FlatKnots, const Standard_Integer Degree, const Standard_Real Tolerance3D, Standard_Real& UTolerance) { Standard_Real *PolesArray ; PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ; BSplCLib::Resolution(PolesArray[0], Dimension_gen, NumPoles, Weights, FlatKnots, Degree, Tolerance3D, UTolerance) ; } //======================================================================= //function : FunctionMultiply //purpose : //======================================================================= void BSplCLib::FunctionMultiply (const BSplCLib_EvaluatorFunction & FunctionPtr, const Standard_Integer BSplineDegree, const TColStd_Array1OfReal & BSplineFlatKnots, const Array1OfPoints & Poles, const TColStd_Array1OfReal & FlatKnots, const Standard_Integer NewDegree, Array1OfPoints & NewPoles, Standard_Integer & theStatus) { Standard_Integer num_bspline_poles = BSplineFlatKnots.Length() - BSplineDegree - 1 ; Standard_Integer num_new_poles = FlatKnots.Length() - NewDegree - 1 ; if (Poles.Length() != num_bspline_poles || NewPoles.Length() != num_new_poles) { throw Standard_ConstructionError(); } Standard_Real * array_of_poles = (Standard_Real *) &Poles(Poles.Lower()) ; Standard_Real * array_of_new_poles = (Standard_Real *) &NewPoles(NewPoles.Lower()) ; BSplCLib::FunctionMultiply(FunctionPtr, BSplineDegree, BSplineFlatKnots, Dimension_gen, array_of_poles[0], FlatKnots, NewDegree, array_of_new_poles[0], theStatus); } //======================================================================= //function : FunctionReparameterise //purpose : //======================================================================= void BSplCLib::FunctionReparameterise (const BSplCLib_EvaluatorFunction & FunctionPtr, const Standard_Integer BSplineDegree, const TColStd_Array1OfReal & BSplineFlatKnots, const Array1OfPoints & Poles, const TColStd_Array1OfReal & FlatKnots, const Standard_Integer NewDegree, Array1OfPoints & NewPoles, Standard_Integer & theStatus) { Standard_Integer num_bspline_poles = BSplineFlatKnots.Length() - BSplineDegree - 1 ; Standard_Integer num_new_poles = FlatKnots.Length() - NewDegree - 1 ; if (Poles.Length() != num_bspline_poles || NewPoles.Length() != num_new_poles) { throw Standard_ConstructionError(); } Standard_Real * array_of_poles = (Standard_Real *) &Poles(Poles.Lower()) ; Standard_Real * array_of_new_poles = (Standard_Real *) &NewPoles(NewPoles.Lower()) ; BSplCLib::FunctionReparameterise(FunctionPtr, BSplineDegree, BSplineFlatKnots, Dimension_gen, array_of_poles[0], FlatKnots, NewDegree, array_of_new_poles[0], theStatus); }