From 652b07f92ce2c30dbec7b64c29490156b738d5be Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Wed, 26 Nov 2025 16:43:09 +0000 Subject: [PATCH] Foundation Classes - Enhance B-Spline Curve Computation (#855) - Introduces `BSplCLib_KnotArrays` template for efficient stack-based knot/multiplicity management - Refactors `BSplCLib_Reverse` to use `std::reverse` for improved performance - Replaces heap-allocated arrays with stack-based `NCollection_LocalArray` in critical code paths - Modernizes copy constructor/assignment operator prevention using `= delete` --- .../TKMath/BSplCLib/BSplCLib.cxx | 431 ++--- .../TKMath/BSplCLib/BSplCLib_1.cxx | 2 +- .../TKMath/BSplCLib/BSplCLib_2.cxx | 376 ++--- .../TKMath/BSplCLib/BSplCLib_3.cxx | 2 +- .../TKMath/BSplCLib/BSplCLib_BzSyntaxes.cxx | 157 +- .../TKMath/BSplCLib/BSplCLib_Cache.cxx | 64 +- .../TKMath/BSplCLib/BSplCLib_Cache.hxx | 18 +- .../TKMath/BSplCLib/BSplCLib_CacheParams.hxx | 49 +- .../BSplCLib/BSplCLib_CurveComputation.gxx | 1483 ----------------- .../BSplCLib/BSplCLib_CurveComputation.pxx | 159 +- .../BSplCLib/BSplCLib_EvaluatorFunction.hxx | 9 +- 11 files changed, 586 insertions(+), 2164 deletions(-) delete mode 100644 src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.gxx diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx index 01c6b3b462..ac3309fe69 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx @@ -35,8 +35,11 @@ #include #include #include +#include #include +#include + typedef gp_Pnt Pnt; typedef gp_Vec Vec; typedef TColgp_Array1OfPnt Array1OfPnt; @@ -159,11 +162,7 @@ Standard_Integer BSplCLib::FlatIndex(const Standard_Integer Degree, return index; } -//======================================================================= -// function : LocateParameter -// purpose : Processing of nodes with multiplicities -// pmn 28-01-97 -> compute eventual of the period. -//======================================================================= +//================================================================================================= void BSplCLib::LocateParameter(const Standard_Integer, // Degree, const Array1OfReal& Knots, @@ -184,12 +183,7 @@ void BSplCLib::LocateParameter(const Standard_Integer, // Degree, BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2, KnotIndex, NewU, uf, ul); } -//======================================================================= -// function : LocateParameter -// purpose : For plane nodes -// pmn 28-01-97 -> There is a need of the degre to calculate -// the eventual period -//======================================================================= +//================================================================================================= void BSplCLib::LocateParameter(const Standard_Integer Degree, const Array1OfReal& Knots, @@ -214,12 +208,7 @@ void BSplCLib::LocateParameter(const Standard_Integer Degree, BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2, KnotIndex, NewU, 0., 1.); } -//======================================================================= -// function : LocateParameter -// purpose : Effective computation -// pmn 28-01-97 : Add limits of the period as input argument, -// as it is impossible to produce them at this level. -//======================================================================= +//================================================================================================= void BSplCLib::LocateParameter(const TColStd_Array1OfReal& Knots, const Standard_Real U, @@ -306,11 +295,7 @@ void BSplCLib::LocateParameter(const TColStd_Array1OfReal& Knots, } } -//======================================================================= -// function : LocateParameter -// purpose : the index is recomputed only if out of range -// pmn 28-01-97 -> eventual computation of the period. -//======================================================================= +//================================================================================================= void BSplCLib::LocateParameter(const Standard_Integer Degree, const TColStd_Array1OfReal& Knots, @@ -359,17 +344,12 @@ Standard_Integer BSplCLib::MaxKnotMult(const Array1OfInteger& Mults, const Standard_Integer FromK1, const Standard_Integer ToK2) { - Standard_Integer MLower = Mults.Lower(); - const Standard_Integer* pmu = &Mults(MLower); - pmu -= MLower; - Standard_Integer MaxMult = pmu[FromK1]; - - for (Standard_Integer i = FromK1; i <= ToK2; i++) + int aMaxMult = Mults(FromK1); + for (int anIndex = FromK1 + 1; anIndex <= ToK2; ++anIndex) { - if (MaxMult < pmu[i]) - MaxMult = pmu[i]; + aMaxMult = std::max(aMaxMult, Mults(anIndex)); } - return MaxMult; + return aMaxMult; } //================================================================================================= @@ -378,17 +358,12 @@ Standard_Integer BSplCLib::MinKnotMult(const Array1OfInteger& Mults, const Standard_Integer FromK1, const Standard_Integer ToK2) { - Standard_Integer MLower = Mults.Lower(); - const Standard_Integer* pmu = &Mults(MLower); - pmu -= MLower; - Standard_Integer MinMult = pmu[FromK1]; - - for (Standard_Integer i = FromK1; i <= ToK2; i++) + int aMinMult = Mults(FromK1); + for (int anIndex = FromK1 + 1; anIndex <= ToK2; ++anIndex) { - if (MinMult > pmu[i]) - MinMult = pmu[i]; + aMinMult = std::min(aMinMult, Mults(anIndex)); } - return MinMult; + return aMinMult; } //================================================================================================= @@ -588,47 +563,32 @@ BSplCLib_KnotDistribution BSplCLib::KnotForm(const Array1OfReal& Knots, const Standard_Integer FromK1, const Standard_Integer ToK2) { - Standard_Real DU0, DU1, Ui, Uj, Eps0, val; - BSplCLib_KnotDistribution KForm = BSplCLib_Uniform; - if (FromK1 + 1 > Knots.Upper()) { return BSplCLib_Uniform; } - Ui = Knots(FromK1); - if (Ui < 0) - Ui = -Ui; - Uj = Knots(FromK1 + 1); - if (Uj < 0) - Uj = -Uj; - DU0 = Uj - Ui; - if (DU0 < 0) - DU0 = -DU0; - Eps0 = Epsilon(Ui) + Epsilon(Uj) + Epsilon(DU0); - Standard_Integer i = FromK1 + 1; + double aUi = std::abs(Knots(FromK1)); + double aUj = std::abs(Knots(FromK1 + 1)); + double aDU0 = std::abs(aUj - aUi); + double anEps = Epsilon(aUi) + Epsilon(aUj) + Epsilon(aDU0); - while (KForm != BSplCLib_NonUniform && i < ToK2) + for (int i = FromK1 + 1; i < ToK2; i++) { - Ui = Knots(i); - if (Ui < 0) - Ui = -Ui; - i++; - Uj = Knots(i); - if (Uj < 0) - Uj = -Uj; - DU1 = Uj - Ui; - if (DU1 < 0) - DU1 = -DU1; - val = DU1 - DU0; - if (val < 0) - val = -val; - if (val > Eps0) - KForm = BSplCLib_NonUniform; - DU0 = DU1; - Eps0 = Epsilon(Ui) + Epsilon(Uj) + Epsilon(DU0); + aUi = std::abs(Knots(i)); + aUj = std::abs(Knots(i + 1)); + const double aDU1 = std::abs(aUj - aUi); + + if (std::abs(aDU1 - aDU0) > anEps) + { + return BSplCLib_NonUniform; + } + + aDU0 = aDU1; + anEps = Epsilon(aUi) + Epsilon(aUj) + Epsilon(aDU0); } - return KForm; + + return BSplCLib_Uniform; } //================================================================================================= @@ -637,57 +597,54 @@ BSplCLib_MultDistribution BSplCLib::MultForm(const Array1OfInteger& Mults, const Standard_Integer FromK1, const Standard_Integer ToK2) { - Standard_Integer First, Last; - if (FromK1 < ToK2) - { - First = FromK1; - Last = ToK2; - } - else - { - First = ToK2; - Last = FromK1; - } - if (First + 1 > Mults.Upper()) + const auto [aFirst, aLast] = std::minmax(FromK1, ToK2); + + if (aFirst + 1 > Mults.Upper()) { return BSplCLib_Constant; } - Standard_Integer FirstMult = Mults(First); - BSplCLib_MultDistribution MForm = BSplCLib_Constant; - Standard_Integer i = First + 1; - Standard_Integer Mult = Mults(i); + const int aFirstMult = Mults(aFirst); + BSplCLib_MultDistribution aForm = BSplCLib_Constant; + int aMult = Mults(aFirst + 1); - // while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR???????? - while (MForm != BSplCLib_NonConstant && i <= Last) + for (int i = aFirst + 1; i <= aLast && aForm != BSplCLib_NonConstant; i++) { - if (i == First + 1) + if (i == aFirst + 1) { - if (Mult != FirstMult) - MForm = BSplCLib_QuasiConstant; + if (aMult != aFirstMult) + { + aForm = BSplCLib_QuasiConstant; + } } - else if (i == Last) + else if (i == aLast) { - if (MForm == BSplCLib_QuasiConstant) + if (aForm == BSplCLib_QuasiConstant) { - if (FirstMult != Mults(i)) - MForm = BSplCLib_NonConstant; + if (aFirstMult != Mults(i)) + { + aForm = BSplCLib_NonConstant; + } } else { - if (Mult != Mults(i)) - MForm = BSplCLib_NonConstant; + if (aMult != Mults(i)) + { + aForm = BSplCLib_NonConstant; + } } } else { - if (Mult != Mults(i)) - MForm = BSplCLib_NonConstant; - Mult = Mults(i); + if (aMult != Mults(i)) + { + aForm = BSplCLib_NonConstant; + } + aMult = Mults(i); } - i++; } - return MForm; + + return aForm; } //================================================================================================= @@ -822,37 +779,14 @@ void BSplCLib::Reverse(TColStd_Array1OfReal& Knots) void BSplCLib::Reverse(TColStd_Array1OfInteger& Mults) { - Standard_Integer first = Mults.Lower(); - Standard_Integer last = Mults.Upper(); - Standard_Integer temp; - - while (first < last) - { - temp = Mults(first); - Mults(first) = Mults(last); - Mults(last) = temp; - first++; - last--; - } + std::reverse(&Mults(Mults.Lower()), &Mults(Mults.Upper()) + 1); } //================================================================================================= void BSplCLib::Reverse(TColStd_Array1OfReal& Weights, const Standard_Integer L) { - Standard_Integer i, l = L; - l = Weights.Lower() + (l - Weights.Lower()) % (Weights.Upper() - Weights.Lower() + 1); - - TColStd_Array1OfReal temp(0, Weights.Length() - 1); - - for (i = Weights.Lower(); i <= l; i++) - temp(l - i) = Weights(i); - - for (i = l + 1; i <= Weights.Upper(); i++) - temp(l - Weights.Lower() + Weights.Upper() - i + 1) = Weights(i); - - for (i = Weights.Lower(); i <= Weights.Upper(); i++) - Weights(i) = temp(i - Weights.Lower()); + BSplCLib_Reverse(Weights, L); } //================================================================================================= @@ -876,10 +810,7 @@ Standard_Boolean BSplCLib::IsRational(const TColStd_Array1OfReal& Weights, return Standard_False; } -//======================================================================= -// function : Eval -// purpose : evaluate point and derivatives -//======================================================================= +//================================================================================================= void BSplCLib::Eval(const Standard_Real U, const Standard_Integer Degree, @@ -1775,10 +1706,7 @@ Standard_Integer BSplCLib::PoleIndex(const Standard_Integer Degree, return pindex; } -//======================================================================= -// function : BuildBoor -// purpose : builds the local array for boor -//======================================================================= +//================================================================================================= void BSplCLib::BuildBoor(const Standard_Integer Index, const Standard_Integer Length, @@ -1880,7 +1808,7 @@ Standard_Boolean BSplCLib::PrepareInsertKnots(const Standard_Integer Deg { // gka for case when segments was produced on full period only one knot // was added in the end of curve - if (fabs(adeltaK1) <= gp::Resolution() && fabs(adeltaK2) <= gp::Resolution()) + if (std::abs(adeltaK1) <= gp::Resolution() && std::abs(adeltaK2) <= gp::Resolution()) ak++; } @@ -1995,18 +1923,7 @@ Standard_Boolean BSplCLib::PrepareInsertKnots(const Standard_Integer Deg return Standard_True; } -//======================================================================= -// function : Copy -// purpose : copy reals from an array to an other -// -// NbValues are copied from OldPoles(OldFirst) -// to NewPoles(NewFirst) -// -// Periodicity is handled. -// OldFirst and NewFirst are updated -// to the position after the last copied pole. -// -//======================================================================= +//================================================================================================= static void Copy(const Standard_Integer NbPoles, Standard_Integer& OldFirst, @@ -2037,10 +1954,7 @@ static void Copy(const Standard_Integer NbPoles, } } -//======================================================================= -// function : InsertKnots -// purpose : insert an array of knots and multiplicities -//======================================================================= +//================================================================================================= void BSplCLib::InsertKnots(const Standard_Integer Degree, const Standard_Boolean Periodic, @@ -2068,9 +1982,9 @@ void BSplCLib::InsertKnots(const Standard_Integer Degree, // ------------------- // create local arrays // ------------------- - - Standard_Real* knots = new Standard_Real[2 * Degree]; - Standard_Real* poles = new Standard_Real[(2 * Degree + 1) * Dimension]; + // Use stack-based allocation for small arrays (MaxDegree=25, typical Dimension=4) + NCollection_LocalArray knots(2 * Degree); + NCollection_LocalArray poles((2 * Degree + 1) * Dimension); //---------------------------- // loop on the knots to insert @@ -2306,9 +2220,6 @@ void BSplCLib::InsertKnots(const Standard_Integer Degree, NewMults(NewMults.Upper()) += depth; } } - // free local arrays - delete[] knots; - delete[] poles; } //================================================================================================= @@ -2356,9 +2267,9 @@ Standard_Boolean BSplCLib::RemoveKnot(const Standard_Integer Index, // ------------------- // create local arrays // ------------------- - - Standard_Real* knots = new Standard_Real[4 * Degree]; - Standard_Real* poles = new Standard_Real[(2 * Degree + 1) * Dimension]; + // Use stack-based allocation for small arrays (MaxDegree=25, typical Dimension=4) + NCollection_LocalArray knots(4 * Degree); + NCollection_LocalArray poles((2 * Degree + 1) * Dimension); // ------------------------------------ // build the knots for anti Boor Scheme @@ -2481,10 +2392,6 @@ Standard_Boolean BSplCLib::RemoveKnot(const Standard_Integer Index, } } - // free local arrays - delete[] knots; - delete[] poles; - return result; } @@ -2674,10 +2581,11 @@ void BSplCLib::IncreaseDegree(const Standard_Integer Degree, // Declare the temporary arrays used in degree incrementation //----------------------------------------------------------- - Standard_Integer nbwp = nbwpoles + (nbwknots - 1) * (NewDegree - Degree); + // Maximum number of poles after all degree elevations + const Standard_Integer nbwpMax = nbwpoles + (nbwknots - 1) * (NewDegree - Degree); // Arrays for storing the temporary curves - TColStd_Array1OfReal tempc1(1, nbwp * Dimension); - TColStd_Array1OfReal tempc2(1, nbwp * Dimension); + TColStd_Array1OfReal tempc1(1, nbwpMax * Dimension); + TColStd_Array1OfReal tempc2(1, nbwpMax * Dimension); // Array for storing the knots to insert TColStd_Array1OfReal iknots(1, nbwknots); @@ -2685,13 +2593,16 @@ void BSplCLib::IncreaseDegree(const Standard_Integer Degree, // Arrays for receiving the knots after insertion TColStd_Array1OfReal nknots(1, nbwknots); + // Pre-allocate nwpoles with maximum size to avoid repeated heap allocations in loop + TColStd_Array1OfReal nwpoles(1, nbwpMax * Dimension); + //------------------------------ // Loop on degree incrementation //------------------------------ Standard_Integer step, curDeg; - Standard_Integer nbp = nbwpoles; - nbwp = nbp; + Standard_Integer nbp = nbwpoles; + Standard_Integer nbwp = nbp; for (curDeg = Degree; curDeg < NewDegree; curDeg++) { @@ -2699,9 +2610,9 @@ void BSplCLib::IncreaseDegree(const Standard_Integer Degree, nbp = nbwp; // current number of poles nbwp = nbp + nbwknots - 1; // new number of poles - // For the averaging - TColStd_Array1OfReal nwpoles(1, nbwp * Dimension); - nwpoles.Init(0.0e0); + // Initialize the portion of nwpoles that will be used this iteration + for (Standard_Integer idx = 1; idx <= nbwp * Dimension; idx++) + nwpoles(idx) = 0.0; for (step = 0; step <= curDeg; step++) { @@ -3095,25 +3006,21 @@ Standard_Integer BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, const Standard_Integer ArrayDimension, Standard_Real& Array) { - Standard_Integer ii, jj, kk, MinIndex, MaxIndex, ReturnCode = 0; - - Standard_Real* PolesArray = &Array; - Standard_Real Inverse; - if (Matrix.LowerCol() != 1 || Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) { - ReturnCode = 1; - goto FINISH; + return 1; } - for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow(); ii++) + Standard_Real* PolesArray = &Array; + + for (int ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow(); ii++) { - MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ? ii - LowerBandWidth : Matrix.LowerRow()); + const int aMinIndex = + (ii - LowerBandWidth >= Matrix.LowerRow() ? ii - LowerBandWidth : Matrix.LowerRow()); - for (jj = MinIndex; jj < ii; jj++) + for (int jj = aMinIndex; jj < ii; jj++) { - - for (kk = 0; kk < ArrayDimension; kk++) + for (int kk = 0; kk < ArrayDimension; kk++) { PolesArray[(ii - 1) * ArrayDimension + kk] += PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1); @@ -3121,39 +3028,36 @@ Standard_Integer BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, } } - for (ii = Matrix.UpperRow(); ii >= Matrix.LowerRow(); ii--) + for (int ii = Matrix.UpperRow(); ii >= Matrix.LowerRow(); ii--) { - MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ? ii + UpperBandWidth : Matrix.UpperRow()); + const int aMaxIndex = + (ii + UpperBandWidth <= Matrix.UpperRow() ? ii + UpperBandWidth : Matrix.UpperRow()); - for (jj = MaxIndex; jj > ii; jj--) + for (int jj = aMaxIndex; jj > ii; jj--) { - - for (kk = 0; kk < ArrayDimension; kk++) + for (int kk = 0; kk < ArrayDimension; kk++) { PolesArray[(ii - 1) * ArrayDimension + kk] -= PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1); } } - // fixing a bug PRO18577 to avoid divizion by zero - - Standard_Real divizor = Matrix(ii, LowerBandWidth + 1); - Standard_Real Toler = 1.0e-16; - if (std::abs(divizor) > Toler) - Inverse = 1.0e0 / divizor; - else + // Fixing a bug PRO18577 to avoid division by zero + const double aDivisor = Matrix(ii, LowerBandWidth + 1); + constexpr double THE_TOLERANCE = 1.0e-16; + if (std::abs(aDivisor) <= THE_TOLERANCE) { - ReturnCode = 1; - goto FINISH; + return 1; } + const double anInverse = 1.0 / aDivisor; - for (kk = 0; kk < ArrayDimension; kk++) + for (int kk = 0; kk < ArrayDimension; kk++) { - PolesArray[(ii - 1) * ArrayDimension + kk] *= Inverse; + PolesArray[(ii - 1) * ArrayDimension + kk] *= anInverse; } } -FINISH: - return (ReturnCode); + + return 0; } //================================================================================================= @@ -3166,55 +3070,52 @@ Standard_Integer BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, Standard_Real& Poles, Standard_Real& Weights) { - Standard_Integer ii, kk, ErrorCode = 0, ReturnCode = 0; - - Standard_Real Inverse, *PolesArray = &Poles, *WeightsArray = &Weights; - if (Matrix.LowerCol() != 1 || Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) { - ReturnCode = 1; - goto FINISH; + return 1; } - if (HomogeneousFlag == Standard_False) - { - for (ii = 0; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) - { + double* aPolesArray = &Poles; + double* aWeightsArray = &Weights; + const int aNumRows = Matrix.UpperRow() - Matrix.LowerRow() + 1; - for (kk = 0; kk < ArrayDimension; kk++) + if (!HomogeneousFlag) + { + for (int ii = 0; ii < aNumRows; ii++) + { + for (int kk = 0; kk < ArrayDimension; kk++) { - PolesArray[ii * ArrayDimension + kk] *= WeightsArray[ii]; + aPolesArray[ii * ArrayDimension + kk] *= aWeightsArray[ii]; } } } - ErrorCode = + + int anErrorCode = BSplCLib::SolveBandedSystem(Matrix, UpperBandWidth, LowerBandWidth, ArrayDimension, Poles); - if (ErrorCode != 0) + if (anErrorCode != 0) { - ReturnCode = 2; - goto FINISH; + return 2; } - ErrorCode = BSplCLib::SolveBandedSystem(Matrix, UpperBandWidth, LowerBandWidth, 1, Weights); - if (ErrorCode != 0) + + anErrorCode = BSplCLib::SolveBandedSystem(Matrix, UpperBandWidth, LowerBandWidth, 1, Weights); + if (anErrorCode != 0) { - ReturnCode = 3; - goto FINISH; + return 3; } - if (HomogeneousFlag == Standard_False) - { - for (ii = 0; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) + if (!HomogeneousFlag) + { + for (int ii = 0; ii < aNumRows; ii++) { - Inverse = 1.0e0 / WeightsArray[ii]; - - for (kk = 0; kk < ArrayDimension; kk++) + const double anInverse = 1.0 / aWeightsArray[ii]; + for (int kk = 0; kk < ArrayDimension; kk++) { - PolesArray[ii * ArrayDimension + kk] *= Inverse; + aPolesArray[ii * ArrayDimension + kk] *= anInverse; } } } -FINISH: - return (ReturnCode); + + return 0; } //================================================================================================= @@ -3320,14 +3221,7 @@ void BSplCLib::Interpolate(const Standard_Integer Degree, throw Standard_OutOfRange("BSplCLib::Interpolate"); } -//======================================================================= -// 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 -// this evaluates the numerator by multiplying by the weights -// and evaluating it but does not call RationalDerivatives after -//======================================================================= +//================================================================================================= void BSplCLib::Eval(const Standard_Real Parameter, const Standard_Boolean PeriodicFlag, @@ -3407,7 +3301,7 @@ void BSplCLib::Eval(const Standard_Real Parameter, BsplineBasis); if (ErrorCode != 0) { - goto FINISH; + return; } if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) { @@ -3503,24 +3397,15 @@ void BSplCLib::Eval(const Standard_Real Parameter, Index1 = Index1 % Modulus; Index1 += 1; } - LocalRealArray[Index + kk] *= Inverse; + LocalRealArray[Index] *= Inverse; Index += 1; Inverse /= (Standard_Real)ii; } PLib::EvalPolynomial(Delta, NewRequest, Degree, 1, LocalRealArray[0], WeightsResults); } -FINISH:; } -//======================================================================= -// 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 -// WARNING : the array Results is supposed to have at least -// (DerivativeRequest + 1) * ArrayDimension slots and the -// -//======================================================================= +//================================================================================================= void BSplCLib::Eval(const Standard_Real Parameter, const Standard_Boolean PeriodicFlag, @@ -3599,7 +3484,7 @@ void BSplCLib::Eval(const Standard_Real Parameter, BsplineBasis); if (ErrorCode != 0) { - goto FINISH; + return; } if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) { @@ -3673,16 +3558,9 @@ void BSplCLib::Eval(const Standard_Real Parameter, } PLib::EvalPolynomial(Delta, NewRequest, Degree, ArrayDimension, LocalRealArray[0], Results); } -FINISH:; } -//======================================================================= -// function : TangExtendToConstraint -// purpose : Extends a Bspline function using the tangency map -// WARNING : -// -// -//======================================================================= +//================================================================================================= void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, const Standard_Real C1Coefficient, @@ -3747,8 +3625,9 @@ void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, Standard_Boolean periodic_flag = Standard_False; Standard_Integer ipos, extrap_mode[2], derivative_request = std::max(Continuity, 1); extrap_mode[0] = extrap_mode[1] = CDegree; - TColStd_Array1OfReal EvalBS(1, CDimension * (derivative_request + 1)); - Standard_Real* Eadr = (Standard_Real*)&EvalBS(1); + // Use math_Vector for stack allocation + math_Vector EvalBS(1, CDimension * (derivative_request + 1)); + Standard_Real* Eadr = &EvalBS(1); BSplCLib::Eval(Tbord, periodic_flag, derivative_request, @@ -3770,23 +3649,24 @@ void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, // matrix of constraints math_Matrix Contraintes(1, Csize, 1, CDimension); + // Precompute powers of C1Coefficient for efficiency (avoid std::pow overhead) + const Standard_Real C1Coeff2 = C1Coefficient * C1Coefficient; + const Standard_Real C1Coeff3 = C1Coeff2 * C1Coefficient; if (After) { - for (ipos = 1; ipos <= CDimension; ipos++) { Contraintes(1, ipos) = EvalBS(ipos); Contraintes(2, ipos) = C1Coefficient * EvalBS(ipos + CDimension); if (Continuity >= 2) - Contraintes(3, ipos) = EvalBS(ipos + 2 * CDimension) * std::pow(C1Coefficient, 2); + Contraintes(3, ipos) = EvalBS(ipos + 2 * CDimension) * C1Coeff2; if (Continuity >= 3) - Contraintes(4, ipos) = EvalBS(ipos + 3 * CDimension) * std::pow(C1Coefficient, 3); + Contraintes(4, ipos) = EvalBS(ipos + 3 * CDimension) * C1Coeff3; Contraintes(Continuity + 2, ipos) = ConstraintPoint(ipos); } } else { - for (ipos = 1; ipos <= CDimension; ipos++) { Contraintes(1, ipos) = ConstraintPoint(ipos); @@ -3794,9 +3674,9 @@ void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, if (Continuity >= 1) Contraintes(3, ipos) = C1Coefficient * EvalBS(ipos + CDimension); if (Continuity >= 2) - Contraintes(4, ipos) = EvalBS(ipos + 2 * CDimension) * std::pow(C1Coefficient, 2); + Contraintes(4, ipos) = EvalBS(ipos + 2 * CDimension) * C1Coeff2; if (Continuity >= 3) - Contraintes(5, ipos) = EvalBS(ipos + 3 * CDimension) * std::pow(C1Coefficient, 3); + Contraintes(5, ipos) = EvalBS(ipos + 3 * CDimension) * C1Coeff3; } } @@ -4100,9 +3980,8 @@ void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, } } -//======================================================================= -// function : Resolution -// purpose : +//================================================================================================= +// Mathematical basis: // d // Let C(t) = SUM Ci Bi(t) a Bspline curve of degree d // i = 1,n @@ -4175,7 +4054,7 @@ void BSplCLib::TangExtendToConstraint(const TColStd_Array1OfReal& FlatKnots, // i =1,n // // -//======================================================================= +//================================================================================================= void BSplCLib::Resolution(Standard_Real& Poles, const Standard_Integer ArrayDimension, diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_1.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_1.cxx index b03cb8eca3..3552268eec 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_1.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_1.cxx @@ -30,7 +30,7 @@ void BSplCLib::Reverse(TColgp_Array1OfPnt2d& Poles, const Standard_Integer L) { - BSplCLib_Reverse(Poles, L); + BSplCLib_Reverse(Poles, L); } //================================================================================================== diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx index 1e457055e3..373a841a6b 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx @@ -30,6 +30,7 @@ #include #include +#include #include "BSplCLib_CurveComputation.pxx" @@ -38,10 +39,7 @@ using BSplCLib_DataContainer = BSplCLib_DataContainer_T<1>; // methods for 1 dimensional BSplines -//======================================================================= -// function : BuildEval -// purpose : builds the local array for evaluation -//======================================================================= +//================================================================================================= void BSplCLib::BuildEval(const Standard_Integer Degree, const Standard_Integer Index, @@ -81,11 +79,7 @@ void BSplCLib::BuildEval(const Standard_Integer Degree, } } -//======================================================================= -// 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, @@ -290,10 +284,7 @@ void BSplCLib::DN(const Standard_Real U, } } -//======================================================================= -// function : Build BSpline Matrix -// purpose : Builds the Bspline Matrix -//======================================================================= +//================================================================================================= Standard_Integer BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Parameters, const TColStd_Array1OfInteger& ContactOrderArray, @@ -303,109 +294,92 @@ Standard_Integer BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Paramet Standard_Integer& UpperBandWidth, Standard_Integer& LowerBandWidth) { - Standard_Integer ii, jj, Index, ErrorCode, ReturnCode = 0, FirstNonZeroBsplineIndex, BandWidth, - MaxOrder = BSplCLib::MaxDegree() + 1, Order; + constexpr int aMaxOrder = BSplCLib::MaxDegree() + 1; + const int anOrder = Degree + 1; + UpperBandWidth = Degree; + LowerBandWidth = Degree; + const int aBandWidth = UpperBandWidth + LowerBandWidth + 1; - math_Matrix BSplineBasis(1, MaxOrder, 1, MaxOrder); - - Order = Degree + 1; - UpperBandWidth = Degree; - LowerBandWidth = Degree; - BandWidth = UpperBandWidth + LowerBandWidth + 1; if (Matrix.LowerRow() != Parameters.Lower() || Matrix.UpperRow() != Parameters.Upper() - || Matrix.LowerCol() != 1 || Matrix.UpperCol() != BandWidth) + || Matrix.LowerCol() != 1 || Matrix.UpperCol() != aBandWidth) { - ReturnCode = 1; - goto FINISH; + return 1; } - for (ii = Parameters.Lower(); ii <= Parameters.Upper(); ii++) + math_Matrix aBSplineBasis(1, aMaxOrder, 1, aMaxOrder); + + for (int i = Parameters.Lower(); i <= Parameters.Upper(); i++) { - ErrorCode = BSplCLib::EvalBsplineBasis(ContactOrderArray(ii), - Order, - FlatKnots, - Parameters(ii), - - FirstNonZeroBsplineIndex, - BSplineBasis); - if (ErrorCode != 0) + int aFirstNonZeroIndex = 0; + const int anErrorCode = BSplCLib::EvalBsplineBasis(ContactOrderArray(i), + anOrder, + FlatKnots, + Parameters(i), + aFirstNonZeroIndex, + aBSplineBasis); + if (anErrorCode != 0) { - ReturnCode = 2; - goto FINISH; + return 2; } - Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii; - for (jj = 1; jj < Index; jj++) + int anIndex = LowerBandWidth + 1 + aFirstNonZeroIndex - i; + for (int j = 1; j < anIndex; j++) { - Matrix.Value(ii, jj) = 0.0e0; + Matrix.Value(i, j) = 0.0; } - - for (jj = 1; jj <= Order; jj++) + for (int j = 1; j <= anOrder; j++) { - Matrix.Value(ii, Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj); - Index += 1; + Matrix.Value(i, anIndex) = aBSplineBasis(ContactOrderArray(i) + 1, j); + anIndex += 1; } - - for (jj = Index; jj <= BandWidth; jj++) + for (int j = anIndex; j <= aBandWidth; j++) { - Matrix.Value(ii, jj) = 0.0e0; + Matrix.Value(i, j) = 0.0; } } -FINISH:; - return (ReturnCode); + + return 0; } -//======================================================================= -// function : Makes LU decompositiomn without Pivoting -// purpose : Builds the Bspline Matrix -//======================================================================= +//================================================================================================= Standard_Integer BSplCLib::FactorBandedMatrix(math_Matrix& Matrix, const Standard_Integer UpperBandWidth, const Standard_Integer LowerBandWidth, Standard_Integer& PivotIndexProblem) { - Standard_Integer ii, jj, kk, Index, MinIndex, MaxIndex, - ReturnCode = 0, BandWidth = UpperBandWidth + LowerBandWidth + 1; + const int aBandWidth = UpperBandWidth + LowerBandWidth + 1; + PivotIndexProblem = 0; - Standard_Real Inverse; - PivotIndexProblem = 0; - - for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow(); ii++) + for (int i = Matrix.LowerRow() + 1; i <= Matrix.UpperRow(); i++) { - MinIndex = (LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1); + const int aMinIndex = (LowerBandWidth - i + 2 >= 1 ? LowerBandWidth - i + 2 : 1); - for (jj = MinIndex; jj <= LowerBandWidth; jj++) + for (int j = aMinIndex; j <= LowerBandWidth; j++) { - Index = ii - LowerBandWidth + jj - 1; - Inverse = Matrix(Index, LowerBandWidth + 1); - if (std::abs(Inverse) > RealSmall()) - { - Inverse = -1.0e0 / Inverse; - } - else + const int anIndex = i - LowerBandWidth + j - 1; + const double aPivot = Matrix(anIndex, LowerBandWidth + 1); + if (std::abs(aPivot) <= RealSmall()) { - ReturnCode = 1; - PivotIndexProblem = Index; - goto FINISH; + PivotIndexProblem = anIndex; + return 1; } - Matrix(ii, jj) = Matrix(ii, jj) * Inverse; - MaxIndex = BandWidth + Index - ii; - for (kk = jj + 1; kk <= MaxIndex; kk++) + const double anInverse = -1.0 / aPivot; + Matrix(i, j) = Matrix(i, j) * anInverse; + const int aMaxIndex = aBandWidth + anIndex - i; + + for (int k = j + 1; k <= aMaxIndex; k++) { - Matrix(ii, kk) += Matrix(ii, jj) * Matrix(Index, kk + ii - Index); + Matrix(i, k) += Matrix(i, j) * Matrix(anIndex, k + i - anIndex); } } } -FINISH: - return (ReturnCode); + + return 0; } -//======================================================================= -// function : Build BSpline Matrix -// purpose : Builds the Bspline Matrix -//======================================================================= +//================================================================================================= Standard_Integer BSplCLib::EvalBsplineBasis(const Standard_Integer DerivativeRequest, const Standard_Integer Order, @@ -437,109 +411,100 @@ Standard_Integer BSplCLib::EvalBsplineBasis(const Standard_Integer Derivati // B (t) B (t) B (t) // i i+1 i+k-1 // - Standard_Integer ReturnCode, ii, pp, qq, ss, NumPoles, LocalRequest; - // ,Index ; - - Standard_Real NewParameter, Inverse, Factor, LocalInverse, Saved; - // , *FlatKnotsArray ; - - ReturnCode = 0; FirstNonZeroBsplineIndex = 0; - LocalRequest = DerivativeRequest; + int aLocalRequest = DerivativeRequest; if (DerivativeRequest >= Order) { - LocalRequest = Order - 1; + aLocalRequest = Order - 1; } if (BsplineBasis.LowerCol() != 1 || BsplineBasis.UpperCol() < Order - || BsplineBasis.LowerRow() != 1 || BsplineBasis.UpperRow() <= LocalRequest) + || BsplineBasis.LowerRow() != 1 || BsplineBasis.UpperRow() <= aLocalRequest) { - ReturnCode = 1; - goto FINISH; + return 1; } - NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order; + + const int aNumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order; + int ii = 0; + double aNewParam = 0.0; BSplCLib::LocateParameter(Order - 1, FlatKnots, Parameter, isPeriodic, Order, - NumPoles + 1, + aNumPoles + 1, ii, - NewParameter); + aNewParam); FirstNonZeroBsplineIndex = ii - Order + 1; - BsplineBasis(1, 1) = 1.0e0; - LocalRequest = DerivativeRequest; + BsplineBasis(1, 1) = 1.0; + aLocalRequest = DerivativeRequest; if (DerivativeRequest >= Order) { - LocalRequest = Order - 1; + aLocalRequest = Order - 1; } - for (qq = 2; qq <= Order - LocalRequest; qq++) + for (int qq = 2; qq <= Order - aLocalRequest; qq++) { - BsplineBasis(1, qq) = 0.0e0; + BsplineBasis(1, qq) = 0.0; - for (pp = 1; pp <= qq - 1; pp++) + for (int pp = 1; pp <= qq - 1; pp++) { - // - // this should be always invertible if ii is correctly computed - // - const Standard_Real aScale = (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)); + const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1); if (std::abs(aScale) < gp::Resolution()) { return 2; } - Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) / aScale; - Saved = Factor * BsplineBasis(1, pp); - BsplineBasis(1, pp) *= (1.0e0 - Factor); + const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) / aScale; + const double aSaved = aFactor * BsplineBasis(1, pp); + BsplineBasis(1, pp) *= (1.0 - aFactor); BsplineBasis(1, pp) += BsplineBasis(1, qq); - BsplineBasis(1, qq) = Saved; + BsplineBasis(1, qq) = aSaved; } } - for (qq = Order - LocalRequest + 1; qq <= Order; qq++) + for (int qq = Order - aLocalRequest + 1; qq <= Order; qq++) { - - for (pp = 1; pp <= qq - 1; pp++) + for (int pp = 1; pp <= qq - 1; pp++) { BsplineBasis(Order - qq + 2, pp) = BsplineBasis(1, pp); } - BsplineBasis(1, qq) = 0.0e0; + BsplineBasis(1, qq) = 0.0; - for (ss = Order - LocalRequest + 1; ss <= qq; ss++) + for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++) { - BsplineBasis(Order - ss + 2, qq) = 0.0e0; + BsplineBasis(Order - ss + 2, qq) = 0.0; } - for (pp = 1; pp <= qq - 1; pp++) + for (int pp = 1; pp <= qq - 1; pp++) { - const Standard_Real aScale = (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)); + const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1); if (std::abs(aScale) < gp::Resolution()) { return 2; } - Inverse = 1.0e0 / aScale; - Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse; - Saved = Factor * BsplineBasis(1, pp); - BsplineBasis(1, pp) *= (1.0e0 - Factor); + const double anInverse = 1.0 / aScale; + const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) * anInverse; + double aSaved = aFactor * BsplineBasis(1, pp); + BsplineBasis(1, pp) *= (1.0 - aFactor); BsplineBasis(1, pp) += BsplineBasis(1, qq); - BsplineBasis(1, qq) = Saved; - LocalInverse = (Standard_Real)(qq - 1) * Inverse; + BsplineBasis(1, qq) = aSaved; + const double aLocalInverse = static_cast(qq - 1) * anInverse; - for (ss = Order - LocalRequest + 1; ss <= qq; ss++) + for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++) { - Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp); - BsplineBasis(Order - ss + 2, pp) *= -LocalInverse; + aSaved = aLocalInverse * BsplineBasis(Order - ss + 2, pp); + BsplineBasis(Order - ss + 2, pp) *= -aLocalInverse; BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2, qq); - BsplineBasis(Order - ss + 2, qq) = Saved; + BsplineBasis(Order - ss + 2, qq) = aSaved; } } } -FINISH: - return (ReturnCode); + + return 0; } //================================================================================================= @@ -816,8 +781,9 @@ void BSplCLib::MovePointAndTangent(const Standard_Real U, } } a_matrix.Invert(); - TColStd_Array1OfReal the_a_vector(0, ArrayDimension - 1); - TColStd_Array1OfReal the_b_vector(0, ArrayDimension - 1); + // Use math_Vector for stack allocation (ArrayDimension is typically 2-4) + math_Vector the_a_vector(0, ArrayDimension - 1); + math_Vector the_b_vector(0, ArrayDimension - 1); for (ii = 0; ii < ArrayDimension; ii++) { @@ -863,74 +829,69 @@ void BSplCLib::FunctionMultiply(const BSplCLib_EvaluatorFunction& FunctionPtr, Standard_Real& NewPoles, Standard_Integer& theStatus) { - Standard_Integer ii, jj, index; - Standard_Integer extrap_mode[2], error_code, num_new_poles, derivative_request = 0; - Standard_Boolean periodic_flag = Standard_False; - Standard_Real result, start_end[2], *array_of_poles, *array_of_new_poles; - - array_of_poles = (Standard_Real*)&NewPoles; - extrap_mode[0] = extrap_mode[1] = BSplineDegree; - num_new_poles = FlatKnots.Length() - NewDegree - 1; - start_end[0] = FlatKnots(NewDegree + 1); - start_end[1] = FlatKnots(num_new_poles + 1); - TColStd_Array1OfReal parameters(1, num_new_poles); - TColStd_Array1OfInteger contact_order_array(1, num_new_poles); - TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension); - - array_of_new_poles = (Standard_Real*)&new_poles_array(1); - BuildSchoenbergPoints(NewDegree, FlatKnots, parameters); - // - // on recadre sur les bornes - // - if (parameters(1) < start_end[0]) + double* anArrayOfPoles = &NewPoles; + const int aNumNewPoles = FlatKnots.Length() - NewDegree - 1; + double aStartEnd[2] = {FlatKnots(NewDegree + 1), FlatKnots(aNumNewPoles + 1)}; + + TColStd_Array1OfReal aParameters(1, aNumNewPoles); + TColStd_Array1OfInteger aContactOrderArray(1, aNumNewPoles); + TColStd_Array1OfReal aNewPolesArray(1, aNumNewPoles * PolesDimension); + + double* anArrayOfNewPoles = &aNewPolesArray(1); + BuildSchoenbergPoints(NewDegree, FlatKnots, aParameters); + + if (aParameters(1) < aStartEnd[0]) { - parameters(1) = start_end[0]; + aParameters(1) = aStartEnd[0]; } - if (parameters(num_new_poles) > start_end[1]) + if (aParameters(aNumNewPoles) > aStartEnd[1]) { - parameters(num_new_poles) = start_end[1]; + aParameters(aNumNewPoles) = aStartEnd[1]; } - index = 0; - for (ii = 1; ii <= num_new_poles; ii++) + int anExtrapMode = BSplineDegree; + int anIndex = 0; + for (int i = 1; i <= aNumNewPoles; i++) { - contact_order_array(ii) = 0; - FunctionPtr.Evaluate(contact_order_array(ii), start_end, parameters(ii), result, error_code); - if (error_code) + aContactOrderArray(i) = 0; + double aResult = 0.0; + int anErrorCode = 0; + FunctionPtr.Evaluate(aContactOrderArray(i), aStartEnd, aParameters(i), aResult, anErrorCode); + if (anErrorCode) { theStatus = 1; - goto FINISH; + return; } - Eval(parameters(ii), - periodic_flag, - derivative_request, - extrap_mode[0], + Eval(aParameters(i), + Standard_False, + 0, + anExtrapMode, BSplineDegree, BSplineFlatKnots, PolesDimension, Poles, - array_of_new_poles[index]); + anArrayOfNewPoles[anIndex]); - for (jj = 0; jj < PolesDimension; jj++) + for (int j = 0; j < PolesDimension; j++) { - array_of_new_poles[index] *= result; - index += 1; + anArrayOfNewPoles[anIndex] *= aResult; + anIndex += 1; } } + Interpolate(NewDegree, FlatKnots, - parameters, - contact_order_array, + aParameters, + aContactOrderArray, PolesDimension, - array_of_new_poles[0], + anArrayOfNewPoles[0], theStatus); - for (ii = 0; ii < num_new_poles * PolesDimension; ii++) + for (int i = 0; i < aNumNewPoles * PolesDimension; i++) { - array_of_poles[ii] = array_of_new_poles[ii]; + anArrayOfPoles[i] = anArrayOfNewPoles[i]; } -FINISH:; } //================================================================================================= @@ -945,60 +906,55 @@ void BSplCLib::FunctionReparameterise(const BSplCLib_EvaluatorFunction& Function Standard_Real& NewPoles, Standard_Integer& theStatus) { - Standard_Integer ii, - // jj, - index; - Standard_Integer extrap_mode[2], error_code, num_new_poles, derivative_request = 0; - Standard_Boolean periodic_flag = Standard_False; - Standard_Real result, start_end[2], *array_of_poles, *array_of_new_poles; - - array_of_poles = (Standard_Real*)&NewPoles; - extrap_mode[0] = extrap_mode[1] = BSplineDegree; - num_new_poles = FlatKnots.Length() - NewDegree - 1; - start_end[0] = FlatKnots(NewDegree + 1); - start_end[1] = FlatKnots(num_new_poles + 1); - TColStd_Array1OfReal parameters(1, num_new_poles); - TColStd_Array1OfInteger contact_order_array(1, num_new_poles); - TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension); - - array_of_new_poles = (Standard_Real*)&new_poles_array(1); - BuildSchoenbergPoints(NewDegree, FlatKnots, parameters); - index = 0; - - for (ii = 1; ii <= num_new_poles; ii++) + double* anArrayOfPoles = &NewPoles; + const int aNumNewPoles = FlatKnots.Length() - NewDegree - 1; + double aStartEnd[2] = {FlatKnots(NewDegree + 1), FlatKnots(aNumNewPoles + 1)}; + + TColStd_Array1OfReal aParameters(1, aNumNewPoles); + TColStd_Array1OfInteger aContactOrderArray(1, aNumNewPoles); + TColStd_Array1OfReal aNewPolesArray(1, aNumNewPoles * PolesDimension); + + double* anArrayOfNewPoles = &aNewPolesArray(1); + BuildSchoenbergPoints(NewDegree, FlatKnots, aParameters); + + int anExtrapMode = BSplineDegree; + int anIndex = 0; + for (int i = 1; i <= aNumNewPoles; i++) { - contact_order_array(ii) = 0; - FunctionPtr.Evaluate(contact_order_array(ii), start_end, parameters(ii), result, error_code); - if (error_code) + aContactOrderArray(i) = 0; + double aResult = 0.0; + int anErrorCode = 0; + FunctionPtr.Evaluate(aContactOrderArray(i), aStartEnd, aParameters(i), aResult, anErrorCode); + if (anErrorCode) { theStatus = 1; - goto FINISH; + return; } - Eval(result, - periodic_flag, - derivative_request, - extrap_mode[0], + Eval(aResult, + Standard_False, + 0, + anExtrapMode, BSplineDegree, BSplineFlatKnots, PolesDimension, Poles, - array_of_new_poles[index]); - index += PolesDimension; + anArrayOfNewPoles[anIndex]); + anIndex += PolesDimension; } + Interpolate(NewDegree, FlatKnots, - parameters, - contact_order_array, + aParameters, + aContactOrderArray, PolesDimension, - array_of_new_poles[0], + anArrayOfNewPoles[0], theStatus); - for (ii = 0; ii < num_new_poles * PolesDimension; ii++) + for (int i = 0; i < aNumNewPoles * PolesDimension; i++) { - array_of_poles[ii] = array_of_new_poles[ii]; + anArrayOfPoles[i] = anArrayOfNewPoles[i]; } -FINISH:; } //================================================================================================= diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_3.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_3.cxx index 2fbcd690a5..c97821b4b9 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_3.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_3.cxx @@ -35,7 +35,7 @@ void BSplCLib::Reverse(TColgp_Array1OfPnt& Poles, const Standard_Integer L) { - BSplCLib_Reverse(Poles, L); + BSplCLib_Reverse(Poles, L); } //================================================================================================== diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_BzSyntaxes.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_BzSyntaxes.cxx index a50f166a97..22e3a63c81 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_BzSyntaxes.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_BzSyntaxes.cxx @@ -15,37 +15,18 @@ // commercial license or contractual agreement. #include -#include -#include -#define No_Standard_RangeError -#define No_Standard_OutOfRange +#include +#include +#include +#include +#include +#include -//======================================================================= -// struct : BSplCLib_BezierArrays -// purpose: Auxiliary structure providing standard definitions of bspline -// knots for bezier (using stack allocation) -//======================================================================= +#include "BSplCLib_CurveComputation.pxx" -class BSplCLib_BezierArrays -{ -public: - BSplCLib_BezierArrays(Standard_Integer Degree) - : aKnots{0., 1.}, - aMults{Degree + 1, Degree + 1}, - knots(aKnots[0], 1, 2), - mults(aMults[0], 1, 2) - { - } - -private: - Standard_Real aKnots[2]; - Standard_Integer aMults[2]; - -public: - TColStd_Array1OfReal knots; - TColStd_Array1OfInteger mults; -}; +#define No_Standard_RangeError +#define No_Standard_OutOfRange //================================================================================================= @@ -55,19 +36,11 @@ void BSplCLib::IncreaseDegree(const Standard_Integer NewDegree, TColgp_Array1OfPnt& NewPoles, TColStd_Array1OfReal* NewWeights) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::IncreaseDegree(deg, - NewDegree, - 0, - Poles, - Weights, - bzarr.knots, - bzarr.mults, - NewPoles, - NewWeights, - bzarr.knots, - bzarr.mults); + BSplCLib_IncreaseDegree_Bezier(NewDegree, + Poles, + Weights, + NewPoles, + NewWeights); } //================================================================================================= @@ -78,19 +51,11 @@ void BSplCLib::IncreaseDegree(const Standard_Integer NewDegree, TColgp_Array1OfPnt2d& NewPoles, TColStd_Array1OfReal* NewWeights) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::IncreaseDegree(deg, - NewDegree, - 0, - Poles, - Weights, - bzarr.knots, - bzarr.mults, - NewPoles, - NewWeights, - bzarr.knots, - bzarr.mults); + BSplCLib_IncreaseDegree_Bezier(NewDegree, + Poles, + Weights, + NewPoles, + NewWeights); } //================================================================================================= @@ -100,9 +65,10 @@ void BSplCLib::PolesCoefficients(const TColgp_Array1OfPnt& Poles, TColgp_Array1OfPnt& CachePoles, TColStd_Array1OfReal* CacheWeights) { - Standard_Integer deg = Poles.Length() - 1; - TColStd_Array1OfReal bidflatknots(FlatBezierKnots(deg), 1, 2 * (deg + 1)); - BSplCLib::BuildCache(0., 1., 0, deg, bidflatknots, Poles, Weights, CachePoles, CacheWeights); + BSplCLib_PolesCoefficients_Bezier(Poles, + Weights, + CachePoles, + CacheWeights); } //================================================================================================= @@ -112,9 +78,10 @@ void BSplCLib::PolesCoefficients(const TColgp_Array1OfPnt2d& Poles, TColgp_Array1OfPnt2d& CachePoles, TColStd_Array1OfReal* CacheWeights) { - Standard_Integer deg = Poles.Length() - 1; - TColStd_Array1OfReal bidflatknots(FlatBezierKnots(deg), 1, 2 * (deg + 1)); - BSplCLib::BuildCache(0., 1., 0, deg, bidflatknots, Poles, Weights, CachePoles, CacheWeights); + BSplCLib_PolesCoefficients_Bezier(Poles, + Weights, + CachePoles, + CacheWeights); } //================================================================================================= @@ -124,9 +91,9 @@ void BSplCLib::D0(const Standard_Real U, const TColStd_Array1OfReal* Weights, gp_Pnt& P) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D0(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D0(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P); } //================================================================================================= @@ -136,9 +103,9 @@ void BSplCLib::D0(const Standard_Real U, const TColStd_Array1OfReal* Weights, gp_Pnt2d& P) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D0(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D0(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P); } //================================================================================================= @@ -149,9 +116,9 @@ void BSplCLib::D1(const Standard_Real U, gp_Pnt& P, gp_Vec& V) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D1(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D1(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P, V); } //================================================================================================= @@ -162,9 +129,9 @@ void BSplCLib::D1(const Standard_Real U, gp_Pnt2d& P, gp_Vec2d& V) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D1(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D1(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P, V); } //================================================================================================= @@ -176,9 +143,9 @@ void BSplCLib::D2(const Standard_Real U, gp_Vec& V1, gp_Vec& V2) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D2(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V1, V2); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D2(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P, V1, V2); } //================================================================================================= @@ -190,9 +157,9 @@ void BSplCLib::D2(const Standard_Real U, gp_Vec2d& V1, gp_Vec2d& V2) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D2(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V1, V2); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D2(U, 1, aDegree, 0, Poles, Weights, aBezierKnots.Knot, &aBezierKnots.Mult, P, V1, V2); } //================================================================================================= @@ -205,9 +172,20 @@ void BSplCLib::D3(const Standard_Real U, gp_Vec& V2, gp_Vec& V3) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D3(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V1, V2, V3); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D3(U, + 1, + aDegree, + 0, + Poles, + Weights, + aBezierKnots.Knot, + &aBezierKnots.Mult, + P, + V1, + V2, + V3); } //================================================================================================= @@ -220,7 +198,18 @@ void BSplCLib::D3(const Standard_Real U, gp_Vec2d& V2, gp_Vec2d& V3) { - Standard_Integer deg = Poles.Length() - 1; - BSplCLib_BezierArrays bzarr(deg); - BSplCLib::D3(U, 1, deg, 0, Poles, Weights, bzarr.knots, &bzarr.mults, P, V1, V2, V3); + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib::D3(U, + 1, + aDegree, + 0, + Poles, + Weights, + aBezierKnots.Knot, + &aBezierKnots.Mult, + P, + V1, + V2, + V3); } diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.cxx index b44e8231da..c267e9a722 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.cxx @@ -18,17 +18,9 @@ #include #include -#include IMPLEMENT_STANDARD_RTTIEXT(BSplCLib_Cache, Standard_Transient) -//! Converts handle of array of Standard_Real into the pointer to Standard_Real -static Standard_Real* ConvertArray(const Handle(TColStd_HArray2OfReal)& theHArray) -{ - const TColStd_Array2OfReal& anArray = theHArray->Array2(); - return (Standard_Real*)&(anArray(anArray.LowerRow(), anArray.LowerCol())); -} - BSplCLib_Cache::BSplCLib_Cache( const Standard_Integer& theDegree, const Standard_Boolean& thePeriodic, @@ -36,10 +28,9 @@ BSplCLib_Cache::BSplCLib_Cache( const TColgp_Array1OfPnt2d& /* only used to distinguish from 3d variant */, const TColStd_Array1OfReal* theWeights) : myIsRational(theWeights != NULL), - myParams(theDegree, thePeriodic, theFlatKnots) + myParams(theDegree, thePeriodic, theFlatKnots), + myRowLength(myIsRational ? 3 : 2) { - Standard_Integer aPWColNumber = (myIsRational ? 3 : 2); - myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber); } BSplCLib_Cache::BSplCLib_Cache( @@ -49,10 +40,9 @@ BSplCLib_Cache::BSplCLib_Cache( const TColgp_Array1OfPnt& /* only used to distinguish from 2d variant */, const TColStd_Array1OfReal* theWeights) : myIsRational(theWeights != NULL), - myParams(theDegree, thePeriodic, theFlatKnots) + myParams(theDegree, thePeriodic, theFlatKnots), + myRowLength(myIsRational ? 4 : 3) { - Standard_Integer aPWColNumber = (myIsRational ? 4 : 3); - myPolesWeights = new TColStd_HArray2OfReal(1, theDegree + 1, 1, aPWColNumber); } Standard_Boolean BSplCLib_Cache::IsCacheValid(Standard_Real theParameter) const @@ -69,6 +59,13 @@ void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter, Standard_Real aNewParam = myParams.PeriodicNormalization(theParameter); myParams.LocateParameter(aNewParam, theFlatKnots); + // Create array wrapper referencing the stack buffer + TColStd_Array2OfReal aPolesWeights(myPolesWeightsBuffer[0], + 1, + myParams.Degree + 1, + 1, + myRowLength); + // Calculate new cache data BSplCLib::BuildCache(myParams.SpanStart, myParams.SpanLength, @@ -78,7 +75,7 @@ void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter, theFlatKnots, thePoles2d, theWeights, - myPolesWeights->ChangeArray2()); + aPolesWeights); } void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter, @@ -90,6 +87,13 @@ void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter, Standard_Real aNewParam = myParams.PeriodicNormalization(theParameter); myParams.LocateParameter(aNewParam, theFlatKnots); + // Create array wrapper referencing the stack buffer + TColStd_Array2OfReal aPolesWeights(myPolesWeightsBuffer[0], + 1, + myParams.Degree + 1, + 1, + myRowLength); + // Calculate new cache data BSplCLib::BuildCache(myParams.SpanStart, myParams.SpanLength, @@ -99,7 +103,7 @@ void BSplCLib_Cache::BuildCache(const Standard_Real& theParameter, theFlatKnots, thePoles, theWeights, - myPolesWeights->ChangeArray2()); + aPolesWeights); } void BSplCLib_Cache::CalculateDerivative(const Standard_Real& theParameter, @@ -109,8 +113,8 @@ void BSplCLib_Cache::CalculateDerivative(const Standard_Real& theParameter, Standard_Real aNewParameter = myParams.PeriodicNormalization(theParameter); aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength; - Standard_Real* aPolesArray = ConvertArray(myPolesWeights); - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Real* aPolesArray = const_cast(myPolesWeightsBuffer); + const Standard_Integer aDimension = myRowLength; // Temporary container. The maximal size of this container is defined by: // 1) maximal derivative for cache evaluation, which is 3, plus one row for function values, @@ -161,9 +165,9 @@ void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt2d& thePoint) c Standard_Real aNewParameter = myParams.PeriodicNormalization(theParameter); aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength; - Standard_Real* aPolesArray = ConvertArray(myPolesWeights); - Standard_Real aPoint[4]; - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Real* aPolesArray = const_cast(myPolesWeightsBuffer); + Standard_Real aPoint[4]; + const Standard_Integer aDimension = myRowLength; PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree, @@ -182,9 +186,9 @@ void BSplCLib_Cache::D0(const Standard_Real& theParameter, gp_Pnt& thePoint) con Standard_Real aNewParameter = myParams.PeriodicNormalization(theParameter); aNewParameter = (aNewParameter - myParams.SpanStart) / myParams.SpanLength; - Standard_Real* aPolesArray = ConvertArray(myPolesWeights); - Standard_Real aPoint[4]; - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Real* aPolesArray = const_cast(myPolesWeightsBuffer); + Standard_Real aPoint[4]; + const Standard_Integer aDimension = myRowLength; PLib::NoDerivativeEvalPolynomial(aNewParameter, myParams.Degree, @@ -202,7 +206,7 @@ void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt2d& thePoint, gp_Vec2d& theTangent) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates) this->CalculateDerivative(theParameter, 1, aPntDeriv[0]); @@ -217,7 +221,7 @@ void BSplCLib_Cache::D1(const Standard_Real& theParameter, gp_Pnt& thePoint, gp_Vec& theTangent) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[8]; // result storage (point and derivative coordinates) this->CalculateDerivative(theParameter, 1, aPntDeriv[0]); @@ -233,7 +237,7 @@ void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Vec2d& theTangent, gp_Vec2d& theCurvature) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates) this->CalculateDerivative(theParameter, 2, aPntDeriv[0]); @@ -250,7 +254,7 @@ void BSplCLib_Cache::D2(const Standard_Real& theParameter, gp_Vec& theTangent, gp_Vec& theCurvature) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[12]; // result storage (point and derivatives coordinates) this->CalculateDerivative(theParameter, 2, aPntDeriv[0]); @@ -270,7 +274,7 @@ void BSplCLib_Cache::D3(const Standard_Real& theParameter, gp_Vec2d& theCurvature, gp_Vec2d& theTorsion) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates) this->CalculateDerivative(theParameter, 3, aPntDeriv[0]); @@ -291,7 +295,7 @@ void BSplCLib_Cache::D3(const Standard_Real& theParameter, gp_Vec& theCurvature, gp_Vec& theTorsion) const { - Standard_Integer aDimension = myPolesWeights->RowLength(); // number of columns + Standard_Integer aDimension = myRowLength; Standard_Real aPntDeriv[16]; // result storage (point and derivatives coordinates) this->CalculateDerivative(theParameter, 3, aPntDeriv[0]); diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.hxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.hxx index d2d4f64aff..0471090f5d 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.hxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_Cache.hxx @@ -15,7 +15,6 @@ #define _BSplCLib_Cache_Headerfile #include -#include //! \brief A cache class for Bezier and B-spline curves. //! @@ -142,13 +141,16 @@ protected: private: // clang-format off - Standard_Boolean myIsRational; //!< identifies the rationality of Bezier/B-spline curve - BSplCLib_CacheParams myParams; //!< cache parameters - Handle(TColStd_HArray2OfReal) myPolesWeights; //!< array of poles and weights of calculated cache - // the array has following structure: - // x1 y1 [z1] [w1] - // x2 y2 [z2] [w2] etc - // for 2D-curves there is no z conponent, for non-rational curves there is no weight + Standard_Boolean myIsRational; //!< identifies the rationality of Bezier/B-spline curve + BSplCLib_CacheParams myParams; //!< cache parameters + Standard_Integer myRowLength; //!< number of columns in the cache array + Standard_Real myPolesWeightsBuffer[128]; //!< stack buffer for poles/weights cache + //!< 128 elements covers degree 31 for 3D rational (32*4) + //!< Array structure: + //!< x1 y1 [z1] [w1] + //!< x2 y2 [z2] [w2] etc + //!< For 2D-curves: no z component + //!< For non-rational curves: no weight // clang-format on }; diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CacheParams.hxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CacheParams.hxx index 5b74077323..d888bd3ab5 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CacheParams.hxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CacheParams.hxx @@ -15,7 +15,6 @@ #define _BSplCLib_CacheParams_Headerfile #include -#include #include #include @@ -25,53 +24,50 @@ //! and data of the current span for its caching struct BSplCLib_CacheParams { - const Standard_Integer Degree; ///< degree of Bezier/B-spline - const Standard_Boolean IsPeriodic; ///< true of the B-spline is periodic - const Standard_Real FirstParameter; ///< first valid parameter - const Standard_Real LastParameter; ///< last valid parameter + const int Degree; ///< degree of Bezier/B-spline + const bool IsPeriodic; ///< true of the B-spline is periodic + const double FirstParameter; ///< first valid parameter + const double LastParameter; ///< last valid parameter - const Standard_Integer SpanIndexMin; ///< minimal index of span - const Standard_Integer SpanIndexMax; ///< maximal index of span + const int SpanIndexMin; ///< minimal index of span + const int SpanIndexMax; ///< maximal index of span - Standard_Real SpanStart; ///< parameter for the frst point of the span - Standard_Real SpanLength; ///< length of the span - Standard_Integer SpanIndex; ///< index of the span + double SpanStart; ///< parameter for the frst point of the span + double SpanLength; ///< length of the span + int SpanIndex; ///< index of the span //! Constructor, prepares data structures for caching. //! \param theDegree degree of the B-spline (or Bezier) //! \param thePeriodic identify whether the B-spline is periodic //! \param theFlatKnots knots of Bezier / B-spline parameterization - BSplCLib_CacheParams(Standard_Integer theDegree, - Standard_Boolean thePeriodic, - const TColStd_Array1OfReal& theFlatKnots) + BSplCLib_CacheParams(int theDegree, bool thePeriodic, const TColStd_Array1OfReal& theFlatKnots) : Degree(theDegree), IsPeriodic(thePeriodic), FirstParameter(theFlatKnots.Value(theFlatKnots.Lower() + theDegree)), LastParameter(theFlatKnots.Value(theFlatKnots.Upper() - theDegree)), SpanIndexMin(theFlatKnots.Lower() + theDegree), SpanIndexMax(theFlatKnots.Upper() - theDegree - 1), - SpanStart(0.), - SpanLength(0.), + SpanStart(0.0), + SpanLength(0.0), SpanIndex(0) { } //! Normalizes the parameter for periodic B-splines //! \param theParameter the value to be normalized into the knots array - Standard_Real PeriodicNormalization(Standard_Real theParameter) const + double PeriodicNormalization(double theParameter) const noexcept { if (IsPeriodic) { + const double aPeriod = LastParameter - FirstParameter; if (theParameter < FirstParameter) { - Standard_Real aPeriod = LastParameter - FirstParameter; - Standard_Real aScale = std::trunc((FirstParameter - theParameter) / aPeriod); + const double aScale = std::trunc((FirstParameter - theParameter) / aPeriod); return theParameter + aPeriod * (aScale + 1.0); } if (theParameter > LastParameter) { - Standard_Real aPeriod = LastParameter - FirstParameter; - Standard_Real aScale = std::trunc((theParameter - LastParameter) / aPeriod); + const double aScale = std::trunc((theParameter - LastParameter) / aPeriod); return theParameter - aPeriod * (aScale + 1.0); } } @@ -80,10 +76,10 @@ struct BSplCLib_CacheParams //! Verifies validity of the cache using flat parameter of the point //! \param theParameter parameter of the point placed in the span - Standard_Boolean IsCacheValid(Standard_Real theParameter) const + bool IsCacheValid(double theParameter) const noexcept { - Standard_Real aNewParam = PeriodicNormalization(theParameter); - Standard_Real aDelta = aNewParam - SpanStart; + const double aNewParam = PeriodicNormalization(theParameter); + const double aDelta = aNewParam - SpanStart; if (!((aDelta >= 0.0 || SpanIndex == SpanIndexMin) && (aDelta < SpanLength || SpanIndex == SpanIndexMax))) { @@ -106,7 +102,7 @@ struct BSplCLib_CacheParams //! Computes span for the specified parameter //! \param theParameter parameter of the point placed in the span //! \param theFlatKnots knots of Bezier / B-spline parameterization - void LocateParameter(Standard_Real& theParameter, const TColStd_Array1OfReal& theFlatKnots) + void LocateParameter(double& theParameter, const TColStd_Array1OfReal& theFlatKnots) { SpanIndex = 0; BSplCLib::LocateParameter(Degree, @@ -120,10 +116,9 @@ struct BSplCLib_CacheParams SpanLength = theFlatKnots.Value(SpanIndex + 1) - SpanStart; } -private: // copying is prohibited - BSplCLib_CacheParams(const BSplCLib_CacheParams&); - void operator=(const BSplCLib_CacheParams&); + BSplCLib_CacheParams(const BSplCLib_CacheParams&) = delete; + BSplCLib_CacheParams& operator=(const BSplCLib_CacheParams&) = delete; }; #endif diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.gxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.gxx deleted file mode 100644 index 33b71588c2..0000000000 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.gxx +++ /dev/null @@ -1,1483 +0,0 @@ -// 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]; -}; - -//================================================================================================= - -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 -// - -//================================================================================================= - -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); -} - -//================================================================================================= - -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)); -} - -//================================================================================================= - -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))); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -// -// 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); - } -} - -//================================================================================================= - -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, ); -} - -//================================================================================================= - -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, ); -} - -//================================================================================================= - -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, ); -} - -//================================================================================================= - -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, ); -} - -//================================================================================================= - -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, ); - } - } -} - -//================================================================================================= - -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 weights 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 = std::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 = std::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); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -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 symmetry - kk2 = kk1; - i = kk1 - FirstNonZeroBsplineIndex + 2; - if ((kk1 + 1) <= LastIndex) - { - if (std::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 -//======================================================================= - -//================================================================================================= - -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_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); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -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); -} - -//================================================================================================= - -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); -} diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.pxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.pxx index d37217e7c9..0534d0afc9 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.pxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_CurveComputation.pxx @@ -26,6 +26,9 @@ #include #include +#include +#include + // Template traits for 2D/3D point and vector operations template struct BSplCLib_CurveTraits @@ -203,6 +206,44 @@ struct BSplCLib_CurveTraits } }; +//! Stack-allocated arrays for knot/multiplicity pairs. +//! Size=1: single knot constructor (knotValue, multiplicity) +//! Size=2: Bezier constructor (degree) with knots {0,1} and mults {degree+1, degree+1} +//! @tparam Size number of knots in the array (1 or 2) +template +struct BSplCLib_KnotArrays +{ + static_assert(Size == 1 || Size == 2, "BSplCLib_KnotArrays: Size must be 1 or 2"); + + //! Constructor for single knot (Size=1) + template ::type> + BSplCLib_KnotArrays(double theKnotValue, int theMultiplicity) + : myKnotBuffer{theKnotValue}, + myMultBuffer{theMultiplicity}, + Knot(myKnotBuffer[0], 1, Size), + Mult(myMultBuffer[0], 1, Size) + { + } + + //! Constructor for Bezier knots (Size=2) + template ::type> + BSplCLib_KnotArrays(int theDegree) + : myKnotBuffer{0.0, 1.0}, + myMultBuffer{theDegree + 1, theDegree + 1}, + Knot(myKnotBuffer[0], 1, Size), + Mult(myMultBuffer[0], 1, Size) + { + } + +private: + double myKnotBuffer[Size]; + int myMultBuffer[Size]; + +public: + TColStd_Array1OfReal Knot; + TColStd_Array1OfInteger Mult; +}; + // Auxiliary structure providing buffers for poles and knots used in evaluation of bspline // (allocated in the stack) template @@ -220,25 +261,19 @@ struct BSplCLib_DataContainer_T Standard_Real ders[Dimension * 4]; }; -// Reverses the order of poles in the array -// @param[in,out] Poles - array of poles to be reversed -// @param[in] L - length parameter for reversal -template -void BSplCLib_Reverse(Array1OfPoints& Poles, const Standard_Integer L) +// Reverses the order of elements in the array using std::reverse. +// Works with any NCollection_Array1-like container. +// @param[in,out] theArray - array to be reversed +// @param[in] theL - length parameter for reversal +template +void BSplCLib_Reverse(Array& theArray, const Standard_Integer theL) { - 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); + const int aLower = theArray.Lower(); + const int aUpper = theArray.Upper(); + const int aL = aLower + (theL - aLower) % theArray.Length(); - 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()); + std::reverse(&theArray(aLower), &theArray(aL) + 1); + std::reverse(&theArray(aL + 1), &theArray(aUpper) + 1); } // Removes a knot from the B-spline curve @@ -395,24 +430,21 @@ void BSplCLib_InsertKnot(const Standard_Integer, 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); + BSplCLib_KnotArrays<1> aSingleKnot(U, UMult); + TColStd_Array1OfReal aNewKnots(1, Knots.Length() + 1); + TColStd_Array1OfInteger aNewMults(1, Knots.Length() + 1); BSplCLib_InsertKnots(Degree, Periodic, Poles, Weights, Knots, Mults, - k, - &m, + aSingleKnot.Knot, + &aSingleKnot.Mult, NewPoles, NewWeights, - nk, - nm, + aNewKnots, + aNewMults, Epsilon(U), Standard_True); } @@ -440,25 +472,22 @@ void BSplCLib_RaiseMultiplicity(const Standard_Integer KnotIndex, 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()); + BSplCLib_KnotArrays<1> aSingleKnot(Knots(KnotIndex), Mult - Mults(KnotIndex)); + TColStd_Array1OfReal aNewKnots(1, Knots.Length()); + TColStd_Array1OfInteger aNewMults(1, Knots.Length()); BSplCLib_InsertKnots(Degree, Periodic, Poles, Weights, Knots, Mults, - k, - &m, + aSingleKnot.Knot, + &aSingleKnot.Mult, NewPoles, NewWeights, - nk, - nm, - Epsilon(k(1)), + aNewKnots, + aNewMults, + Epsilon(Knots(KnotIndex)), Standard_True); } @@ -1910,4 +1939,58 @@ void BSplCLib_FunctionReparameterise(const BSplCLib_EvaluatorFunction& FunctionP theStatus); } +// Computes coefficients of a Bezier curve from poles (Bezier syntax wrapper) +// Uses stack-allocated flat knots for efficiency. +// @param[in] Poles - Bezier poles array +// @param[in] Weights - optional weights for rational curves +// @param[out] CachePoles - cached poles array +// @param[out] CacheWeights - cached weights array +template +void BSplCLib_PolesCoefficients_Bezier(const Array1OfPoints& Poles, + const TColStd_Array1OfReal* Weights, + Array1OfPoints& CachePoles, + TColStd_Array1OfReal* CacheWeights) +{ + const int aDegree = Poles.Length() - 1; + TColStd_Array1OfReal aBidFlatKnots(BSplCLib::FlatBezierKnots(aDegree), 1, 2 * (aDegree + 1)); + BSplCLib_BuildCache(0., + 1., + Standard_False, + aDegree, + aBidFlatKnots, + Poles, + Weights, + CachePoles, + CacheWeights); +} + +// Increases the degree of a Bezier curve (Bezier syntax wrapper) +// Uses stack-allocated knot arrays for efficiency. +// @param[in] NewDegree - new degree +// @param[in] Poles - Bezier poles array +// @param[in] Weights - optional weights for rational curves +// @param[out] NewPoles - resulting poles array +// @param[out] NewWeights - resulting weights array +template +void BSplCLib_IncreaseDegree_Bezier(const Standard_Integer NewDegree, + const Array1OfPoints& Poles, + const TColStd_Array1OfReal* Weights, + Array1OfPoints& NewPoles, + TColStd_Array1OfReal* NewWeights) +{ + const int aDegree = Poles.Length() - 1; + BSplCLib_KnotArrays<2> aBezierKnots(aDegree); + BSplCLib_IncreaseDegree(aDegree, + NewDegree, + Standard_False, + Poles, + Weights, + aBezierKnots.Knot, + aBezierKnots.Mult, + NewPoles, + NewWeights, + aBezierKnots.Knot, + aBezierKnots.Mult); +} + #endif // _BSplCLib_CurveComputation_pxx_HeaderFile diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_EvaluatorFunction.hxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_EvaluatorFunction.hxx index 3755c09505..77ebb2dde8 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_EvaluatorFunction.hxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_EvaluatorFunction.hxx @@ -51,12 +51,9 @@ public: Evaluate(theDerivativeRequest, theStartEnd, theParameter, theResult, theErrorCode); } -private: - //! Copy constructor is declared private to forbid copying - BSplCLib_EvaluatorFunction(const BSplCLib_EvaluatorFunction&) {} - - //! Assignment operator is declared private to forbid copying - void operator=(const BSplCLib_EvaluatorFunction&) {} + // copying is prohibited + BSplCLib_EvaluatorFunction(const BSplCLib_EvaluatorFunction&) = delete; + BSplCLib_EvaluatorFunction& operator=(const BSplCLib_EvaluatorFunction&) = delete; }; #endif -- 2.39.5