2 // Created : Mon Aug 26 07:39:13 1991
5 // Modifed RLE Aug 93 - Complete rewrite
6 // xab 21-Mar-95 implemented cache mecanism
7 // pmn 25-09-96 Interpolation
8 // jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
9 // pmn 07-10-96 : Correction de DN dans le cas rationnal.
10 // pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
12 #include <BSplSLib.ixx>
14 #include <PLib_LocalArray.hxx>
15 #include <BSplCLib.hxx>
16 #include <TColgp_Array2OfXYZ.hxx>
17 #include <TColgp_Array1OfXYZ.hxx>
18 #include <TColStd_HArray1OfInteger.hxx>
19 #include <Standard_NotImplemented.hxx>
20 #include <Standard_ConstructionError.hxx>
21 #include <math_Matrix.hxx>
23 // for null derivatives
24 static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
26 #define M_SQRT2 1.4142135623730950488016887
29 //=======================================================================
30 //struct : BSplCLib_DataContainer
31 //purpose: Auxiliary structure providing buffers for poles and knots used in
32 // evaluation of bspline (allocated in the stack)
33 //=======================================================================
35 struct BSplSLib_DataContainer
37 BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
39 Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
40 VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
41 "BSplSLib: bspline degree is greater than maximum supported");
44 Standard_Real poles[4*(25+1)*(25+1)];
45 Standard_Real knots1[2*25];
46 Standard_Real knots2[2*25];
47 Standard_Real ders[48];
50 typedef PLib_LocalArray BSplSLib_LocalArray;
52 //**************************************************************************
54 //**************************************************************************
56 //=======================================================================
57 //function : RationalDerivative
58 //purpose : computes the rational derivatives when whe have the
59 // the derivatives of the homogeneous numerator and the
60 // the derivatives of the denominator
61 //=======================================================================
63 void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
64 const Standard_Integer VDeg,
65 const Standard_Integer N,
66 const Standard_Integer M,
67 Standard_Real& HDerivatives,
68 Standard_Real& RDerivatives,
69 const Standard_Boolean All)
72 // if All is True all derivatives are computed. if Not only
73 // the requested N, M is computed
76 // let f(u,v) be a rational function = ------------------
80 // Let (N,M) the order of the derivatives we want : then since
83 // Numerator = f * Denominator
87 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
88 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
89 // (0,0) ( p<N q<M p q )
99 // HDerivatives is an array where derivatives are stored in the following form
100 // Numerator is assumee to have 3 functions that is a vector of dimension
103 // (0,0) (0,0) (0, DegV) (0, DegV)
104 // Numerator Denominator ... Numerator Denominator
106 // (1,0) (1,0) (1, DegV) (1, DegV)
107 // Numerator Denominator ... Numerator Denominator
109 // ...........................................................
112 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
113 // Numerator Denominator ... Numerator Denominator
116 Standard_Integer ii,jj,pp,qq,index,index1,index2;
117 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
118 Standard_Integer MinN,MinN1,MinM,MinM1;
119 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
125 M4 = (VDeg + 1) << 2;
127 BSplSLib_LocalArray StoreDerivatives (All ? 0 : ii * 3);
128 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
129 BSplSLib_LocalArray StoreW (ii);
130 Standard_Real *HomogeneousArray = &HDerivatives;
131 Standard_Real denominator,Pii,Pip,Pjq;
133 denominator = 1.0e0 / HomogeneousArray[3];
136 if (UDeg < N) MinN = UDeg;
138 if (VDeg < M) MinM = VDeg;
144 for (ii = 0 ; ii < MinN1 ; ii++) {
150 for (jj = 0 ; jj < MinM1 ; jj++) {
151 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
152 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
153 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
154 StoreW[index_w] = HomogeneousArray[index_v1]; index_w++; index_v1++;
157 for (jj = MinM1 ; jj < M1 ; jj++) {
158 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
159 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
160 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
161 StoreW[index_w] = HomogeneousArray[index_v1]; index_w++; index_v1++;
166 index_v = MinN1 * M3;
167 index_w = MinN1 * M1;
169 for (ii = MinN1 ; ii < N1 ; ii++) {
171 for (jj = 0 ; jj < M1 ; jj++) {
172 RArray[index_v] = 0.0e0; index_v++;
173 RArray[index_v] = 0.0e0; index_v++;
174 RArray[index_v] = 0.0e0; index_v++;
175 StoreW[index_w] = 0.0e0; index_w++;
179 // --------------- Calculation ----------------
184 for (ii = 0 ; ii <= N ; ii++) {
190 for (jj = 0 ; jj <= M ; jj++) {
196 for (pp = 0 ; pp < ii ; pp++) {
200 index2 = jjM1 - ppM1;
201 Pip = PLib::Bin(ii,pp);
203 for (qq = 0 ; qq <= jj ; qq++) {
205 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
206 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
207 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
208 RArray[index1] -= Pjq * RArray[index]; index++;
214 Pii = PLib::Bin(ii,ii);
216 for (qq = 0 ; qq < jj ; qq++) {
218 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
219 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
220 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
221 RArray[index1] -= Pjq * RArray[index]; index++;
224 RArray[index1] *= denominator; index1++;
225 RArray[index1] *= denominator; index1++;
226 RArray[index1] *= denominator;
231 RArray = &RDerivatives;
233 index = (index << 1) + index;
234 RArray[0] = StoreDerivatives[index]; index++;
235 RArray[1] = StoreDerivatives[index]; index++;
236 RArray[2] = StoreDerivatives[index];
240 //=======================================================================
241 //function : PrepareEval
243 //=======================================================================
248 // Prepare all data for computing points :
249 // local arrays of knots
250 // local array of poles (multiplied by the weights if rational)
252 // The first direction to compute (smaller degree) is returned
253 // and the poles are stored according to this direction.
255 static Standard_Boolean PrepareEval
256 (const Standard_Real U,
257 const Standard_Real V,
258 const Standard_Integer Uindex,
259 const Standard_Integer Vindex,
260 const Standard_Integer UDegree,
261 const Standard_Integer VDegree,
262 const Standard_Boolean URat,
263 const Standard_Boolean VRat,
264 const Standard_Boolean UPer,
265 const Standard_Boolean VPer,
266 const TColgp_Array2OfPnt& Poles,
267 const TColStd_Array2OfReal& Weights,
268 const TColStd_Array1OfReal& UKnots,
269 const TColStd_Array1OfReal& VKnots,
270 const TColStd_Array1OfInteger& UMults,
271 const TColStd_Array1OfInteger& VMults,
272 Standard_Real& u1, // first parameter to use
273 Standard_Real& u2, // second parameter to use
274 Standard_Integer& d1, // first degree
275 Standard_Integer& d2, // second degree
276 Standard_Boolean& rational,
277 BSplSLib_DataContainer& dc)
279 rational = URat || VRat;
280 Standard_Integer uindex = Uindex;
281 Standard_Integer vindex = Vindex;
282 Standard_Integer UKLower = UKnots.Lower();
283 Standard_Integer UKUpper = UKnots.Upper();
284 Standard_Integer VKLower = VKnots.Lower();
285 Standard_Integer VKUpper = VKnots.Upper();
286 if (UDegree <= VDegree) {
287 // compute the indices
288 if (uindex < UKLower || uindex > UKUpper)
289 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
291 if (vindex < VKLower || vindex > VKUpper)
292 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
297 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
298 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
299 if (&UMults == NULL) uindex -= UKLower + UDegree;
300 else uindex = BSplCLib::PoleIndex
301 (UDegree,uindex,UPer,UMults);
302 if (&VMults == NULL) vindex -= VKLower + VDegree;
303 else vindex = BSplCLib::PoleIndex
304 (VDegree,vindex,VPer,VMults);
306 // Standard_Integer i,j,k,ip,jp;
307 Standard_Integer i,j,ip,jp;
308 Standard_Real w, *pole = dc.poles;
311 Standard_Integer PLowerRow = Poles.LowerRow();
312 Standard_Integer PUpperRow = Poles.UpperRow();
313 Standard_Integer PLowerCol = Poles.LowerCol();
314 Standard_Integer PUpperCol = Poles.UpperCol();
315 if (rational) { // verify if locally non rational
316 rational = Standard_False;
317 ip = PLowerRow + uindex;
318 jp = PLowerCol + vindex;
319 w = Weights.Value(ip,jp);
320 Standard_Real eps = Epsilon(w);
323 for (i = 0; i <= UDegree && !rational; i++) {
324 jp = PLowerCol + vindex;
326 for (j = 0; j <= VDegree && !rational; j++) {
327 dw = Weights.Value(ip,jp) - w;
328 if (dw < 0) dw = - dw;
331 if (jp > PUpperCol) jp = PLowerCol;
334 if (ip > PUpperRow) ip = PLowerRow;
338 ip = PLowerRow + uindex;
341 for (i = 0; i <= d1; i++) {
342 jp = PLowerCol + vindex;
344 for (j = 0; j <= d2; j++) {
345 const gp_Pnt& P = Poles .Value(ip,jp);
346 pole[3] = w = Weights.Value(ip,jp);
352 if (jp > PUpperCol) jp = PLowerCol;
355 if (ip > PUpperRow) ip = PLowerRow;
360 for (i = 0; i <= d1; i++) {
361 jp = PLowerCol + vindex;
363 for (j = 0; j <= d2; j++) {
364 const gp_Pnt& P = Poles.Value(ip,jp);
370 if (jp > PUpperCol) jp = PLowerCol;
373 if (ip > PUpperRow) ip = PLowerRow;
376 return Standard_True;
379 // compute the indices
380 if (uindex < UKLower || uindex > UKUpper)
381 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
383 if (vindex < VKLower || vindex > VKUpper)
384 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
389 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
390 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
391 if (&UMults == NULL) uindex -= UKLower + UDegree;
392 else uindex = BSplCLib::PoleIndex
393 (UDegree,uindex,UPer,UMults);
394 if (&VMults == NULL) vindex -= VKLower + VDegree;
395 else vindex = BSplCLib::PoleIndex
396 (VDegree,vindex,VPer,VMults);
398 // Standard_Integer i,j,k,ip,jp;
399 Standard_Integer i,j,ip,jp;
400 Standard_Real w, *pole = dc.poles;
403 Standard_Integer PLowerRow = Poles.LowerRow();
404 Standard_Integer PUpperRow = Poles.UpperRow();
405 Standard_Integer PLowerCol = Poles.LowerCol();
406 Standard_Integer PUpperCol = Poles.UpperCol();
407 if (rational) { // verify if locally non rational
408 rational = Standard_False;
409 ip = PLowerRow + uindex;
410 jp = PLowerCol + vindex;
411 w = Weights.Value(ip,jp);
412 Standard_Real eps = Epsilon(w);
415 for (i = 0; i <= UDegree && !rational; i++) {
416 jp = PLowerCol + vindex;
418 for (j = 0; j <= VDegree && !rational; j++) {
419 dw = Weights.Value(ip,jp) - w;
420 if (dw < 0) dw = - dw;
423 if (jp > PUpperCol) jp = PLowerCol;
426 if (ip > PUpperRow) ip = PLowerRow;
430 jp = PLowerCol + vindex;
433 for (i = 0; i <= d1; i++) {
434 ip = PLowerRow + uindex;
436 for (j = 0; j <= d2; j++) {
437 const gp_Pnt& P = Poles .Value(ip,jp);
438 pole[3] = w = Weights.Value(ip,jp);
444 if (ip > PUpperRow) ip = PLowerRow;
447 if (jp > PUpperCol) jp = PLowerCol;
452 for (i = 0; i <= d1; i++) {
453 ip = PLowerRow + uindex;
455 for (j = 0; j <= d2; j++) {
456 const gp_Pnt& P = Poles.Value(ip,jp);
462 if (ip > PUpperRow) ip = PLowerRow;
465 if (jp > PUpperCol) jp = PLowerCol;
468 return Standard_False;
472 //=======================================================================
475 //=======================================================================
478 (const Standard_Real U,
479 const Standard_Real V,
480 const Standard_Integer UIndex,
481 const Standard_Integer VIndex,
482 const TColgp_Array2OfPnt& Poles,
483 const TColStd_Array2OfReal& Weights,
484 const TColStd_Array1OfReal& UKnots,
485 const TColStd_Array1OfReal& VKnots,
486 const TColStd_Array1OfInteger& UMults,
487 const TColStd_Array1OfInteger& VMults,
488 const Standard_Integer UDegree,
489 const Standard_Integer VDegree,
490 const Standard_Boolean URat,
491 const Standard_Boolean VRat,
492 const Standard_Boolean UPer,
493 const Standard_Boolean VPer,
496 // Standard_Integer k ;
521 //=======================================================================
524 //=======================================================================
526 void BSplSLib::HomogeneousD0
527 (const Standard_Real U,
528 const Standard_Real V,
529 const Standard_Integer UIndex,
530 const Standard_Integer VIndex,
531 const TColgp_Array2OfPnt& Poles,
532 const TColStd_Array2OfReal& Weights,
533 const TColStd_Array1OfReal& UKnots,
534 const TColStd_Array1OfReal& VKnots,
535 const TColStd_Array1OfInteger& UMults,
536 const TColStd_Array1OfInteger& VMults,
537 const Standard_Integer UDegree,
538 const Standard_Integer VDegree,
539 const Standard_Boolean URat,
540 const Standard_Boolean VRat,
541 const Standard_Boolean UPer,
542 const Standard_Boolean VPer,
546 Standard_Boolean rational;
547 // Standard_Integer k,dim;
548 Standard_Integer dim;
550 Standard_Integer d1,d2;
553 BSplSLib_DataContainer dc (UDegree, VDegree);
554 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
555 Poles,Weights,UKnots,VKnots,UMults,VMults,
556 u1,u2,d1,d2,rational,dc);
559 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
560 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
568 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
569 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
576 //=======================================================================
579 //=======================================================================
582 (const Standard_Real U,
583 const Standard_Real V,
584 const Standard_Integer UIndex,
585 const Standard_Integer VIndex,
586 const TColgp_Array2OfPnt& Poles,
587 const TColStd_Array2OfReal& Weights,
588 const TColStd_Array1OfReal& UKnots,
589 const TColStd_Array1OfReal& VKnots,
590 const TColStd_Array1OfInteger& UMults,
591 const TColStd_Array1OfInteger& VMults,
592 const Standard_Integer UDegree,
593 const Standard_Integer VDegree,
594 const Standard_Boolean URat,
595 const Standard_Boolean VRat,
596 const Standard_Boolean UPer,
597 const Standard_Boolean VPer,
602 Standard_Boolean rational;
603 // Standard_Integer k,dim,dim2;
604 Standard_Integer dim,dim2;
606 Standard_Integer d1,d2;
607 Standard_Real *result, *resVu, *resVv;
608 BSplSLib_DataContainer dc (UDegree, VDegree);
610 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
611 Poles,Weights,UKnots,VKnots,UMults,VMults,
612 u1,u2,d1,d2,rational,dc)) {
615 dim2 = (d2 + 1) << 2;
616 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
617 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
618 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
619 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
627 dim2 = (dim2 << 1) + dim2;
628 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
629 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
630 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
632 resVu = result + dim2;
639 dim2 = (d2 + 1) << 2;
640 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
641 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
642 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
643 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
651 dim2 = (dim2 << 1) + dim2;
652 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
653 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
654 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
657 resVv = result + dim2;
674 //=======================================================================
677 //=======================================================================
679 void BSplSLib::HomogeneousD1
680 (const Standard_Real U,
681 const Standard_Real V,
682 const Standard_Integer UIndex,
683 const Standard_Integer VIndex,
684 const TColgp_Array2OfPnt& Poles,
685 const TColStd_Array2OfReal& Weights,
686 const TColStd_Array1OfReal& UKnots,
687 const TColStd_Array1OfReal& VKnots,
688 const TColStd_Array1OfInteger& UMults,
689 const TColStd_Array1OfInteger& VMults,
690 const Standard_Integer UDegree,
691 const Standard_Integer VDegree,
692 const Standard_Boolean URat,
693 const Standard_Boolean VRat,
694 const Standard_Boolean UPer,
695 const Standard_Boolean VPer,
703 Standard_Boolean rational;
704 // Standard_Integer k,dim;
705 Standard_Integer dim;
707 Standard_Integer d1,d2;
712 BSplSLib_DataContainer dc (UDegree, VDegree);
713 Standard_Boolean ufirst = PrepareEval
714 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
715 Poles,Weights,UKnots,VKnots,UMults,VMults,
716 u1,u2,d1,d2,rational,dc);
717 dim = rational ? 4 : 3;
719 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
720 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
721 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
723 Standard_Real *result, *resVu, *resVv;
725 resVu = result + (ufirst ? dim*(d2+1) : dim);
726 resVv = result + (ufirst ? dim : dim*(d2+1));
747 //=======================================================================
750 //=======================================================================
753 (const Standard_Real U,
754 const Standard_Real V,
755 const Standard_Integer UIndex,
756 const Standard_Integer VIndex,
757 const TColgp_Array2OfPnt& Poles,
758 const TColStd_Array2OfReal& Weights,
759 const TColStd_Array1OfReal& UKnots,
760 const TColStd_Array1OfReal& VKnots,
761 const TColStd_Array1OfInteger& UMults,
762 const TColStd_Array1OfInteger& VMults,
763 const Standard_Integer UDegree,
764 const Standard_Integer VDegree,
765 const Standard_Boolean URat,
766 const Standard_Boolean VRat,
767 const Standard_Boolean UPer,
768 const Standard_Boolean VPer,
776 Standard_Boolean rational;
777 // Standard_Integer k,dim,dim2;
778 Standard_Integer dim,dim2;
780 Standard_Integer d1,d2;
781 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
782 BSplSLib_DataContainer dc (UDegree, VDegree);
784 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
785 Poles,Weights,UKnots,VKnots,UMults,VMults,
786 u1,u2,d1,d2,rational,dc)) {
789 dim2 = (d2 + 1) << 2;
790 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
791 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
792 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
794 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
795 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
799 resVuu = result + 18;
801 resVuv = result + 12;
806 dim2 = (dim2 << 1) + dim2;
807 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
808 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
809 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
811 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
813 resVu = result + dim2;
815 if (UDegree <= 1) resVuu = BSplSLib_zero;
816 else resVuu = result + (dim2 << 1);
817 if (VDegree <= 1) resVvv = BSplSLib_zero;
818 else resVvv = result + 6;
819 resVuv = result + (d2 << 1) + d2 + 6;
825 dim2 = (d2 + 1) << 2;
826 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
827 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
828 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
830 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
831 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
836 resVvv = result + 18;
837 resVuv = result + 12;
842 dim2 = (dim2 << 1) + dim2;
843 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
844 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
845 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
847 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
850 resVv = result + dim2;
851 if (UDegree <= 1) resVuu = BSplSLib_zero;
852 else resVuu = result + 6;
853 if (VDegree <= 1) resVvv = BSplSLib_zero;
854 else resVvv = result + (dim2 << 1);
855 resVuv = result + (d2 << 1) + d2 + 6;
881 //=======================================================================
884 //=======================================================================
887 (const Standard_Real U,
888 const Standard_Real V,
889 const Standard_Integer UIndex,
890 const Standard_Integer VIndex,
891 const TColgp_Array2OfPnt& Poles,
892 const TColStd_Array2OfReal& Weights,
893 const TColStd_Array1OfReal& UKnots,
894 const TColStd_Array1OfReal& VKnots,
895 const TColStd_Array1OfInteger& UMults,
896 const TColStd_Array1OfInteger& VMults,
897 const Standard_Integer UDegree,
898 const Standard_Integer VDegree,
899 const Standard_Boolean URat,
900 const Standard_Boolean VRat,
901 const Standard_Boolean UPer,
902 const Standard_Boolean VPer,
914 Standard_Boolean rational;
915 // Standard_Integer k,dim,dim2;
916 Standard_Integer dim,dim2;
918 Standard_Integer d1,d2;
919 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
920 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
921 BSplSLib_DataContainer dc (UDegree, VDegree);
923 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
924 Poles,Weights,UKnots,VKnots,UMults,VMults,
925 u1,u2,d1,d2,rational,dc)) {
928 dim2 = (d2 + 1) << 2;
929 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
930 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
931 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
933 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
935 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
936 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
940 resVuu = result + 24;
942 resVuv = result + 15;
943 resVuuu = result + 36;
944 resVvvv = result + 9;
945 resVuuv = result + 27;
946 resVuvv = result + 18;
951 dim2 = (dim2 << 1) + dim2;
952 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
953 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
954 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
956 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
958 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
960 resVu = result + dim2;
963 resVuu = BSplSLib_zero;
964 resVuuv = BSplSLib_zero;
967 resVuu = result + (dim2 << 1);
968 resVuuv = result + (dim2 << 1) + 3;
971 resVvv = BSplSLib_zero;
972 resVuvv = BSplSLib_zero;
976 resVuvv = result + dim2 + 6;
978 resVuv = result + (d2 << 1) + d2 + 6;
979 if (UDegree <= 2) resVuuu = BSplSLib_zero;
980 else resVuuu = result + (dim2 << 1) + dim2;
981 if (VDegree <= 2) resVvvv = BSplSLib_zero;
982 else resVvvv = result + 9;
988 dim2 = (d2 + 1) << 2;
989 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
990 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
991 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
993 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
995 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
996 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1000 resVuu = result + 6;
1001 resVvv = result + 24;
1002 resVuv = result + 15;
1003 resVuuu = result + 9;
1004 resVvvv = result + 36;
1005 resVuuv = result + 18;
1006 resVuvv = result + 27;
1011 dim2 = (dim2 << 1) + dim2;
1012 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1013 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1014 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1016 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1018 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1021 resVv = result + dim2;
1023 resVuu = BSplSLib_zero;
1024 resVuuv = BSplSLib_zero;
1027 resVuu = result + 6;
1028 resVuuv = result + dim2 + 6;
1031 resVvv = BSplSLib_zero;
1032 resVuvv = BSplSLib_zero;
1035 resVvv = result + (dim2 << 1);
1036 resVuvv = result + (dim2 << 1) + 3;
1038 resVuv = result + (d2 << 1) + d2 + 6;
1039 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1040 else resVuuu = result + 9;
1041 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1042 else resVvvv = result + (dim2 << 1) + dim2;
1046 P .SetX(result [0]);
1047 Vu .SetX(resVu [0]);
1048 Vv .SetX(resVv [0]);
1049 Vuu .SetX(resVuu [0]);
1050 Vvv .SetX(resVvv [0]);
1051 Vuv .SetX(resVuv [0]);
1052 Vuuu.SetX(resVuuu[0]);
1053 Vvvv.SetX(resVvvv[0]);
1054 Vuuv.SetX(resVuuv[0]);
1055 Vuvv.SetX(resVuvv[0]);
1057 P .SetY(result [1]);
1058 Vu .SetY(resVu [1]);
1059 Vv .SetY(resVv [1]);
1060 Vuu .SetY(resVuu [1]);
1061 Vvv .SetY(resVvv [1]);
1062 Vuv .SetY(resVuv [1]);
1063 Vuuu.SetY(resVuuu[1]);
1064 Vvvv.SetY(resVvvv[1]);
1065 Vuuv.SetY(resVuuv[1]);
1066 Vuvv.SetY(resVuvv[1]);
1068 P .SetZ(result [2]);
1069 Vu .SetZ(resVu [2]);
1070 Vv .SetZ(resVv [2]);
1071 Vuu .SetZ(resVuu [2]);
1072 Vvv .SetZ(resVvv [2]);
1073 Vuv .SetZ(resVuv [2]);
1074 Vuuu.SetZ(resVuuu[2]);
1075 Vvvv.SetZ(resVvvv[2]);
1076 Vuuv.SetZ(resVuuv[2]);
1077 Vuvv.SetZ(resVuvv[2]);
1080 //=======================================================================
1083 //=======================================================================
1086 (const Standard_Real U,
1087 const Standard_Real V,
1088 const Standard_Integer Nu,
1089 const Standard_Integer Nv,
1090 const Standard_Integer UIndex,
1091 const Standard_Integer VIndex,
1092 const TColgp_Array2OfPnt& Poles,
1093 const TColStd_Array2OfReal& Weights,
1094 const TColStd_Array1OfReal& UKnots,
1095 const TColStd_Array1OfReal& VKnots,
1096 const TColStd_Array1OfInteger& UMults,
1097 const TColStd_Array1OfInteger& VMults,
1098 const Standard_Integer UDegree,
1099 const Standard_Integer VDegree,
1100 const Standard_Boolean URat,
1101 const Standard_Boolean VRat,
1102 const Standard_Boolean UPer,
1103 const Standard_Boolean VPer,
1106 Standard_Boolean rational;
1107 Standard_Integer k,dim;
1108 Standard_Real u1,u2;
1109 Standard_Integer d1,d2;
1111 BSplSLib_DataContainer dc (UDegree, VDegree);
1112 Standard_Boolean ufirst = PrepareEval
1113 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1114 Poles,Weights,UKnots,VKnots,UMults,VMults,
1115 u1,u2,d1,d2,rational,dc);
1116 dim = rational ? 4 : 3;
1119 if ((Nu > UDegree) || (Nv > VDegree)) {
1127 Standard_Integer n1 = ufirst ? Nu : Nv;
1128 Standard_Integer n2 = ufirst ? Nv : Nu;
1130 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1132 for (k = 0; k <= Min(n1,d1); k++)
1133 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1135 Standard_Real *result;
1137 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1138 result = dc.ders; // because Standard_False ci-dessus.
1142 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1151 // Surface modifications
1153 // a surface is processed as a curve of curves.
1154 // i.e : a curve of parameter u where the current point is the set of poles
1158 //=======================================================================
1161 //=======================================================================
1163 void BSplSLib::Iso(const Standard_Real Param,
1164 const Standard_Boolean IsU,
1165 const TColgp_Array2OfPnt& Poles,
1166 const TColStd_Array2OfReal& Weights,
1167 const TColStd_Array1OfReal& Knots,
1168 const TColStd_Array1OfInteger& Mults,
1169 const Standard_Integer Degree,
1170 const Standard_Boolean Periodic,
1171 TColgp_Array1OfPnt& CPoles,
1172 TColStd_Array1OfReal& CWeights)
1174 Standard_Integer index = 0;
1175 Standard_Real u = Param;
1176 Standard_Boolean rational = &Weights != NULL;
1177 Standard_Integer dim = rational ? 4 : 3;
1179 // compute local knots
1181 BSplSLib_LocalArray locknots1 (2*Degree);
1182 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1183 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1185 index -= Knots.Lower() + Degree;
1187 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1190 // copy the local poles
1192 // Standard_Integer f1,l1,f2,l2,i,j,k;
1193 Standard_Integer f1,l1,f2,l2,i,j;
1196 f1 = Poles.LowerRow();
1197 l1 = Poles.UpperRow();
1198 f2 = Poles.LowerCol();
1199 l2 = Poles.UpperCol();
1202 f1 = Poles.LowerCol();
1203 l1 = Poles.UpperCol();
1204 f2 = Poles.LowerRow();
1205 l2 = Poles.UpperRow();
1208 BSplSLib_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim);
1210 Standard_Real w, *pole = locpoles;
1213 for (i = 0; i <= Degree; i++) {
1215 for (j = f2; j <= l2; j++) {
1217 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1219 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1220 pole[0] = P.X() * w;
1221 pole[1] = P.Y() * w;
1222 pole[2] = P.Z() * w;
1232 if (index > l1) index = f1;
1236 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1241 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1242 gp_Pnt& P = CPoles(i);
1244 CWeights(i) = w = pole[3];
1245 P.SetX( pole[0] / w);
1246 P.SetY( pole[1] / w);
1247 P.SetZ( pole[2] / w);
1257 // if the input is not rational but weights are wanted
1258 if (!rational && (&CWeights != NULL)) {
1260 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1265 //=======================================================================
1266 //function : Reverse
1268 //=======================================================================
1270 void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1271 const Standard_Integer Last,
1272 const Standard_Boolean UDirection)
1274 Standard_Integer i,j, l = Last;
1276 l = Poles.LowerRow() +
1277 (l - Poles.LowerRow())%(Poles.ColLength());
1278 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1279 Poles.LowerCol(), Poles.UpperCol());
1281 for (i = Poles.LowerRow(); i <= l; i++) {
1283 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1284 temp(l-i,j) = Poles(i,j);
1288 for (i = l+1; i <= Poles.UpperRow(); i++) {
1290 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1291 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1295 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1297 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1298 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1303 l = Poles.LowerCol() +
1304 (l - Poles.LowerCol())%(Poles.RowLength());
1305 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1306 0, Poles.RowLength()-1);
1308 for (j = Poles.LowerCol(); j <= l; j++) {
1310 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1311 temp(i,l-j) = Poles(i,j);
1315 for (j = l+1; j <= Poles.UpperCol(); j++) {
1317 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1318 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1322 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1324 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1325 Poles(i,j) = temp (i,j-Poles.LowerCol());
1331 //=======================================================================
1332 //function : Reverse
1334 //=======================================================================
1336 void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1337 const Standard_Integer Last,
1338 const Standard_Boolean UDirection)
1340 Standard_Integer i,j, l = Last;
1342 l = Weights.LowerRow() +
1343 (l - Weights.LowerRow())%(Weights.ColLength());
1344 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1345 Weights.LowerCol(), Weights.UpperCol());
1347 for (i = Weights.LowerRow(); i <= l; i++) {
1349 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1350 temp(l-i,j) = Weights(i,j);
1354 for (i = l+1; i <= Weights.UpperRow(); i++) {
1356 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1357 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1361 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1363 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1364 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1369 l = Weights.LowerCol() +
1370 (l - Weights.LowerCol())%(Weights.RowLength());
1371 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1372 0, Weights.RowLength()-1);
1374 for (j = Weights.LowerCol(); j <= l; j++) {
1376 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1377 temp(i,l-j) = Weights(i,j);
1381 for (j = l+1; j <= Weights.UpperCol(); j++) {
1383 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1384 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1388 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1390 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1391 Weights(i,j) = temp (i,j-Weights.LowerCol());
1397 //=======================================================================
1398 //function : IsRational
1400 //=======================================================================
1402 Standard_Boolean BSplSLib::IsRational
1403 (const TColStd_Array2OfReal& Weights,
1404 const Standard_Integer I1,
1405 const Standard_Integer I2,
1406 const Standard_Integer J1,
1407 const Standard_Integer J2,
1408 const Standard_Real Epsi)
1410 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1411 Standard_Integer i,j;
1412 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1413 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1415 for (i = I1 - fi; i < I2 - fi; i++) {
1417 for (j = J1 - fj; j < J2 - fj; j++) {
1418 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1419 return Standard_True;
1422 return Standard_False;
1425 //=======================================================================
1426 //function : SetPoles
1428 //=======================================================================
1430 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1431 TColStd_Array1OfReal& FP,
1432 const Standard_Boolean UDirection)
1434 Standard_Integer i, j, l = FP.Lower();
1435 Standard_Integer PLowerRow = Poles.LowerRow();
1436 Standard_Integer PUpperRow = Poles.UpperRow();
1437 Standard_Integer PLowerCol = Poles.LowerCol();
1438 Standard_Integer PUpperCol = Poles.UpperCol();
1441 for ( i = PLowerRow; i <= PUpperRow; i++) {
1443 for ( j = PLowerCol; j <= PUpperCol; j++) {
1444 const gp_Pnt& P = Poles.Value(i,j);
1453 for ( j = PLowerCol; j <= PUpperCol; j++) {
1455 for ( i = PLowerRow; i <= PUpperRow; i++) {
1456 const gp_Pnt& P = Poles.Value(i,j);
1465 //=======================================================================
1466 //function : SetPoles
1468 //=======================================================================
1470 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1471 const TColStd_Array2OfReal& Weights,
1472 TColStd_Array1OfReal& FP,
1473 const Standard_Boolean UDirection)
1475 Standard_Integer i, j, l = FP.Lower();
1476 Standard_Integer PLowerRow = Poles.LowerRow();
1477 Standard_Integer PUpperRow = Poles.UpperRow();
1478 Standard_Integer PLowerCol = Poles.LowerCol();
1479 Standard_Integer PUpperCol = Poles.UpperCol();
1482 for ( i = PLowerRow; i <= PUpperRow; i++) {
1484 for ( j = PLowerCol; j <= PUpperCol; j++) {
1485 const gp_Pnt& P = Poles .Value(i,j);
1486 Standard_Real w = Weights.Value(i,j);
1487 FP(l) = P.X() * w; l++;
1488 FP(l) = P.Y() * w; l++;
1489 FP(l) = P.Z() * w; l++;
1496 for ( j = PLowerCol; j <= PUpperCol; j++) {
1498 for ( i = PLowerRow; i <= PUpperRow; i++) {
1499 const gp_Pnt& P = Poles .Value(i,j);
1500 Standard_Real w = Weights.Value(i,j);
1501 FP(l) = P.X() * w; l++;
1502 FP(l) = P.Y() * w; l++;
1503 FP(l) = P.Z() * w; l++;
1510 //=======================================================================
1511 //function : GetPoles
1513 //=======================================================================
1515 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1516 TColgp_Array2OfPnt& Poles,
1517 const Standard_Boolean UDirection)
1519 Standard_Integer i, j, l = FP.Lower();
1520 Standard_Integer PLowerRow = Poles.LowerRow();
1521 Standard_Integer PUpperRow = Poles.UpperRow();
1522 Standard_Integer PLowerCol = Poles.LowerCol();
1523 Standard_Integer PUpperCol = Poles.UpperCol();
1526 for ( i = PLowerRow; i <= PUpperRow; i++) {
1528 for ( j = PLowerCol; j <= PUpperCol; j++) {
1529 gp_Pnt& P = Poles.ChangeValue(i,j);
1538 for ( j = PLowerCol; j <= PUpperCol; j++) {
1540 for ( i = PLowerRow; i <= PUpperRow; i++) {
1541 gp_Pnt& P = Poles.ChangeValue(i,j);
1550 //=======================================================================
1551 //function : GetPoles
1553 //=======================================================================
1555 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1556 TColgp_Array2OfPnt& Poles,
1557 TColStd_Array2OfReal& Weights,
1558 const Standard_Boolean UDirection)
1560 Standard_Integer i, j, l = FP.Lower();
1561 Standard_Integer PLowerRow = Poles.LowerRow();
1562 Standard_Integer PUpperRow = Poles.UpperRow();
1563 Standard_Integer PLowerCol = Poles.LowerCol();
1564 Standard_Integer PUpperCol = Poles.UpperCol();
1567 for ( i = PLowerRow; i <= PUpperRow; i++) {
1569 for ( j = PLowerCol; j <= PUpperCol; j++) {
1570 Standard_Real w = FP( l + 3);
1572 gp_Pnt& P = Poles.ChangeValue(i,j);
1573 P.SetX(FP(l) / w); l++;
1574 P.SetY(FP(l) / w); l++;
1575 P.SetZ(FP(l) / w); l++;
1582 for ( j = PLowerCol; j <= PUpperCol; j++) {
1584 for ( i = PLowerRow; i <= PUpperRow; i++) {
1585 Standard_Real w = FP( l + 3);
1587 gp_Pnt& P = Poles.ChangeValue(i,j);
1588 P.SetX(FP(l) / w); l++;
1589 P.SetY(FP(l) / w); l++;
1590 P.SetZ(FP(l) / w); l++;
1597 //=======================================================================
1598 //function : InsertKnots
1600 //=======================================================================
1602 void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1603 const Standard_Integer Degree,
1604 const Standard_Boolean Periodic,
1605 const TColgp_Array2OfPnt& Poles,
1606 const TColStd_Array2OfReal& Weights,
1607 const TColStd_Array1OfReal& Knots,
1608 const TColStd_Array1OfInteger& Mults,
1609 const TColStd_Array1OfReal& AddKnots,
1610 const TColStd_Array1OfInteger& AddMults,
1611 TColgp_Array2OfPnt& NewPoles,
1612 TColStd_Array2OfReal& NewWeights,
1613 TColStd_Array1OfReal& NewKnots,
1614 TColStd_Array1OfInteger& NewMults,
1615 const Standard_Real Epsilon,
1616 const Standard_Boolean Add )
1618 Standard_Boolean rational = &Weights != NULL;
1619 Standard_Integer dim = 3;
1620 if (rational) dim++;
1622 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1623 TColStd_Array1OfReal
1624 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1626 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1627 else SetPoles(Poles,poles,UDirection);
1630 dim *= Poles.RowLength();
1633 dim *= Poles.ColLength();
1635 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1636 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1639 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1640 else GetPoles(newpoles,NewPoles,UDirection);
1643 //=======================================================================
1644 //function : RemoveKnot
1646 //=======================================================================
1648 Standard_Boolean BSplSLib::RemoveKnot
1649 (const Standard_Boolean UDirection,
1650 const Standard_Integer Index,
1651 const Standard_Integer Mult,
1652 const Standard_Integer Degree,
1653 const Standard_Boolean Periodic,
1654 const TColgp_Array2OfPnt& Poles,
1655 const TColStd_Array2OfReal& Weights,
1656 const TColStd_Array1OfReal& Knots,
1657 const TColStd_Array1OfInteger& Mults,
1658 TColgp_Array2OfPnt& NewPoles,
1659 TColStd_Array2OfReal& NewWeights,
1660 TColStd_Array1OfReal& NewKnots,
1661 TColStd_Array1OfInteger& NewMults,
1662 const Standard_Real Tolerance)
1664 Standard_Boolean rational = &Weights != NULL;
1665 Standard_Integer dim = 3;
1666 if (rational) dim++;
1668 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1669 TColStd_Array1OfReal
1670 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1672 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1673 else SetPoles(Poles,poles,UDirection);
1676 dim *= Poles.RowLength();
1679 dim *= Poles.ColLength();
1682 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1683 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1685 return Standard_False;
1687 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1688 else GetPoles(newpoles,NewPoles,UDirection);
1689 return Standard_True;
1692 //=======================================================================
1693 //function : IncreaseDegree
1695 //=======================================================================
1697 void BSplSLib::IncreaseDegree
1698 (const Standard_Boolean UDirection,
1699 const Standard_Integer Degree,
1700 const Standard_Integer NewDegree,
1701 const Standard_Boolean Periodic,
1702 const TColgp_Array2OfPnt& Poles,
1703 const TColStd_Array2OfReal& Weights,
1704 const TColStd_Array1OfReal& Knots,
1705 const TColStd_Array1OfInteger& Mults,
1706 TColgp_Array2OfPnt& NewPoles,
1707 TColStd_Array2OfReal& NewWeights,
1708 TColStd_Array1OfReal& NewKnots,
1709 TColStd_Array1OfInteger& NewMults)
1711 Standard_Boolean rational = &Weights != NULL;
1712 Standard_Integer dim = 3;
1713 if (rational) dim++;
1715 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1716 TColStd_Array1OfReal
1717 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1719 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1720 else SetPoles(Poles,poles,UDirection);
1723 dim *= Poles.RowLength();
1726 dim *= Poles.ColLength();
1729 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1730 newpoles,NewKnots,NewMults);
1732 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1733 else GetPoles(newpoles,NewPoles,UDirection);
1736 //=======================================================================
1737 //function : Unperiodize
1739 //=======================================================================
1741 void BSplSLib::Unperiodize
1742 (const Standard_Boolean UDirection,
1743 const Standard_Integer Degree,
1744 const TColStd_Array1OfInteger& Mults,
1745 const TColStd_Array1OfReal& Knots,
1746 const TColgp_Array2OfPnt& Poles,
1747 const TColStd_Array2OfReal& Weights,
1748 TColStd_Array1OfInteger& NewMults,
1749 TColStd_Array1OfReal& NewKnots,
1750 TColgp_Array2OfPnt& NewPoles,
1751 TColStd_Array2OfReal& NewWeights)
1753 Standard_Boolean rational = &Weights != NULL;
1754 Standard_Integer dim = 3;
1755 if (rational) dim++;
1757 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1758 TColStd_Array1OfReal
1759 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1761 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1762 else SetPoles(Poles,poles,UDirection);
1765 dim *= Poles.RowLength();
1768 dim *= Poles.ColLength();
1771 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1772 NewMults, NewKnots, newpoles);
1774 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1775 else GetPoles(newpoles,NewPoles,UDirection);
1778 //=======================================================================
1779 //function : BuildCache
1780 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1781 // Cache : in case of a rational function the Poles are
1782 // stored in homogeneous form
1783 //=======================================================================
1785 void BSplSLib::BuildCache
1786 (const Standard_Real U,
1787 const Standard_Real V,
1788 const Standard_Real USpanDomain,
1789 const Standard_Real VSpanDomain,
1790 const Standard_Boolean UPeriodic,
1791 const Standard_Boolean VPeriodic,
1792 const Standard_Integer UDegree,
1793 const Standard_Integer VDegree,
1794 const Standard_Integer UIndex,
1795 const Standard_Integer VIndex,
1796 const TColStd_Array1OfReal& UFlatKnots,
1797 const TColStd_Array1OfReal& VFlatKnots,
1798 const TColgp_Array2OfPnt& Poles,
1799 const TColStd_Array2OfReal& Weights,
1800 TColgp_Array2OfPnt& CachePoles,
1801 TColStd_Array2OfReal& CacheWeights)
1803 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1804 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1805 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1806 if (&Weights != NULL)
1807 rational_u = rational_v = Standard_True;
1809 rational_u = rational_v = Standard_False;
1810 BSplSLib_DataContainer dc (UDegree, VDegree);
1826 (BSplCLib::NoMults()),
1827 (BSplCLib::NoMults()),
1837 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1839 for (kk = 0; kk <= d1 ; kk++)
1840 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1842 min_degree_domain = USpanDomain ;
1843 max_degree_domain = VSpanDomain ;
1846 min_degree_domain = VSpanDomain ;
1847 max_degree_domain = USpanDomain ;
1851 for (ii = 0 ; ii <= d2 ; ii++) {
1855 for (jj = 0 ; jj <= d1 ; jj++) {
1857 Index = jj * d2p1 + ii ;
1859 gp_Pnt& P = CachePoles(iii,jjj);
1860 f = factor[0] * factor[1];
1861 P.SetX( f * dc.poles[Index]); Index++;
1862 P.SetY( f * dc.poles[Index]); Index++;
1863 P.SetZ( f * dc.poles[Index]); Index++;
1864 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1865 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1867 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1871 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1873 for (kk = 0; kk <= d1 ; kk++)
1874 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
1876 min_degree_domain = USpanDomain ;
1877 max_degree_domain = VSpanDomain ;
1880 min_degree_domain = VSpanDomain ;
1881 max_degree_domain = USpanDomain ;
1885 for (ii = 0 ; ii <= d2 ; ii++) {
1889 for (jj = 0 ; jj <= d1 ; jj++) {
1891 Index = jj * d2p1 + ii ;
1892 Index = (Index << 1) + Index;
1893 gp_Pnt& P = CachePoles(iii,jjj);
1894 f = factor[0] * factor[1];
1895 P.SetX( f * dc.poles[Index]); Index++;
1896 P.SetY( f * dc.poles[Index]); Index++;
1897 P.SetZ( f * dc.poles[Index]);
1898 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1900 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1902 if (&Weights != NULL) {
1904 // means that PrepareEval did found out that the surface was
1905 // locally polynomial but since the surface is constructed
1906 // with some weights we need to set the weight polynomial to constant
1909 for (ii = 1 ; ii <= d2p1 ; ii++) {
1911 for (jj = 1 ; jj <= d1p1 ; jj++) {
1912 CacheWeights(ii,jj) = 0.0e0 ;
1915 CacheWeights(1,1) = 1.0e0 ;
1920 //=======================================================================
1921 //function : CacheD0
1922 //purpose : Evaluates the polynomial cache of the Bspline Curve
1924 //=======================================================================
1926 void BSplSLib::CacheD0(const Standard_Real UParameter,
1927 const Standard_Real VParameter,
1928 const Standard_Integer UDegree,
1929 const Standard_Integer VDegree,
1930 const Standard_Real UCacheParameter,
1931 const Standard_Real VCacheParameter,
1932 const Standard_Real USpanLenght,
1933 const Standard_Real VSpanLenght,
1934 const TColgp_Array2OfPnt& PolesArray,
1935 const TColStd_Array2OfReal& WeightsArray,
1939 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
1941 // the SpanLenght is the normalizing factor so that the polynomial is between
1954 PArray = (Standard_Real *)
1955 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
1957 myPoint = (Standard_Real *) &aPoint ;
1958 if (UDegree <= VDegree) {
1959 min_degree = UDegree ;
1960 max_degree = VDegree ;
1961 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
1962 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
1963 dimension = 3 * (UDegree + 1) ;
1966 min_degree = VDegree ;
1967 max_degree = UDegree ;
1968 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
1969 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
1970 dimension = 3 * (VDegree + 1) ;
1972 BSplSLib_LocalArray locpoles(dimension);
1974 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
1977 max_degree*dimension,
1981 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
1984 (min_degree << 1) + min_degree,
1987 if (&WeightsArray != NULL) {
1988 dimension = min_degree + 1 ;
1990 WArray = (Standard_Real *)
1991 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
1992 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
1995 max_degree*dimension,
1999 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2005 inverse = 1.0e0/ inverse ;
2007 myPoint[0] *= inverse ;
2008 myPoint[1] *= inverse ;
2009 myPoint[2] *= inverse ;
2013 //=======================================================================
2014 //function : CacheD1
2015 //purpose : Evaluates the polynomial cache of the Bspline Curve
2017 //=======================================================================
2019 void BSplSLib::CacheD1(const Standard_Real UParameter,
2020 const Standard_Real VParameter,
2021 const Standard_Integer UDegree,
2022 const Standard_Integer VDegree,
2023 const Standard_Real UCacheParameter,
2024 const Standard_Real VCacheParameter,
2025 const Standard_Real USpanLenght,
2026 const Standard_Real VSpanLenght,
2027 const TColgp_Array2OfPnt& PolesArray,
2028 const TColStd_Array2OfReal& WeightsArray,
2034 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2036 // the SpanLenght is the normalizing factor so that the polynomial is between
2052 PArray = (Standard_Real *)
2053 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2054 Standard_Real local_poles_array[2][2][3],
2055 local_poles_and_weights_array[2][2][4],
2056 local_weights_array[2][2] ;
2057 Standard_Real * my_vec_min,
2060 my_point = (Standard_Real *) &aPoint ;
2062 // initialize in case of rational evaluation
2063 // because RationalDerivative will use all
2067 if (&WeightsArray != NULL) {
2069 local_poles_array [0][0][0] = 0.0e0 ;
2070 local_poles_array [0][0][1] = 0.0e0 ;
2071 local_poles_array [0][0][2] = 0.0e0 ;
2072 local_weights_array [0][0] = 0.0e0 ;
2073 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2074 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2075 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2076 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2078 local_poles_array [0][1][0] = 0.0e0 ;
2079 local_poles_array [0][1][1] = 0.0e0 ;
2080 local_poles_array [0][1][2] = 0.0e0 ;
2081 local_weights_array [0][1] = 0.0e0 ;
2082 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2083 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2084 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2085 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2087 local_poles_array [1][0][0] = 0.0e0 ;
2088 local_poles_array [1][0][1] = 0.0e0 ;
2089 local_poles_array [1][0][2] = 0.0e0 ;
2090 local_weights_array [1][0] = 0.0e0 ;
2091 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2092 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2093 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2094 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2096 local_poles_array [1][1][0] = 0.0e0 ;
2097 local_poles_array [1][1][1] = 0.0e0 ;
2098 local_poles_array [1][1][2] = 0.0e0 ;
2099 local_weights_array [1][1] = 0.0e0 ;
2100 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2101 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2102 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2103 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2106 if (UDegree <= VDegree) {
2107 min_degree = UDegree ;
2108 max_degree = VDegree ;
2109 inverse_min = 1.0e0/USpanLenght ;
2110 inverse_max = 1.0e0/VSpanLenght ;
2111 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2112 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2114 dimension = 3 * (UDegree + 1) ;
2115 my_vec_min = (Standard_Real *) &aVecU ;
2116 my_vec_max = (Standard_Real *) &aVecV ;
2119 min_degree = VDegree ;
2120 max_degree = UDegree ;
2121 inverse_min = 1.0e0/VSpanLenght ;
2122 inverse_max = 1.0e0/USpanLenght ;
2123 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2124 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2125 dimension = 3 * (VDegree + 1) ;
2126 my_vec_min = (Standard_Real *) &aVecV ;
2127 my_vec_max = (Standard_Real *) &aVecU ;
2130 BSplSLib_LocalArray locpoles (2 * dimension);
2132 PLib::EvalPolynomial(new_parameter[0],
2139 PLib::EvalPolynomial(new_parameter[1],
2144 local_poles_array[0][0][0]) ;
2145 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2148 (min_degree << 1) + min_degree,
2149 locpoles[dimension],
2150 local_poles_array[1][0][0]) ;
2152 if (&WeightsArray != NULL) {
2153 dimension = min_degree + 1 ;
2155 WArray = (Standard_Real *)
2156 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2157 PLib::EvalPolynomial(new_parameter[0],
2164 PLib::EvalPolynomial(new_parameter[1],
2169 local_weights_array[0][0]) ;
2170 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2174 locpoles[dimension],
2175 local_weights_array[1][0]) ;
2177 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2178 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2179 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2180 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2182 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2183 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2184 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2185 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2187 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2188 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2189 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2190 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2192 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2193 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2194 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2195 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2197 BSplSLib::RationalDerivative(1,
2201 local_poles_and_weights_array[0][0][0],
2202 local_poles_array[0][0][0]) ;
2205 my_point [0] = local_poles_array [0][0][0] ;
2206 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2207 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2209 my_point [1] = local_poles_array [0][0][1] ;
2210 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2211 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2213 my_point [2] = local_poles_array [0][0][2] ;
2214 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2215 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2218 //=======================================================================
2219 //function : CacheD2
2220 //purpose : Evaluates the polynomial cache of the Bspline Curve
2222 //=======================================================================
2224 void BSplSLib::CacheD2(const Standard_Real UParameter,
2225 const Standard_Real VParameter,
2226 const Standard_Integer UDegree,
2227 const Standard_Integer VDegree,
2228 const Standard_Real UCacheParameter,
2229 const Standard_Real VCacheParameter,
2230 const Standard_Real USpanLenght,
2231 const Standard_Real VSpanLenght,
2232 const TColgp_Array2OfPnt& PolesArray,
2233 const TColStd_Array2OfReal& WeightsArray,
2242 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2244 // the SpanLenght is the normalizing factor so that the polynomial is between
2261 PArray = (Standard_Real *)
2262 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2263 Standard_Real local_poles_array[3][3][3],
2264 local_poles_and_weights_array[3][3][4],
2265 local_weights_array[3][3] ;
2266 Standard_Real * my_vec_min,
2272 my_point = (Standard_Real *) &aPoint ;
2275 // initialize in case the min and max degree are less than 2
2277 local_poles_array[0][0][0] = 0.0e0 ;
2278 local_poles_array[0][0][1] = 0.0e0 ;
2279 local_poles_array[0][0][2] = 0.0e0 ;
2280 local_poles_array[0][1][0] = 0.0e0 ;
2281 local_poles_array[0][1][1] = 0.0e0 ;
2282 local_poles_array[0][1][2] = 0.0e0 ;
2283 local_poles_array[0][2][0] = 0.0e0 ;
2284 local_poles_array[0][2][1] = 0.0e0 ;
2285 local_poles_array[0][2][2] = 0.0e0 ;
2287 local_poles_array[1][0][0] = 0.0e0 ;
2288 local_poles_array[1][0][1] = 0.0e0 ;
2289 local_poles_array[1][0][2] = 0.0e0 ;
2290 local_poles_array[1][1][0] = 0.0e0 ;
2291 local_poles_array[1][1][1] = 0.0e0 ;
2292 local_poles_array[1][1][2] = 0.0e0 ;
2293 local_poles_array[1][2][0] = 0.0e0 ;
2294 local_poles_array[1][2][1] = 0.0e0 ;
2295 local_poles_array[1][2][2] = 0.0e0 ;
2297 local_poles_array[2][0][0] = 0.0e0 ;
2298 local_poles_array[2][0][1] = 0.0e0 ;
2299 local_poles_array[2][0][2] = 0.0e0 ;
2300 local_poles_array[2][1][0] = 0.0e0 ;
2301 local_poles_array[2][1][1] = 0.0e0 ;
2302 local_poles_array[2][1][2] = 0.0e0 ;
2303 local_poles_array[2][2][0] = 0.0e0 ;
2304 local_poles_array[2][2][1] = 0.0e0 ;
2305 local_poles_array[2][2][2] = 0.0e0 ;
2307 // initialize in case of rational evaluation
2308 // because RationalDerivative will use all
2312 if (&WeightsArray != NULL) {
2314 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2315 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2316 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2317 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2318 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2319 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2320 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2321 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2322 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2324 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2325 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2326 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2327 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2328 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2329 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2330 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2331 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2332 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2334 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2335 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2336 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2337 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2338 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2339 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2340 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2341 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2342 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2344 local_poles_and_weights_array[0][0][3] =
2345 local_weights_array[0][0] = 0.0e0 ;
2346 local_poles_and_weights_array[0][1][3] =
2347 local_weights_array[0][1] = 0.0e0 ;
2348 local_poles_and_weights_array[0][2][3] =
2349 local_weights_array[0][2] = 0.0e0 ;
2350 local_poles_and_weights_array[1][0][3] =
2351 local_weights_array[1][0] = 0.0e0 ;
2352 local_poles_and_weights_array[1][1][3] =
2353 local_weights_array[1][1] = 0.0e0 ;
2354 local_poles_and_weights_array[1][2][3] =
2355 local_weights_array[1][2] = 0.0e0 ;
2356 local_poles_and_weights_array[2][0][3] =
2357 local_weights_array[2][0] = 0.0e0 ;
2358 local_poles_and_weights_array[2][1][3] =
2359 local_weights_array[2][1] = 0.0e0 ;
2360 local_poles_and_weights_array[2][2][3] =
2361 local_weights_array[2][2] = 0.0e0 ;
2364 if (UDegree <= VDegree) {
2365 min_degree = UDegree ;
2366 max_degree = VDegree ;
2367 inverse_min = 1.0e0/USpanLenght ;
2368 inverse_max = 1.0e0/VSpanLenght ;
2369 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2370 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2372 dimension = 3 * (UDegree + 1) ;
2373 my_vec_min = (Standard_Real *) &aVecU ;
2374 my_vec_max = (Standard_Real *) &aVecV ;
2375 my_vec_min_min = (Standard_Real *) &aVecUU ;
2376 my_vec_min_max = (Standard_Real *) &aVecUV ;
2377 my_vec_max_max = (Standard_Real *) &aVecVV ;
2380 min_degree = VDegree ;
2381 max_degree = UDegree ;
2382 inverse_min = 1.0e0/VSpanLenght ;
2383 inverse_max = 1.0e0/USpanLenght ;
2384 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2385 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2386 dimension = 3 * (VDegree + 1) ;
2387 my_vec_min = (Standard_Real *) &aVecV ;
2388 my_vec_max = (Standard_Real *) &aVecU ;
2389 my_vec_min_min = (Standard_Real *) &aVecVV ;
2390 my_vec_min_max = (Standard_Real *) &aVecUV ;
2391 my_vec_max_max = (Standard_Real *) &aVecUU ;
2394 BSplSLib_LocalArray locpoles (3 * dimension);
2397 // initialize in case min or max degree are less than 2
2399 Standard_Integer MinIndMax = 2;
2400 if ( max_degree < 2) MinIndMax = max_degree;
2401 Standard_Integer MinIndMin = 2;
2402 if ( min_degree < 2) MinIndMin = min_degree;
2404 index = MinIndMax * dimension ;
2406 for (ii = MinIndMax ; ii < 3 ; ii++) {
2408 for (kk = 0 ; kk < dimension ; kk++) {
2409 locpoles[index] = 0.0e0 ;
2414 PLib::EvalPolynomial(new_parameter[0],
2421 PLib::EvalPolynomial(new_parameter[1],
2426 local_poles_array[0][0][0]) ;
2427 PLib::EvalPolynomial(new_parameter[1],
2431 locpoles[dimension],
2432 local_poles_array[1][0][0]) ;
2434 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2437 (min_degree << 1) + min_degree,
2438 locpoles[dimension + dimension],
2439 local_poles_array[2][0][0]) ;
2441 if (&WeightsArray != NULL) {
2442 dimension = min_degree + 1 ;
2444 WArray = (Standard_Real *)
2445 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2446 PLib::EvalPolynomial(new_parameter[0],
2453 PLib::EvalPolynomial(new_parameter[1],
2458 local_weights_array[0][0]) ;
2459 PLib::EvalPolynomial(new_parameter[1],
2463 locpoles[dimension],
2464 local_weights_array[1][0]) ;
2465 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2469 locpoles[dimension + dimension],
2470 local_weights_array[2][0]) ;
2473 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2474 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2475 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2476 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2477 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2478 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2479 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2480 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2481 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2483 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2484 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2485 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2486 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2487 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2488 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2489 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2490 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2491 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2493 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2494 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2495 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2496 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2497 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2498 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2499 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2500 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2501 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2504 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2505 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2506 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2507 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2508 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2509 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2510 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2511 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2512 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2514 BSplSLib::RationalDerivative(2,
2518 local_poles_and_weights_array[0][0][0],
2519 local_poles_array[0][0][0]) ;
2523 Standard_Real minmin = inverse_min * inverse_min;
2524 Standard_Real minmax = inverse_min * inverse_max;
2525 Standard_Real maxmax = inverse_max * inverse_max;
2527 my_point [0] = local_poles_array [0][0][0] ;
2528 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2529 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2530 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2531 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2532 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2534 my_point [1] = local_poles_array [0][0][1] ;
2535 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2536 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2537 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2538 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2539 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2541 my_point [2] = local_poles_array [0][0][2] ;
2542 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2543 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2544 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2545 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2546 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2549 //=======================================================================
2550 //function : MovePoint
2551 //purpose : Find the new poles which allows an old point (with a
2552 // given u and v as parameters) to reach a new position
2553 //=======================================================================
2555 void BSplSLib::MovePoint (const Standard_Real U,
2556 const Standard_Real V,
2557 const gp_Vec& Displ,
2558 const Standard_Integer UIndex1,
2559 const Standard_Integer UIndex2,
2560 const Standard_Integer VIndex1,
2561 const Standard_Integer VIndex2,
2562 const Standard_Integer UDegree,
2563 const Standard_Integer VDegree,
2564 const Standard_Boolean Rational,
2565 const TColgp_Array2OfPnt& Poles,
2566 const TColStd_Array2OfReal& Weights,
2567 const TColStd_Array1OfReal& UFlatKnots,
2568 const TColStd_Array1OfReal& VFlatKnots,
2569 Standard_Integer& UFirstIndex,
2570 Standard_Integer& ULastIndex,
2571 Standard_Integer& VFirstIndex,
2572 Standard_Integer& VLastIndex,
2573 TColgp_Array2OfPnt& NewPoles)
2575 // calculate the UBSplineBasis in the parameter U
2576 Standard_Integer UFirstNonZeroBsplineIndex;
2577 math_Matrix UBSplineBasis(1, 1,
2579 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2584 UFirstNonZeroBsplineIndex,
2586 // calculate the VBSplineBasis in the parameter V
2587 Standard_Integer VFirstNonZeroBsplineIndex;
2588 math_Matrix VBSplineBasis(1, 1,
2590 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2595 VFirstNonZeroBsplineIndex,
2597 if (ErrorCod1 || ErrorCod2) {
2605 // find the span which is predominant for parameter U
2606 UFirstIndex = UFirstNonZeroBsplineIndex;
2607 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2608 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2609 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2611 Standard_Real maxValue = 0.0;
2612 Standard_Integer i, ukk1=0, ukk2;
2614 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2615 if (UBSplineBasis(1,i) > maxValue) {
2616 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2617 maxValue = UBSplineBasis(1,i);
2621 // find a ukk2 if symetriy
2623 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2624 if ((ukk1+1) <= ULastIndex) {
2625 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2630 // find the span which is predominant for parameter V
2631 VFirstIndex = VFirstNonZeroBsplineIndex;
2632 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2634 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2635 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2638 Standard_Integer j, vkk1=0, vkk2;
2640 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2641 if (VBSplineBasis(1,j) > maxValue) {
2642 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2643 maxValue = VBSplineBasis(1,j);
2647 // find a vkk2 if symetriy
2649 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2650 if ((vkk1+1) <= VLastIndex) {
2651 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2656 // compute the vector of displacement
2657 Standard_Real D1 = 0.0;
2658 Standard_Real D2 = 0.0;
2659 Standard_Real hN, Coef, DvalU, DvalV;
2661 Standard_Integer ii, jj;
2663 for (i = 1; i <= UDegree+1; i++) {
2664 ii = i + UFirstNonZeroBsplineIndex - 1;
2668 else if (ii > ukk2) {
2675 for (j = 1; j <= VDegree+1; j++) {
2676 jj = j + VFirstNonZeroBsplineIndex - 1;
2678 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2682 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2684 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2688 else if (jj > vkk2) {
2694 D1 += 1./(DvalU + DvalV + 1.) * hN;
2706 // compute the new poles
2708 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2712 else if (i > ukk2) {
2719 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2720 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2724 else if (j > vkk2) {
2730 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2733 NewPoles(i,j) = Poles(i,j);
2739 //=======================================================================
2740 // function : Resolution
2741 // purpose : this computes an estimate for the maximum of the
2742 // partial derivatives both in U and in V
2745 // The calculation resembles at the calculation of curves with
2746 // additional index for the control point. Let Si,j be the
2747 // control points for ls surface and Di,j the weights.
2748 // The checking of upper bounds for the partial derivatives
2749 // will be omitted and Su is the next upper bound in the polynomial case :
2753 // | Si,j - Si-1,j |
2754 // d * Max | ------------- |
2755 // i = 2,n | ti+d - ti |
2759 // and in the rational case :
2763 // Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2764 // Max Max d * -----------------------------------------------
2765 // k=1,n i dans Rj ti+d - ti
2767 // ----------------------------------------------------------------------
2775 // with Rj = {j-d, ...., j+d+d+1}.
2778 //=======================================================================
2780 void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
2781 const TColStd_Array2OfReal& Weights,
2782 const TColStd_Array1OfReal& UKnots,
2783 const TColStd_Array1OfReal& VKnots,
2784 const TColStd_Array1OfInteger& UMults,
2785 const TColStd_Array1OfInteger& VMults,
2786 const Standard_Integer UDegree,
2787 const Standard_Integer VDegree,
2788 const Standard_Boolean URational,
2789 const Standard_Boolean VRational,
2790 const Standard_Boolean UPeriodic,
2791 const Standard_Boolean VPeriodic,
2792 const Standard_Real Tolerance3D,
2793 Standard_Real& UTolerance,
2794 Standard_Real& VTolerance)
2796 Standard_Real Wij,Wmj,Wji,Wjm;
2797 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2798 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2799 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2800 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2802 max_derivative[0] = max_derivative[1] = 0.0e0 ;
2804 Standard_Integer PRowLength, PColLength;
2805 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2806 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2807 Standard_Integer num_poles[2],num_flat_knots[2];
2810 BSplCLib::KnotSequenceLength(UMults,
2814 BSplCLib::KnotSequenceLength(VMults,
2817 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
2818 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
2819 BSplCLib::KnotSequence(UKnots,
2824 BSplCLib::KnotSequence(VKnots,
2829 PRowLength = Poles.RowLength();
2830 PColLength = Poles.ColLength();
2831 if (URational || VRational) {
2832 Standard_Integer Wsize = PRowLength * PColLength;
2833 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2834 min_weights = WG[0];
2836 for (ii = 1 ; ii < Wsize ; ii++) {
2838 if (min_weights > min) min_weights = min;
2841 Standard_Integer UD1 = UDegree + 1;
2842 Standard_Integer VD1 = VDegree + 1;
2843 num_poles[0] = num_flat_knots[0] - UD1;
2844 num_poles[1] = num_flat_knots[1] - VD1;
2845 poles_length[0] = PColLength;
2846 poles_length[1] = PRowLength;
2848 Standard_Integer UD2 = UDegree << 1;
2849 Standard_Integer VD2 = VDegree << 1;
2851 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2852 ii_index = (ii - 1) % poles_length[0] + 1 ;
2853 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2854 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2855 inverse = 1.0e0 / inverse ;
2856 lower[0] = ii - UD1;
2857 if (lower[0] < 1) lower[0] = 1;
2858 upper[0] = ii + UD2 + 1;
2859 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2861 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2862 jj_index = (jj - 1) % poles_length[1] + 1 ;
2863 lower[1] = jj - VD1;
2864 if (lower[1] < 1) lower[1] = 1;
2865 upper[1] = jj + VD2 + 1;
2866 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
2867 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
2868 Wij = Weights.Value(ii_index,jj_index);
2869 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
2870 Wmj = Weights.Value(ii_minus,jj_index);
2878 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
2879 pp_index = (pp - 1) % poles_length[0] + 1 ;
2881 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
2883 qq_index = (qq - 1) % poles_length[1] + 1 ;
2884 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
2888 factor = (Xpq - Xij) * Wij;
2889 factor -= (Xpq - Xmj) * Wmj;
2890 if (factor < 0) factor = - factor;
2892 factor = (Ypq - Yij) * Wij;
2893 factor -= (Ypq - Ymj) * Wmj;
2894 if (factor < 0) factor = - factor;
2896 factor = (Zpq - Zij) * Wij;
2897 factor -= (Zpq - Zmj) * Wmj;
2898 if (factor < 0) factor = - factor;
2901 if (max_derivative[0] < value) max_derivative[0] = value ;
2906 max_derivative[0] /= min_weights ;
2910 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2911 ii_index = (ii - 1) % poles_length[0] + 1 ;
2912 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2913 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2914 inverse = 1.0e0 / inverse ;
2916 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2917 jj_index = (jj - 1) % poles_length[1] + 1 ;
2919 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
2920 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
2921 factor = Pij.X() - Pmj.X();
2922 if (factor < 0) factor = - factor;
2924 factor = Pij.Y() - Pmj.Y();
2925 if (factor < 0) factor = - factor;
2927 factor = Pij.Z() - Pmj.Z();
2928 if (factor < 0) factor = - factor;
2931 if (max_derivative[0] < value) max_derivative[0] = value ;
2935 max_derivative[0] *= UDegree ;
2937 Standard_Integer UD2 = UDegree << 1;
2938 Standard_Integer VD2 = VDegree << 1;
2940 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
2941 ii_index = (ii - 1) % poles_length[1] + 1 ;
2942 ii_minus = (ii - 2) % poles_length[1] + 1 ;
2943 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
2944 inverse = 1.0e0 / inverse ;
2945 lower[0] = ii - VD1;
2946 if (lower[0] < 1) lower[0] = 1;
2947 upper[0] = ii + VD2 + 1;
2948 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
2950 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
2951 jj_index = (jj - 1) % poles_length[0] + 1 ;
2952 lower[1] = jj - UD1;
2953 if (lower[1] < 1) lower[1] = 1;
2954 upper[1] = jj + UD2 + 1;
2955 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
2956 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
2957 Wji = Weights.Value(jj_index,ii_index);
2958 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
2959 Wjm = Weights.Value(jj_index,ii_minus);
2967 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
2968 pp_index = (pp - 1) % poles_length[1] + 1 ;
2970 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
2972 qq_index = (qq - 1) % poles_length[0] + 1 ;
2973 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
2977 factor = (Xqp - Xji) * Wji;
2978 factor -= (Xqp - Xjm) * Wjm;
2979 if (factor < 0) factor = - factor;
2981 factor = (Yqp - Yji) * Wji;
2982 factor -= (Yqp - Yjm) * Wjm;
2983 if (factor < 0) factor = - factor;
2985 factor = (Zqp - Zji) * Wji;
2986 factor -= (Zqp - Zjm) * Wjm;
2987 if (factor < 0) factor = - factor;
2990 if (max_derivative[1] < value) max_derivative[1] = value ;
2995 max_derivative[1] /= min_weights ;
2999 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3000 ii_index = (ii - 1) % poles_length[1] + 1 ;
3001 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3002 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3003 inverse = 1.0e0 / inverse ;
3005 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3006 jj_index = (jj - 1) % poles_length[0] + 1 ;
3008 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3009 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3010 factor = Pji.X() - Pjm.X() ;
3011 if (factor < 0) factor = - factor;
3013 factor = Pji.Y() - Pjm.Y() ;
3014 if (factor < 0) factor = - factor;
3016 factor = Pji.Z() - Pjm.Z() ;
3017 if (factor < 0) factor = - factor;
3020 if (max_derivative[1] < value) max_derivative[1] = value ;
3024 max_derivative[1] *= VDegree ;
3025 max_derivative[0] *= M_SQRT2 ;
3026 max_derivative[1] *= M_SQRT2 ;
3027 if(max_derivative[0] && max_derivative[1]) {
3028 UTolerance = Tolerance3D / max_derivative[0] ;
3029 VTolerance = Tolerance3D / max_derivative[1] ;
3032 UTolerance=VTolerance=0.0;
3034 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3039 //=======================================================================
3040 //function : Interpolate
3042 //=======================================================================
3044 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3045 const Standard_Integer VDegree,
3046 const TColStd_Array1OfReal& UFlatKnots,
3047 const TColStd_Array1OfReal& VFlatKnots,
3048 const TColStd_Array1OfReal& UParameters,
3049 const TColStd_Array1OfReal& VParameters,
3050 TColgp_Array2OfPnt& Poles,
3051 TColStd_Array2OfReal& Weights,
3052 Standard_Integer& InversionProblem)
3054 Standard_Integer ii, jj, ll, kk, dimension;
3055 Standard_Integer ULength = UParameters.Length();
3056 Standard_Integer VLength = VParameters.Length();
3057 Standard_Real * poles_array;
3059 // extraction of iso u
3060 dimension = 4*ULength;
3061 TColStd_Array2OfReal Points(1, VLength,
3064 Handle(TColStd_HArray1OfInteger) ContactOrder =
3065 new (TColStd_HArray1OfInteger)(1, VLength);
3066 ContactOrder->Init(0);
3068 for (ii=1; ii <= VLength; ii++) {
3070 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3071 Points(ii,ll) = Poles(jj, ii).X();
3072 Points(ii,ll+1) = Poles(jj, ii).Y();
3073 Points(ii,ll+2) = Poles(jj, ii).Z();
3074 Points(ii,ll+3) = Weights(jj,ii) ;
3078 // interpolation of iso u
3079 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3080 BSplCLib::Interpolate(VDegree,
3083 ContactOrder->Array1(),
3087 if (InversionProblem != 0) return;
3089 // extraction of iso v
3091 dimension = VLength*4;
3092 TColStd_Array2OfReal IsoPoles(1, ULength,
3095 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3096 ContactOrder->Init(0);
3097 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3099 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3101 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3102 IsoPoles (ii,ll) = Points(jj, kk);
3103 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3104 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3105 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3108 // interpolation of iso v
3109 BSplCLib::Interpolate(UDegree,
3112 ContactOrder->Array1(),
3119 for (ii=1; ii <= ULength; ii++) {
3121 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3122 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3123 Poles.SetValue(ii, jj, Pnt);
3124 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3129 //=======================================================================
3130 //function : Interpolate
3132 //=======================================================================
3134 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3135 const Standard_Integer VDegree,
3136 const TColStd_Array1OfReal& UFlatKnots,
3137 const TColStd_Array1OfReal& VFlatKnots,
3138 const TColStd_Array1OfReal& UParameters,
3139 const TColStd_Array1OfReal& VParameters,
3140 TColgp_Array2OfPnt& Poles,
3141 Standard_Integer& InversionProblem)
3143 Standard_Integer ii, jj, ll, kk, dimension;
3144 Standard_Integer ULength = UParameters.Length();
3145 Standard_Integer VLength = VParameters.Length();
3146 Standard_Real * poles_array;
3148 // extraction of iso u
3149 dimension = 3*ULength;
3150 TColStd_Array2OfReal Points(1, VLength,
3153 Handle(TColStd_HArray1OfInteger) ContactOrder =
3154 new (TColStd_HArray1OfInteger)(1, VLength);
3155 ContactOrder->Init(0);
3157 for (ii=1; ii <= VLength; ii++) {
3159 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3160 Points(ii,ll) = Poles(jj, ii).X();
3161 Points(ii,ll+1) = Poles(jj, ii).Y();
3162 Points(ii,ll+2) = Poles(jj, ii).Z();
3166 // interpolation of iso u
3167 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3168 BSplCLib::Interpolate(VDegree,
3171 ContactOrder->Array1(),
3175 if (InversionProblem != 0) return;
3177 // extraction of iso v
3179 dimension = VLength*3;
3180 TColStd_Array2OfReal IsoPoles(1, ULength,
3183 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3184 ContactOrder->Init(0);
3185 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3187 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3189 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3190 IsoPoles (ii,ll) = Points(jj, kk);
3191 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3192 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3195 // interpolation of iso v
3196 BSplCLib::Interpolate(UDegree,
3199 ContactOrder->Array1(),
3206 for (ii=1; ii <= ULength; ii++) {
3208 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3209 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3210 Poles.SetValue(ii, jj, Pnt);
3215 //=======================================================================
3216 //function : FunctionMultiply
3218 //=======================================================================
3220 void BSplSLib::FunctionMultiply
3221 (const BSplSLib_EvaluatorFunction& Function,
3222 const Standard_Integer UBSplineDegree,
3223 const Standard_Integer VBSplineDegree,
3224 const TColStd_Array1OfReal& UBSplineKnots,
3225 const TColStd_Array1OfReal& VBSplineKnots,
3226 const TColStd_Array1OfInteger & UMults,
3227 const TColStd_Array1OfInteger & VMults,
3228 const TColgp_Array2OfPnt& Poles,
3229 const TColStd_Array2OfReal& Weights,
3230 const TColStd_Array1OfReal& UFlatKnots,
3231 const TColStd_Array1OfReal& VFlatKnots,
3232 const Standard_Integer UNewDegree,
3233 const Standard_Integer VNewDegree,
3234 TColgp_Array2OfPnt& NewNumerator,
3235 TColStd_Array2OfReal& NewDenominator,
3236 Standard_Integer& Status)
3238 Standard_Integer num_uparameters,
3243 Standard_Real result ;
3245 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3246 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3247 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3248 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3250 if ((NewNumerator.ColLength() == num_uparameters) &&
3251 (NewNumerator.RowLength() == num_vparameters) &&
3252 (NewDenominator.ColLength() == num_uparameters) &&
3253 (NewDenominator.RowLength() == num_vparameters)) {
3256 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3260 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3264 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3266 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3267 HomogeneousD0(UParameters(ii),
3283 NewDenominator(ii,jj),
3284 NewNumerator(ii,jj)) ;
3286 Function.Evaluate (0,
3292 Standard_ConstructionError::Raise();
3294 gp_Pnt& P = NewNumerator(ii,jj);
3295 P.SetX(P.X() * result);
3296 P.SetY(P.Y() * result);
3297 P.SetZ(P.Z() * result);
3298 NewDenominator(ii,jj) *= result ;
3301 Interpolate(UNewDegree,
3312 Standard_ConstructionError::Raise();