1 // Created on: 1991-08-26
3 // Copyright (c) 1991-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // Modifed RLE Aug 93 - Complete rewrite
18 // xab 21-Mar-95 implemented cache mecanism
19 // pmn 25-09-96 Interpolation
20 // jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
21 // pmn 07-10-96 : Correction de DN dans le cas rationnal.
22 // pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
24 #include <BSplSLib.ixx>
26 #include <NCollection_LocalArray.hxx>
27 #include <BSplCLib.hxx>
28 #include <TColgp_Array2OfXYZ.hxx>
29 #include <TColgp_Array1OfXYZ.hxx>
30 #include <TColStd_HArray1OfInteger.hxx>
31 #include <Standard_NotImplemented.hxx>
32 #include <Standard_ConstructionError.hxx>
33 #include <math_Matrix.hxx>
35 // for null derivatives
36 static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
38 //=======================================================================
39 //struct : BSplCLib_DataContainer
40 //purpose: Auxiliary structure providing buffers for poles and knots used in
41 // evaluation of bspline (allocated in the stack)
42 //=======================================================================
44 struct BSplSLib_DataContainer
46 BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
48 (void)UDegree; (void)VDegree; // just to avoid compiler warning in Release mode
49 Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
50 VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
51 "BSplSLib: bspline degree is greater than maximum supported");
53 Standard_Real poles[4*(25+1)*(25+1)];
54 Standard_Real knots1[2*25];
55 Standard_Real knots2[2*25];
56 Standard_Real ders[48];
59 //**************************************************************************
61 //**************************************************************************
63 //=======================================================================
64 //function : RationalDerivative
65 //purpose : computes the rational derivatives when whe have the
66 // the derivatives of the homogeneous numerator and the
67 // the derivatives of the denominator
68 //=======================================================================
70 void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
71 const Standard_Integer VDeg,
72 const Standard_Integer N,
73 const Standard_Integer M,
74 Standard_Real& HDerivatives,
75 Standard_Real& RDerivatives,
76 const Standard_Boolean All)
79 // if All is True all derivatives are computed. if Not only
80 // the requested N, M is computed
83 // let f(u,v) be a rational function = ------------------
87 // Let (N,M) the order of the derivatives we want : then since
90 // Numerator = f * Denominator
94 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
95 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
96 // (0,0) ( p<N q<M p q )
106 // HDerivatives is an array where derivatives are stored in the following form
107 // Numerator is assumee to have 3 functions that is a vector of dimension
110 // (0,0) (0,0) (0, DegV) (0, DegV)
111 // Numerator Denominator ... Numerator Denominator
113 // (1,0) (1,0) (1, DegV) (1, DegV)
114 // Numerator Denominator ... Numerator Denominator
116 // ...........................................................
119 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
120 // Numerator Denominator ... Numerator Denominator
123 Standard_Integer ii,jj,pp,qq,index,index1,index2;
124 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
125 Standard_Integer MinN,MinN1,MinM,MinM1;
126 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
132 M4 = (VDeg + 1) << 2;
134 NCollection_LocalArray<Standard_Real> StoreDerivatives (All ? 0 : ii * 3);
135 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
136 NCollection_LocalArray<Standard_Real> StoreW (ii);
137 Standard_Real *HomogeneousArray = &HDerivatives;
138 Standard_Real denominator,Pii,Pip,Pjq;
140 denominator = 1.0e0 / HomogeneousArray[3];
143 if (UDeg < N) MinN = UDeg;
145 if (VDeg < M) MinM = VDeg;
151 for (ii = 0 ; ii < MinN1 ; ii++) {
157 for (jj = 0 ; jj < MinM1 ; jj++) {
158 RArray[index_v++] = HomogeneousArray[index_v1++];
159 RArray[index_v++] = HomogeneousArray[index_v1++];
160 RArray[index_v++] = HomogeneousArray[index_v1++];
161 StoreW[index_w++] = HomogeneousArray[index_v1++];
164 for (jj = MinM1 ; jj < M1 ; jj++) {
165 RArray[index_v++] = 0.;
166 RArray[index_v++] = 0.;
167 RArray[index_v++] = 0.;
168 StoreW[index_w++] = 0.;
173 index_v = MinN1 * M3;
174 index_w = MinN1 * M1;
176 for (ii = MinN1 ; ii < N1 ; ii++) {
178 for (jj = 0 ; jj < M1 ; jj++) {
179 RArray[index_v++] = 0.0e0;
180 RArray[index_v++] = 0.0e0;
181 RArray[index_v++] = 0.0e0;
182 StoreW[index_w++] = 0.0e0;
186 // --------------- Calculation ----------------
191 for (ii = 0 ; ii <= N ; ii++) {
197 for (jj = 0 ; jj <= M ; jj++) {
203 for (pp = 0 ; pp < ii ; pp++) {
207 index2 = jjM1 - ppM1;
208 Pip = PLib::Bin(ii,pp);
210 for (qq = 0 ; qq <= jj ; qq++) {
212 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
213 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
214 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
215 RArray[index1] -= Pjq * RArray[index]; index++;
221 Pii = PLib::Bin(ii,ii);
223 for (qq = 0 ; qq < jj ; qq++) {
225 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
226 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
227 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
228 RArray[index1] -= Pjq * RArray[index]; index++;
231 RArray[index1] *= denominator; index1++;
232 RArray[index1] *= denominator; index1++;
233 RArray[index1] *= denominator;
238 RArray = &RDerivatives;
240 index = (index << 1) + index;
241 RArray[0] = StoreDerivatives[index]; index++;
242 RArray[1] = StoreDerivatives[index]; index++;
243 RArray[2] = StoreDerivatives[index];
247 //=======================================================================
248 //function : PrepareEval
250 //=======================================================================
255 // Prepare all data for computing points :
256 // local arrays of knots
257 // local array of poles (multiplied by the weights if rational)
259 // The first direction to compute (smaller degree) is returned
260 // and the poles are stored according to this direction.
262 static Standard_Boolean PrepareEval (const Standard_Real U,
263 const Standard_Real V,
264 const Standard_Integer Uindex,
265 const Standard_Integer Vindex,
266 const Standard_Integer UDegree,
267 const Standard_Integer VDegree,
268 const Standard_Boolean URat,
269 const Standard_Boolean VRat,
270 const Standard_Boolean UPer,
271 const Standard_Boolean VPer,
272 const TColgp_Array2OfPnt& Poles,
273 const TColStd_Array2OfReal& Weights,
274 const TColStd_Array1OfReal& UKnots,
275 const TColStd_Array1OfReal& VKnots,
276 const TColStd_Array1OfInteger& UMults,
277 const TColStd_Array1OfInteger& VMults,
278 Standard_Real& u1, // first parameter to use
279 Standard_Real& u2, // second parameter to use
280 Standard_Integer& d1, // first degree
281 Standard_Integer& d2, // second degree
282 Standard_Boolean& rational,
283 BSplSLib_DataContainer& dc)
285 rational = URat || VRat;
286 Standard_Integer uindex = Uindex;
287 Standard_Integer vindex = Vindex;
288 Standard_Integer UKLower = UKnots.Lower();
289 Standard_Integer UKUpper = UKnots.Upper();
290 Standard_Integer VKLower = VKnots.Lower();
291 Standard_Integer VKUpper = VKnots.Upper();
293 if (UDegree <= VDegree)
295 // compute the indices
296 if (uindex < UKLower || uindex > UKUpper)
297 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
301 if (vindex < VKLower || vindex > VKUpper)
302 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
309 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
310 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
313 uindex -= UKLower + UDegree;
315 uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
318 vindex -= VKLower + VDegree;
320 vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
323 Standard_Integer i,j,ip,jp;
324 Standard_Real w, *pole = dc.poles;
327 Standard_Integer PLowerRow = Poles.LowerRow();
328 Standard_Integer PUpperRow = Poles.UpperRow();
329 Standard_Integer PLowerCol = Poles.LowerCol();
330 Standard_Integer PUpperCol = Poles.UpperCol();
332 // verify if locally non rational
335 rational = Standard_False;
336 ip = PLowerRow + uindex;
337 jp = PLowerCol + vindex;
339 if(ip < PLowerRow) ip = PUpperRow;
340 if(jp < PLowerCol) jp = PUpperCol;
342 w = Weights.Value(ip,jp);
343 Standard_Real eps = Epsilon(w);
346 for (i = 0; i <= UDegree && !rational; i++)
348 jp = PLowerCol + vindex;
353 for (j = 0; j <= VDegree && !rational; j++)
355 dw = Weights.Value(ip,jp) - w;
359 rational = (dw > eps);
376 ip = PLowerRow + uindex;
383 for (i = 0; i <= d1; i++)
385 jp = PLowerCol + vindex;
390 for (j = 0; j <= d2; j++)
392 const gp_Pnt& P = Poles .Value(ip,jp);
393 pole[3] = w = Weights.Value(ip,jp);
413 for (i = 0; i <= d1; i++)
415 jp = PLowerCol + vindex;
420 for (j = 0; j <= d2; j++)
422 const gp_Pnt& P = Poles.Value(ip,jp);
440 return Standard_True;
444 // compute the indices
445 if (uindex < UKLower || uindex > UKUpper)
446 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
450 if (vindex < VKLower || vindex > VKUpper)
451 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
460 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
461 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
464 uindex -= UKLower + UDegree;
466 uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
469 vindex -= VKLower + VDegree;
471 vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
474 Standard_Integer i,j,ip,jp;
475 Standard_Real w, *pole = dc.poles;
478 Standard_Integer PLowerRow = Poles.LowerRow();
479 Standard_Integer PUpperRow = Poles.UpperRow();
480 Standard_Integer PLowerCol = Poles.LowerCol();
481 Standard_Integer PUpperCol = Poles.UpperCol();
483 // verify if locally non rational
486 rational = Standard_False;
487 ip = PLowerRow + uindex;
488 jp = PLowerCol + vindex;
496 w = Weights.Value(ip,jp);
497 Standard_Real eps = Epsilon(w);
500 for (i = 0; i <= UDegree && !rational; i++)
502 jp = PLowerCol + vindex;
507 for (j = 0; j <= VDegree && !rational; j++)
509 dw = Weights.Value(ip,jp) - w;
510 if (dw < 0) dw = - dw;
528 jp = PLowerCol + vindex;
535 for (i = 0; i <= d1; i++)
537 ip = PLowerRow + uindex;
542 for (j = 0; j <= d2; j++)
544 const gp_Pnt& P = Poles.Value(ip,jp);
545 pole[3] = w = Weights.Value(ip,jp);
566 for (i = 0; i <= d1; i++)
568 ip = PLowerRow + uindex;
573 for (j = 0; j <= d2; j++)
575 const gp_Pnt& P = Poles.Value(ip,jp);
594 return Standard_False;
598 //=======================================================================
601 //=======================================================================
604 (const Standard_Real U,
605 const Standard_Real V,
606 const Standard_Integer UIndex,
607 const Standard_Integer VIndex,
608 const TColgp_Array2OfPnt& Poles,
609 const TColStd_Array2OfReal& Weights,
610 const TColStd_Array1OfReal& UKnots,
611 const TColStd_Array1OfReal& VKnots,
612 const TColStd_Array1OfInteger& UMults,
613 const TColStd_Array1OfInteger& VMults,
614 const Standard_Integer UDegree,
615 const Standard_Integer VDegree,
616 const Standard_Boolean URat,
617 const Standard_Boolean VRat,
618 const Standard_Boolean UPer,
619 const Standard_Boolean VPer,
622 // Standard_Integer k ;
647 //=======================================================================
650 //=======================================================================
652 void BSplSLib::HomogeneousD0
653 (const Standard_Real U,
654 const Standard_Real V,
655 const Standard_Integer UIndex,
656 const Standard_Integer VIndex,
657 const TColgp_Array2OfPnt& Poles,
658 const TColStd_Array2OfReal& Weights,
659 const TColStd_Array1OfReal& UKnots,
660 const TColStd_Array1OfReal& VKnots,
661 const TColStd_Array1OfInteger& UMults,
662 const TColStd_Array1OfInteger& VMults,
663 const Standard_Integer UDegree,
664 const Standard_Integer VDegree,
665 const Standard_Boolean URat,
666 const Standard_Boolean VRat,
667 const Standard_Boolean UPer,
668 const Standard_Boolean VPer,
672 Standard_Boolean rational;
673 // Standard_Integer k,dim;
674 Standard_Integer dim;
676 Standard_Integer d1,d2;
679 BSplSLib_DataContainer dc (UDegree, VDegree);
680 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
681 Poles,Weights,UKnots,VKnots,UMults,VMults,
682 u1,u2,d1,d2,rational,dc);
685 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
686 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
694 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
695 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
702 //=======================================================================
705 //=======================================================================
708 (const Standard_Real U,
709 const Standard_Real V,
710 const Standard_Integer UIndex,
711 const Standard_Integer VIndex,
712 const TColgp_Array2OfPnt& Poles,
713 const TColStd_Array2OfReal& Weights,
714 const TColStd_Array1OfReal& UKnots,
715 const TColStd_Array1OfReal& VKnots,
716 const TColStd_Array1OfInteger& UMults,
717 const TColStd_Array1OfInteger& VMults,
718 const Standard_Integer UDegree,
719 const Standard_Integer VDegree,
720 const Standard_Boolean URat,
721 const Standard_Boolean VRat,
722 const Standard_Boolean UPer,
723 const Standard_Boolean VPer,
728 Standard_Boolean rational;
729 // Standard_Integer k,dim,dim2;
730 Standard_Integer dim,dim2;
732 Standard_Integer d1,d2;
733 Standard_Real *result, *resVu, *resVv;
734 BSplSLib_DataContainer dc (UDegree, VDegree);
736 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
737 Poles,Weights,UKnots,VKnots,UMults,VMults,
738 u1,u2,d1,d2,rational,dc)) {
741 dim2 = (d2 + 1) << 2;
742 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
743 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
744 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
745 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
753 dim2 = (dim2 << 1) + dim2;
754 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
755 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
756 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
758 resVu = result + dim2;
765 dim2 = (d2 + 1) << 2;
766 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
767 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
768 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
769 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
777 dim2 = (dim2 << 1) + dim2;
778 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
779 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
780 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
783 resVv = result + dim2;
800 //=======================================================================
803 //=======================================================================
805 void BSplSLib::HomogeneousD1
806 (const Standard_Real U,
807 const Standard_Real V,
808 const Standard_Integer UIndex,
809 const Standard_Integer VIndex,
810 const TColgp_Array2OfPnt& Poles,
811 const TColStd_Array2OfReal& Weights,
812 const TColStd_Array1OfReal& UKnots,
813 const TColStd_Array1OfReal& VKnots,
814 const TColStd_Array1OfInteger& UMults,
815 const TColStd_Array1OfInteger& VMults,
816 const Standard_Integer UDegree,
817 const Standard_Integer VDegree,
818 const Standard_Boolean URat,
819 const Standard_Boolean VRat,
820 const Standard_Boolean UPer,
821 const Standard_Boolean VPer,
829 Standard_Boolean rational;
830 // Standard_Integer k,dim;
831 Standard_Integer dim;
833 Standard_Integer d1,d2;
838 BSplSLib_DataContainer dc (UDegree, VDegree);
839 Standard_Boolean ufirst = PrepareEval
840 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
841 Poles,Weights,UKnots,VKnots,UMults,VMults,
842 u1,u2,d1,d2,rational,dc);
843 dim = rational ? 4 : 3;
845 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
846 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
847 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
849 Standard_Real *result, *resVu, *resVv;
851 resVu = result + (ufirst ? dim*(d2+1) : dim);
852 resVv = result + (ufirst ? dim : dim*(d2+1));
873 //=======================================================================
876 //=======================================================================
879 (const Standard_Real U,
880 const Standard_Real V,
881 const Standard_Integer UIndex,
882 const Standard_Integer VIndex,
883 const TColgp_Array2OfPnt& Poles,
884 const TColStd_Array2OfReal& Weights,
885 const TColStd_Array1OfReal& UKnots,
886 const TColStd_Array1OfReal& VKnots,
887 const TColStd_Array1OfInteger& UMults,
888 const TColStd_Array1OfInteger& VMults,
889 const Standard_Integer UDegree,
890 const Standard_Integer VDegree,
891 const Standard_Boolean URat,
892 const Standard_Boolean VRat,
893 const Standard_Boolean UPer,
894 const Standard_Boolean VPer,
902 Standard_Boolean rational;
903 // Standard_Integer k,dim,dim2;
904 Standard_Integer dim,dim2;
906 Standard_Integer d1,d2;
907 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
908 BSplSLib_DataContainer dc (UDegree, VDegree);
910 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
911 Poles,Weights,UKnots,VKnots,UMults,VMults,
912 u1,u2,d1,d2,rational,dc)) {
915 dim2 = (d2 + 1) << 2;
916 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
917 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
918 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
920 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
921 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
925 resVuu = result + 18;
927 resVuv = result + 12;
932 dim2 = (dim2 << 1) + dim2;
933 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
934 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
935 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
937 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
939 resVu = result + dim2;
941 if (UDegree <= 1) resVuu = BSplSLib_zero;
942 else resVuu = result + (dim2 << 1);
943 if (VDegree <= 1) resVvv = BSplSLib_zero;
944 else resVvv = result + 6;
945 resVuv = result + (d2 << 1) + d2 + 6;
951 dim2 = (d2 + 1) << 2;
952 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
953 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
954 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
956 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
957 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
962 resVvv = result + 18;
963 resVuv = result + 12;
968 dim2 = (dim2 << 1) + dim2;
969 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
970 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
971 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
973 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
976 resVv = result + dim2;
977 if (UDegree <= 1) resVuu = BSplSLib_zero;
978 else resVuu = result + 6;
979 if (VDegree <= 1) resVvv = BSplSLib_zero;
980 else resVvv = result + (dim2 << 1);
981 resVuv = result + (d2 << 1) + d2 + 6;
1000 Vu .SetZ(resVu [2]);
1001 Vv .SetZ(resVv [2]);
1002 Vuu.SetZ(resVuu[2]);
1003 Vvv.SetZ(resVvv[2]);
1004 Vuv.SetZ(resVuv[2]);
1007 //=======================================================================
1010 //=======================================================================
1013 (const Standard_Real U,
1014 const Standard_Real V,
1015 const Standard_Integer UIndex,
1016 const Standard_Integer VIndex,
1017 const TColgp_Array2OfPnt& Poles,
1018 const TColStd_Array2OfReal& Weights,
1019 const TColStd_Array1OfReal& UKnots,
1020 const TColStd_Array1OfReal& VKnots,
1021 const TColStd_Array1OfInteger& UMults,
1022 const TColStd_Array1OfInteger& VMults,
1023 const Standard_Integer UDegree,
1024 const Standard_Integer VDegree,
1025 const Standard_Boolean URat,
1026 const Standard_Boolean VRat,
1027 const Standard_Boolean UPer,
1028 const Standard_Boolean VPer,
1040 Standard_Boolean rational;
1041 // Standard_Integer k,dim,dim2;
1042 Standard_Integer dim,dim2;
1043 Standard_Real u1,u2;
1044 Standard_Integer d1,d2;
1045 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
1046 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
1047 BSplSLib_DataContainer dc (UDegree, VDegree);
1049 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1050 Poles,Weights,UKnots,VKnots,UMults,VMults,
1051 u1,u2,d1,d2,rational,dc)) {
1054 dim2 = (d2 + 1) << 2;
1055 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1056 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1057 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1059 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1061 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1062 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1064 resVu = result + 12;
1066 resVuu = result + 24;
1067 resVvv = result + 6;
1068 resVuv = result + 15;
1069 resVuuu = result + 36;
1070 resVvvv = result + 9;
1071 resVuuv = result + 27;
1072 resVuvv = result + 18;
1077 dim2 = (dim2 << 1) + dim2;
1078 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1079 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1080 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1082 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1084 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1086 resVu = result + dim2;
1089 resVuu = BSplSLib_zero;
1090 resVuuv = BSplSLib_zero;
1093 resVuu = result + (dim2 << 1);
1094 resVuuv = result + (dim2 << 1) + 3;
1097 resVvv = BSplSLib_zero;
1098 resVuvv = BSplSLib_zero;
1101 resVvv = result + 6;
1102 resVuvv = result + dim2 + 6;
1104 resVuv = result + (d2 << 1) + d2 + 6;
1105 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1106 else resVuuu = result + (dim2 << 1) + dim2;
1107 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1108 else resVvvv = result + 9;
1114 dim2 = (d2 + 1) << 2;
1115 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1116 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1117 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1119 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1121 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1122 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1125 resVv = result + 12;
1126 resVuu = result + 6;
1127 resVvv = result + 24;
1128 resVuv = result + 15;
1129 resVuuu = result + 9;
1130 resVvvv = result + 36;
1131 resVuuv = result + 18;
1132 resVuvv = result + 27;
1137 dim2 = (dim2 << 1) + dim2;
1138 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1139 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1140 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1142 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1144 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1147 resVv = result + dim2;
1149 resVuu = BSplSLib_zero;
1150 resVuuv = BSplSLib_zero;
1153 resVuu = result + 6;
1154 resVuuv = result + dim2 + 6;
1157 resVvv = BSplSLib_zero;
1158 resVuvv = BSplSLib_zero;
1161 resVvv = result + (dim2 << 1);
1162 resVuvv = result + (dim2 << 1) + 3;
1164 resVuv = result + (d2 << 1) + d2 + 6;
1165 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1166 else resVuuu = result + 9;
1167 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1168 else resVvvv = result + (dim2 << 1) + dim2;
1172 P .SetX(result [0]);
1173 Vu .SetX(resVu [0]);
1174 Vv .SetX(resVv [0]);
1175 Vuu .SetX(resVuu [0]);
1176 Vvv .SetX(resVvv [0]);
1177 Vuv .SetX(resVuv [0]);
1178 Vuuu.SetX(resVuuu[0]);
1179 Vvvv.SetX(resVvvv[0]);
1180 Vuuv.SetX(resVuuv[0]);
1181 Vuvv.SetX(resVuvv[0]);
1183 P .SetY(result [1]);
1184 Vu .SetY(resVu [1]);
1185 Vv .SetY(resVv [1]);
1186 Vuu .SetY(resVuu [1]);
1187 Vvv .SetY(resVvv [1]);
1188 Vuv .SetY(resVuv [1]);
1189 Vuuu.SetY(resVuuu[1]);
1190 Vvvv.SetY(resVvvv[1]);
1191 Vuuv.SetY(resVuuv[1]);
1192 Vuvv.SetY(resVuvv[1]);
1194 P .SetZ(result [2]);
1195 Vu .SetZ(resVu [2]);
1196 Vv .SetZ(resVv [2]);
1197 Vuu .SetZ(resVuu [2]);
1198 Vvv .SetZ(resVvv [2]);
1199 Vuv .SetZ(resVuv [2]);
1200 Vuuu.SetZ(resVuuu[2]);
1201 Vvvv.SetZ(resVvvv[2]);
1202 Vuuv.SetZ(resVuuv[2]);
1203 Vuvv.SetZ(resVuvv[2]);
1206 //=======================================================================
1209 //=======================================================================
1212 (const Standard_Real U,
1213 const Standard_Real V,
1214 const Standard_Integer Nu,
1215 const Standard_Integer Nv,
1216 const Standard_Integer UIndex,
1217 const Standard_Integer VIndex,
1218 const TColgp_Array2OfPnt& Poles,
1219 const TColStd_Array2OfReal& Weights,
1220 const TColStd_Array1OfReal& UKnots,
1221 const TColStd_Array1OfReal& VKnots,
1222 const TColStd_Array1OfInteger& UMults,
1223 const TColStd_Array1OfInteger& VMults,
1224 const Standard_Integer UDegree,
1225 const Standard_Integer VDegree,
1226 const Standard_Boolean URat,
1227 const Standard_Boolean VRat,
1228 const Standard_Boolean UPer,
1229 const Standard_Boolean VPer,
1232 Standard_Boolean rational;
1233 Standard_Integer k,dim;
1234 Standard_Real u1,u2;
1235 Standard_Integer d1,d2;
1237 BSplSLib_DataContainer dc (UDegree, VDegree);
1238 Standard_Boolean ufirst = PrepareEval
1239 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1240 Poles,Weights,UKnots,VKnots,UMults,VMults,
1241 u1,u2,d1,d2,rational,dc);
1242 dim = rational ? 4 : 3;
1245 if ((Nu > UDegree) || (Nv > VDegree)) {
1253 Standard_Integer n1 = ufirst ? Nu : Nv;
1254 Standard_Integer n2 = ufirst ? Nv : Nu;
1256 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1258 for (k = 0; k <= Min(n1,d1); k++)
1259 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1261 Standard_Real *result;
1263 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1264 result = dc.ders; // because Standard_False ci-dessus.
1268 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1277 // Surface modifications
1279 // a surface is processed as a curve of curves.
1280 // i.e : a curve of parameter u where the current point is the set of poles
1284 //=======================================================================
1287 //=======================================================================
1289 void BSplSLib::Iso(const Standard_Real Param,
1290 const Standard_Boolean IsU,
1291 const TColgp_Array2OfPnt& Poles,
1292 const TColStd_Array2OfReal& Weights,
1293 const TColStd_Array1OfReal& Knots,
1294 const TColStd_Array1OfInteger& Mults,
1295 const Standard_Integer Degree,
1296 const Standard_Boolean Periodic,
1297 TColgp_Array1OfPnt& CPoles,
1298 TColStd_Array1OfReal& CWeights)
1300 Standard_Integer index = 0;
1301 Standard_Real u = Param;
1302 Standard_Boolean rational = &Weights != NULL;
1303 Standard_Integer dim = rational ? 4 : 3;
1305 // compute local knots
1307 NCollection_LocalArray<Standard_Real> locknots1 (2*Degree);
1308 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1309 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1311 index -= Knots.Lower() + Degree;
1313 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1316 // copy the local poles
1318 // Standard_Integer f1,l1,f2,l2,i,j,k;
1319 Standard_Integer f1,l1,f2,l2,i,j;
1322 f1 = Poles.LowerRow();
1323 l1 = Poles.UpperRow();
1324 f2 = Poles.LowerCol();
1325 l2 = Poles.UpperCol();
1328 f1 = Poles.LowerCol();
1329 l1 = Poles.UpperCol();
1330 f2 = Poles.LowerRow();
1331 l2 = Poles.UpperRow();
1334 NCollection_LocalArray<Standard_Real> locpoles ((Degree+1) * (l2-f2+1) * dim);
1336 Standard_Real w, *pole = locpoles;
1339 for (i = 0; i <= Degree; i++) {
1341 for (j = f2; j <= l2; j++) {
1343 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1345 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1346 pole[0] = P.X() * w;
1347 pole[1] = P.Y() * w;
1348 pole[2] = P.Z() * w;
1358 if (index > l1) index = f1;
1362 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1367 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1368 gp_Pnt& P = CPoles(i);
1370 CWeights(i) = w = pole[3];
1371 P.SetX( pole[0] / w);
1372 P.SetY( pole[1] / w);
1373 P.SetZ( pole[2] / w);
1383 // if the input is not rational but weights are wanted
1384 if (!rational && (&CWeights != NULL)) {
1386 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1391 //=======================================================================
1392 //function : Reverse
1394 //=======================================================================
1396 void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1397 const Standard_Integer Last,
1398 const Standard_Boolean UDirection)
1400 Standard_Integer i,j, l = Last;
1402 l = Poles.LowerRow() +
1403 (l - Poles.LowerRow())%(Poles.ColLength());
1404 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1405 Poles.LowerCol(), Poles.UpperCol());
1407 for (i = Poles.LowerRow(); i <= l; i++) {
1409 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1410 temp(l-i,j) = Poles(i,j);
1414 for (i = l+1; i <= Poles.UpperRow(); i++) {
1416 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1417 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1421 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1423 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1424 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1429 l = Poles.LowerCol() +
1430 (l - Poles.LowerCol())%(Poles.RowLength());
1431 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1432 0, Poles.RowLength()-1);
1434 for (j = Poles.LowerCol(); j <= l; j++) {
1436 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1437 temp(i,l-j) = Poles(i,j);
1441 for (j = l+1; j <= Poles.UpperCol(); j++) {
1443 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1444 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1448 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1450 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1451 Poles(i,j) = temp (i,j-Poles.LowerCol());
1457 //=======================================================================
1458 //function : Reverse
1460 //=======================================================================
1462 void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1463 const Standard_Integer Last,
1464 const Standard_Boolean UDirection)
1466 Standard_Integer i,j, l = Last;
1468 l = Weights.LowerRow() +
1469 (l - Weights.LowerRow())%(Weights.ColLength());
1470 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1471 Weights.LowerCol(), Weights.UpperCol());
1473 for (i = Weights.LowerRow(); i <= l; i++) {
1475 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1476 temp(l-i,j) = Weights(i,j);
1480 for (i = l+1; i <= Weights.UpperRow(); i++) {
1482 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1483 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1487 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1489 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1490 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1495 l = Weights.LowerCol() +
1496 (l - Weights.LowerCol())%(Weights.RowLength());
1497 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1498 0, Weights.RowLength()-1);
1500 for (j = Weights.LowerCol(); j <= l; j++) {
1502 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1503 temp(i,l-j) = Weights(i,j);
1507 for (j = l+1; j <= Weights.UpperCol(); j++) {
1509 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1510 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1514 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1516 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1517 Weights(i,j) = temp (i,j-Weights.LowerCol());
1523 //=======================================================================
1524 //function : IsRational
1526 //=======================================================================
1528 Standard_Boolean BSplSLib::IsRational
1529 (const TColStd_Array2OfReal& Weights,
1530 const Standard_Integer I1,
1531 const Standard_Integer I2,
1532 const Standard_Integer J1,
1533 const Standard_Integer J2,
1534 const Standard_Real Epsi)
1536 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1537 Standard_Integer i,j;
1538 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1539 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1541 for (i = I1 - fi; i < I2 - fi; i++) {
1543 for (j = J1 - fj; j < J2 - fj; j++) {
1544 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1545 return Standard_True;
1548 return Standard_False;
1551 //=======================================================================
1552 //function : SetPoles
1554 //=======================================================================
1556 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1557 TColStd_Array1OfReal& FP,
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 const gp_Pnt& P = Poles.Value(i,j);
1579 for ( j = PLowerCol; j <= PUpperCol; j++) {
1581 for ( i = PLowerRow; i <= PUpperRow; i++) {
1582 const gp_Pnt& P = Poles.Value(i,j);
1591 //=======================================================================
1592 //function : SetPoles
1594 //=======================================================================
1596 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1597 const TColStd_Array2OfReal& Weights,
1598 TColStd_Array1OfReal& FP,
1599 const Standard_Boolean UDirection)
1601 Standard_Integer i, j, l = FP.Lower();
1602 Standard_Integer PLowerRow = Poles.LowerRow();
1603 Standard_Integer PUpperRow = Poles.UpperRow();
1604 Standard_Integer PLowerCol = Poles.LowerCol();
1605 Standard_Integer PUpperCol = Poles.UpperCol();
1608 for ( i = PLowerRow; i <= PUpperRow; i++) {
1610 for ( j = PLowerCol; j <= PUpperCol; j++) {
1611 const gp_Pnt& P = Poles .Value(i,j);
1612 Standard_Real w = Weights.Value(i,j);
1613 FP(l) = P.X() * w; l++;
1614 FP(l) = P.Y() * w; l++;
1615 FP(l) = P.Z() * w; l++;
1622 for ( j = PLowerCol; j <= PUpperCol; j++) {
1624 for ( i = PLowerRow; i <= PUpperRow; i++) {
1625 const gp_Pnt& P = Poles .Value(i,j);
1626 Standard_Real w = Weights.Value(i,j);
1627 FP(l) = P.X() * w; l++;
1628 FP(l) = P.Y() * w; l++;
1629 FP(l) = P.Z() * w; l++;
1636 //=======================================================================
1637 //function : GetPoles
1639 //=======================================================================
1641 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1642 TColgp_Array2OfPnt& Poles,
1643 const Standard_Boolean UDirection)
1645 Standard_Integer i, j, l = FP.Lower();
1646 Standard_Integer PLowerRow = Poles.LowerRow();
1647 Standard_Integer PUpperRow = Poles.UpperRow();
1648 Standard_Integer PLowerCol = Poles.LowerCol();
1649 Standard_Integer PUpperCol = Poles.UpperCol();
1652 for ( i = PLowerRow; i <= PUpperRow; i++) {
1654 for ( j = PLowerCol; j <= PUpperCol; j++) {
1655 gp_Pnt& P = Poles.ChangeValue(i,j);
1664 for ( j = PLowerCol; j <= PUpperCol; j++) {
1666 for ( i = PLowerRow; i <= PUpperRow; i++) {
1667 gp_Pnt& P = Poles.ChangeValue(i,j);
1676 //=======================================================================
1677 //function : GetPoles
1679 //=======================================================================
1681 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1682 TColgp_Array2OfPnt& Poles,
1683 TColStd_Array2OfReal& Weights,
1684 const Standard_Boolean UDirection)
1686 Standard_Integer i, j, l = FP.Lower();
1687 Standard_Integer PLowerRow = Poles.LowerRow();
1688 Standard_Integer PUpperRow = Poles.UpperRow();
1689 Standard_Integer PLowerCol = Poles.LowerCol();
1690 Standard_Integer PUpperCol = Poles.UpperCol();
1693 for ( i = PLowerRow; i <= PUpperRow; i++) {
1695 for ( j = PLowerCol; j <= PUpperCol; j++) {
1696 Standard_Real w = FP( l + 3);
1698 gp_Pnt& P = Poles.ChangeValue(i,j);
1699 P.SetX(FP(l) / w); l++;
1700 P.SetY(FP(l) / w); l++;
1701 P.SetZ(FP(l) / w); l++;
1708 for ( j = PLowerCol; j <= PUpperCol; j++) {
1710 for ( i = PLowerRow; i <= PUpperRow; i++) {
1711 Standard_Real w = FP( l + 3);
1713 gp_Pnt& P = Poles.ChangeValue(i,j);
1714 P.SetX(FP(l) / w); l++;
1715 P.SetY(FP(l) / w); l++;
1716 P.SetZ(FP(l) / w); l++;
1723 //=======================================================================
1724 //function : InsertKnots
1726 //=======================================================================
1728 void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1729 const Standard_Integer Degree,
1730 const Standard_Boolean Periodic,
1731 const TColgp_Array2OfPnt& Poles,
1732 const TColStd_Array2OfReal& Weights,
1733 const TColStd_Array1OfReal& Knots,
1734 const TColStd_Array1OfInteger& Mults,
1735 const TColStd_Array1OfReal& AddKnots,
1736 const TColStd_Array1OfInteger& AddMults,
1737 TColgp_Array2OfPnt& NewPoles,
1738 TColStd_Array2OfReal& NewWeights,
1739 TColStd_Array1OfReal& NewKnots,
1740 TColStd_Array1OfInteger& NewMults,
1741 const Standard_Real Epsilon,
1742 const Standard_Boolean Add )
1744 Standard_Boolean rational = &Weights != NULL;
1745 Standard_Integer dim = 3;
1746 if (rational) dim++;
1748 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1749 TColStd_Array1OfReal
1750 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1752 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1753 else SetPoles(Poles,poles,UDirection);
1756 dim *= Poles.RowLength();
1759 dim *= Poles.ColLength();
1761 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1762 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1765 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1766 else GetPoles(newpoles,NewPoles,UDirection);
1769 //=======================================================================
1770 //function : RemoveKnot
1772 //=======================================================================
1774 Standard_Boolean BSplSLib::RemoveKnot
1775 (const Standard_Boolean UDirection,
1776 const Standard_Integer Index,
1777 const Standard_Integer Mult,
1778 const Standard_Integer Degree,
1779 const Standard_Boolean Periodic,
1780 const TColgp_Array2OfPnt& Poles,
1781 const TColStd_Array2OfReal& Weights,
1782 const TColStd_Array1OfReal& Knots,
1783 const TColStd_Array1OfInteger& Mults,
1784 TColgp_Array2OfPnt& NewPoles,
1785 TColStd_Array2OfReal& NewWeights,
1786 TColStd_Array1OfReal& NewKnots,
1787 TColStd_Array1OfInteger& NewMults,
1788 const Standard_Real Tolerance)
1790 Standard_Boolean rational = &Weights != NULL;
1791 Standard_Integer dim = 3;
1792 if (rational) dim++;
1794 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1795 TColStd_Array1OfReal
1796 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1798 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1799 else SetPoles(Poles,poles,UDirection);
1802 dim *= Poles.RowLength();
1805 dim *= Poles.ColLength();
1808 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1809 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1811 return Standard_False;
1813 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1814 else GetPoles(newpoles,NewPoles,UDirection);
1815 return Standard_True;
1818 //=======================================================================
1819 //function : IncreaseDegree
1821 //=======================================================================
1823 void BSplSLib::IncreaseDegree
1824 (const Standard_Boolean UDirection,
1825 const Standard_Integer Degree,
1826 const Standard_Integer NewDegree,
1827 const Standard_Boolean Periodic,
1828 const TColgp_Array2OfPnt& Poles,
1829 const TColStd_Array2OfReal& Weights,
1830 const TColStd_Array1OfReal& Knots,
1831 const TColStd_Array1OfInteger& Mults,
1832 TColgp_Array2OfPnt& NewPoles,
1833 TColStd_Array2OfReal& NewWeights,
1834 TColStd_Array1OfReal& NewKnots,
1835 TColStd_Array1OfInteger& NewMults)
1837 Standard_Boolean rational = &Weights != NULL;
1838 Standard_Integer dim = 3;
1839 if (rational) dim++;
1841 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1842 TColStd_Array1OfReal
1843 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1845 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1846 else SetPoles(Poles,poles,UDirection);
1849 dim *= Poles.RowLength();
1852 dim *= Poles.ColLength();
1855 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1856 newpoles,NewKnots,NewMults);
1858 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1859 else GetPoles(newpoles,NewPoles,UDirection);
1862 //=======================================================================
1863 //function : Unperiodize
1865 //=======================================================================
1867 void BSplSLib::Unperiodize
1868 (const Standard_Boolean UDirection,
1869 const Standard_Integer Degree,
1870 const TColStd_Array1OfInteger& Mults,
1871 const TColStd_Array1OfReal& Knots,
1872 const TColgp_Array2OfPnt& Poles,
1873 const TColStd_Array2OfReal& Weights,
1874 TColStd_Array1OfInteger& NewMults,
1875 TColStd_Array1OfReal& NewKnots,
1876 TColgp_Array2OfPnt& NewPoles,
1877 TColStd_Array2OfReal& NewWeights)
1879 Standard_Boolean rational = &Weights != NULL;
1880 Standard_Integer dim = 3;
1881 if (rational) dim++;
1883 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1884 TColStd_Array1OfReal
1885 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1887 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1888 else SetPoles(Poles,poles,UDirection);
1891 dim *= Poles.RowLength();
1894 dim *= Poles.ColLength();
1897 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1898 NewMults, NewKnots, newpoles);
1900 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1901 else GetPoles(newpoles,NewPoles,UDirection);
1904 //=======================================================================
1905 //function : BuildCache
1906 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1907 // Cache : in case of a rational function the Poles are
1908 // stored in homogeneous form
1909 //=======================================================================
1911 void BSplSLib::BuildCache
1912 (const Standard_Real U,
1913 const Standard_Real V,
1914 const Standard_Real USpanDomain,
1915 const Standard_Real VSpanDomain,
1916 const Standard_Boolean UPeriodic,
1917 const Standard_Boolean VPeriodic,
1918 const Standard_Integer UDegree,
1919 const Standard_Integer VDegree,
1920 const Standard_Integer UIndex,
1921 const Standard_Integer VIndex,
1922 const TColStd_Array1OfReal& UFlatKnots,
1923 const TColStd_Array1OfReal& VFlatKnots,
1924 const TColgp_Array2OfPnt& Poles,
1925 const TColStd_Array2OfReal& Weights,
1926 TColgp_Array2OfPnt& CachePoles,
1927 TColStd_Array2OfReal& CacheWeights)
1929 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1930 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1931 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1932 if (&Weights != NULL)
1933 rational_u = rational_v = Standard_True;
1935 rational_u = rational_v = Standard_False;
1936 BSplSLib_DataContainer dc (UDegree, VDegree);
1952 (BSplCLib::NoMults()),
1953 (BSplCLib::NoMults()),
1963 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1965 for (kk = 0; kk <= d1 ; kk++)
1966 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1968 min_degree_domain = USpanDomain ;
1969 max_degree_domain = VSpanDomain ;
1972 min_degree_domain = VSpanDomain ;
1973 max_degree_domain = USpanDomain ;
1977 for (ii = 0 ; ii <= d2 ; ii++) {
1981 for (jj = 0 ; jj <= d1 ; jj++) {
1983 Index = jj * d2p1 + ii ;
1985 gp_Pnt& P = CachePoles(iii,jjj);
1986 f = factor[0] * factor[1];
1987 P.SetX( f * dc.poles[Index]); Index++;
1988 P.SetY( f * dc.poles[Index]); Index++;
1989 P.SetZ( f * dc.poles[Index]); Index++;
1990 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1991 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1993 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1997 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1999 for (kk = 0; kk <= d1 ; kk++)
2000 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
2002 min_degree_domain = USpanDomain ;
2003 max_degree_domain = VSpanDomain ;
2006 min_degree_domain = VSpanDomain ;
2007 max_degree_domain = USpanDomain ;
2011 for (ii = 0 ; ii <= d2 ; ii++) {
2015 for (jj = 0 ; jj <= d1 ; jj++) {
2017 Index = jj * d2p1 + ii ;
2018 Index = (Index << 1) + Index;
2019 gp_Pnt& P = CachePoles(iii,jjj);
2020 f = factor[0] * factor[1];
2021 P.SetX( f * dc.poles[Index]); Index++;
2022 P.SetY( f * dc.poles[Index]); Index++;
2023 P.SetZ( f * dc.poles[Index]);
2024 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
2026 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
2028 if (&Weights != NULL) {
2030 // means that PrepareEval did found out that the surface was
2031 // locally polynomial but since the surface is constructed
2032 // with some weights we need to set the weight polynomial to constant
2035 for (ii = 1 ; ii <= d2p1 ; ii++) {
2037 for (jj = 1 ; jj <= d1p1 ; jj++) {
2038 CacheWeights(ii,jj) = 0.0e0 ;
2041 CacheWeights(1,1) = 1.0e0 ;
2046 void BSplSLib::BuildCache(const Standard_Real theU,
2047 const Standard_Real theV,
2048 const Standard_Real theUSpanDomain,
2049 const Standard_Real theVSpanDomain,
2050 const Standard_Boolean theUPeriodicFlag,
2051 const Standard_Boolean theVPeriodicFlag,
2052 const Standard_Integer theUDegree,
2053 const Standard_Integer theVDegree,
2054 const Standard_Integer theUIndex,
2055 const Standard_Integer theVIndex,
2056 const TColStd_Array1OfReal& theUFlatKnots,
2057 const TColStd_Array1OfReal& theVFlatKnots,
2058 const TColgp_Array2OfPnt& thePoles,
2059 const TColStd_Array2OfReal& theWeights,
2060 TColStd_Array2OfReal& theCacheArray)
2062 Standard_Boolean flag_u_or_v;
2063 Standard_Integer d1, d2;
2064 Standard_Real u1, u2;
2065 Standard_Boolean isRationalOnParam = (&theWeights != NULL);
2066 Standard_Boolean isRational;
2068 BSplSLib_DataContainer dc(theUDegree, theVDegree);
2070 PrepareEval(theU, theV, theUIndex, theVIndex, theUDegree, theVDegree,
2071 isRationalOnParam, isRationalOnParam,
2072 theUPeriodicFlag, theVPeriodicFlag,
2073 thePoles, theWeights,
2074 theUFlatKnots, theVFlatKnots,
2075 (BSplCLib::NoMults()), (BSplCLib::NoMults()),
2076 u1, u2, d1, d2, isRational, dc);
2078 Standard_Integer d2p1 = d2 + 1;
2079 Standard_Integer aDimension = isRational ? 4 : 3;
2080 Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the surface is locally polynomial
2081 (isRationalOnParam && !isRational) ? aDimension + 1 : aDimension;
2083 Standard_Real aDomains[2];
2084 // aDomains[0] corresponds to variable with minimal degree
2085 // aDomains[1] corresponds to variable with maximal degree
2088 aDomains[0] = theUSpanDomain;
2089 aDomains[1] = theVSpanDomain;
2093 aDomains[0] = theVSpanDomain;
2094 aDomains[1] = theUSpanDomain;
2097 BSplCLib::Bohm(u1, d1, d1, *dc.knots1, aDimension * d2p1, *dc.poles);
2098 for (Standard_Integer kk = 0; kk <= d1 ; kk++)
2099 BSplCLib::Bohm(u2, d2, d2, *dc.knots2, aDimension, *(dc.poles + kk * aDimension * d2p1));
2101 Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
2102 Standard_Real* aPolyCoeffs = dc.poles;
2104 Standard_Real aFactors[2];
2105 // aFactors[0] corresponds to variable with minimal degree
2106 // aFactors[1] corresponds to variable with maximal degree
2108 Standard_Integer aRow, aCol, i;
2109 Standard_Real aCoeff;
2110 for (aRow = 0; aRow <= d2; aRow++)
2113 for (aCol = 0; aCol <= d1; aCol++)
2115 aPolyCoeffs = dc.poles + (aCol * d2p1 + aRow) * aDimension;
2116 aCoeff = aFactors[0] * aFactors[1];
2117 for (i = 0; i < aDimension; i++)
2118 aCache[i] = aPolyCoeffs[i] * aCoeff;
2119 aCache += aCacheShift;
2120 aFactors[0] *= aDomains[0] / (aCol + 1);
2122 aFactors[1] *= aDomains[1] / (aRow + 1);
2125 // Fill the weights for the surface which is not locally polynomial
2126 if (aCacheShift > aDimension)
2128 aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
2129 aCache += aCacheShift - 1;
2130 for (aRow = 0; aRow <= d2; aRow++)
2131 for (aCol = 0; aCol <= d1; aCol++)
2134 aCache += aCacheShift;
2136 theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
2141 //=======================================================================
2142 //function : CacheD0
2143 //purpose : Evaluates the polynomial cache of the Bspline Curve
2145 //=======================================================================
2147 void BSplSLib::CacheD0(const Standard_Real UParameter,
2148 const Standard_Real VParameter,
2149 const Standard_Integer UDegree,
2150 const Standard_Integer VDegree,
2151 const Standard_Real UCacheParameter,
2152 const Standard_Real VCacheParameter,
2153 const Standard_Real USpanLenght,
2154 const Standard_Real VSpanLenght,
2155 const TColgp_Array2OfPnt& PolesArray,
2156 const TColStd_Array2OfReal& WeightsArray,
2160 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2162 // the SpanLenght is the normalizing factor so that the polynomial is between
2175 PArray = (Standard_Real *)
2176 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2178 myPoint = (Standard_Real *) &aPoint ;
2179 if (UDegree <= VDegree) {
2180 min_degree = UDegree ;
2181 max_degree = VDegree ;
2182 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
2183 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
2184 dimension = 3 * (UDegree + 1) ;
2187 min_degree = VDegree ;
2188 max_degree = UDegree ;
2189 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
2190 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
2191 dimension = 3 * (VDegree + 1) ;
2193 NCollection_LocalArray<Standard_Real> locpoles(dimension);
2195 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2198 max_degree*dimension,
2202 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2205 (min_degree << 1) + min_degree,
2208 if (&WeightsArray != NULL) {
2209 dimension = min_degree + 1 ;
2211 WArray = (Standard_Real *)
2212 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2213 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2216 max_degree*dimension,
2220 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2226 inverse = 1.0e0/ inverse ;
2228 myPoint[0] *= inverse ;
2229 myPoint[1] *= inverse ;
2230 myPoint[2] *= inverse ;
2234 //=======================================================================
2235 //function : CacheD1
2236 //purpose : Evaluates the polynomial cache of the Bspline Curve
2238 //=======================================================================
2240 void BSplSLib::CacheD1(const Standard_Real UParameter,
2241 const Standard_Real VParameter,
2242 const Standard_Integer UDegree,
2243 const Standard_Integer VDegree,
2244 const Standard_Real UCacheParameter,
2245 const Standard_Real VCacheParameter,
2246 const Standard_Real USpanLenght,
2247 const Standard_Real VSpanLenght,
2248 const TColgp_Array2OfPnt& PolesArray,
2249 const TColStd_Array2OfReal& WeightsArray,
2255 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2257 // the SpanLenght is the normalizing factor so that the polynomial is between
2273 PArray = (Standard_Real *)
2274 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2275 Standard_Real local_poles_array[2][2][3],
2276 local_poles_and_weights_array[2][2][4],
2277 local_weights_array[2][2] ;
2278 Standard_Real * my_vec_min,
2281 my_point = (Standard_Real *) &aPoint ;
2283 // initialize in case of rational evaluation
2284 // because RationalDerivative will use all
2288 if (&WeightsArray != NULL) {
2290 local_poles_array [0][0][0] = 0.0e0 ;
2291 local_poles_array [0][0][1] = 0.0e0 ;
2292 local_poles_array [0][0][2] = 0.0e0 ;
2293 local_weights_array [0][0] = 0.0e0 ;
2294 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2295 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2296 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2297 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2299 local_poles_array [0][1][0] = 0.0e0 ;
2300 local_poles_array [0][1][1] = 0.0e0 ;
2301 local_poles_array [0][1][2] = 0.0e0 ;
2302 local_weights_array [0][1] = 0.0e0 ;
2303 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2304 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2305 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2306 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2308 local_poles_array [1][0][0] = 0.0e0 ;
2309 local_poles_array [1][0][1] = 0.0e0 ;
2310 local_poles_array [1][0][2] = 0.0e0 ;
2311 local_weights_array [1][0] = 0.0e0 ;
2312 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2313 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2314 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2315 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2317 local_poles_array [1][1][0] = 0.0e0 ;
2318 local_poles_array [1][1][1] = 0.0e0 ;
2319 local_poles_array [1][1][2] = 0.0e0 ;
2320 local_weights_array [1][1] = 0.0e0 ;
2321 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2322 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2323 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2324 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2327 if (UDegree <= VDegree) {
2328 min_degree = UDegree ;
2329 max_degree = VDegree ;
2330 inverse_min = 1.0e0/USpanLenght ;
2331 inverse_max = 1.0e0/VSpanLenght ;
2332 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2333 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2335 dimension = 3 * (UDegree + 1) ;
2336 my_vec_min = (Standard_Real *) &aVecU ;
2337 my_vec_max = (Standard_Real *) &aVecV ;
2340 min_degree = VDegree ;
2341 max_degree = UDegree ;
2342 inverse_min = 1.0e0/VSpanLenght ;
2343 inverse_max = 1.0e0/USpanLenght ;
2344 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2345 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2346 dimension = 3 * (VDegree + 1) ;
2347 my_vec_min = (Standard_Real *) &aVecV ;
2348 my_vec_max = (Standard_Real *) &aVecU ;
2351 NCollection_LocalArray<Standard_Real> locpoles (2 * dimension);
2353 PLib::EvalPolynomial(new_parameter[0],
2360 PLib::EvalPolynomial(new_parameter[1],
2365 local_poles_array[0][0][0]) ;
2366 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2369 (min_degree << 1) + min_degree,
2370 locpoles[dimension],
2371 local_poles_array[1][0][0]) ;
2373 if (&WeightsArray != NULL) {
2374 dimension = min_degree + 1 ;
2376 WArray = (Standard_Real *)
2377 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2378 PLib::EvalPolynomial(new_parameter[0],
2385 PLib::EvalPolynomial(new_parameter[1],
2390 local_weights_array[0][0]) ;
2391 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2395 locpoles[dimension],
2396 local_weights_array[1][0]) ;
2398 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2399 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2400 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2401 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2403 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2404 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2405 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2406 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2408 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2409 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2410 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2411 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2413 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2414 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2415 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2416 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2418 BSplSLib::RationalDerivative(1,
2422 local_poles_and_weights_array[0][0][0],
2423 local_poles_array[0][0][0]) ;
2426 my_point [0] = local_poles_array [0][0][0] ;
2427 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2428 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2430 my_point [1] = local_poles_array [0][0][1] ;
2431 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2432 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2434 my_point [2] = local_poles_array [0][0][2] ;
2435 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2436 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2439 //=======================================================================
2440 //function : CacheD2
2441 //purpose : Evaluates the polynomial cache of the Bspline Curve
2443 //=======================================================================
2445 void BSplSLib::CacheD2(const Standard_Real UParameter,
2446 const Standard_Real VParameter,
2447 const Standard_Integer UDegree,
2448 const Standard_Integer VDegree,
2449 const Standard_Real UCacheParameter,
2450 const Standard_Real VCacheParameter,
2451 const Standard_Real USpanLenght,
2452 const Standard_Real VSpanLenght,
2453 const TColgp_Array2OfPnt& PolesArray,
2454 const TColStd_Array2OfReal& WeightsArray,
2463 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2465 // the SpanLenght is the normalizing factor so that the polynomial is between
2482 PArray = (Standard_Real *)
2483 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2484 Standard_Real local_poles_array[3][3][3],
2485 local_poles_and_weights_array[3][3][4],
2486 local_weights_array[3][3] ;
2487 Standard_Real * my_vec_min,
2493 my_point = (Standard_Real *) &aPoint ;
2496 // initialize in case the min and max degree are less than 2
2498 local_poles_array[0][0][0] = 0.0e0 ;
2499 local_poles_array[0][0][1] = 0.0e0 ;
2500 local_poles_array[0][0][2] = 0.0e0 ;
2501 local_poles_array[0][1][0] = 0.0e0 ;
2502 local_poles_array[0][1][1] = 0.0e0 ;
2503 local_poles_array[0][1][2] = 0.0e0 ;
2504 local_poles_array[0][2][0] = 0.0e0 ;
2505 local_poles_array[0][2][1] = 0.0e0 ;
2506 local_poles_array[0][2][2] = 0.0e0 ;
2508 local_poles_array[1][0][0] = 0.0e0 ;
2509 local_poles_array[1][0][1] = 0.0e0 ;
2510 local_poles_array[1][0][2] = 0.0e0 ;
2511 local_poles_array[1][1][0] = 0.0e0 ;
2512 local_poles_array[1][1][1] = 0.0e0 ;
2513 local_poles_array[1][1][2] = 0.0e0 ;
2514 local_poles_array[1][2][0] = 0.0e0 ;
2515 local_poles_array[1][2][1] = 0.0e0 ;
2516 local_poles_array[1][2][2] = 0.0e0 ;
2518 local_poles_array[2][0][0] = 0.0e0 ;
2519 local_poles_array[2][0][1] = 0.0e0 ;
2520 local_poles_array[2][0][2] = 0.0e0 ;
2521 local_poles_array[2][1][0] = 0.0e0 ;
2522 local_poles_array[2][1][1] = 0.0e0 ;
2523 local_poles_array[2][1][2] = 0.0e0 ;
2524 local_poles_array[2][2][0] = 0.0e0 ;
2525 local_poles_array[2][2][1] = 0.0e0 ;
2526 local_poles_array[2][2][2] = 0.0e0 ;
2528 // initialize in case of rational evaluation
2529 // because RationalDerivative will use all
2533 if (&WeightsArray != NULL) {
2535 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2536 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2537 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2538 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2539 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2540 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2541 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2542 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2543 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2545 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2546 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2547 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2548 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2549 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2550 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2551 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2552 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2553 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2555 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2556 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2557 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2558 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2559 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2560 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2561 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2562 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2563 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2565 local_poles_and_weights_array[0][0][3] =
2566 local_weights_array[0][0] = 0.0e0 ;
2567 local_poles_and_weights_array[0][1][3] =
2568 local_weights_array[0][1] = 0.0e0 ;
2569 local_poles_and_weights_array[0][2][3] =
2570 local_weights_array[0][2] = 0.0e0 ;
2571 local_poles_and_weights_array[1][0][3] =
2572 local_weights_array[1][0] = 0.0e0 ;
2573 local_poles_and_weights_array[1][1][3] =
2574 local_weights_array[1][1] = 0.0e0 ;
2575 local_poles_and_weights_array[1][2][3] =
2576 local_weights_array[1][2] = 0.0e0 ;
2577 local_poles_and_weights_array[2][0][3] =
2578 local_weights_array[2][0] = 0.0e0 ;
2579 local_poles_and_weights_array[2][1][3] =
2580 local_weights_array[2][1] = 0.0e0 ;
2581 local_poles_and_weights_array[2][2][3] =
2582 local_weights_array[2][2] = 0.0e0 ;
2585 if (UDegree <= VDegree) {
2586 min_degree = UDegree ;
2587 max_degree = VDegree ;
2588 inverse_min = 1.0e0/USpanLenght ;
2589 inverse_max = 1.0e0/VSpanLenght ;
2590 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2591 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2593 dimension = 3 * (UDegree + 1) ;
2594 my_vec_min = (Standard_Real *) &aVecU ;
2595 my_vec_max = (Standard_Real *) &aVecV ;
2596 my_vec_min_min = (Standard_Real *) &aVecUU ;
2597 my_vec_min_max = (Standard_Real *) &aVecUV ;
2598 my_vec_max_max = (Standard_Real *) &aVecVV ;
2601 min_degree = VDegree ;
2602 max_degree = UDegree ;
2603 inverse_min = 1.0e0/VSpanLenght ;
2604 inverse_max = 1.0e0/USpanLenght ;
2605 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2606 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2607 dimension = 3 * (VDegree + 1) ;
2608 my_vec_min = (Standard_Real *) &aVecV ;
2609 my_vec_max = (Standard_Real *) &aVecU ;
2610 my_vec_min_min = (Standard_Real *) &aVecVV ;
2611 my_vec_min_max = (Standard_Real *) &aVecUV ;
2612 my_vec_max_max = (Standard_Real *) &aVecUU ;
2615 NCollection_LocalArray<Standard_Real> locpoles (3 * dimension);
2618 // initialize in case min or max degree are less than 2
2620 Standard_Integer MinIndMax = 2;
2621 if ( max_degree < 2) MinIndMax = max_degree;
2622 Standard_Integer MinIndMin = 2;
2623 if ( min_degree < 2) MinIndMin = min_degree;
2625 index = MinIndMax * dimension ;
2627 for (ii = MinIndMax ; ii < 3 ; ii++) {
2629 for (kk = 0 ; kk < dimension ; kk++) {
2630 locpoles[index] = 0.0e0 ;
2635 PLib::EvalPolynomial(new_parameter[0],
2642 PLib::EvalPolynomial(new_parameter[1],
2647 local_poles_array[0][0][0]) ;
2648 PLib::EvalPolynomial(new_parameter[1],
2652 locpoles[dimension],
2653 local_poles_array[1][0][0]) ;
2655 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2658 (min_degree << 1) + min_degree,
2659 locpoles[dimension + dimension],
2660 local_poles_array[2][0][0]) ;
2662 if (&WeightsArray != NULL) {
2663 dimension = min_degree + 1 ;
2665 WArray = (Standard_Real *)
2666 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2667 PLib::EvalPolynomial(new_parameter[0],
2674 PLib::EvalPolynomial(new_parameter[1],
2679 local_weights_array[0][0]) ;
2680 PLib::EvalPolynomial(new_parameter[1],
2684 locpoles[dimension],
2685 local_weights_array[1][0]) ;
2686 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2690 locpoles[dimension + dimension],
2691 local_weights_array[2][0]) ;
2694 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2695 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2696 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2697 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2698 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2699 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2700 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2701 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2702 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2704 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2705 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2706 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2707 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2708 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2709 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2710 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2711 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2712 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2714 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2715 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2716 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2717 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2718 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2719 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2720 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2721 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2722 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2725 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2726 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2727 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2728 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2729 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2730 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2731 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2732 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2733 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2735 BSplSLib::RationalDerivative(2,
2739 local_poles_and_weights_array[0][0][0],
2740 local_poles_array[0][0][0]) ;
2744 Standard_Real minmin = inverse_min * inverse_min;
2745 Standard_Real minmax = inverse_min * inverse_max;
2746 Standard_Real maxmax = inverse_max * inverse_max;
2748 my_point [0] = local_poles_array [0][0][0] ;
2749 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2750 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2751 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2752 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2753 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2755 my_point [1] = local_poles_array [0][0][1] ;
2756 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2757 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2758 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2759 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2760 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2762 my_point [2] = local_poles_array [0][0][2] ;
2763 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2764 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2765 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2766 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2767 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2770 //=======================================================================
2771 //function : MovePoint
2772 //purpose : Find the new poles which allows an old point (with a
2773 // given u and v as parameters) to reach a new position
2774 //=======================================================================
2776 void BSplSLib::MovePoint (const Standard_Real U,
2777 const Standard_Real V,
2778 const gp_Vec& Displ,
2779 const Standard_Integer UIndex1,
2780 const Standard_Integer UIndex2,
2781 const Standard_Integer VIndex1,
2782 const Standard_Integer VIndex2,
2783 const Standard_Integer UDegree,
2784 const Standard_Integer VDegree,
2785 const Standard_Boolean Rational,
2786 const TColgp_Array2OfPnt& Poles,
2787 const TColStd_Array2OfReal& Weights,
2788 const TColStd_Array1OfReal& UFlatKnots,
2789 const TColStd_Array1OfReal& VFlatKnots,
2790 Standard_Integer& UFirstIndex,
2791 Standard_Integer& ULastIndex,
2792 Standard_Integer& VFirstIndex,
2793 Standard_Integer& VLastIndex,
2794 TColgp_Array2OfPnt& NewPoles)
2796 // calculate the UBSplineBasis in the parameter U
2797 Standard_Integer UFirstNonZeroBsplineIndex;
2798 math_Matrix UBSplineBasis(1, 1,
2800 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2805 UFirstNonZeroBsplineIndex,
2807 // calculate the VBSplineBasis in the parameter V
2808 Standard_Integer VFirstNonZeroBsplineIndex;
2809 math_Matrix VBSplineBasis(1, 1,
2811 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2816 VFirstNonZeroBsplineIndex,
2818 if (ErrorCod1 || ErrorCod2) {
2826 // find the span which is predominant for parameter U
2827 UFirstIndex = UFirstNonZeroBsplineIndex;
2828 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2829 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2830 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2832 Standard_Real maxValue = 0.0;
2833 Standard_Integer i, ukk1=0, ukk2;
2835 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2836 if (UBSplineBasis(1,i) > maxValue) {
2837 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2838 maxValue = UBSplineBasis(1,i);
2842 // find a ukk2 if symetriy
2844 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2845 if ((ukk1+1) <= ULastIndex) {
2846 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2851 // find the span which is predominant for parameter V
2852 VFirstIndex = VFirstNonZeroBsplineIndex;
2853 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2855 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2856 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2859 Standard_Integer j, vkk1=0, vkk2;
2861 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2862 if (VBSplineBasis(1,j) > maxValue) {
2863 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2864 maxValue = VBSplineBasis(1,j);
2868 // find a vkk2 if symetriy
2870 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2871 if ((vkk1+1) <= VLastIndex) {
2872 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2877 // compute the vector of displacement
2878 Standard_Real D1 = 0.0;
2879 Standard_Real D2 = 0.0;
2880 Standard_Real hN, Coef, DvalU, DvalV;
2882 Standard_Integer ii, jj;
2884 for (i = 1; i <= UDegree+1; i++) {
2885 ii = i + UFirstNonZeroBsplineIndex - 1;
2889 else if (ii > ukk2) {
2896 for (j = 1; j <= VDegree+1; j++) {
2897 jj = j + VFirstNonZeroBsplineIndex - 1;
2899 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2903 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2905 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2909 else if (jj > vkk2) {
2915 D1 += 1./(DvalU + DvalV + 1.) * hN;
2927 // compute the new poles
2929 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2933 else if (i > ukk2) {
2940 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2941 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2945 else if (j > vkk2) {
2951 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2954 NewPoles(i,j) = Poles(i,j);
2960 //=======================================================================
2961 // function : Resolution
2962 // purpose : this computes an estimate for the maximum of the
2963 // partial derivatives both in U and in V
2966 // The calculation resembles at the calculation of curves with
2967 // additional index for the control point. Let Si,j be the
2968 // control points for ls surface and Di,j the weights.
2969 // The checking of upper bounds for the partial derivatives
2970 // will be omitted and Su is the next upper bound in the polynomial case :
2974 // | Si,j - Si-1,j |
2975 // d * Max | ------------- |
2976 // i = 2,n | ti+d - ti |
2980 // and in the rational case :
2984 // Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2985 // Max Max d * -----------------------------------------------
2986 // k=1,n i dans Rj ti+d - ti
2988 // ----------------------------------------------------------------------
2996 // with Rj = {j-d, ...., j+d+d+1}.
2999 //=======================================================================
3001 void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
3002 const TColStd_Array2OfReal& Weights,
3003 const TColStd_Array1OfReal& UKnots,
3004 const TColStd_Array1OfReal& VKnots,
3005 const TColStd_Array1OfInteger& UMults,
3006 const TColStd_Array1OfInteger& VMults,
3007 const Standard_Integer UDegree,
3008 const Standard_Integer VDegree,
3009 const Standard_Boolean URational,
3010 const Standard_Boolean VRational,
3011 const Standard_Boolean UPeriodic,
3012 const Standard_Boolean VPeriodic,
3013 const Standard_Real Tolerance3D,
3014 Standard_Real& UTolerance,
3015 Standard_Real& VTolerance)
3017 Standard_Real Wij,Wmj,Wji,Wjm;
3018 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
3019 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
3020 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
3021 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
3023 max_derivative[0] = max_derivative[1] = 0.0e0 ;
3025 Standard_Integer PRowLength, PColLength;
3026 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
3027 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
3028 Standard_Integer num_poles[2],num_flat_knots[2];
3031 BSplCLib::KnotSequenceLength(UMults,
3035 BSplCLib::KnotSequenceLength(VMults,
3038 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
3039 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
3040 BSplCLib::KnotSequence(UKnots,
3045 BSplCLib::KnotSequence(VKnots,
3050 PRowLength = Poles.RowLength();
3051 PColLength = Poles.ColLength();
3052 if (URational || VRational) {
3053 Standard_Integer Wsize = PRowLength * PColLength;
3054 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
3055 min_weights = WG[0];
3057 for (ii = 1 ; ii < Wsize ; ii++) {
3059 if (min_weights > min) min_weights = min;
3062 Standard_Integer UD1 = UDegree + 1;
3063 Standard_Integer VD1 = VDegree + 1;
3064 num_poles[0] = num_flat_knots[0] - UD1;
3065 num_poles[1] = num_flat_knots[1] - VD1;
3066 poles_length[0] = PColLength;
3067 poles_length[1] = PRowLength;
3069 Standard_Integer UD2 = UDegree << 1;
3070 Standard_Integer VD2 = VDegree << 1;
3072 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3073 ii_index = (ii - 1) % poles_length[0] + 1 ;
3074 ii_minus = (ii - 2) % poles_length[0] + 1 ;
3075 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3076 inverse = 1.0e0 / inverse ;
3077 lower[0] = ii - UD1;
3078 if (lower[0] < 1) lower[0] = 1;
3079 upper[0] = ii + UD2 + 1;
3080 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
3082 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3083 jj_index = (jj - 1) % poles_length[1] + 1 ;
3084 lower[1] = jj - VD1;
3085 if (lower[1] < 1) lower[1] = 1;
3086 upper[1] = jj + VD2 + 1;
3087 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
3088 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
3089 Wij = Weights.Value(ii_index,jj_index);
3090 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
3091 Wmj = Weights.Value(ii_minus,jj_index);
3099 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
3100 pp_index = (pp - 1) % poles_length[0] + 1 ;
3102 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
3104 qq_index = (qq - 1) % poles_length[1] + 1 ;
3105 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
3109 factor = (Xpq - Xij) * Wij;
3110 factor -= (Xpq - Xmj) * Wmj;
3111 if (factor < 0) factor = - factor;
3113 factor = (Ypq - Yij) * Wij;
3114 factor -= (Ypq - Ymj) * Wmj;
3115 if (factor < 0) factor = - factor;
3117 factor = (Zpq - Zij) * Wij;
3118 factor -= (Zpq - Zmj) * Wmj;
3119 if (factor < 0) factor = - factor;
3122 if (max_derivative[0] < value) max_derivative[0] = value ;
3127 max_derivative[0] /= min_weights ;
3131 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3132 ii_index = (ii - 1) % poles_length[0] + 1 ;
3133 ii_minus = (ii - 2) % poles_length[0] + 1 ;
3134 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3135 inverse = 1.0e0 / inverse ;
3137 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3138 jj_index = (jj - 1) % poles_length[1] + 1 ;
3140 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
3141 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
3142 factor = Pij.X() - Pmj.X();
3143 if (factor < 0) factor = - factor;
3145 factor = Pij.Y() - Pmj.Y();
3146 if (factor < 0) factor = - factor;
3148 factor = Pij.Z() - Pmj.Z();
3149 if (factor < 0) factor = - factor;
3152 if (max_derivative[0] < value) max_derivative[0] = value ;
3156 max_derivative[0] *= UDegree ;
3158 Standard_Integer UD2 = UDegree << 1;
3159 Standard_Integer VD2 = VDegree << 1;
3161 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3162 ii_index = (ii - 1) % poles_length[1] + 1 ;
3163 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3164 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3165 inverse = 1.0e0 / inverse ;
3166 lower[0] = ii - VD1;
3167 if (lower[0] < 1) lower[0] = 1;
3168 upper[0] = ii + VD2 + 1;
3169 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
3171 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3172 jj_index = (jj - 1) % poles_length[0] + 1 ;
3173 lower[1] = jj - UD1;
3174 if (lower[1] < 1) lower[1] = 1;
3175 upper[1] = jj + UD2 + 1;
3176 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
3177 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
3178 Wji = Weights.Value(jj_index,ii_index);
3179 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
3180 Wjm = Weights.Value(jj_index,ii_minus);
3188 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
3189 pp_index = (pp - 1) % poles_length[1] + 1 ;
3191 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
3193 qq_index = (qq - 1) % poles_length[0] + 1 ;
3194 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
3198 factor = (Xqp - Xji) * Wji;
3199 factor -= (Xqp - Xjm) * Wjm;
3200 if (factor < 0) factor = - factor;
3202 factor = (Yqp - Yji) * Wji;
3203 factor -= (Yqp - Yjm) * Wjm;
3204 if (factor < 0) factor = - factor;
3206 factor = (Zqp - Zji) * Wji;
3207 factor -= (Zqp - Zjm) * Wjm;
3208 if (factor < 0) factor = - factor;
3211 if (max_derivative[1] < value) max_derivative[1] = value ;
3216 max_derivative[1] /= min_weights ;
3220 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3221 ii_index = (ii - 1) % poles_length[1] + 1 ;
3222 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3223 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3224 inverse = 1.0e0 / inverse ;
3226 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3227 jj_index = (jj - 1) % poles_length[0] + 1 ;
3229 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3230 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3231 factor = Pji.X() - Pjm.X() ;
3232 if (factor < 0) factor = - factor;
3234 factor = Pji.Y() - Pjm.Y() ;
3235 if (factor < 0) factor = - factor;
3237 factor = Pji.Z() - Pjm.Z() ;
3238 if (factor < 0) factor = - factor;
3241 if (max_derivative[1] < value) max_derivative[1] = value ;
3245 max_derivative[1] *= VDegree ;
3246 max_derivative[0] *= M_SQRT2 ;
3247 max_derivative[1] *= M_SQRT2 ;
3248 if(max_derivative[0] && max_derivative[1]) {
3249 UTolerance = Tolerance3D / max_derivative[0] ;
3250 VTolerance = Tolerance3D / max_derivative[1] ;
3253 UTolerance=VTolerance=0.0;
3255 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3260 //=======================================================================
3261 //function : Interpolate
3263 //=======================================================================
3265 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3266 const Standard_Integer VDegree,
3267 const TColStd_Array1OfReal& UFlatKnots,
3268 const TColStd_Array1OfReal& VFlatKnots,
3269 const TColStd_Array1OfReal& UParameters,
3270 const TColStd_Array1OfReal& VParameters,
3271 TColgp_Array2OfPnt& Poles,
3272 TColStd_Array2OfReal& Weights,
3273 Standard_Integer& InversionProblem)
3275 Standard_Integer ii, jj, ll, kk, dimension;
3276 Standard_Integer ULength = UParameters.Length();
3277 Standard_Integer VLength = VParameters.Length();
3278 Standard_Real * poles_array;
3280 // extraction of iso u
3281 dimension = 4*ULength;
3282 TColStd_Array2OfReal Points(1, VLength,
3285 Handle(TColStd_HArray1OfInteger) ContactOrder =
3286 new (TColStd_HArray1OfInteger)(1, VLength);
3287 ContactOrder->Init(0);
3289 for (ii=1; ii <= VLength; ii++) {
3291 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3292 Points(ii,ll) = Poles(jj, ii).X();
3293 Points(ii,ll+1) = Poles(jj, ii).Y();
3294 Points(ii,ll+2) = Poles(jj, ii).Z();
3295 Points(ii,ll+3) = Weights(jj,ii) ;
3299 // interpolation of iso u
3300 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3301 BSplCLib::Interpolate(VDegree,
3304 ContactOrder->Array1(),
3308 if (InversionProblem != 0) return;
3310 // extraction of iso v
3312 dimension = VLength*4;
3313 TColStd_Array2OfReal IsoPoles(1, ULength,
3316 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3317 ContactOrder->Init(0);
3318 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3320 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3322 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3323 IsoPoles (ii,ll) = Points(jj, kk);
3324 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3325 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3326 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3329 // interpolation of iso v
3330 BSplCLib::Interpolate(UDegree,
3333 ContactOrder->Array1(),
3340 for (ii=1; ii <= ULength; ii++) {
3342 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3343 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3344 Poles.SetValue(ii, jj, Pnt);
3345 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3350 //=======================================================================
3351 //function : Interpolate
3353 //=======================================================================
3355 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3356 const Standard_Integer VDegree,
3357 const TColStd_Array1OfReal& UFlatKnots,
3358 const TColStd_Array1OfReal& VFlatKnots,
3359 const TColStd_Array1OfReal& UParameters,
3360 const TColStd_Array1OfReal& VParameters,
3361 TColgp_Array2OfPnt& Poles,
3362 Standard_Integer& InversionProblem)
3364 Standard_Integer ii, jj, ll, kk, dimension;
3365 Standard_Integer ULength = UParameters.Length();
3366 Standard_Integer VLength = VParameters.Length();
3367 Standard_Real * poles_array;
3369 // extraction of iso u
3370 dimension = 3*ULength;
3371 TColStd_Array2OfReal Points(1, VLength,
3374 Handle(TColStd_HArray1OfInteger) ContactOrder =
3375 new (TColStd_HArray1OfInteger)(1, VLength);
3376 ContactOrder->Init(0);
3378 for (ii=1; ii <= VLength; ii++) {
3380 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3381 Points(ii,ll) = Poles(jj, ii).X();
3382 Points(ii,ll+1) = Poles(jj, ii).Y();
3383 Points(ii,ll+2) = Poles(jj, ii).Z();
3387 // interpolation of iso u
3388 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3389 BSplCLib::Interpolate(VDegree,
3392 ContactOrder->Array1(),
3396 if (InversionProblem != 0) return;
3398 // extraction of iso v
3400 dimension = VLength*3;
3401 TColStd_Array2OfReal IsoPoles(1, ULength,
3404 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3405 ContactOrder->Init(0);
3406 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3408 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3410 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3411 IsoPoles (ii,ll) = Points(jj, kk);
3412 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3413 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3416 // interpolation of iso v
3417 BSplCLib::Interpolate(UDegree,
3420 ContactOrder->Array1(),
3427 for (ii=1; ii <= ULength; ii++) {
3429 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3430 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3431 Poles.SetValue(ii, jj, Pnt);
3436 //=======================================================================
3437 //function : FunctionMultiply
3439 //=======================================================================
3441 void BSplSLib::FunctionMultiply
3442 (const BSplSLib_EvaluatorFunction& Function,
3443 const Standard_Integer UBSplineDegree,
3444 const Standard_Integer VBSplineDegree,
3445 const TColStd_Array1OfReal& UBSplineKnots,
3446 const TColStd_Array1OfReal& VBSplineKnots,
3447 const TColStd_Array1OfInteger & UMults,
3448 const TColStd_Array1OfInteger & VMults,
3449 const TColgp_Array2OfPnt& Poles,
3450 const TColStd_Array2OfReal& Weights,
3451 const TColStd_Array1OfReal& UFlatKnots,
3452 const TColStd_Array1OfReal& VFlatKnots,
3453 const Standard_Integer UNewDegree,
3454 const Standard_Integer VNewDegree,
3455 TColgp_Array2OfPnt& NewNumerator,
3456 TColStd_Array2OfReal& NewDenominator,
3457 Standard_Integer& Status)
3459 Standard_Integer num_uparameters,
3464 Standard_Real result ;
3466 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3467 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3468 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3469 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3471 if ((NewNumerator.ColLength() == num_uparameters) &&
3472 (NewNumerator.RowLength() == num_vparameters) &&
3473 (NewDenominator.ColLength() == num_uparameters) &&
3474 (NewDenominator.RowLength() == num_vparameters)) {
3477 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3481 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3485 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3487 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3488 HomogeneousD0(UParameters(ii),
3504 NewDenominator(ii,jj),
3505 NewNumerator(ii,jj)) ;
3507 Function.Evaluate (0,
3513 Standard_ConstructionError::Raise();
3515 gp_Pnt& P = NewNumerator(ii,jj);
3516 P.SetX(P.X() * result);
3517 P.SetY(P.Y() * result);
3518 P.SetZ(P.Z() * result);
3519 NewDenominator(ii,jj) *= result ;
3522 Interpolate(UNewDegree,
3533 Standard_ConstructionError::Raise();