1 // Created on: 1995-08-28
2 // Created by: Laurent BOURESCHE
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // Modified: 28/02/1996 by PMN : HermiteCoefficients added
18 // Modified: 18/06/1996 by PMN : NULL reference.
19 // Modified: 19/02/1997 by JCT : EvalPoly2Var added
21 #include <GeomAbs_Shape.hxx>
23 #include <math_Gauss.hxx>
24 #include <math_Matrix.hxx>
25 #include <NCollection_LocalArray.hxx>
27 #include <Standard_ConstructionError.hxx>
29 // To convert points array into Real ..
30 // *********************************
31 //=======================================================================
34 //=======================================================================
35 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
36 TColStd_Array1OfReal& FP)
38 Standard_Integer j = FP .Lower();
39 Standard_Integer PLower = Poles.Lower();
40 Standard_Integer PUpper = Poles.Upper();
42 for (Standard_Integer i = PLower; i <= PUpper; i++) {
43 const gp_Pnt2d& P = Poles(i);
44 FP(j) = P.Coord(1); j++;
45 FP(j) = P.Coord(2); j++;
49 //=======================================================================
52 //=======================================================================
54 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
55 const TColStd_Array1OfReal& Weights,
56 TColStd_Array1OfReal& FP)
58 Standard_Integer j = FP .Lower();
59 Standard_Integer PLower = Poles.Lower();
60 Standard_Integer PUpper = Poles.Upper();
62 for (Standard_Integer i = PLower; i <= PUpper; i++) {
63 Standard_Real w = Weights(i);
64 const gp_Pnt2d& P = Poles(i);
65 FP(j) = P.Coord(1) * w; j++;
66 FP(j) = P.Coord(2) * w; j++;
71 //=======================================================================
74 //=======================================================================
76 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
77 TColgp_Array1OfPnt2d& Poles)
79 Standard_Integer j = FP .Lower();
80 Standard_Integer PLower = Poles.Lower();
81 Standard_Integer PUpper = Poles.Upper();
83 for (Standard_Integer i = PLower; i <= PUpper; i++) {
84 gp_Pnt2d& P = Poles(i);
85 P.SetCoord(1,FP(j)); j++;
86 P.SetCoord(2,FP(j)); j++;
90 //=======================================================================
93 //=======================================================================
95 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
96 TColgp_Array1OfPnt2d& Poles,
97 TColStd_Array1OfReal& Weights)
99 Standard_Integer j = FP .Lower();
100 Standard_Integer PLower = Poles.Lower();
101 Standard_Integer PUpper = Poles.Upper();
103 for (Standard_Integer i = PLower; i <= PUpper; i++) {
104 Standard_Real w = FP(j + 2);
106 gp_Pnt2d& P = Poles(i);
107 P.SetCoord(1,FP(j) / w); j++;
108 P.SetCoord(2,FP(j) / w); j++;
113 //=======================================================================
114 //function : SetPoles
116 //=======================================================================
118 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
119 TColStd_Array1OfReal& FP)
121 Standard_Integer j = FP .Lower();
122 Standard_Integer PLower = Poles.Lower();
123 Standard_Integer PUpper = Poles.Upper();
125 for (Standard_Integer i = PLower; i <= PUpper; i++) {
126 const gp_Pnt& P = Poles(i);
127 FP(j) = P.Coord(1); j++;
128 FP(j) = P.Coord(2); j++;
129 FP(j) = P.Coord(3); j++;
133 //=======================================================================
134 //function : SetPoles
136 //=======================================================================
138 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
139 const TColStd_Array1OfReal& Weights,
140 TColStd_Array1OfReal& FP)
142 Standard_Integer j = FP .Lower();
143 Standard_Integer PLower = Poles.Lower();
144 Standard_Integer PUpper = Poles.Upper();
146 for (Standard_Integer i = PLower; i <= PUpper; i++) {
147 Standard_Real w = Weights(i);
148 const gp_Pnt& P = Poles(i);
149 FP(j) = P.Coord(1) * w; j++;
150 FP(j) = P.Coord(2) * w; j++;
151 FP(j) = P.Coord(3) * w; j++;
156 //=======================================================================
157 //function : GetPoles
159 //=======================================================================
161 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
162 TColgp_Array1OfPnt& Poles)
164 Standard_Integer j = FP .Lower();
165 Standard_Integer PLower = Poles.Lower();
166 Standard_Integer PUpper = Poles.Upper();
168 for (Standard_Integer i = PLower; i <= PUpper; i++) {
169 gp_Pnt& P = Poles(i);
170 P.SetCoord(1,FP(j)); j++;
171 P.SetCoord(2,FP(j)); j++;
172 P.SetCoord(3,FP(j)); j++;
176 //=======================================================================
177 //function : GetPoles
179 //=======================================================================
181 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
182 TColgp_Array1OfPnt& Poles,
183 TColStd_Array1OfReal& Weights)
185 Standard_Integer j = FP .Lower();
186 Standard_Integer PLower = Poles.Lower();
187 Standard_Integer PUpper = Poles.Upper();
189 for (Standard_Integer i = PLower; i <= PUpper; i++) {
190 Standard_Real w = FP(j + 3);
192 gp_Pnt& P = Poles(i);
193 P.SetCoord(1,FP(j) / w); j++;
194 P.SetCoord(2,FP(j) / w); j++;
195 P.SetCoord(3,FP(j) / w); j++;
200 // specialized allocator
209 BinomAllocator (const Standard_Integer theMaxBinom)
211 myMaxBinom (theMaxBinom)
213 Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
214 Standard_Integer np1 = myMaxBinom + 1;
215 myBinom = new Standard_Integer*[np1];
216 myBinom[0] = new Standard_Integer[1];
218 for (i = 1; i < np1; ++i)
226 myBinom[i] = new Standard_Integer[ip1];
228 for (j = 0; j < id2; ++j)
230 myBinom[i][j] = k + myBinom[im1][j];
234 if (j > md2) j = im1 - j;
235 myBinom[i][id2] = k + myBinom[im1][j];
237 for (j = ip1 - md3; j < ip1; j++)
239 myBinom[i][j] = myBinom[i][i - j];
248 for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
255 Standard_Real Value (const Standard_Integer N,
256 const Standard_Integer P) const
258 Standard_OutOfRange_Raise_if (N > myMaxBinom,
259 "PLib, BinomAllocator: requested degree is greater than maximum supported");
260 return Standard_Real (myBinom[N][P]);
264 BinomAllocator (const BinomAllocator&);
265 BinomAllocator& operator= (const BinomAllocator&);
268 Standard_Integer** myBinom;
269 Standard_Integer myMaxBinom;
273 // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
274 //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
275 static BinomAllocator THE_BINOM (25 + 1);
278 //=======================================================================
281 //=======================================================================
283 Standard_Real PLib::Bin(const Standard_Integer N,
284 const Standard_Integer P)
286 return THE_BINOM.Value (N, P);
289 //=======================================================================
290 //function : RationalDerivative
292 //=======================================================================
294 void PLib::RationalDerivative(const Standard_Integer Degree,
295 const Standard_Integer DerivativeRequest,
296 const Standard_Integer Dimension,
298 Standard_Real& RDers,
299 const Standard_Boolean All)
302 // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
305 // Let C(N,P) be the binomial
310 // u = SUM C (q,p) f v
317 // (q) ( (q) (p) (q-p) )
318 // f = (1/v) ( u - SUM C (q,p) f v )
322 // make arrays for the binomial since computing it each time could raise a performance
324 // As oppose to the method below the <Der> array is organized in the following
327 // u (1) u (2) .... u (Dimension) v (1)
330 // u (1) u (2) .... u (Dimension) v (1)
332 // ............................................
334 // (Degree) (Degree) (Degree) (Degree)
335 // u (1) u (2) .... u (Dimension) v (1)
338 Standard_Real Inverse;
339 Standard_Real *PolesArray = &Ders;
340 Standard_Real *RationalArray = &RDers;
341 Standard_Real Factor ;
342 Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
343 NCollection_LocalArray<Standard_Real> binomial_array;
344 NCollection_LocalArray<Standard_Real> derivative_storage;
345 if (Dimension == 3) {
346 Standard_Integer DeRequest1 = DerivativeRequest + 1;
347 Standard_Integer MinDegRequ = DerivativeRequest;
348 if (MinDegRequ > Degree) MinDegRequ = Degree;
349 binomial_array.Allocate (DeRequest1);
351 for (ii = 0 ; ii < DeRequest1 ; ii++) {
352 binomial_array[ii] = 1.0e0 ;
355 Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
356 derivative_storage.Allocate (DimDeRequ1);
357 RationalArray = derivative_storage ;
360 Inverse = 1.0e0 / PolesArray[3] ;
365 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
368 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
369 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
370 RationalArray[Index] = PolesArray[OtherIndex];
374 for (jj = ii - 1 ; jj >= 0 ; jj--) {
375 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
376 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
377 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
378 RationalArray[Index] -= Factor * RationalArray[Index1];
383 for (jj = ii ; jj >= 1 ; jj--) {
384 binomial_array[jj] += binomial_array[jj - 1] ;
386 RationalArray[Index] *= Inverse ; Index++;
387 RationalArray[Index] *= Inverse ; Index++;
388 RationalArray[Index] *= Inverse ; Index++;
391 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
394 RationalArray[Index] = 0.0e0; Index++;
395 RationalArray[Index] = 0.0e0; Index++;
396 RationalArray[Index] = 0.0e0;
399 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
400 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
401 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
402 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
403 RationalArray[Index] -= Factor * RationalArray[Index1];
408 for (jj = ii ; jj >= 1 ; jj--) {
409 binomial_array[jj] += binomial_array[jj - 1] ;
411 RationalArray[Index] *= Inverse; Index++;
412 RationalArray[Index] *= Inverse; Index++;
413 RationalArray[Index] *= Inverse; Index++;
417 RationalArray = &RDers ;
418 Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
419 RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
420 RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
421 RationalArray[2] = derivative_storage[DimDeRequ];
426 Standard_Integer Dimension1 = Dimension + 1;
427 Standard_Integer Dimension2 = Dimension << 1;
428 Standard_Integer DeRequest1 = DerivativeRequest + 1;
429 Standard_Integer MinDegRequ = DerivativeRequest;
430 if (MinDegRequ > Degree) MinDegRequ = Degree;
431 binomial_array.Allocate (DeRequest1);
433 for (ii = 0 ; ii < DeRequest1 ; ii++) {
434 binomial_array[ii] = 1.0e0 ;
437 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
438 derivative_storage.Allocate (DimDeRequ1);
439 RationalArray = derivative_storage ;
442 Inverse = 1.0e0 / PolesArray[Dimension] ;
444 Index2 = - Dimension2;
447 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
451 for (kk = 0 ; kk < Dimension ; kk++) {
452 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
457 for (jj = ii - 1 ; jj >= 0 ; jj--) {
458 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
460 for (kk = 0 ; kk < Dimension ; kk++) {
461 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
464 Index1 -= Dimension2 ;
467 for (jj = ii ; jj >= 1 ; jj--) {
468 binomial_array[jj] += binomial_array[jj - 1] ;
471 for (kk = 0 ; kk < Dimension ; kk++) {
472 RationalArray[Index] *= Inverse ; Index++;
476 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
480 for (kk = 0 ; kk < Dimension ; kk++) {
481 RationalArray[Index] = 0.0e0 ; Index++;
485 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
486 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
488 for (kk = 0 ; kk < Dimension ; kk++) {
489 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
492 Index1 -= Dimension2 ;
495 for (jj = ii ; jj >= 1 ; jj--) {
496 binomial_array[jj] += binomial_array[jj - 1] ;
499 for (kk = 0 ; kk < Dimension ; kk++) {
500 RationalArray[Index] *= Inverse; Index++;
505 RationalArray = &RDers ;
506 Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
508 for (kk = 0 ; kk < Dimension ; kk++) {
509 RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
515 //=======================================================================
516 //function : RationalDerivatives
517 //purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights
518 //=======================================================================
520 void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
521 const Standard_Integer Dimension,
522 Standard_Real& PolesDerivates,
523 // must be an array with
524 // (DerivativeRequest + 1) * Dimension slots
525 Standard_Real& WeightsDerivates,
526 // must be an array with
527 // (DerivativeRequest + 1) slots
528 Standard_Real& RationalDerivates)
531 // Our purpose is to compute f = (u/v) derivated N times
534 // Let C(N,P) be the binomial
539 // u = SUM C (q,p) f v
546 // (q) ( (q) (p) (q-p) )
547 // f = (1/v) ( u - SUM C (q,p) f v )
551 // make arrays for the binomial since computing it each time could
552 // raize a performance issue
554 Standard_Real Inverse;
555 Standard_Real *PolesArray = &PolesDerivates;
556 Standard_Real *WeightsArray = &WeightsDerivates;
557 Standard_Real *RationalArray = &RationalDerivates;
558 Standard_Real Factor ;
560 Standard_Integer ii, Index, Index1, Index2, jj;
561 Standard_Integer DeRequest1 = DerivativeRequest + 1;
563 NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
564 NCollection_LocalArray<Standard_Real> derivative_storage;
566 for (ii = 0 ; ii < DeRequest1 ; ii++) {
567 binomial_array[ii] = 1.0e0 ;
569 Inverse = 1.0e0 / WeightsArray[0] ;
570 if (Dimension == 3) {
574 for (ii = 0 ; ii < DeRequest1 ; ii++) {
577 RationalArray[Index] = PolesArray[Index] ; Index++;
578 RationalArray[Index] = PolesArray[Index] ; Index++;
579 RationalArray[Index] = PolesArray[Index] ;
582 for (jj = ii - 1 ; jj >= 0 ; jj--) {
583 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
584 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
585 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
586 RationalArray[Index] -= Factor * RationalArray[Index1];
591 for (jj = ii ; jj >= 1 ; jj--) {
592 binomial_array[jj] += binomial_array[jj - 1] ;
594 RationalArray[Index] *= Inverse ; Index++;
595 RationalArray[Index] *= Inverse ; Index++;
596 RationalArray[Index] *= Inverse ; Index++;
601 Standard_Integer Dimension2 = Dimension << 1;
603 Index2 = - Dimension2;
605 for (ii = 0 ; ii < DeRequest1 ; ii++) {
609 for (kk = 0 ; kk < Dimension ; kk++) {
610 RationalArray[Index] = PolesArray[Index]; Index++;
614 for (jj = ii - 1 ; jj >= 0 ; jj--) {
615 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
617 for (kk = 0 ; kk < Dimension ; kk++) {
618 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
621 Index1 -= Dimension2;
624 for (jj = ii ; jj >= 1 ; jj--) {
625 binomial_array[jj] += binomial_array[jj - 1] ;
628 for (kk = 0 ; kk < Dimension ; kk++) {
629 RationalArray[Index] *= Inverse ; Index++;
635 //=======================================================================
636 // Auxiliary template functions used for optimized evaluation of polynome
637 // and its derivatives for smaller dimensions of the polynome
638 //=======================================================================
641 // recursive template for evaluating value or first derivative
643 inline void eval_step1 (double* poly, double par, double* coef)
645 eval_step1<dim - 1> (poly, par, coef);
646 poly[dim] = poly[dim] * par + coef[dim];
651 inline void eval_step1<-1> (double*, double, double*)
655 // recursive template for evaluating second derivative
657 inline void eval_step2 (double* poly, double par, double* coef)
659 eval_step2<dim - 1> (poly, par, coef);
660 poly[dim] = poly[dim] * par + coef[dim] * 2.;
665 inline void eval_step2<-1> (double*, double, double*)
669 // evaluation of only value
671 inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
673 Standard_Real* aRes0 = aRes;
674 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
676 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
679 // Calculating the value of the polynomial
680 eval_step1<dim-1> (aRes0, Par, aCoeffs);
684 // evaluation of value and first derivative
686 inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
688 Standard_Real* aRes0 = aRes;
689 Standard_Real* aRes1 = aRes + dim;
691 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
692 memset(aRes1, 0, sizeof(Standard_Real) * dim);
694 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
697 // Calculating derivatives of the polynomial
698 eval_step1<dim-1> (aRes1, Par, aRes0);
699 // Calculating the value of the polynomial
700 eval_step1<dim-1> (aRes0, Par, aCoeffs);
704 // evaluation of value and first and second derivatives
706 inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
708 Standard_Real* aRes0 = aRes;
709 Standard_Real* aRes1 = aRes + dim;
710 Standard_Real* aRes2 = aRes + 2 * dim;
712 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
713 memset(aRes1, 0, sizeof(Standard_Real) * dim);
714 memset(aRes2, 0, sizeof(Standard_Real) * dim);
716 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
719 // Calculating second derivatives of the polynomial
720 eval_step2<dim-1> (aRes2, Par, aRes1);
721 // Calculating derivatives of the polynomial
722 eval_step1<dim-1> (aRes1, Par, aRes0);
723 // Calculating the value of the polynomial
724 eval_step1<dim-1> (aRes0, Par, aCoeffs);
729 //=======================================================================
730 //function : This evaluates a polynomial and its derivatives
731 //purpose : up to the requested order
732 //=======================================================================
734 void PLib::EvalPolynomial(const Standard_Real Par,
735 const Standard_Integer DerivativeRequest,
736 const Standard_Integer Degree,
737 const Standard_Integer Dimension,
738 Standard_Real& PolynomialCoeff,
739 Standard_Real& Results)
741 // the polynomial coefficients are assumed to be stored as follows :
743 // [0] [Dimension -1] X coefficient
745 // [Dimension] [Dimension + Dimension -1] X coefficient
747 // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient
749 // ...................................................
753 // [d * Dimension] [d * Dimension + Dimension-1] X coefficient
755 // where d is the Degree
758 Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
759 Standard_Real* aRes = &Results;
760 Standard_Real* anOriginal;
761 Standard_Integer ind = 0;
762 switch (DerivativeRequest)
768 case 1: eval_poly1<1> (aRes, aCoeffs, Degree, Par); break;
769 case 2: eval_poly1<2> (aRes, aCoeffs, Degree, Par); break;
770 case 3: eval_poly1<3> (aRes, aCoeffs, Degree, Par); break;
771 case 4: eval_poly1<4> (aRes, aCoeffs, Degree, Par); break;
772 case 5: eval_poly1<5> (aRes, aCoeffs, Degree, Par); break;
773 case 6: eval_poly1<6> (aRes, aCoeffs, Degree, Par); break;
774 case 7: eval_poly1<7> (aRes, aCoeffs, Degree, Par); break;
775 case 8: eval_poly1<8> (aRes, aCoeffs, Degree, Par); break;
776 case 9: eval_poly1<9> (aRes, aCoeffs, Degree, Par); break;
777 case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
778 case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
779 case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
780 case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
781 case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
782 case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
785 Standard_Real* aRes0 = aRes;
786 Standard_Real* aRes1 = aRes + Dimension;
788 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
789 memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
791 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
793 aCoeffs -= Dimension;
794 // Calculating derivatives of the polynomial
795 for (ind = 0; ind < Dimension; ind++)
796 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
797 // Calculating the value of the polynomial
798 for (ind = 0; ind < Dimension; ind++)
799 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
809 case 1: eval_poly2<1> (aRes, aCoeffs, Degree, Par); break;
810 case 2: eval_poly2<2> (aRes, aCoeffs, Degree, Par); break;
811 case 3: eval_poly2<3> (aRes, aCoeffs, Degree, Par); break;
812 case 4: eval_poly2<4> (aRes, aCoeffs, Degree, Par); break;
813 case 5: eval_poly2<5> (aRes, aCoeffs, Degree, Par); break;
814 case 6: eval_poly2<6> (aRes, aCoeffs, Degree, Par); break;
815 case 7: eval_poly2<7> (aRes, aCoeffs, Degree, Par); break;
816 case 8: eval_poly2<8> (aRes, aCoeffs, Degree, Par); break;
817 case 9: eval_poly2<9> (aRes, aCoeffs, Degree, Par); break;
818 case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
819 case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
820 case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
821 case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
822 case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
823 case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
826 Standard_Real* aRes0 = aRes;
827 Standard_Real* aRes1 = aRes + Dimension;
828 Standard_Real* aRes2 = aRes1 + Dimension;
830 // Nullify the results
831 Standard_Integer aSize = 2 * Dimension;
832 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
833 memset(aRes1, 0, sizeof(Standard_Real) * aSize);
835 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
837 aCoeffs -= Dimension;
838 // Calculating derivatives of the polynomial
839 for (ind = 0; ind < Dimension; ind++)
840 aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
841 for (ind = 0; ind < Dimension; ind++)
842 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
843 // Calculating the value of the polynomial
844 for (ind = 0; ind < Dimension; ind++)
845 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
854 // Nullify the results
855 Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
856 memset(aRes, 0, sizeof(Standard_Real) * aResSize);
858 for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
860 aRes = &Results + aResSize - Dimension;
861 // Calculating derivatives of the polynomial
862 for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
864 anOriginal = aRes - Dimension; // pointer to the derivative minus 1
865 for (ind = 0; ind < Dimension; ind++)
866 aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
869 // Calculating the value of the polynomial
870 for (ind = 0; ind < Dimension; ind++)
871 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
872 aCoeffs -= Dimension;
878 //=======================================================================
879 //function : This evaluates a polynomial without derivative
881 //=======================================================================
883 void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par,
884 const Standard_Integer Degree,
885 const Standard_Integer Dimension,
886 const Standard_Integer DegreeDimension,
887 Standard_Real& PolynomialCoeff,
888 Standard_Real& Results)
890 Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
891 Standard_Real* aRes = &Results;
895 case 1: eval_poly0<1> (aRes, aCoeffs, Degree, Par); break;
896 case 2: eval_poly0<2> (aRes, aCoeffs, Degree, Par); break;
897 case 3: eval_poly0<3> (aRes, aCoeffs, Degree, Par); break;
898 case 4: eval_poly0<4> (aRes, aCoeffs, Degree, Par); break;
899 case 5: eval_poly0<5> (aRes, aCoeffs, Degree, Par); break;
900 case 6: eval_poly0<6> (aRes, aCoeffs, Degree, Par); break;
901 case 7: eval_poly0<7> (aRes, aCoeffs, Degree, Par); break;
902 case 8: eval_poly0<8> (aRes, aCoeffs, Degree, Par); break;
903 case 9: eval_poly0<9> (aRes, aCoeffs, Degree, Par); break;
904 case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
905 case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
906 case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
907 case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
908 case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
909 case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
912 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
913 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
915 aCoeffs -= Dimension;
916 for (Standard_Integer ind = 0; ind < Dimension; ind++)
917 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
923 //=======================================================================
924 //function : This evaluates a polynomial of 2 variables
925 //purpose : or its derivative at the requested orders
926 //=======================================================================
928 void PLib::EvalPoly2Var(const Standard_Real UParameter,
929 const Standard_Real VParameter,
930 const Standard_Integer UDerivativeRequest,
931 const Standard_Integer VDerivativeRequest,
932 const Standard_Integer UDegree,
933 const Standard_Integer VDegree,
934 const Standard_Integer Dimension,
935 Standard_Real& PolynomialCoeff,
936 Standard_Real& Results)
938 // the polynomial coefficients are assumed to be stored as follows :
940 // [0] [Dimension -1] U V coefficient
942 // [Dimension] [Dimension + Dimension -1] U V coefficient
944 // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient
946 // ...................................................
950 // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient
955 // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient
957 // ...................................................
960 // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient
962 // ...................................................
965 // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient
969 Standard_Integer Udim = (VDegree+1)*Dimension,
970 index = Udim*UDerivativeRequest;
971 TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
972 TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
973 Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1);
974 Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1);
975 Standard_Real * ResultArray ;
976 ResultArray = &Results ;
978 PLib::EvalPolynomial(UParameter,
985 PLib::EvalPolynomial(VParameter,
992 index = Dimension*VDerivativeRequest;
994 for (Standard_Integer i=0;i<Dimension;i++) {
995 ResultArray[i] = Digit[index+i];
1001 //=======================================================================
1002 //function : This evaluates the lagrange polynomial and its derivatives
1003 //purpose : up to the requested order that interpolates a series of
1004 //points of dimension <Dimension> with given assigned parameters
1005 //=======================================================================
1008 PLib::EvalLagrange(const Standard_Real Parameter,
1009 const Standard_Integer DerivativeRequest,
1010 const Standard_Integer Degree,
1011 const Standard_Integer Dimension,
1012 Standard_Real& Values,
1013 Standard_Real& Parameters,
1014 Standard_Real& Results)
1017 // the points are assumed to be stored as follows in the Values array :
1019 // [0] [Dimension -1] first point coefficients
1021 // [Dimension] [Dimension + Dimension -1] second point coefficients
1023 // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients
1025 // ...................................................
1029 // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients
1031 // where d is the Degree
1033 // The ParameterArray stores the parameter value assign to each point in
1034 // order described above, that is
1035 // [0] is assign to first point
1036 // [1] is assign to second point
1038 Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1039 Standard_Integer local_request = DerivativeRequest;
1040 Standard_Real *ParameterArray;
1041 Standard_Real difference;
1042 Standard_Real *PointsArray;
1043 Standard_Real *ResultArray ;
1045 PointsArray = &Values ;
1046 ParameterArray = &Parameters ;
1047 ResultArray = &Results ;
1048 if (local_request >= Degree) {
1049 local_request = Degree ;
1051 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1053 // Build the divided differences array
1056 for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) {
1057 divided_differences_array[ii] = PointsArray[ii] ;
1060 for (ii = Degree ; ii >= 0 ; ii--) {
1062 for (jj = Degree ; jj > Degree - ii ; jj--) {
1063 Index = jj * Dimension ;
1064 Index1 = Index - Dimension ;
1066 for (kk = 0 ; kk < Dimension ; kk++) {
1067 divided_differences_array[Index + kk] -=
1068 divided_differences_array[Index1 + kk] ;
1071 ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ;
1072 if (Abs(difference) < RealSmall()) {
1076 difference = 1.0e0 / difference ;
1078 for (kk = 0 ; kk < Dimension ; kk++) {
1079 divided_differences_array[Index + kk] *= difference ;
1085 // Evaluate the divided difference array polynomial which expresses as
1087 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1088 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1090 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1093 Index = Degree * Dimension ;
1095 for (kk = 0 ; kk < Dimension ; kk++) {
1096 ResultArray[kk] = divided_differences_array[Index + kk] ;
1099 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1100 ResultArray[ii] = 0.0e0 ;
1103 for (ii = Degree ; ii >= 1 ; ii--) {
1104 difference = Parameter - ParameterArray[ii - 1] ;
1106 for (jj = local_request ; jj > 0 ; jj--) {
1107 Index = jj * Dimension ;
1108 Index1 = Index - Dimension ;
1110 for (kk = 0 ; kk < Dimension ; kk++) {
1111 ResultArray[Index + kk] *= difference ;
1112 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
1115 Index = (ii -1) * Dimension ;
1117 for (kk = 0 ; kk < Dimension ; kk++) {
1118 ResultArray[kk] *= difference ;
1119 ResultArray[kk] += divided_differences_array[Index+kk] ;
1123 return (ReturnCode) ;
1126 //=======================================================================
1127 //function : This evaluates the hermite polynomial and its derivatives
1128 //purpose : up to the requested order that interpolates a series of
1129 //points of dimension <Dimension> with given assigned parameters
1130 //=======================================================================
1132 Standard_Integer PLib::EvalCubicHermite
1133 (const Standard_Real Parameter,
1134 const Standard_Integer DerivativeRequest,
1135 const Standard_Integer Dimension,
1136 Standard_Real& Values,
1137 Standard_Real& Derivatives,
1138 Standard_Real& theParameters,
1139 Standard_Real& Results)
1142 // the points are assumed to be stored as follows in the Values array :
1144 // [0] [Dimension -1] first point coefficients
1146 // [Dimension] [Dimension + Dimension -1] last point coefficients
1149 // the derivatives are assumed to be stored as follows in
1150 // the Derivatives array :
1152 // [0] [Dimension -1] first point coefficients
1154 // [Dimension] [Dimension + Dimension -1] last point coefficients
1156 // The ParameterArray stores the parameter value assign to each point in
1157 // order described above, that is
1158 // [0] is assign to first point
1159 // [1] is assign to last point
1161 Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
1162 Standard_Integer local_request = DerivativeRequest ;
1166 Standard_Real ParametersArray[4];
1167 Standard_Real difference;
1168 Standard_Real inverse;
1169 Standard_Real *FirstLast;
1170 Standard_Real *PointsArray;
1171 Standard_Real *DerivativesArray;
1172 Standard_Real *ResultArray ;
1174 DerivativesArray = &Derivatives ;
1175 PointsArray = &Values ;
1176 FirstLast = &theParameters ;
1177 ResultArray = &Results ;
1178 if (local_request >= Degree) {
1179 local_request = Degree ;
1181 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1183 for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) {
1184 ParametersArray[jj] =
1185 ParametersArray[jj+1] = FirstLast[ii] ;
1188 // Build the divided differences array
1191 // initialise it at the stage 2 of the building algorithm
1192 // for devided differences
1194 inverse = FirstLast[1] - FirstLast[0] ;
1195 inverse = 1.0e0 / inverse ;
1197 for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ;
1199 ii++, jj++, kk++, pp++) {
1200 divided_differences_array[ii] = PointsArray[ii] ;
1201 divided_differences_array[kk] = inverse *
1202 (PointsArray[jj] - PointsArray[ii]) ;
1203 divided_differences_array[jj] = DerivativesArray[ii] ;
1204 divided_differences_array[pp] = DerivativesArray[jj] ;
1207 for (ii = 1 ; ii <= Degree ; ii++) {
1209 for (jj = Degree ; jj >= ii+1 ; jj--) {
1210 Index = jj * Dimension ;
1211 Index1 = Index - Dimension ;
1213 for (kk = 0 ; kk < Dimension ; kk++) {
1214 divided_differences_array[Index + kk] -=
1215 divided_differences_array[Index1 + kk] ;
1218 for (kk = 0 ; kk < Dimension ; kk++) {
1219 divided_differences_array[Index + kk] *= inverse ;
1225 // Evaluate the divided difference array polynomial which expresses as
1227 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1228 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1230 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1233 Index = Degree * Dimension ;
1235 for (kk = 0 ; kk < Dimension ; kk++) {
1236 ResultArray[kk] = divided_differences_array[Index + kk] ;
1239 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1240 ResultArray[ii] = 0.0e0 ;
1243 for (ii = Degree ; ii >= 1 ; ii--) {
1244 difference = Parameter - ParametersArray[ii - 1] ;
1246 for (jj = local_request ; jj > 0 ; jj--) {
1247 Index = jj * Dimension ;
1248 Index1 = Index - Dimension ;
1250 for (kk = 0 ; kk < Dimension ; kk++) {
1251 ResultArray[Index + kk] *= difference ;
1252 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
1255 Index = (ii -1) * Dimension ;
1257 for (kk = 0 ; kk < Dimension ; kk++) {
1258 ResultArray[kk] *= difference ;
1259 ResultArray[kk] += divided_differences_array[Index+kk] ;
1263 return (ReturnCode) ;
1266 //=======================================================================
1267 //function : HermiteCoefficients
1268 //purpose : calcul des polynomes d'Hermite
1269 //=======================================================================
1271 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
1272 const Standard_Real LastParameter,
1273 const Standard_Integer FirstOrder,
1274 const Standard_Integer LastOrder,
1275 math_Matrix& MatrixCoefs)
1277 Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2];
1278 Standard_Integer ii, jj, pp, cote, iof=0;
1279 Standard_Real Prod, TBorne = FirstParameter;
1280 math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
1281 math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
1283 // Test de validites
1285 if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
1286 Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
1287 if (D1 > 100 || D2 > 100) return Standard_False;
1289 if (D2 < 0.01) return Standard_False;
1290 if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
1292 // Calcul de la matrice a inverser (MAT)
1294 Ordre[0] = FirstOrder+1;
1295 Ordre[1] = LastOrder+1;
1297 for (cote=0; cote<=1; cote++) {
1300 for (pp=1; pp<=Ordre[cote]; pp++) {
1304 for (jj=pp; jj<=NbCoeff; jj++) {
1305 // tout se passe dans les 3 lignes suivantes
1306 MAT(ii, jj) = Coeff(jj) * Prod;
1307 Coeff(jj) *= jj - pp;
1311 TBorne = LastParameter;
1315 // resolution du systemes
1316 math_Gauss ResolCoeff(MAT, 1.0e-10);
1317 if (!ResolCoeff.IsDone()) return Standard_False;
1319 for (ii=1; ii<=NbCoeff; ii++) {
1321 ResolCoeff.Solve(B, Coeff);
1322 MatrixCoefs.SetRow( ii, Coeff);
1325 return Standard_True;
1328 //=======================================================================
1329 //function : CoefficientsPoles
1331 //=======================================================================
1333 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs,
1334 const TColStd_Array1OfReal* WCoefs,
1335 TColgp_Array1OfPnt& Poles,
1336 TColStd_Array1OfReal* Weights)
1338 TColStd_Array1OfReal tempC(1,3*Coefs.Length());
1339 PLib::SetPoles(Coefs,tempC);
1340 TColStd_Array1OfReal tempP(1,3*Poles.Length());
1341 PLib::SetPoles(Coefs,tempP);
1342 PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
1343 PLib::GetPoles(tempP,Poles);
1346 //=======================================================================
1347 //function : CoefficientsPoles
1349 //=======================================================================
1351 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
1352 const TColStd_Array1OfReal* WCoefs,
1353 TColgp_Array1OfPnt2d& Poles,
1354 TColStd_Array1OfReal* Weights)
1356 TColStd_Array1OfReal tempC(1,2*Coefs.Length());
1357 PLib::SetPoles(Coefs,tempC);
1358 TColStd_Array1OfReal tempP(1,2*Poles.Length());
1359 PLib::SetPoles(Coefs,tempP);
1360 PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
1361 PLib::GetPoles(tempP,Poles);
1364 //=======================================================================
1365 //function : CoefficientsPoles
1367 //=======================================================================
1369 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
1370 const TColStd_Array1OfReal* WCoefs,
1371 TColStd_Array1OfReal& Poles,
1372 TColStd_Array1OfReal* Weights)
1374 PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
1377 //=======================================================================
1378 //function : CoefficientsPoles
1380 //=======================================================================
1382 void PLib::CoefficientsPoles (const Standard_Integer dim,
1383 const TColStd_Array1OfReal& Coefs,
1384 const TColStd_Array1OfReal* WCoefs,
1385 TColStd_Array1OfReal& Poles,
1386 TColStd_Array1OfReal* Weights)
1388 Standard_Boolean rat = WCoefs != NULL;
1389 Standard_Integer loc = Coefs.Lower();
1390 Standard_Integer lop = Poles.Lower();
1391 Standard_Integer lowc=0;
1392 Standard_Integer lowp=0;
1393 Standard_Integer upc = Coefs.Upper();
1394 Standard_Integer upp = Poles.Upper();
1395 Standard_Integer upwc=0;
1396 Standard_Integer upwp=0;
1397 Standard_Integer reflen = Coefs.Length()/dim;
1398 Standard_Integer i,j,k;
1401 lowc = WCoefs->Lower(); lowp = Weights->Lower();
1402 upwc = WCoefs->Upper(); upwp = Weights->Upper();
1405 for (i = 0; i < dim; i++){
1406 Poles (lop + i) = Coefs (loc + i);
1407 Poles (upp - i) = Coefs (upc - i);
1410 (*Weights) (lowp) = (*WCoefs) (lowc);
1411 (*Weights) (upwp) = (*WCoefs) (upwc);
1415 for (i = 2; i < reflen; i++ ) {
1416 Cnp = PLib::Bin(reflen - 1, i - 1);
1417 if (rat) (*Weights)(lowp + i - 1) = (*WCoefs)(lowc + i - 1) / Cnp;
1419 for(j = 0; j < dim; j++){
1420 Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
1424 for (i = 1; i <= reflen - 1; i++) {
1426 for (j = reflen - 1; j >= i; j--) {
1427 if (rat) (*Weights)(lowp + j) += (*Weights)(lowp + j - 1);
1429 for(k = 0; k < dim; k++){
1430 Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
1436 for (i = 1; i <= reflen; i++) {
1438 for(j = 0; j < dim; j++){
1439 Poles(lop + dim * (i-1) + j) /= (*Weights)(lowp + i -1);
1445 //=======================================================================
1446 //function : Trimming
1448 //=======================================================================
1450 void PLib::Trimming(const Standard_Real U1,
1451 const Standard_Real U2,
1452 TColgp_Array1OfPnt& Coefs,
1453 TColStd_Array1OfReal* WCoefs)
1455 TColStd_Array1OfReal temp(1,3*Coefs.Length());
1456 PLib::SetPoles(Coefs,temp);
1457 PLib::Trimming(U1,U2,3,temp,WCoefs);
1458 PLib::GetPoles(temp,Coefs);
1461 //=======================================================================
1462 //function : Trimming
1464 //=======================================================================
1466 void PLib::Trimming(const Standard_Real U1,
1467 const Standard_Real U2,
1468 TColgp_Array1OfPnt2d& Coefs,
1469 TColStd_Array1OfReal* WCoefs)
1471 TColStd_Array1OfReal temp(1,2*Coefs.Length());
1472 PLib::SetPoles(Coefs,temp);
1473 PLib::Trimming(U1,U2,2,temp,WCoefs);
1474 PLib::GetPoles(temp,Coefs);
1477 //=======================================================================
1478 //function : Trimming
1480 //=======================================================================
1482 void PLib::Trimming(const Standard_Real U1,
1483 const Standard_Real U2,
1484 TColStd_Array1OfReal& Coefs,
1485 TColStd_Array1OfReal* WCoefs)
1487 PLib::Trimming(U1,U2,1,Coefs,WCoefs);
1490 //=======================================================================
1491 //function : Trimming
1493 //=======================================================================
1495 void PLib::Trimming(const Standard_Real U1,
1496 const Standard_Real U2,
1497 const Standard_Integer dim,
1498 TColStd_Array1OfReal& Coefs,
1499 TColStd_Array1OfReal* WCoefs)
1503 // on fait le changement de variable v = (u-U1) / (U2-U1)
1504 // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
1505 // decomposee sous la forme du schema iteratif de horner.
1507 Standard_Real lsp = U2 - U1;
1508 Standard_Integer indc, indw=0;
1509 Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
1510 Standard_Integer len = Coefs.Length()/dim;
1511 Standard_Boolean rat = WCoefs != NULL;
1514 if(len != WCoefs->Length())
1515 throw Standard_Failure("PLib::Trimming : nbcoefs/dim != nbweights !!!");
1516 upw = WCoefs->Upper();
1520 for (Standard_Integer i = 1; i <= len; i++) {
1521 Standard_Integer j ;
1522 indc = upc - dim*(i-1);
1523 if (rat) indw = upw - i + 1;
1524 //calcul du coefficient de degre le plus faible a l'iteration i
1526 for( j = 0; j < dim; j++){
1527 Coefs(indc - dim + j) += U1 * Coefs(indc + j);
1529 if (rat) (*WCoefs)(indw - 1) += U1 * (*WCoefs)(indw);
1531 //calcul des coefficients intermediaires :
1536 for(Standard_Integer k = 0; k < dim; k++){
1537 Coefs(indc - dim + k) =
1538 U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
1542 (*WCoefs)(indw - 1) = U1 * (*WCoefs)(indw) + lsp * (*WCoefs)(indw - 1);
1546 //calcul du coefficient de degre le plus eleve :
1548 for(j = 0; j < dim; j++){
1549 Coefs(upc + j) *= lsp;
1551 if (rat) (*WCoefs)(upw) *= lsp;
1555 //=======================================================================
1556 //function : CoefficientsPoles
1558 // Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852).
1559 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
1560 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
1561 //=======================================================================
1563 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs,
1564 const TColStd_Array2OfReal* WCoefs,
1565 TColgp_Array2OfPnt& Poles,
1566 TColStd_Array2OfReal* Weights)
1568 Standard_Boolean rat = (WCoefs != NULL);
1569 Standard_Integer LowerRow = Poles.LowerRow();
1570 Standard_Integer UpperRow = Poles.UpperRow();
1571 Standard_Integer LowerCol = Poles.LowerCol();
1572 Standard_Integer UpperCol = Poles.UpperCol();
1573 Standard_Integer ColLength = Poles.ColLength();
1574 Standard_Integer RowLength = Poles.RowLength();
1576 // Bidouille pour retablir u et v pour les coefs calcules
1578 // Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
1580 Standard_Integer Row, Col;
1581 Standard_Real W, Cnp;
1583 Standard_Integer I1, I2;
1584 Standard_Integer NPoleu , NPolev;
1587 for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
1588 Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol);
1590 (*Weights) (NPoleu, LowerCol) = (*WCoefs) (NPoleu, LowerCol);
1593 for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
1594 Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
1595 Temp = Coefs (NPoleu, Col).XYZ();
1597 Poles (NPoleu, Col).SetXYZ (Temp);
1599 (*Weights) (NPoleu, Col) = (*WCoefs) (NPoleu, Col) / Cnp;
1602 Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
1604 (*Weights) (NPoleu, UpperCol) = (*WCoefs) (NPoleu, UpperCol);
1607 for (I1 = 1; I1 <= RowLength - 1; I1++) {
1609 for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
1611 (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
1612 Poles (NPoleu, I2).SetXYZ (Temp);
1613 if (rat) (*Weights)(NPoleu, I2) += (*Weights)(NPoleu, I2-1);
1618 for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
1620 for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
1621 Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
1622 Temp = Poles (Row, NPolev).XYZ();
1624 Poles (Row, NPolev).SetXYZ (Temp);
1625 if (rat) (*Weights)(Row, NPolev) /= Cnp;
1628 for (I1 = 1; I1 <= ColLength - 1; I1++) {
1630 for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
1632 (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
1633 Poles (I2, NPolev).SetXYZ (Temp);
1634 if (rat) (*Weights)(I2, NPolev) += (*Weights)(I2-1, NPolev);
1640 for (Row = LowerRow; Row <= UpperRow; Row++) {
1642 for (Col = LowerCol; Col <= UpperCol; Col++) {
1643 W = (*Weights) (Row, Col);
1644 Temp = Poles(Row, Col).XYZ();
1646 Poles(Row, Col).SetXYZ (Temp);
1652 //=======================================================================
1653 //function : UTrimming
1655 //=======================================================================
1657 void PLib::UTrimming(const Standard_Real U1,
1658 const Standard_Real U2,
1659 TColgp_Array2OfPnt& Coeffs,
1660 TColStd_Array2OfReal* WCoeffs)
1662 Standard_Boolean rat = WCoeffs != NULL;
1663 Standard_Integer lr = Coeffs.LowerRow();
1664 Standard_Integer ur = Coeffs.UpperRow();
1665 Standard_Integer lc = Coeffs.LowerCol();
1666 Standard_Integer uc = Coeffs.UpperCol();
1667 TColgp_Array1OfPnt Temp (lr,ur);
1668 TColStd_Array1OfReal Temw (lr,ur);
1670 for (Standard_Integer icol = lc; icol <= uc; icol++) {
1671 Standard_Integer irow ;
1672 for ( irow = lr; irow <= ur; irow++) {
1673 Temp (irow) = Coeffs (irow, icol);
1674 if (rat) Temw (irow) = (*WCoeffs) (irow, icol);
1676 if (rat) PLib::Trimming (U1, U2, Temp, &Temw);
1677 else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
1679 for (irow = lr; irow <= ur; irow++) {
1680 Coeffs (irow, icol) = Temp (irow);
1681 if (rat) (*WCoeffs) (irow, icol) = Temw (irow);
1686 //=======================================================================
1687 //function : VTrimming
1689 //=======================================================================
1691 void PLib::VTrimming(const Standard_Real V1,
1692 const Standard_Real V2,
1693 TColgp_Array2OfPnt& Coeffs,
1694 TColStd_Array2OfReal* WCoeffs)
1696 Standard_Boolean rat = WCoeffs != NULL;
1697 Standard_Integer lr = Coeffs.LowerRow();
1698 Standard_Integer ur = Coeffs.UpperRow();
1699 Standard_Integer lc = Coeffs.LowerCol();
1700 Standard_Integer uc = Coeffs.UpperCol();
1701 TColgp_Array1OfPnt Temp (lc,uc);
1702 TColStd_Array1OfReal Temw (lc,uc);
1704 for (Standard_Integer irow = lr; irow <= ur; irow++) {
1705 Standard_Integer icol ;
1706 for ( icol = lc; icol <= uc; icol++) {
1707 Temp (icol) = Coeffs (irow, icol);
1708 if (rat) Temw (icol) = (*WCoeffs) (irow, icol);
1710 if (rat) PLib::Trimming (V1, V2, Temp, &Temw);
1711 else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
1713 for (icol = lc; icol <= uc; icol++) {
1714 Coeffs (irow, icol) = Temp (icol);
1715 if (rat) (*WCoeffs) (irow, icol) = Temw (icol);
1720 //=======================================================================
1721 //function : HermiteInterpolate
1723 //=======================================================================
1725 Standard_Boolean PLib::HermiteInterpolate
1726 (const Standard_Integer Dimension,
1727 const Standard_Real FirstParameter,
1728 const Standard_Real LastParameter,
1729 const Standard_Integer FirstOrder,
1730 const Standard_Integer LastOrder,
1731 const TColStd_Array2OfReal& FirstConstr,
1732 const TColStd_Array2OfReal& LastConstr,
1733 TColStd_Array1OfReal& Coefficients)
1735 Standard_Real Pattern[3][6];
1737 // portage HP : il faut les initialiser 1 par 1
1758 math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
1759 // The initialisation of the matrix A
1760 Standard_Integer irow ;
1761 for ( irow=0; irow<=FirstOrder; irow++) {
1762 Standard_Real FirstVal = 1.;
1764 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1765 A(irow,icol) = Pattern[irow][icol]*FirstVal;
1766 if (irow <= icol) FirstVal *= FirstParameter;
1770 for (irow=0; irow<=LastOrder; irow++) {
1771 Standard_Real LastVal = 1.;
1773 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1774 A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
1775 if (irow <= icol) LastVal *= LastParameter;
1779 // The filled matrix A for FirstOrder=LastOrder=2 is:
1781 // 1 FP FP**2 FP**3 FP**4 FP**5
1782 // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter
1783 // 0 0 2 6*FP 12*FP**2 20*FP**3
1784 // 1 LP LP**2 LP**3 LP**4 LP**5
1785 // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter
1786 // 0 0 2 6*LP 12*LP**2 20*LP**3
1788 // If FirstOrder or LastOrder <=2 then some rows and columns are missing.
1790 // If FirstOrder=1 then 3th row and 6th column are missing
1791 // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
1793 math_Gauss Equations(A);
1794 // cout << "A=" << A << endl;
1796 for (Standard_Integer idim=1; idim<=Dimension; idim++) {
1797 // cout << "idim=" << idim << endl;
1799 math_Vector B(0,FirstOrder+LastOrder+1);
1800 Standard_Integer icol ;
1801 for ( icol=0; icol<=FirstOrder; icol++)
1802 B(icol) = FirstConstr(idim,icol);
1804 for (icol=0; icol<=LastOrder; icol++)
1805 B(FirstOrder+1+icol) = LastConstr(idim,icol);
1806 // cout << "B=" << B << endl;
1808 // The solving of equations system A * X = B. Then B = X
1810 // cout << "After Solving" << endl << "B=" << B << endl;
1812 if (Equations.IsDone()==Standard_False) return Standard_False;
1814 // the filling of the Coefficients
1816 for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
1817 Coefficients(Dimension*icol+idim-1) = B(icol);
1819 return Standard_True;
1822 //=======================================================================
1823 //function : JacobiParameters
1825 //=======================================================================
1827 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
1828 const Standard_Integer MaxDegree,
1829 const Standard_Integer Code,
1830 Standard_Integer& NbGaussPoints,
1831 Standard_Integer& WorkDegree)
1833 // ConstraintOrder: Ordre de contrainte aux extremites :
1834 // C0 = contraintes de passage aux bornes;
1835 // C1 = C0 + contraintes de derivees 1eres;
1836 // C2 = C1 + contraintes de derivees 2ndes.
1837 // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
1838 // d' approximation (doit etre superieur ou egal a
1839 // 2*NivConstr+2 et inferieur ou egal a 50).
1840 // Code: Code d' init. des parametres de discretisation.
1841 // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
1842 // = -5 Calcul tres rapide mais peu precis (8pts)
1843 // = -4 ' ' ' ' ' ' (10pts)
1844 // = -3 ' ' ' ' ' ' (15pts)
1845 // = -2 ' ' ' ' ' ' (20pts)
1846 // = -1 ' ' ' ' ' ' (25pts)
1847 // = 1 calcul rapide avec precision moyenne (30pts).
1848 // = 2 calcul rapide avec meilleure precision (40pts).
1849 // = 3 calcul un peu plus lent avec bonne precision (50 pts).
1850 // = 4 calcul lent avec la meilleure precision possible
1853 // The possible values of NbGaussPoints
1855 const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
1856 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
1858 Standard_Integer NivConstr=0;
1859 switch (ConstraintOrder) {
1860 case GeomAbs_C0: NivConstr = 0; break;
1861 case GeomAbs_C1: NivConstr = 1; break;
1862 case GeomAbs_C2: NivConstr = 2; break;
1864 throw Standard_ConstructionError("Invalid ConstraintOrder");
1866 if (MaxDegree < 2*NivConstr+1)
1867 throw Standard_ConstructionError("Invalid MaxDegree");
1870 WorkDegree = MaxDegree + 9;
1872 WorkDegree = MaxDegree + 6;
1874 //---> Nbre mini de points necessaires.
1875 Standard_Integer IPMIN=0;
1876 if (WorkDegree < NDEG8)
1878 else if (WorkDegree < NDEG10)
1880 else if (WorkDegree < NDEG15)
1882 else if (WorkDegree < NDEG20)
1884 else if (WorkDegree < NDEG25)
1886 else if (WorkDegree < NDEG30)
1888 else if (WorkDegree < NDEG40)
1890 else if (WorkDegree < NDEG50)
1892 else if (WorkDegree < NDEG61)
1895 throw Standard_ConstructionError("Invalid MaxDegree");
1896 // ---> Nbre de points voulus.
1897 Standard_Integer IWANT=0;
1899 case -5: IWANT=NDEG8; break;
1900 case -4: IWANT=NDEG10; break;
1901 case -3: IWANT=NDEG15; break;
1902 case -2: IWANT=NDEG20; break;
1903 case -1: IWANT=NDEG25; break;
1904 case 1: IWANT=NDEG30; break;
1905 case 2: IWANT=NDEG40; break;
1906 case 3: IWANT=NDEG50; break;
1907 case 4: IWANT=NDEG61; break;
1909 throw Standard_ConstructionError("Invalid Code");
1911 //--> NbGaussPoints est le nombre de points de discretisation de la fonction,
1912 // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
1913 // NbGaussPoints doit etre superieur strictement a WorkDegree.
1914 NbGaussPoints = Max(IPMIN,IWANT);
1915 // NbGaussPoints +=2;
1918 //=======================================================================
1919 //function : NivConstr
1920 //purpose : translates from GeomAbs_Shape to Integer
1921 //=======================================================================
1923 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
1925 Standard_Integer NivConstr=0;
1926 switch (ConstraintOrder) {
1927 case GeomAbs_C0: NivConstr = 0; break;
1928 case GeomAbs_C1: NivConstr = 1; break;
1929 case GeomAbs_C2: NivConstr = 2; break;
1931 throw Standard_ConstructionError("Invalid ConstraintOrder");
1936 //=======================================================================
1937 //function : ConstraintOrder
1938 //purpose : translates from Integer to GeomAbs_Shape
1939 //=======================================================================
1941 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
1943 GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
1944 switch (NivConstr) {
1945 case 0: ConstraintOrder = GeomAbs_C0; break;
1946 case 1: ConstraintOrder = GeomAbs_C1; break;
1947 case 2: ConstraintOrder = GeomAbs_C2; break;
1949 throw Standard_ConstructionError("Invalid NivConstr");
1951 return ConstraintOrder;
1954 //=======================================================================
1955 //function : EvalLength
1957 //=======================================================================
1959 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
1960 Standard_Real& PolynomialCoeff,
1961 const Standard_Real U1, const Standard_Real U2,
1962 Standard_Real& Length)
1964 Standard_Integer i,j,idim, degdim;
1965 Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
1967 Standard_Real *PolynomialArray = &PolynomialCoeff ;
1969 Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
1971 math_Vector GaussPoints(1,NbGaussPoints);
1972 math::GaussPoints(NbGaussPoints,GaussPoints);
1974 math_Vector GaussWeights(1,NbGaussPoints);
1975 math::GaussWeights(NbGaussPoints,GaussWeights);
1977 C1 = (U2 + U1) / 2.;
1978 C2 = (U2 - U1) / 2.;
1980 //-----------------------------------------------------------
1981 //****** Integration - Boucle sur les intervalles de GAUSS **
1982 //-----------------------------------------------------------
1986 for (j=1; j<=NbGaussPoints/2; j++) {
1987 // Integration en tenant compte de la symetrie
1988 Tran = C2 * GaussPoints(j);
1992 //****** Derivation sur la dimension de l'espace **
1994 degdim = Degree*Dimension;
1996 for (idim=0; idim<Dimension; idim++) {
1997 D1 = D2 = Degree * PolynomialArray [idim + degdim];
1998 for (i=Degree-1; i>=1; i--) {
1999 DD = i * PolynomialArray [idim + i*Dimension];
2007 //****** Integration **
2009 Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2011 //****** Fin de boucle dur les intervalles de GAUSS **
2017 //=======================================================================
2018 //function : EvalLength
2020 //=======================================================================
2022 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2023 Standard_Real& PolynomialCoeff,
2024 const Standard_Real U1, const Standard_Real U2,
2025 const Standard_Real Tol,
2026 Standard_Real& Length, Standard_Real& Error)
2029 Standard_Integer NbSubInt = 1, // Current number of subintervals
2030 MaxNbIter = 13, // Max number of iterations
2031 NbIter = 1; // Current number of iterations
2032 Standard_Real dU,OldLen,LenI;
2034 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2040 dU = (U2-U1)/NbSubInt;
2041 for (i=1; i<=NbSubInt; i++) {
2042 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2046 Error = Abs(OldLen-Length);
2048 while (Error > Tol && NbIter <= MaxNbIter);