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.
19 #include <GeomAbs_Shape.hxx>
21 #include <math_Gauss.hxx>
22 #include <math_Matrix.hxx>
23 #include <NCollection_LocalArray.hxx>
24 #include <Standard_ConstructionError.hxx>
26 // To convert points array into Real ..
27 // *********************************
28 //=======================================================================
31 //=======================================================================
32 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
33 TColStd_Array1OfReal& FP)
35 Standard_Integer j = FP .Lower();
36 Standard_Integer PLower = Poles.Lower();
37 Standard_Integer PUpper = Poles.Upper();
39 for (Standard_Integer i = PLower; i <= PUpper; i++) {
40 const gp_Pnt2d& P = Poles(i);
41 FP(j) = P.Coord(1); j++;
42 FP(j) = P.Coord(2); j++;
46 //=======================================================================
49 //=======================================================================
51 void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
52 const TColStd_Array1OfReal& Weights,
53 TColStd_Array1OfReal& FP)
55 Standard_Integer j = FP .Lower();
56 Standard_Integer PLower = Poles.Lower();
57 Standard_Integer PUpper = Poles.Upper();
59 for (Standard_Integer i = PLower; i <= PUpper; i++) {
60 Standard_Real w = Weights(i);
61 const gp_Pnt2d& P = Poles(i);
62 FP(j) = P.Coord(1) * w; j++;
63 FP(j) = P.Coord(2) * w; j++;
68 //=======================================================================
71 //=======================================================================
73 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
74 TColgp_Array1OfPnt2d& Poles)
76 Standard_Integer j = FP .Lower();
77 Standard_Integer PLower = Poles.Lower();
78 Standard_Integer PUpper = Poles.Upper();
80 for (Standard_Integer i = PLower; i <= PUpper; i++) {
81 gp_Pnt2d& P = Poles(i);
82 P.SetCoord(1,FP(j)); j++;
83 P.SetCoord(2,FP(j)); j++;
87 //=======================================================================
90 //=======================================================================
92 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
93 TColgp_Array1OfPnt2d& Poles,
94 TColStd_Array1OfReal& Weights)
96 Standard_Integer j = FP .Lower();
97 Standard_Integer PLower = Poles.Lower();
98 Standard_Integer PUpper = Poles.Upper();
100 for (Standard_Integer i = PLower; i <= PUpper; i++) {
101 Standard_Real w = FP(j + 2);
103 gp_Pnt2d& P = Poles(i);
104 P.SetCoord(1,FP(j) / w); j++;
105 P.SetCoord(2,FP(j) / w); j++;
110 //=======================================================================
111 //function : SetPoles
113 //=======================================================================
115 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
116 TColStd_Array1OfReal& FP)
118 Standard_Integer j = FP .Lower();
119 Standard_Integer PLower = Poles.Lower();
120 Standard_Integer PUpper = Poles.Upper();
122 for (Standard_Integer i = PLower; i <= PUpper; i++) {
123 const gp_Pnt& P = Poles(i);
124 FP(j) = P.Coord(1); j++;
125 FP(j) = P.Coord(2); j++;
126 FP(j) = P.Coord(3); j++;
130 //=======================================================================
131 //function : SetPoles
133 //=======================================================================
135 void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
136 const TColStd_Array1OfReal& Weights,
137 TColStd_Array1OfReal& FP)
139 Standard_Integer j = FP .Lower();
140 Standard_Integer PLower = Poles.Lower();
141 Standard_Integer PUpper = Poles.Upper();
143 for (Standard_Integer i = PLower; i <= PUpper; i++) {
144 Standard_Real w = Weights(i);
145 const gp_Pnt& P = Poles(i);
146 FP(j) = P.Coord(1) * w; j++;
147 FP(j) = P.Coord(2) * w; j++;
148 FP(j) = P.Coord(3) * w; j++;
153 //=======================================================================
154 //function : GetPoles
156 //=======================================================================
158 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
159 TColgp_Array1OfPnt& Poles)
161 Standard_Integer j = FP .Lower();
162 Standard_Integer PLower = Poles.Lower();
163 Standard_Integer PUpper = Poles.Upper();
165 for (Standard_Integer i = PLower; i <= PUpper; i++) {
166 gp_Pnt& P = Poles(i);
167 P.SetCoord(1,FP(j)); j++;
168 P.SetCoord(2,FP(j)); j++;
169 P.SetCoord(3,FP(j)); j++;
173 //=======================================================================
174 //function : GetPoles
176 //=======================================================================
178 void PLib::GetPoles(const TColStd_Array1OfReal& FP,
179 TColgp_Array1OfPnt& Poles,
180 TColStd_Array1OfReal& Weights)
182 Standard_Integer j = FP .Lower();
183 Standard_Integer PLower = Poles.Lower();
184 Standard_Integer PUpper = Poles.Upper();
186 for (Standard_Integer i = PLower; i <= PUpper; i++) {
187 Standard_Real w = FP(j + 3);
189 gp_Pnt& P = Poles(i);
190 P.SetCoord(1,FP(j) / w); j++;
191 P.SetCoord(2,FP(j) / w); j++;
192 P.SetCoord(3,FP(j) / w); j++;
197 // specialized allocator
206 BinomAllocator (const Standard_Integer theMaxBinom)
208 myMaxBinom (theMaxBinom)
210 Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
211 Standard_Integer np1 = myMaxBinom + 1;
212 myBinom = new Standard_Integer*[np1];
213 myBinom[0] = new Standard_Integer[1];
215 for (i = 1; i < np1; ++i)
223 myBinom[i] = new Standard_Integer[ip1];
225 for (j = 0; j < id2; ++j)
227 myBinom[i][j] = k + myBinom[im1][j];
231 if (j > md2) j = im1 - j;
232 myBinom[i][id2] = k + myBinom[im1][j];
234 for (j = ip1 - md3; j < ip1; j++)
236 myBinom[i][j] = myBinom[i][i - j];
245 for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
252 Standard_Real Value (const Standard_Integer N,
253 const Standard_Integer P) const
255 Standard_OutOfRange_Raise_if (N > myMaxBinom,
256 "PLib, BinomAllocator: requested degree is greater than maximum supported");
257 return Standard_Real (myBinom[N][P]);
261 BinomAllocator (const BinomAllocator&);
262 BinomAllocator& operator= (const BinomAllocator&);
265 Standard_Integer** myBinom;
266 Standard_Integer myMaxBinom;
270 // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
271 //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
272 static BinomAllocator THE_BINOM (25 + 1);
275 //=======================================================================
278 //=======================================================================
280 Standard_Real PLib::Bin(const Standard_Integer N,
281 const Standard_Integer P)
283 return THE_BINOM.Value (N, P);
286 //=======================================================================
287 //function : RationalDerivative
289 //=======================================================================
291 void PLib::RationalDerivative(const Standard_Integer Degree,
292 const Standard_Integer DerivativeRequest,
293 const Standard_Integer Dimension,
295 Standard_Real& RDers,
296 const Standard_Boolean All)
299 // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
302 // Let C(N,P) be the binomial
307 // u = SUM C (q,p) f v
314 // (q) ( (q) (p) (q-p) )
315 // f = (1/v) ( u - SUM C (q,p) f v )
319 // make arrays for the binomial since computing it each time could raise a performance
321 // As oppose to the method below the <Der> array is organized in the following
324 // u (1) u (2) .... u (Dimension) v (1)
327 // u (1) u (2) .... u (Dimension) v (1)
329 // ............................................
331 // (Degree) (Degree) (Degree) (Degree)
332 // u (1) u (2) .... u (Dimension) v (1)
335 Standard_Real Inverse;
336 Standard_Real *PolesArray = &Ders;
337 Standard_Real *RationalArray = &RDers;
338 Standard_Real Factor ;
339 Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
340 NCollection_LocalArray<Standard_Real> binomial_array;
341 NCollection_LocalArray<Standard_Real> derivative_storage;
342 if (Dimension == 3) {
343 Standard_Integer DeRequest1 = DerivativeRequest + 1;
344 Standard_Integer MinDegRequ = DerivativeRequest;
345 if (MinDegRequ > Degree) MinDegRequ = Degree;
346 binomial_array.Allocate (DeRequest1);
348 for (ii = 0 ; ii < DeRequest1 ; ii++) {
349 binomial_array[ii] = 1.0e0 ;
352 Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
353 derivative_storage.Allocate (DimDeRequ1);
354 RationalArray = derivative_storage ;
357 Inverse = 1.0e0 / PolesArray[3] ;
362 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
365 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
366 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
367 RationalArray[Index] = PolesArray[OtherIndex];
371 for (jj = ii - 1 ; jj >= 0 ; jj--) {
372 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
373 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
374 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
375 RationalArray[Index] -= Factor * RationalArray[Index1];
380 for (jj = ii ; jj >= 1 ; jj--) {
381 binomial_array[jj] += binomial_array[jj - 1] ;
383 RationalArray[Index] *= Inverse ; Index++;
384 RationalArray[Index] *= Inverse ; Index++;
385 RationalArray[Index] *= Inverse ; Index++;
388 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
391 RationalArray[Index] = 0.0e0; Index++;
392 RationalArray[Index] = 0.0e0; Index++;
393 RationalArray[Index] = 0.0e0;
396 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
397 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
398 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
399 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
400 RationalArray[Index] -= Factor * RationalArray[Index1];
405 for (jj = ii ; jj >= 1 ; jj--) {
406 binomial_array[jj] += binomial_array[jj - 1] ;
408 RationalArray[Index] *= Inverse; Index++;
409 RationalArray[Index] *= Inverse; Index++;
410 RationalArray[Index] *= Inverse; Index++;
414 RationalArray = &RDers ;
415 Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
416 RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
417 RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
418 RationalArray[2] = derivative_storage[DimDeRequ];
423 Standard_Integer Dimension1 = Dimension + 1;
424 Standard_Integer Dimension2 = Dimension << 1;
425 Standard_Integer DeRequest1 = DerivativeRequest + 1;
426 Standard_Integer MinDegRequ = DerivativeRequest;
427 if (MinDegRequ > Degree) MinDegRequ = Degree;
428 binomial_array.Allocate (DeRequest1);
430 for (ii = 0 ; ii < DeRequest1 ; ii++) {
431 binomial_array[ii] = 1.0e0 ;
434 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
435 derivative_storage.Allocate (DimDeRequ1);
436 RationalArray = derivative_storage ;
439 Inverse = 1.0e0 / PolesArray[Dimension] ;
441 Index2 = - Dimension2;
444 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
448 for (kk = 0 ; kk < Dimension ; kk++) {
449 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
454 for (jj = ii - 1 ; jj >= 0 ; jj--) {
455 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
457 for (kk = 0 ; kk < Dimension ; kk++) {
458 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
461 Index1 -= Dimension2 ;
464 for (jj = ii ; jj >= 1 ; jj--) {
465 binomial_array[jj] += binomial_array[jj - 1] ;
468 for (kk = 0 ; kk < Dimension ; kk++) {
469 RationalArray[Index] *= Inverse ; Index++;
473 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
477 for (kk = 0 ; kk < Dimension ; kk++) {
478 RationalArray[Index] = 0.0e0 ; Index++;
482 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
483 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
485 for (kk = 0 ; kk < Dimension ; kk++) {
486 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
489 Index1 -= Dimension2 ;
492 for (jj = ii ; jj >= 1 ; jj--) {
493 binomial_array[jj] += binomial_array[jj - 1] ;
496 for (kk = 0 ; kk < Dimension ; kk++) {
497 RationalArray[Index] *= Inverse; Index++;
502 RationalArray = &RDers ;
503 Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
505 for (kk = 0 ; kk < Dimension ; kk++) {
506 RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
512 //=======================================================================
513 //function : RationalDerivatives
514 //purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights
515 //=======================================================================
517 void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
518 const Standard_Integer Dimension,
519 Standard_Real& PolesDerivates,
520 // must be an array with
521 // (DerivativeRequest + 1) * Dimension slots
522 Standard_Real& WeightsDerivates,
523 // must be an array with
524 // (DerivativeRequest + 1) slots
525 Standard_Real& RationalDerivates)
528 // Our purpose is to compute f = (u/v) derivated N times
531 // Let C(N,P) be the binomial
536 // u = SUM C (q,p) f v
543 // (q) ( (q) (p) (q-p) )
544 // f = (1/v) ( u - SUM C (q,p) f v )
548 // make arrays for the binomial since computing it each time could
549 // raize a performance issue
551 Standard_Real Inverse;
552 Standard_Real *PolesArray = &PolesDerivates;
553 Standard_Real *WeightsArray = &WeightsDerivates;
554 Standard_Real *RationalArray = &RationalDerivates;
555 Standard_Real Factor ;
557 Standard_Integer ii, Index, Index1, Index2, jj;
558 Standard_Integer DeRequest1 = DerivativeRequest + 1;
560 NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
561 NCollection_LocalArray<Standard_Real> derivative_storage;
563 for (ii = 0 ; ii < DeRequest1 ; ii++) {
564 binomial_array[ii] = 1.0e0 ;
566 Inverse = 1.0e0 / WeightsArray[0] ;
567 if (Dimension == 3) {
571 for (ii = 0 ; ii < DeRequest1 ; ii++) {
574 RationalArray[Index] = PolesArray[Index] ; Index++;
575 RationalArray[Index] = PolesArray[Index] ; Index++;
576 RationalArray[Index] = PolesArray[Index] ;
579 for (jj = ii - 1 ; jj >= 0 ; jj--) {
580 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
581 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
582 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
583 RationalArray[Index] -= Factor * RationalArray[Index1];
588 for (jj = ii ; jj >= 1 ; jj--) {
589 binomial_array[jj] += binomial_array[jj - 1] ;
591 RationalArray[Index] *= Inverse ; Index++;
592 RationalArray[Index] *= Inverse ; Index++;
593 RationalArray[Index] *= Inverse ; Index++;
598 Standard_Integer Dimension2 = Dimension << 1;
600 Index2 = - Dimension2;
602 for (ii = 0 ; ii < DeRequest1 ; ii++) {
606 for (kk = 0 ; kk < Dimension ; kk++) {
607 RationalArray[Index] = PolesArray[Index]; Index++;
611 for (jj = ii - 1 ; jj >= 0 ; jj--) {
612 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
614 for (kk = 0 ; kk < Dimension ; kk++) {
615 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
618 Index1 -= Dimension2;
621 for (jj = ii ; jj >= 1 ; jj--) {
622 binomial_array[jj] += binomial_array[jj - 1] ;
625 for (kk = 0 ; kk < Dimension ; kk++) {
626 RationalArray[Index] *= Inverse ; Index++;
632 //=======================================================================
633 // Auxiliary template functions used for optimized evaluation of polynome
634 // and its derivatives for smaller dimensions of the polynome
635 //=======================================================================
638 // recursive template for evaluating value or first derivative
640 inline void eval_step1 (double* poly, double par, double* coef)
642 eval_step1<dim - 1> (poly, par, coef);
643 poly[dim] = poly[dim] * par + coef[dim];
648 inline void eval_step1<-1> (double*, double, double*)
652 // recursive template for evaluating second derivative
654 inline void eval_step2 (double* poly, double par, double* coef)
656 eval_step2<dim - 1> (poly, par, coef);
657 poly[dim] = poly[dim] * par + coef[dim] * 2.;
662 inline void eval_step2<-1> (double*, double, double*)
666 // evaluation of only value
668 inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
670 Standard_Real* aRes0 = aRes;
671 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
673 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
676 // Calculating the value of the polynomial
677 eval_step1<dim-1> (aRes0, Par, aCoeffs);
681 // evaluation of value and first derivative
683 inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
685 Standard_Real* aRes0 = aRes;
686 Standard_Real* aRes1 = aRes + dim;
688 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
689 memset(aRes1, 0, sizeof(Standard_Real) * dim);
691 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
694 // Calculating derivatives of the polynomial
695 eval_step1<dim-1> (aRes1, Par, aRes0);
696 // Calculating the value of the polynomial
697 eval_step1<dim-1> (aRes0, Par, aCoeffs);
701 // evaluation of value and first and second derivatives
703 inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
705 Standard_Real* aRes0 = aRes;
706 Standard_Real* aRes1 = aRes + dim;
707 Standard_Real* aRes2 = aRes + 2 * dim;
709 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
710 memset(aRes1, 0, sizeof(Standard_Real) * dim);
711 memset(aRes2, 0, sizeof(Standard_Real) * dim);
713 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
716 // Calculating second derivatives of the polynomial
717 eval_step2<dim-1> (aRes2, Par, aRes1);
718 // Calculating derivatives of the polynomial
719 eval_step1<dim-1> (aRes1, Par, aRes0);
720 // Calculating the value of the polynomial
721 eval_step1<dim-1> (aRes0, Par, aCoeffs);
726 //=======================================================================
727 //function : This evaluates a polynomial and its derivatives
728 //purpose : up to the requested order
729 //=======================================================================
731 void PLib::EvalPolynomial(const Standard_Real Par,
732 const Standard_Integer DerivativeRequest,
733 const Standard_Integer Degree,
734 const Standard_Integer Dimension,
735 Standard_Real& PolynomialCoeff,
736 Standard_Real& Results)
738 // the polynomial coefficients are assumed to be stored as follows :
740 // [0] [Dimension -1] X coefficient
742 // [Dimension] [Dimension + Dimension -1] X coefficient
744 // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient
746 // ...................................................
750 // [d * Dimension] [d * Dimension + Dimension-1] X coefficient
752 // where d is the Degree
755 Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
756 Standard_Real* aRes = &Results;
757 Standard_Real* anOriginal;
758 Standard_Integer ind = 0;
759 switch (DerivativeRequest)
765 case 1: eval_poly1<1> (aRes, aCoeffs, Degree, Par); break;
766 case 2: eval_poly1<2> (aRes, aCoeffs, Degree, Par); break;
767 case 3: eval_poly1<3> (aRes, aCoeffs, Degree, Par); break;
768 case 4: eval_poly1<4> (aRes, aCoeffs, Degree, Par); break;
769 case 5: eval_poly1<5> (aRes, aCoeffs, Degree, Par); break;
770 case 6: eval_poly1<6> (aRes, aCoeffs, Degree, Par); break;
771 case 7: eval_poly1<7> (aRes, aCoeffs, Degree, Par); break;
772 case 8: eval_poly1<8> (aRes, aCoeffs, Degree, Par); break;
773 case 9: eval_poly1<9> (aRes, aCoeffs, Degree, Par); break;
774 case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
775 case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
776 case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
777 case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
778 case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
779 case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
782 Standard_Real* aRes0 = aRes;
783 Standard_Real* aRes1 = aRes + Dimension;
785 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
786 memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
788 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
790 aCoeffs -= Dimension;
791 // Calculating derivatives of the polynomial
792 for (ind = 0; ind < Dimension; ind++)
793 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
794 // Calculating the value of the polynomial
795 for (ind = 0; ind < Dimension; ind++)
796 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
806 case 1: eval_poly2<1> (aRes, aCoeffs, Degree, Par); break;
807 case 2: eval_poly2<2> (aRes, aCoeffs, Degree, Par); break;
808 case 3: eval_poly2<3> (aRes, aCoeffs, Degree, Par); break;
809 case 4: eval_poly2<4> (aRes, aCoeffs, Degree, Par); break;
810 case 5: eval_poly2<5> (aRes, aCoeffs, Degree, Par); break;
811 case 6: eval_poly2<6> (aRes, aCoeffs, Degree, Par); break;
812 case 7: eval_poly2<7> (aRes, aCoeffs, Degree, Par); break;
813 case 8: eval_poly2<8> (aRes, aCoeffs, Degree, Par); break;
814 case 9: eval_poly2<9> (aRes, aCoeffs, Degree, Par); break;
815 case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
816 case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
817 case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
818 case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
819 case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
820 case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
823 Standard_Real* aRes0 = aRes;
824 Standard_Real* aRes1 = aRes + Dimension;
825 Standard_Real* aRes2 = aRes1 + Dimension;
827 // Nullify the results
828 Standard_Integer aSize = 2 * Dimension;
829 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
830 memset(aRes1, 0, sizeof(Standard_Real) * aSize);
832 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
834 aCoeffs -= Dimension;
835 // Calculating derivatives of the polynomial
836 for (ind = 0; ind < Dimension; ind++)
837 aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
838 for (ind = 0; ind < Dimension; ind++)
839 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
840 // Calculating the value of the polynomial
841 for (ind = 0; ind < Dimension; ind++)
842 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
851 // Nullify the results
852 Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
853 memset(aRes, 0, sizeof(Standard_Real) * aResSize);
855 for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
857 aRes = &Results + aResSize - Dimension;
858 // Calculating derivatives of the polynomial
859 for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
861 anOriginal = aRes - Dimension; // pointer to the derivative minus 1
862 for (ind = 0; ind < Dimension; ind++)
863 aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
866 // Calculating the value of the polynomial
867 for (ind = 0; ind < Dimension; ind++)
868 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
869 aCoeffs -= Dimension;
875 //=======================================================================
876 //function : This evaluates a polynomial without derivative
878 //=======================================================================
880 void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par,
881 const Standard_Integer Degree,
882 const Standard_Integer Dimension,
883 const Standard_Integer DegreeDimension,
884 Standard_Real& PolynomialCoeff,
885 Standard_Real& Results)
887 Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
888 Standard_Real* aRes = &Results;
892 case 1: eval_poly0<1> (aRes, aCoeffs, Degree, Par); break;
893 case 2: eval_poly0<2> (aRes, aCoeffs, Degree, Par); break;
894 case 3: eval_poly0<3> (aRes, aCoeffs, Degree, Par); break;
895 case 4: eval_poly0<4> (aRes, aCoeffs, Degree, Par); break;
896 case 5: eval_poly0<5> (aRes, aCoeffs, Degree, Par); break;
897 case 6: eval_poly0<6> (aRes, aCoeffs, Degree, Par); break;
898 case 7: eval_poly0<7> (aRes, aCoeffs, Degree, Par); break;
899 case 8: eval_poly0<8> (aRes, aCoeffs, Degree, Par); break;
900 case 9: eval_poly0<9> (aRes, aCoeffs, Degree, Par); break;
901 case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
902 case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
903 case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
904 case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
905 case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
906 case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
909 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
910 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
912 aCoeffs -= Dimension;
913 for (Standard_Integer ind = 0; ind < Dimension; ind++)
914 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
920 //=======================================================================
921 //function : This evaluates a polynomial of 2 variables
922 //purpose : or its derivative at the requested orders
923 //=======================================================================
925 void PLib::EvalPoly2Var(const Standard_Real UParameter,
926 const Standard_Real VParameter,
927 const Standard_Integer UDerivativeRequest,
928 const Standard_Integer VDerivativeRequest,
929 const Standard_Integer UDegree,
930 const Standard_Integer VDegree,
931 const Standard_Integer Dimension,
932 Standard_Real& PolynomialCoeff,
933 Standard_Real& Results)
935 // the polynomial coefficients are assumed to be stored as follows :
937 // [0] [Dimension -1] U V coefficient
939 // [Dimension] [Dimension + Dimension -1] U V coefficient
941 // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient
943 // ...................................................
947 // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient
952 // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient
954 // ...................................................
957 // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient
959 // ...................................................
962 // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient
966 Standard_Integer Udim = (VDegree+1)*Dimension,
967 index = Udim*UDerivativeRequest;
968 TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
969 TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
970 Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1);
971 Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1);
972 Standard_Real * ResultArray ;
973 ResultArray = &Results ;
975 PLib::EvalPolynomial(UParameter,
982 PLib::EvalPolynomial(VParameter,
989 index = Dimension*VDerivativeRequest;
991 for (Standard_Integer i=0;i<Dimension;i++) {
992 ResultArray[i] = Digit[index+i];
998 //=======================================================================
999 //function : This evaluates the lagrange polynomial and its derivatives
1000 //purpose : up to the requested order that interpolates a series of
1001 //points of dimension <Dimension> with given assigned parameters
1002 //=======================================================================
1005 PLib::EvalLagrange(const Standard_Real Parameter,
1006 const Standard_Integer DerivativeRequest,
1007 const Standard_Integer Degree,
1008 const Standard_Integer Dimension,
1009 Standard_Real& Values,
1010 Standard_Real& Parameters,
1011 Standard_Real& Results)
1014 // the points are assumed to be stored as follows in the Values array :
1016 // [0] [Dimension -1] first point coefficients
1018 // [Dimension] [Dimension + Dimension -1] second point coefficients
1020 // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients
1022 // ...................................................
1026 // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients
1028 // where d is the Degree
1030 // The ParameterArray stores the parameter value assign to each point in
1031 // order described above, that is
1032 // [0] is assign to first point
1033 // [1] is assign to second point
1035 Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1036 Standard_Integer local_request = DerivativeRequest;
1037 Standard_Real *ParameterArray;
1038 Standard_Real difference;
1039 Standard_Real *PointsArray;
1040 Standard_Real *ResultArray ;
1042 PointsArray = &Values ;
1043 ParameterArray = &Parameters ;
1044 ResultArray = &Results ;
1045 if (local_request >= Degree) {
1046 local_request = Degree ;
1048 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1050 // Build the divided differences array
1053 for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) {
1054 divided_differences_array[ii] = PointsArray[ii] ;
1057 for (ii = Degree ; ii >= 0 ; ii--) {
1059 for (jj = Degree ; jj > Degree - ii ; jj--) {
1060 Index = jj * Dimension ;
1061 Index1 = Index - Dimension ;
1063 for (kk = 0 ; kk < Dimension ; kk++) {
1064 divided_differences_array[Index + kk] -=
1065 divided_differences_array[Index1 + kk] ;
1068 ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ;
1069 if (Abs(difference) < RealSmall()) {
1073 difference = 1.0e0 / difference ;
1075 for (kk = 0 ; kk < Dimension ; kk++) {
1076 divided_differences_array[Index + kk] *= difference ;
1082 // Evaluate the divided difference array polynomial which expresses as
1084 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1085 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1087 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1090 Index = Degree * Dimension ;
1092 for (kk = 0 ; kk < Dimension ; kk++) {
1093 ResultArray[kk] = divided_differences_array[Index + kk] ;
1096 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1097 ResultArray[ii] = 0.0e0 ;
1100 for (ii = Degree ; ii >= 1 ; ii--) {
1101 difference = Parameter - ParameterArray[ii - 1] ;
1103 for (jj = local_request ; jj > 0 ; jj--) {
1104 Index = jj * Dimension ;
1105 Index1 = Index - Dimension ;
1107 for (kk = 0 ; kk < Dimension ; kk++) {
1108 ResultArray[Index + kk] *= difference ;
1109 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
1112 Index = (ii -1) * Dimension ;
1114 for (kk = 0 ; kk < Dimension ; kk++) {
1115 ResultArray[kk] *= difference ;
1116 ResultArray[kk] += divided_differences_array[Index+kk] ;
1120 return (ReturnCode) ;
1123 //=======================================================================
1124 //function : This evaluates the hermite polynomial and its derivatives
1125 //purpose : up to the requested order that interpolates a series of
1126 //points of dimension <Dimension> with given assigned parameters
1127 //=======================================================================
1129 Standard_Integer PLib::EvalCubicHermite
1130 (const Standard_Real Parameter,
1131 const Standard_Integer DerivativeRequest,
1132 const Standard_Integer Dimension,
1133 Standard_Real& Values,
1134 Standard_Real& Derivatives,
1135 Standard_Real& theParameters,
1136 Standard_Real& Results)
1139 // the points are assumed to be stored as follows in the Values array :
1141 // [0] [Dimension -1] first point coefficients
1143 // [Dimension] [Dimension + Dimension -1] last point coefficients
1146 // the derivatives are assumed to be stored as follows in
1147 // the Derivatives array :
1149 // [0] [Dimension -1] first point coefficients
1151 // [Dimension] [Dimension + Dimension -1] last point coefficients
1153 // The ParameterArray stores the parameter value assign to each point in
1154 // order described above, that is
1155 // [0] is assign to first point
1156 // [1] is assign to last point
1158 Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
1159 Standard_Integer local_request = DerivativeRequest ;
1163 Standard_Real ParametersArray[4];
1164 Standard_Real difference;
1165 Standard_Real inverse;
1166 Standard_Real *FirstLast;
1167 Standard_Real *PointsArray;
1168 Standard_Real *DerivativesArray;
1169 Standard_Real *ResultArray ;
1171 DerivativesArray = &Derivatives ;
1172 PointsArray = &Values ;
1173 FirstLast = &theParameters ;
1174 ResultArray = &Results ;
1175 if (local_request >= Degree) {
1176 local_request = Degree ;
1178 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
1180 for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) {
1181 ParametersArray[jj] =
1182 ParametersArray[jj+1] = FirstLast[ii] ;
1185 // Build the divided differences array
1188 // initialise it at the stage 2 of the building algorithm
1189 // for divided differences
1191 inverse = FirstLast[1] - FirstLast[0] ;
1192 inverse = 1.0e0 / inverse ;
1194 for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ;
1196 ii++, jj++, kk++, pp++) {
1197 divided_differences_array[ii] = PointsArray[ii] ;
1198 divided_differences_array[kk] = inverse *
1199 (PointsArray[jj] - PointsArray[ii]) ;
1200 divided_differences_array[jj] = DerivativesArray[ii] ;
1201 divided_differences_array[pp] = DerivativesArray[jj] ;
1204 for (ii = 1 ; ii <= Degree ; ii++) {
1206 for (jj = Degree ; jj >= ii+1 ; jj--) {
1207 Index = jj * Dimension ;
1208 Index1 = Index - Dimension ;
1210 for (kk = 0 ; kk < Dimension ; kk++) {
1211 divided_differences_array[Index + kk] -=
1212 divided_differences_array[Index1 + kk] ;
1215 for (kk = 0 ; kk < Dimension ; kk++) {
1216 divided_differences_array[Index + kk] *= inverse ;
1222 // Evaluate the divided difference array polynomial which expresses as
1224 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1225 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1227 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1230 Index = Degree * Dimension ;
1232 for (kk = 0 ; kk < Dimension ; kk++) {
1233 ResultArray[kk] = divided_differences_array[Index + kk] ;
1236 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1237 ResultArray[ii] = 0.0e0 ;
1240 for (ii = Degree ; ii >= 1 ; ii--) {
1241 difference = Parameter - ParametersArray[ii - 1] ;
1243 for (jj = local_request ; jj > 0 ; jj--) {
1244 Index = jj * Dimension ;
1245 Index1 = Index - Dimension ;
1247 for (kk = 0 ; kk < Dimension ; kk++) {
1248 ResultArray[Index + kk] *= difference ;
1249 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
1252 Index = (ii -1) * Dimension ;
1254 for (kk = 0 ; kk < Dimension ; kk++) {
1255 ResultArray[kk] *= difference ;
1256 ResultArray[kk] += divided_differences_array[Index+kk] ;
1260 return (ReturnCode) ;
1263 //=======================================================================
1264 //function : HermiteCoefficients
1265 //purpose : calcul des polynomes d'Hermite
1266 //=======================================================================
1268 Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
1269 const Standard_Real LastParameter,
1270 const Standard_Integer FirstOrder,
1271 const Standard_Integer LastOrder,
1272 math_Matrix& MatrixCoefs)
1274 Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2];
1275 Standard_Integer ii, jj, pp, cote, iof=0;
1276 Standard_Real Prod, TBorne = FirstParameter;
1277 math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
1278 math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
1280 // Test de validites
1282 if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
1283 Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
1284 if (D1 > 100 || D2 > 100) return Standard_False;
1286 if (D2 < 0.01) return Standard_False;
1287 if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
1289 // Calcul de la matrice a inverser (MAT)
1291 Ordre[0] = FirstOrder+1;
1292 Ordre[1] = LastOrder+1;
1294 for (cote=0; cote<=1; cote++) {
1297 for (pp=1; pp<=Ordre[cote]; pp++) {
1301 for (jj=pp; jj<=NbCoeff; jj++) {
1302 // tout se passe dans les 3 lignes suivantes
1303 MAT(ii, jj) = Coeff(jj) * Prod;
1304 Coeff(jj) *= jj - pp;
1308 TBorne = LastParameter;
1312 // resolution du systemes
1313 math_Gauss ResolCoeff(MAT, 1.0e-10);
1314 if (!ResolCoeff.IsDone()) return Standard_False;
1316 for (ii=1; ii<=NbCoeff; ii++) {
1318 ResolCoeff.Solve(B, Coeff);
1319 MatrixCoefs.SetRow( ii, Coeff);
1322 return Standard_True;
1325 //=======================================================================
1326 //function : CoefficientsPoles
1328 //=======================================================================
1330 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs,
1331 const TColStd_Array1OfReal* WCoefs,
1332 TColgp_Array1OfPnt& Poles,
1333 TColStd_Array1OfReal* Weights)
1335 TColStd_Array1OfReal tempC(1,3*Coefs.Length());
1336 PLib::SetPoles(Coefs,tempC);
1337 TColStd_Array1OfReal tempP(1,3*Poles.Length());
1338 PLib::SetPoles(Coefs,tempP);
1339 PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
1340 PLib::GetPoles(tempP,Poles);
1343 //=======================================================================
1344 //function : CoefficientsPoles
1346 //=======================================================================
1348 void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
1349 const TColStd_Array1OfReal* WCoefs,
1350 TColgp_Array1OfPnt2d& Poles,
1351 TColStd_Array1OfReal* Weights)
1353 TColStd_Array1OfReal tempC(1,2*Coefs.Length());
1354 PLib::SetPoles(Coefs,tempC);
1355 TColStd_Array1OfReal tempP(1,2*Poles.Length());
1356 PLib::SetPoles(Coefs,tempP);
1357 PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
1358 PLib::GetPoles(tempP,Poles);
1361 //=======================================================================
1362 //function : CoefficientsPoles
1364 //=======================================================================
1366 void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
1367 const TColStd_Array1OfReal* WCoefs,
1368 TColStd_Array1OfReal& Poles,
1369 TColStd_Array1OfReal* Weights)
1371 PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
1374 //=======================================================================
1375 //function : CoefficientsPoles
1377 //=======================================================================
1379 void PLib::CoefficientsPoles (const Standard_Integer dim,
1380 const TColStd_Array1OfReal& Coefs,
1381 const TColStd_Array1OfReal* WCoefs,
1382 TColStd_Array1OfReal& Poles,
1383 TColStd_Array1OfReal* Weights)
1385 Standard_Boolean rat = WCoefs != NULL;
1386 Standard_Integer loc = Coefs.Lower();
1387 Standard_Integer lop = Poles.Lower();
1388 Standard_Integer lowc=0;
1389 Standard_Integer lowp=0;
1390 Standard_Integer upc = Coefs.Upper();
1391 Standard_Integer upp = Poles.Upper();
1392 Standard_Integer upwc=0;
1393 Standard_Integer upwp=0;
1394 Standard_Integer reflen = Coefs.Length()/dim;
1395 Standard_Integer i,j,k;
1398 lowc = WCoefs->Lower(); lowp = Weights->Lower();
1399 upwc = WCoefs->Upper(); upwp = Weights->Upper();
1402 for (i = 0; i < dim; i++){
1403 Poles (lop + i) = Coefs (loc + i);
1404 Poles (upp - i) = Coefs (upc - i);
1407 (*Weights) (lowp) = (*WCoefs) (lowc);
1408 (*Weights) (upwp) = (*WCoefs) (upwc);
1412 for (i = 2; i < reflen; i++ ) {
1413 Cnp = PLib::Bin(reflen - 1, i - 1);
1414 if (rat) (*Weights)(lowp + i - 1) = (*WCoefs)(lowc + i - 1) / Cnp;
1416 for(j = 0; j < dim; j++){
1417 Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
1421 for (i = 1; i <= reflen - 1; i++) {
1423 for (j = reflen - 1; j >= i; j--) {
1424 if (rat) (*Weights)(lowp + j) += (*Weights)(lowp + j - 1);
1426 for(k = 0; k < dim; k++){
1427 Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
1433 for (i = 1; i <= reflen; i++) {
1435 for(j = 0; j < dim; j++){
1436 Poles(lop + dim * (i-1) + j) /= (*Weights)(lowp + i -1);
1442 //=======================================================================
1443 //function : Trimming
1445 //=======================================================================
1447 void PLib::Trimming(const Standard_Real U1,
1448 const Standard_Real U2,
1449 TColgp_Array1OfPnt& Coefs,
1450 TColStd_Array1OfReal* WCoefs)
1452 TColStd_Array1OfReal temp(1,3*Coefs.Length());
1453 PLib::SetPoles(Coefs,temp);
1454 PLib::Trimming(U1,U2,3,temp,WCoefs);
1455 PLib::GetPoles(temp,Coefs);
1458 //=======================================================================
1459 //function : Trimming
1461 //=======================================================================
1463 void PLib::Trimming(const Standard_Real U1,
1464 const Standard_Real U2,
1465 TColgp_Array1OfPnt2d& Coefs,
1466 TColStd_Array1OfReal* WCoefs)
1468 TColStd_Array1OfReal temp(1,2*Coefs.Length());
1469 PLib::SetPoles(Coefs,temp);
1470 PLib::Trimming(U1,U2,2,temp,WCoefs);
1471 PLib::GetPoles(temp,Coefs);
1474 //=======================================================================
1475 //function : Trimming
1477 //=======================================================================
1479 void PLib::Trimming(const Standard_Real U1,
1480 const Standard_Real U2,
1481 TColStd_Array1OfReal& Coefs,
1482 TColStd_Array1OfReal* WCoefs)
1484 PLib::Trimming(U1,U2,1,Coefs,WCoefs);
1487 //=======================================================================
1488 //function : Trimming
1490 //=======================================================================
1492 void PLib::Trimming(const Standard_Real U1,
1493 const Standard_Real U2,
1494 const Standard_Integer dim,
1495 TColStd_Array1OfReal& Coefs,
1496 TColStd_Array1OfReal* WCoefs)
1500 // on fait le changement de variable v = (u-U1) / (U2-U1)
1501 // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
1502 // decomposee sous la forme du schema iteratif de horner.
1504 Standard_Real lsp = U2 - U1;
1505 Standard_Integer indc, indw=0;
1506 Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
1507 Standard_Integer len = Coefs.Length()/dim;
1508 Standard_Boolean rat = WCoefs != NULL;
1511 if(len != WCoefs->Length())
1512 throw Standard_Failure("PLib::Trimming : nbcoefs/dim != nbweights !!!");
1513 upw = WCoefs->Upper();
1517 for (Standard_Integer i = 1; i <= len; i++) {
1518 Standard_Integer j ;
1519 indc = upc - dim*(i-1);
1520 if (rat) indw = upw - i + 1;
1521 //calcul du coefficient de degre le plus faible a l'iteration i
1523 for( j = 0; j < dim; j++){
1524 Coefs(indc - dim + j) += U1 * Coefs(indc + j);
1526 if (rat) (*WCoefs)(indw - 1) += U1 * (*WCoefs)(indw);
1528 //calcul des coefficients intermediaires :
1533 for(Standard_Integer k = 0; k < dim; k++){
1534 Coefs(indc - dim + k) =
1535 U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
1539 (*WCoefs)(indw - 1) = U1 * (*WCoefs)(indw) + lsp * (*WCoefs)(indw - 1);
1543 //calcul du coefficient de degre le plus eleve :
1545 for(j = 0; j < dim; j++){
1546 Coefs(upc + j) *= lsp;
1548 if (rat) (*WCoefs)(upw) *= lsp;
1552 //=======================================================================
1553 //function : CoefficientsPoles
1555 // Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852).
1556 // on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
1557 // fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
1558 //=======================================================================
1560 void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs,
1561 const TColStd_Array2OfReal* WCoefs,
1562 TColgp_Array2OfPnt& Poles,
1563 TColStd_Array2OfReal* Weights)
1565 Standard_Boolean rat = (WCoefs != NULL);
1566 Standard_Integer LowerRow = Poles.LowerRow();
1567 Standard_Integer UpperRow = Poles.UpperRow();
1568 Standard_Integer LowerCol = Poles.LowerCol();
1569 Standard_Integer UpperCol = Poles.UpperCol();
1570 Standard_Integer ColLength = Poles.ColLength();
1571 Standard_Integer RowLength = Poles.RowLength();
1573 // Bidouille pour retablir u et v pour les coefs calcules
1575 // Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
1577 Standard_Integer Row, Col;
1578 Standard_Real W, Cnp;
1580 Standard_Integer I1, I2;
1581 Standard_Integer NPoleu , NPolev;
1584 for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
1585 Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol);
1587 (*Weights) (NPoleu, LowerCol) = (*WCoefs) (NPoleu, LowerCol);
1590 for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
1591 Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
1592 Temp = Coefs (NPoleu, Col).XYZ();
1594 Poles (NPoleu, Col).SetXYZ (Temp);
1596 (*Weights) (NPoleu, Col) = (*WCoefs) (NPoleu, Col) / Cnp;
1599 Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
1601 (*Weights) (NPoleu, UpperCol) = (*WCoefs) (NPoleu, UpperCol);
1604 for (I1 = 1; I1 <= RowLength - 1; I1++) {
1606 for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
1608 (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
1609 Poles (NPoleu, I2).SetXYZ (Temp);
1610 if (rat) (*Weights)(NPoleu, I2) += (*Weights)(NPoleu, I2-1);
1615 for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
1617 for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
1618 Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
1619 Temp = Poles (Row, NPolev).XYZ();
1621 Poles (Row, NPolev).SetXYZ (Temp);
1622 if (rat) (*Weights)(Row, NPolev) /= Cnp;
1625 for (I1 = 1; I1 <= ColLength - 1; I1++) {
1627 for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
1629 (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
1630 Poles (I2, NPolev).SetXYZ (Temp);
1631 if (rat) (*Weights)(I2, NPolev) += (*Weights)(I2-1, NPolev);
1637 for (Row = LowerRow; Row <= UpperRow; Row++) {
1639 for (Col = LowerCol; Col <= UpperCol; Col++) {
1640 W = (*Weights) (Row, Col);
1641 Temp = Poles(Row, Col).XYZ();
1643 Poles(Row, Col).SetXYZ (Temp);
1649 //=======================================================================
1650 //function : UTrimming
1652 //=======================================================================
1654 void PLib::UTrimming(const Standard_Real U1,
1655 const Standard_Real U2,
1656 TColgp_Array2OfPnt& Coeffs,
1657 TColStd_Array2OfReal* WCoeffs)
1659 Standard_Boolean rat = WCoeffs != NULL;
1660 Standard_Integer lr = Coeffs.LowerRow();
1661 Standard_Integer ur = Coeffs.UpperRow();
1662 Standard_Integer lc = Coeffs.LowerCol();
1663 Standard_Integer uc = Coeffs.UpperCol();
1664 TColgp_Array1OfPnt Temp (lr,ur);
1665 TColStd_Array1OfReal Temw (lr,ur);
1667 for (Standard_Integer icol = lc; icol <= uc; icol++) {
1668 Standard_Integer irow ;
1669 for ( irow = lr; irow <= ur; irow++) {
1670 Temp (irow) = Coeffs (irow, icol);
1671 if (rat) Temw (irow) = (*WCoeffs) (irow, icol);
1673 if (rat) PLib::Trimming (U1, U2, Temp, &Temw);
1674 else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
1676 for (irow = lr; irow <= ur; irow++) {
1677 Coeffs (irow, icol) = Temp (irow);
1678 if (rat) (*WCoeffs) (irow, icol) = Temw (irow);
1683 //=======================================================================
1684 //function : VTrimming
1686 //=======================================================================
1688 void PLib::VTrimming(const Standard_Real V1,
1689 const Standard_Real V2,
1690 TColgp_Array2OfPnt& Coeffs,
1691 TColStd_Array2OfReal* WCoeffs)
1693 Standard_Boolean rat = WCoeffs != NULL;
1694 Standard_Integer lr = Coeffs.LowerRow();
1695 Standard_Integer ur = Coeffs.UpperRow();
1696 Standard_Integer lc = Coeffs.LowerCol();
1697 Standard_Integer uc = Coeffs.UpperCol();
1698 TColgp_Array1OfPnt Temp (lc,uc);
1699 TColStd_Array1OfReal Temw (lc,uc);
1701 for (Standard_Integer irow = lr; irow <= ur; irow++) {
1702 Standard_Integer icol ;
1703 for ( icol = lc; icol <= uc; icol++) {
1704 Temp (icol) = Coeffs (irow, icol);
1705 if (rat) Temw (icol) = (*WCoeffs) (irow, icol);
1707 if (rat) PLib::Trimming (V1, V2, Temp, &Temw);
1708 else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
1710 for (icol = lc; icol <= uc; icol++) {
1711 Coeffs (irow, icol) = Temp (icol);
1712 if (rat) (*WCoeffs) (irow, icol) = Temw (icol);
1717 //=======================================================================
1718 //function : HermiteInterpolate
1720 //=======================================================================
1722 Standard_Boolean PLib::HermiteInterpolate
1723 (const Standard_Integer Dimension,
1724 const Standard_Real FirstParameter,
1725 const Standard_Real LastParameter,
1726 const Standard_Integer FirstOrder,
1727 const Standard_Integer LastOrder,
1728 const TColStd_Array2OfReal& FirstConstr,
1729 const TColStd_Array2OfReal& LastConstr,
1730 TColStd_Array1OfReal& Coefficients)
1732 Standard_Real Pattern[3][6];
1734 // portage HP : il faut les initialiser 1 par 1
1755 math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
1756 // The initialisation of the matrix A
1757 Standard_Integer irow ;
1758 for ( irow=0; irow<=FirstOrder; irow++) {
1759 Standard_Real FirstVal = 1.;
1761 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1762 A(irow,icol) = Pattern[irow][icol]*FirstVal;
1763 if (irow <= icol) FirstVal *= FirstParameter;
1767 for (irow=0; irow<=LastOrder; irow++) {
1768 Standard_Real LastVal = 1.;
1770 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1771 A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
1772 if (irow <= icol) LastVal *= LastParameter;
1776 // The filled matrix A for FirstOrder=LastOrder=2 is:
1778 // 1 FP FP**2 FP**3 FP**4 FP**5
1779 // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter
1780 // 0 0 2 6*FP 12*FP**2 20*FP**3
1781 // 1 LP LP**2 LP**3 LP**4 LP**5
1782 // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter
1783 // 0 0 2 6*LP 12*LP**2 20*LP**3
1785 // If FirstOrder or LastOrder <=2 then some rows and columns are missing.
1787 // If FirstOrder=1 then 3th row and 6th column are missing
1788 // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
1790 math_Gauss Equations(A);
1791 // std::cout << "A=" << A << std::endl;
1793 for (Standard_Integer idim=1; idim<=Dimension; idim++) {
1794 // std::cout << "idim=" << idim << std::endl;
1796 math_Vector B(0,FirstOrder+LastOrder+1);
1797 Standard_Integer icol ;
1798 for ( icol=0; icol<=FirstOrder; icol++)
1799 B(icol) = FirstConstr(idim,icol);
1801 for (icol=0; icol<=LastOrder; icol++)
1802 B(FirstOrder+1+icol) = LastConstr(idim,icol);
1803 // std::cout << "B=" << B << std::endl;
1805 // The solving of equations system A * X = B. Then B = X
1807 // std::cout << "After Solving" << std::endl << "B=" << B << std::endl;
1809 if (Equations.IsDone()==Standard_False) return Standard_False;
1811 // the filling of the Coefficients
1813 for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
1814 Coefficients(Dimension*icol+idim-1) = B(icol);
1816 return Standard_True;
1819 //=======================================================================
1820 //function : JacobiParameters
1822 //=======================================================================
1824 void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
1825 const Standard_Integer MaxDegree,
1826 const Standard_Integer Code,
1827 Standard_Integer& NbGaussPoints,
1828 Standard_Integer& WorkDegree)
1830 // ConstraintOrder: Ordre de contrainte aux extremites :
1831 // C0 = contraintes de passage aux bornes;
1832 // C1 = C0 + contraintes de derivees 1eres;
1833 // C2 = C1 + contraintes de derivees 2ndes.
1834 // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
1835 // d' approximation (doit etre superieur ou egal a
1836 // 2*NivConstr+2 et inferieur ou egal a 50).
1837 // Code: Code d' init. des parametres de discretisation.
1838 // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
1839 // = -5 Calcul tres rapide mais peu precis (8pts)
1840 // = -4 ' ' ' ' ' ' (10pts)
1841 // = -3 ' ' ' ' ' ' (15pts)
1842 // = -2 ' ' ' ' ' ' (20pts)
1843 // = -1 ' ' ' ' ' ' (25pts)
1844 // = 1 calcul rapide avec precision moyenne (30pts).
1845 // = 2 calcul rapide avec meilleure precision (40pts).
1846 // = 3 calcul un peu plus lent avec bonne precision (50 pts).
1847 // = 4 calcul lent avec la meilleure precision possible
1850 // The possible values of NbGaussPoints
1852 const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
1853 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
1855 Standard_Integer NivConstr=0;
1856 switch (ConstraintOrder) {
1857 case GeomAbs_C0: NivConstr = 0; break;
1858 case GeomAbs_C1: NivConstr = 1; break;
1859 case GeomAbs_C2: NivConstr = 2; break;
1861 throw Standard_ConstructionError("Invalid ConstraintOrder");
1863 if (MaxDegree < 2*NivConstr+1)
1864 throw Standard_ConstructionError("Invalid MaxDegree");
1867 WorkDegree = MaxDegree + 9;
1869 WorkDegree = MaxDegree + 6;
1871 //---> Nbre mini de points necessaires.
1872 Standard_Integer IPMIN=0;
1873 if (WorkDegree < NDEG8)
1875 else if (WorkDegree < NDEG10)
1877 else if (WorkDegree < NDEG15)
1879 else if (WorkDegree < NDEG20)
1881 else if (WorkDegree < NDEG25)
1883 else if (WorkDegree < NDEG30)
1885 else if (WorkDegree < NDEG40)
1887 else if (WorkDegree < NDEG50)
1889 else if (WorkDegree < NDEG61)
1892 throw Standard_ConstructionError("Invalid MaxDegree");
1893 // ---> Nbre de points voulus.
1894 Standard_Integer IWANT=0;
1896 case -5: IWANT=NDEG8; break;
1897 case -4: IWANT=NDEG10; break;
1898 case -3: IWANT=NDEG15; break;
1899 case -2: IWANT=NDEG20; break;
1900 case -1: IWANT=NDEG25; break;
1901 case 1: IWANT=NDEG30; break;
1902 case 2: IWANT=NDEG40; break;
1903 case 3: IWANT=NDEG50; break;
1904 case 4: IWANT=NDEG61; break;
1906 throw Standard_ConstructionError("Invalid Code");
1908 //--> NbGaussPoints est le nombre de points de discretisation de la fonction,
1909 // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
1910 // NbGaussPoints doit etre superieur strictement a WorkDegree.
1911 NbGaussPoints = Max(IPMIN,IWANT);
1912 // NbGaussPoints +=2;
1915 //=======================================================================
1916 //function : NivConstr
1917 //purpose : translates from GeomAbs_Shape to Integer
1918 //=======================================================================
1920 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
1922 Standard_Integer NivConstr=0;
1923 switch (ConstraintOrder) {
1924 case GeomAbs_C0: NivConstr = 0; break;
1925 case GeomAbs_C1: NivConstr = 1; break;
1926 case GeomAbs_C2: NivConstr = 2; break;
1928 throw Standard_ConstructionError("Invalid ConstraintOrder");
1933 //=======================================================================
1934 //function : ConstraintOrder
1935 //purpose : translates from Integer to GeomAbs_Shape
1936 //=======================================================================
1938 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
1940 GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
1941 switch (NivConstr) {
1942 case 0: ConstraintOrder = GeomAbs_C0; break;
1943 case 1: ConstraintOrder = GeomAbs_C1; break;
1944 case 2: ConstraintOrder = GeomAbs_C2; break;
1946 throw Standard_ConstructionError("Invalid NivConstr");
1948 return ConstraintOrder;
1951 //=======================================================================
1952 //function : EvalLength
1954 //=======================================================================
1956 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
1957 Standard_Real& PolynomialCoeff,
1958 const Standard_Real U1, const Standard_Real U2,
1959 Standard_Real& Length)
1961 Standard_Integer i,j,idim, degdim;
1962 Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
1964 Standard_Real *PolynomialArray = &PolynomialCoeff ;
1966 Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
1968 math_Vector GaussPoints(1,NbGaussPoints);
1969 math::GaussPoints(NbGaussPoints,GaussPoints);
1971 math_Vector GaussWeights(1,NbGaussPoints);
1972 math::GaussWeights(NbGaussPoints,GaussWeights);
1974 C1 = (U2 + U1) / 2.;
1975 C2 = (U2 - U1) / 2.;
1977 //-----------------------------------------------------------
1978 //****** Integration - Boucle sur les intervalles de GAUSS **
1979 //-----------------------------------------------------------
1983 for (j=1; j<=NbGaussPoints/2; j++) {
1984 // Integration en tenant compte de la symetrie
1985 Tran = C2 * GaussPoints(j);
1989 //****** Derivation sur la dimension de l'espace **
1991 degdim = Degree*Dimension;
1993 for (idim=0; idim<Dimension; idim++) {
1994 D1 = D2 = Degree * PolynomialArray [idim + degdim];
1995 for (i=Degree-1; i>=1; i--) {
1996 DD = i * PolynomialArray [idim + i*Dimension];
2004 //****** Integration **
2006 Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2008 //****** Fin de boucle dur les intervalles de GAUSS **
2014 //=======================================================================
2015 //function : EvalLength
2017 //=======================================================================
2019 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2020 Standard_Real& PolynomialCoeff,
2021 const Standard_Real U1, const Standard_Real U2,
2022 const Standard_Real Tol,
2023 Standard_Real& Length, Standard_Real& Error)
2026 Standard_Integer NbSubInt = 1, // Current number of subintervals
2027 MaxNbIter = 13, // Max number of iterations
2028 NbIter = 1; // Current number of iterations
2029 Standard_Real dU,OldLen,LenI;
2031 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2037 dU = (U2-U1)/NbSubInt;
2038 for (i=1; i<=NbSubInt; i++) {
2039 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2043 Error = Abs(OldLen-Length);
2045 while (Error > Tol && NbIter <= MaxNbIter);