1 // Created on: 1991-08-26
3 // Copyright (c) 1991-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
22 // Modifed RLE Aug 93 - Complete rewrite
23 // xab 21-Mar-95 implemented cache mecanism
24 // pmn 25-09-96 Interpolation
25 // jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
26 // pmn 07-10-96 : Correction de DN dans le cas rationnal.
27 // pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
29 #include <BSplSLib.ixx>
31 #include <NCollection_LocalArray.hxx>
32 #include <BSplCLib.hxx>
33 #include <TColgp_Array2OfXYZ.hxx>
34 #include <TColgp_Array1OfXYZ.hxx>
35 #include <TColStd_HArray1OfInteger.hxx>
36 #include <Standard_NotImplemented.hxx>
37 #include <Standard_ConstructionError.hxx>
38 #include <math_Matrix.hxx>
40 // for null derivatives
41 static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
43 #define M_SQRT2 1.4142135623730950488016887
46 //=======================================================================
47 //struct : BSplCLib_DataContainer
48 //purpose: Auxiliary structure providing buffers for poles and knots used in
49 // evaluation of bspline (allocated in the stack)
50 //=======================================================================
52 struct BSplSLib_DataContainer
54 BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
56 (void)UDegree; (void)VDegree; // just to avoid compiler warning in Release mode
57 Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
58 VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
59 "BSplSLib: bspline degree is greater than maximum supported");
61 Standard_Real poles[4*(25+1)*(25+1)];
62 Standard_Real knots1[2*25];
63 Standard_Real knots2[2*25];
64 Standard_Real ders[48];
67 //**************************************************************************
69 //**************************************************************************
71 //=======================================================================
72 //function : RationalDerivative
73 //purpose : computes the rational derivatives when whe have the
74 // the derivatives of the homogeneous numerator and the
75 // the derivatives of the denominator
76 //=======================================================================
78 void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
79 const Standard_Integer VDeg,
80 const Standard_Integer N,
81 const Standard_Integer M,
82 Standard_Real& HDerivatives,
83 Standard_Real& RDerivatives,
84 const Standard_Boolean All)
87 // if All is True all derivatives are computed. if Not only
88 // the requested N, M is computed
91 // let f(u,v) be a rational function = ------------------
95 // Let (N,M) the order of the derivatives we want : then since
98 // Numerator = f * Denominator
102 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
103 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
104 // (0,0) ( p<N q<M p q )
114 // HDerivatives is an array where derivatives are stored in the following form
115 // Numerator is assumee to have 3 functions that is a vector of dimension
118 // (0,0) (0,0) (0, DegV) (0, DegV)
119 // Numerator Denominator ... Numerator Denominator
121 // (1,0) (1,0) (1, DegV) (1, DegV)
122 // Numerator Denominator ... Numerator Denominator
124 // ...........................................................
127 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
128 // Numerator Denominator ... Numerator Denominator
131 Standard_Integer ii,jj,pp,qq,index,index1,index2;
132 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
133 Standard_Integer MinN,MinN1,MinM,MinM1;
134 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
140 M4 = (VDeg + 1) << 2;
142 NCollection_LocalArray<Standard_Real> StoreDerivatives (All ? 0 : ii * 3);
143 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
144 NCollection_LocalArray<Standard_Real> StoreW (ii);
145 Standard_Real *HomogeneousArray = &HDerivatives;
146 Standard_Real denominator,Pii,Pip,Pjq;
148 denominator = 1.0e0 / HomogeneousArray[3];
151 if (UDeg < N) MinN = UDeg;
153 if (VDeg < M) MinM = VDeg;
159 for (ii = 0 ; ii < MinN1 ; ii++) {
165 for (jj = 0 ; jj < MinM1 ; jj++) {
166 RArray[index_v++] = HomogeneousArray[index_v1++];
167 RArray[index_v++] = HomogeneousArray[index_v1++];
168 RArray[index_v++] = HomogeneousArray[index_v1++];
169 StoreW[index_w++] = HomogeneousArray[index_v1++];
172 for (jj = MinM1 ; jj < M1 ; jj++) {
173 RArray[index_v++] = 0.;
174 RArray[index_v++] = 0.;
175 RArray[index_v++] = 0.;
176 StoreW[index_w++] = 0.;
181 index_v = MinN1 * M3;
182 index_w = MinN1 * M1;
184 for (ii = MinN1 ; ii < N1 ; ii++) {
186 for (jj = 0 ; jj < M1 ; jj++) {
187 RArray[index_v++] = 0.0e0;
188 RArray[index_v++] = 0.0e0;
189 RArray[index_v++] = 0.0e0;
190 StoreW[index_w++] = 0.0e0;
194 // --------------- Calculation ----------------
199 for (ii = 0 ; ii <= N ; ii++) {
205 for (jj = 0 ; jj <= M ; jj++) {
211 for (pp = 0 ; pp < ii ; pp++) {
215 index2 = jjM1 - ppM1;
216 Pip = PLib::Bin(ii,pp);
218 for (qq = 0 ; qq <= jj ; qq++) {
220 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
221 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
222 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
223 RArray[index1] -= Pjq * RArray[index]; index++;
229 Pii = PLib::Bin(ii,ii);
231 for (qq = 0 ; qq < jj ; qq++) {
233 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
234 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
235 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
236 RArray[index1] -= Pjq * RArray[index]; index++;
239 RArray[index1] *= denominator; index1++;
240 RArray[index1] *= denominator; index1++;
241 RArray[index1] *= denominator;
246 RArray = &RDerivatives;
248 index = (index << 1) + index;
249 RArray[0] = StoreDerivatives[index]; index++;
250 RArray[1] = StoreDerivatives[index]; index++;
251 RArray[2] = StoreDerivatives[index];
255 //=======================================================================
256 //function : PrepareEval
258 //=======================================================================
263 // Prepare all data for computing points :
264 // local arrays of knots
265 // local array of poles (multiplied by the weights if rational)
267 // The first direction to compute (smaller degree) is returned
268 // and the poles are stored according to this direction.
270 static Standard_Boolean PrepareEval (const Standard_Real U,
271 const Standard_Real V,
272 const Standard_Integer Uindex,
273 const Standard_Integer Vindex,
274 const Standard_Integer UDegree,
275 const Standard_Integer VDegree,
276 const Standard_Boolean URat,
277 const Standard_Boolean VRat,
278 const Standard_Boolean UPer,
279 const Standard_Boolean VPer,
280 const TColgp_Array2OfPnt& Poles,
281 const TColStd_Array2OfReal& Weights,
282 const TColStd_Array1OfReal& UKnots,
283 const TColStd_Array1OfReal& VKnots,
284 const TColStd_Array1OfInteger& UMults,
285 const TColStd_Array1OfInteger& VMults,
286 Standard_Real& u1, // first parameter to use
287 Standard_Real& u2, // second parameter to use
288 Standard_Integer& d1, // first degree
289 Standard_Integer& d2, // second degree
290 Standard_Boolean& rational,
291 BSplSLib_DataContainer& dc)
293 rational = URat || VRat;
294 Standard_Integer uindex = Uindex;
295 Standard_Integer vindex = Vindex;
296 Standard_Integer UKLower = UKnots.Lower();
297 Standard_Integer UKUpper = UKnots.Upper();
298 Standard_Integer VKLower = VKnots.Lower();
299 Standard_Integer VKUpper = VKnots.Upper();
301 if (UDegree <= VDegree)
303 // compute the indices
304 if (uindex < UKLower || uindex > UKUpper)
305 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
309 if (vindex < VKLower || vindex > VKUpper)
310 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
317 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
318 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
321 uindex -= UKLower + UDegree;
323 uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
326 vindex -= VKLower + VDegree;
328 vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
331 Standard_Integer i,j,ip,jp;
332 Standard_Real w, *pole = dc.poles;
335 Standard_Integer PLowerRow = Poles.LowerRow();
336 Standard_Integer PUpperRow = Poles.UpperRow();
337 Standard_Integer PLowerCol = Poles.LowerCol();
338 Standard_Integer PUpperCol = Poles.UpperCol();
340 // verify if locally non rational
343 rational = Standard_False;
344 ip = PLowerRow + uindex;
345 jp = PLowerCol + vindex;
347 if(ip < PLowerRow) ip = PUpperRow;
348 if(jp < PLowerCol) jp = PUpperCol;
350 w = Weights.Value(ip,jp);
351 Standard_Real eps = Epsilon(w);
354 for (i = 0; i <= UDegree && !rational; i++)
356 jp = PLowerCol + vindex;
361 for (j = 0; j <= VDegree && !rational; j++)
363 dw = Weights.Value(ip,jp) - w;
367 rational = (dw > eps);
384 ip = PLowerRow + uindex;
391 for (i = 0; i <= d1; i++)
393 jp = PLowerCol + vindex;
398 for (j = 0; j <= d2; j++)
400 const gp_Pnt& P = Poles .Value(ip,jp);
401 pole[3] = w = Weights.Value(ip,jp);
421 for (i = 0; i <= d1; i++)
423 jp = PLowerCol + vindex;
428 for (j = 0; j <= d2; j++)
430 const gp_Pnt& P = Poles.Value(ip,jp);
448 return Standard_True;
452 // compute the indices
453 if (uindex < UKLower || uindex > UKUpper)
454 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
458 if (vindex < VKLower || vindex > VKUpper)
459 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
468 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
469 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
472 uindex -= UKLower + UDegree;
474 uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
477 vindex -= VKLower + VDegree;
479 vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
482 Standard_Integer i,j,ip,jp;
483 Standard_Real w, *pole = dc.poles;
486 Standard_Integer PLowerRow = Poles.LowerRow();
487 Standard_Integer PUpperRow = Poles.UpperRow();
488 Standard_Integer PLowerCol = Poles.LowerCol();
489 Standard_Integer PUpperCol = Poles.UpperCol();
491 // verify if locally non rational
494 rational = Standard_False;
495 ip = PLowerRow + uindex;
496 jp = PLowerCol + vindex;
504 w = Weights.Value(ip,jp);
505 Standard_Real eps = Epsilon(w);
508 for (i = 0; i <= UDegree && !rational; i++)
510 jp = PLowerCol + vindex;
515 for (j = 0; j <= VDegree && !rational; j++)
517 dw = Weights.Value(ip,jp) - w;
518 if (dw < 0) dw = - dw;
536 jp = PLowerCol + vindex;
543 for (i = 0; i <= d1; i++)
545 ip = PLowerRow + uindex;
550 for (j = 0; j <= d2; j++)
552 const gp_Pnt& P = Poles.Value(ip,jp);
553 pole[3] = w = Weights.Value(ip,jp);
574 for (i = 0; i <= d1; i++)
576 ip = PLowerRow + uindex;
581 for (j = 0; j <= d2; j++)
583 const gp_Pnt& P = Poles.Value(ip,jp);
602 return Standard_False;
606 //=======================================================================
609 //=======================================================================
612 (const Standard_Real U,
613 const Standard_Real V,
614 const Standard_Integer UIndex,
615 const Standard_Integer VIndex,
616 const TColgp_Array2OfPnt& Poles,
617 const TColStd_Array2OfReal& Weights,
618 const TColStd_Array1OfReal& UKnots,
619 const TColStd_Array1OfReal& VKnots,
620 const TColStd_Array1OfInteger& UMults,
621 const TColStd_Array1OfInteger& VMults,
622 const Standard_Integer UDegree,
623 const Standard_Integer VDegree,
624 const Standard_Boolean URat,
625 const Standard_Boolean VRat,
626 const Standard_Boolean UPer,
627 const Standard_Boolean VPer,
630 // Standard_Integer k ;
655 //=======================================================================
658 //=======================================================================
660 void BSplSLib::HomogeneousD0
661 (const Standard_Real U,
662 const Standard_Real V,
663 const Standard_Integer UIndex,
664 const Standard_Integer VIndex,
665 const TColgp_Array2OfPnt& Poles,
666 const TColStd_Array2OfReal& Weights,
667 const TColStd_Array1OfReal& UKnots,
668 const TColStd_Array1OfReal& VKnots,
669 const TColStd_Array1OfInteger& UMults,
670 const TColStd_Array1OfInteger& VMults,
671 const Standard_Integer UDegree,
672 const Standard_Integer VDegree,
673 const Standard_Boolean URat,
674 const Standard_Boolean VRat,
675 const Standard_Boolean UPer,
676 const Standard_Boolean VPer,
680 Standard_Boolean rational;
681 // Standard_Integer k,dim;
682 Standard_Integer dim;
684 Standard_Integer d1,d2;
687 BSplSLib_DataContainer dc (UDegree, VDegree);
688 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
689 Poles,Weights,UKnots,VKnots,UMults,VMults,
690 u1,u2,d1,d2,rational,dc);
693 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
694 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
702 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
703 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
710 //=======================================================================
713 //=======================================================================
716 (const Standard_Real U,
717 const Standard_Real V,
718 const Standard_Integer UIndex,
719 const Standard_Integer VIndex,
720 const TColgp_Array2OfPnt& Poles,
721 const TColStd_Array2OfReal& Weights,
722 const TColStd_Array1OfReal& UKnots,
723 const TColStd_Array1OfReal& VKnots,
724 const TColStd_Array1OfInteger& UMults,
725 const TColStd_Array1OfInteger& VMults,
726 const Standard_Integer UDegree,
727 const Standard_Integer VDegree,
728 const Standard_Boolean URat,
729 const Standard_Boolean VRat,
730 const Standard_Boolean UPer,
731 const Standard_Boolean VPer,
736 Standard_Boolean rational;
737 // Standard_Integer k,dim,dim2;
738 Standard_Integer dim,dim2;
740 Standard_Integer d1,d2;
741 Standard_Real *result, *resVu, *resVv;
742 BSplSLib_DataContainer dc (UDegree, VDegree);
744 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
745 Poles,Weights,UKnots,VKnots,UMults,VMults,
746 u1,u2,d1,d2,rational,dc)) {
749 dim2 = (d2 + 1) << 2;
750 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
751 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
752 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
753 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
761 dim2 = (dim2 << 1) + dim2;
762 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
763 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
764 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
766 resVu = result + dim2;
773 dim2 = (d2 + 1) << 2;
774 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
775 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
776 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
777 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
785 dim2 = (dim2 << 1) + dim2;
786 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
787 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
788 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
791 resVv = result + dim2;
808 //=======================================================================
811 //=======================================================================
813 void BSplSLib::HomogeneousD1
814 (const Standard_Real U,
815 const Standard_Real V,
816 const Standard_Integer UIndex,
817 const Standard_Integer VIndex,
818 const TColgp_Array2OfPnt& Poles,
819 const TColStd_Array2OfReal& Weights,
820 const TColStd_Array1OfReal& UKnots,
821 const TColStd_Array1OfReal& VKnots,
822 const TColStd_Array1OfInteger& UMults,
823 const TColStd_Array1OfInteger& VMults,
824 const Standard_Integer UDegree,
825 const Standard_Integer VDegree,
826 const Standard_Boolean URat,
827 const Standard_Boolean VRat,
828 const Standard_Boolean UPer,
829 const Standard_Boolean VPer,
837 Standard_Boolean rational;
838 // Standard_Integer k,dim;
839 Standard_Integer dim;
841 Standard_Integer d1,d2;
846 BSplSLib_DataContainer dc (UDegree, VDegree);
847 Standard_Boolean ufirst = PrepareEval
848 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
849 Poles,Weights,UKnots,VKnots,UMults,VMults,
850 u1,u2,d1,d2,rational,dc);
851 dim = rational ? 4 : 3;
853 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
854 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
855 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
857 Standard_Real *result, *resVu, *resVv;
859 resVu = result + (ufirst ? dim*(d2+1) : dim);
860 resVv = result + (ufirst ? dim : dim*(d2+1));
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,
910 Standard_Boolean rational;
911 // Standard_Integer k,dim,dim2;
912 Standard_Integer dim,dim2;
914 Standard_Integer d1,d2;
915 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
916 BSplSLib_DataContainer dc (UDegree, VDegree);
918 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
919 Poles,Weights,UKnots,VKnots,UMults,VMults,
920 u1,u2,d1,d2,rational,dc)) {
923 dim2 = (d2 + 1) << 2;
924 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
925 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
926 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
928 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
929 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
933 resVuu = result + 18;
935 resVuv = result + 12;
940 dim2 = (dim2 << 1) + dim2;
941 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
942 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
943 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
945 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
947 resVu = result + dim2;
949 if (UDegree <= 1) resVuu = BSplSLib_zero;
950 else resVuu = result + (dim2 << 1);
951 if (VDegree <= 1) resVvv = BSplSLib_zero;
952 else resVvv = result + 6;
953 resVuv = result + (d2 << 1) + d2 + 6;
959 dim2 = (d2 + 1) << 2;
960 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
961 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
962 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
964 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
965 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
970 resVvv = result + 18;
971 resVuv = result + 12;
976 dim2 = (dim2 << 1) + dim2;
977 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
978 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
979 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
981 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
984 resVv = result + dim2;
985 if (UDegree <= 1) resVuu = BSplSLib_zero;
986 else resVuu = result + 6;
987 if (VDegree <= 1) resVvv = BSplSLib_zero;
988 else resVvv = result + (dim2 << 1);
989 resVuv = result + (d2 << 1) + d2 + 6;
1001 Vu .SetY(resVu [1]);
1002 Vv .SetY(resVv [1]);
1003 Vuu.SetY(resVuu[1]);
1004 Vvv.SetY(resVvv[1]);
1005 Vuv.SetY(resVuv[1]);
1008 Vu .SetZ(resVu [2]);
1009 Vv .SetZ(resVv [2]);
1010 Vuu.SetZ(resVuu[2]);
1011 Vvv.SetZ(resVvv[2]);
1012 Vuv.SetZ(resVuv[2]);
1015 //=======================================================================
1018 //=======================================================================
1021 (const Standard_Real U,
1022 const Standard_Real V,
1023 const Standard_Integer UIndex,
1024 const Standard_Integer VIndex,
1025 const TColgp_Array2OfPnt& Poles,
1026 const TColStd_Array2OfReal& Weights,
1027 const TColStd_Array1OfReal& UKnots,
1028 const TColStd_Array1OfReal& VKnots,
1029 const TColStd_Array1OfInteger& UMults,
1030 const TColStd_Array1OfInteger& VMults,
1031 const Standard_Integer UDegree,
1032 const Standard_Integer VDegree,
1033 const Standard_Boolean URat,
1034 const Standard_Boolean VRat,
1035 const Standard_Boolean UPer,
1036 const Standard_Boolean VPer,
1048 Standard_Boolean rational;
1049 // Standard_Integer k,dim,dim2;
1050 Standard_Integer dim,dim2;
1051 Standard_Real u1,u2;
1052 Standard_Integer d1,d2;
1053 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
1054 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
1055 BSplSLib_DataContainer dc (UDegree, VDegree);
1057 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1058 Poles,Weights,UKnots,VKnots,UMults,VMults,
1059 u1,u2,d1,d2,rational,dc)) {
1062 dim2 = (d2 + 1) << 2;
1063 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1064 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1065 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1067 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1069 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1070 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1072 resVu = result + 12;
1074 resVuu = result + 24;
1075 resVvv = result + 6;
1076 resVuv = result + 15;
1077 resVuuu = result + 36;
1078 resVvvv = result + 9;
1079 resVuuv = result + 27;
1080 resVuvv = result + 18;
1085 dim2 = (dim2 << 1) + dim2;
1086 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1087 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1088 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1090 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1092 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1094 resVu = result + dim2;
1097 resVuu = BSplSLib_zero;
1098 resVuuv = BSplSLib_zero;
1101 resVuu = result + (dim2 << 1);
1102 resVuuv = result + (dim2 << 1) + 3;
1105 resVvv = BSplSLib_zero;
1106 resVuvv = BSplSLib_zero;
1109 resVvv = result + 6;
1110 resVuvv = result + dim2 + 6;
1112 resVuv = result + (d2 << 1) + d2 + 6;
1113 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1114 else resVuuu = result + (dim2 << 1) + dim2;
1115 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1116 else resVvvv = result + 9;
1122 dim2 = (d2 + 1) << 2;
1123 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1124 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1125 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1127 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1129 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1130 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1133 resVv = result + 12;
1134 resVuu = result + 6;
1135 resVvv = result + 24;
1136 resVuv = result + 15;
1137 resVuuu = result + 9;
1138 resVvvv = result + 36;
1139 resVuuv = result + 18;
1140 resVuvv = result + 27;
1145 dim2 = (dim2 << 1) + dim2;
1146 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1147 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1148 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1150 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1152 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1155 resVv = result + dim2;
1157 resVuu = BSplSLib_zero;
1158 resVuuv = BSplSLib_zero;
1161 resVuu = result + 6;
1162 resVuuv = result + dim2 + 6;
1165 resVvv = BSplSLib_zero;
1166 resVuvv = BSplSLib_zero;
1169 resVvv = result + (dim2 << 1);
1170 resVuvv = result + (dim2 << 1) + 3;
1172 resVuv = result + (d2 << 1) + d2 + 6;
1173 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1174 else resVuuu = result + 9;
1175 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1176 else resVvvv = result + (dim2 << 1) + dim2;
1180 P .SetX(result [0]);
1181 Vu .SetX(resVu [0]);
1182 Vv .SetX(resVv [0]);
1183 Vuu .SetX(resVuu [0]);
1184 Vvv .SetX(resVvv [0]);
1185 Vuv .SetX(resVuv [0]);
1186 Vuuu.SetX(resVuuu[0]);
1187 Vvvv.SetX(resVvvv[0]);
1188 Vuuv.SetX(resVuuv[0]);
1189 Vuvv.SetX(resVuvv[0]);
1191 P .SetY(result [1]);
1192 Vu .SetY(resVu [1]);
1193 Vv .SetY(resVv [1]);
1194 Vuu .SetY(resVuu [1]);
1195 Vvv .SetY(resVvv [1]);
1196 Vuv .SetY(resVuv [1]);
1197 Vuuu.SetY(resVuuu[1]);
1198 Vvvv.SetY(resVvvv[1]);
1199 Vuuv.SetY(resVuuv[1]);
1200 Vuvv.SetY(resVuvv[1]);
1202 P .SetZ(result [2]);
1203 Vu .SetZ(resVu [2]);
1204 Vv .SetZ(resVv [2]);
1205 Vuu .SetZ(resVuu [2]);
1206 Vvv .SetZ(resVvv [2]);
1207 Vuv .SetZ(resVuv [2]);
1208 Vuuu.SetZ(resVuuu[2]);
1209 Vvvv.SetZ(resVvvv[2]);
1210 Vuuv.SetZ(resVuuv[2]);
1211 Vuvv.SetZ(resVuvv[2]);
1214 //=======================================================================
1217 //=======================================================================
1220 (const Standard_Real U,
1221 const Standard_Real V,
1222 const Standard_Integer Nu,
1223 const Standard_Integer Nv,
1224 const Standard_Integer UIndex,
1225 const Standard_Integer VIndex,
1226 const TColgp_Array2OfPnt& Poles,
1227 const TColStd_Array2OfReal& Weights,
1228 const TColStd_Array1OfReal& UKnots,
1229 const TColStd_Array1OfReal& VKnots,
1230 const TColStd_Array1OfInteger& UMults,
1231 const TColStd_Array1OfInteger& VMults,
1232 const Standard_Integer UDegree,
1233 const Standard_Integer VDegree,
1234 const Standard_Boolean URat,
1235 const Standard_Boolean VRat,
1236 const Standard_Boolean UPer,
1237 const Standard_Boolean VPer,
1240 Standard_Boolean rational;
1241 Standard_Integer k,dim;
1242 Standard_Real u1,u2;
1243 Standard_Integer d1,d2;
1245 BSplSLib_DataContainer dc (UDegree, VDegree);
1246 Standard_Boolean ufirst = PrepareEval
1247 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1248 Poles,Weights,UKnots,VKnots,UMults,VMults,
1249 u1,u2,d1,d2,rational,dc);
1250 dim = rational ? 4 : 3;
1253 if ((Nu > UDegree) || (Nv > VDegree)) {
1261 Standard_Integer n1 = ufirst ? Nu : Nv;
1262 Standard_Integer n2 = ufirst ? Nv : Nu;
1264 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1266 for (k = 0; k <= Min(n1,d1); k++)
1267 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1269 Standard_Real *result;
1271 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1272 result = dc.ders; // because Standard_False ci-dessus.
1276 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1285 // Surface modifications
1287 // a surface is processed as a curve of curves.
1288 // i.e : a curve of parameter u where the current point is the set of poles
1292 //=======================================================================
1295 //=======================================================================
1297 void BSplSLib::Iso(const Standard_Real Param,
1298 const Standard_Boolean IsU,
1299 const TColgp_Array2OfPnt& Poles,
1300 const TColStd_Array2OfReal& Weights,
1301 const TColStd_Array1OfReal& Knots,
1302 const TColStd_Array1OfInteger& Mults,
1303 const Standard_Integer Degree,
1304 const Standard_Boolean Periodic,
1305 TColgp_Array1OfPnt& CPoles,
1306 TColStd_Array1OfReal& CWeights)
1308 Standard_Integer index = 0;
1309 Standard_Real u = Param;
1310 Standard_Boolean rational = &Weights != NULL;
1311 Standard_Integer dim = rational ? 4 : 3;
1313 // compute local knots
1315 NCollection_LocalArray<Standard_Real> locknots1 (2*Degree);
1316 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1317 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1319 index -= Knots.Lower() + Degree;
1321 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1324 // copy the local poles
1326 // Standard_Integer f1,l1,f2,l2,i,j,k;
1327 Standard_Integer f1,l1,f2,l2,i,j;
1330 f1 = Poles.LowerRow();
1331 l1 = Poles.UpperRow();
1332 f2 = Poles.LowerCol();
1333 l2 = Poles.UpperCol();
1336 f1 = Poles.LowerCol();
1337 l1 = Poles.UpperCol();
1338 f2 = Poles.LowerRow();
1339 l2 = Poles.UpperRow();
1342 NCollection_LocalArray<Standard_Real> locpoles ((Degree+1) * (l2-f2+1) * dim);
1344 Standard_Real w, *pole = locpoles;
1347 for (i = 0; i <= Degree; i++) {
1349 for (j = f2; j <= l2; j++) {
1351 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1353 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1354 pole[0] = P.X() * w;
1355 pole[1] = P.Y() * w;
1356 pole[2] = P.Z() * w;
1366 if (index > l1) index = f1;
1370 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1375 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1376 gp_Pnt& P = CPoles(i);
1378 CWeights(i) = w = pole[3];
1379 P.SetX( pole[0] / w);
1380 P.SetY( pole[1] / w);
1381 P.SetZ( pole[2] / w);
1391 // if the input is not rational but weights are wanted
1392 if (!rational && (&CWeights != NULL)) {
1394 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1399 //=======================================================================
1400 //function : Reverse
1402 //=======================================================================
1404 void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1405 const Standard_Integer Last,
1406 const Standard_Boolean UDirection)
1408 Standard_Integer i,j, l = Last;
1410 l = Poles.LowerRow() +
1411 (l - Poles.LowerRow())%(Poles.ColLength());
1412 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1413 Poles.LowerCol(), Poles.UpperCol());
1415 for (i = Poles.LowerRow(); i <= l; i++) {
1417 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1418 temp(l-i,j) = Poles(i,j);
1422 for (i = l+1; i <= Poles.UpperRow(); i++) {
1424 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1425 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1429 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1431 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1432 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1437 l = Poles.LowerCol() +
1438 (l - Poles.LowerCol())%(Poles.RowLength());
1439 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1440 0, Poles.RowLength()-1);
1442 for (j = Poles.LowerCol(); j <= l; j++) {
1444 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1445 temp(i,l-j) = Poles(i,j);
1449 for (j = l+1; j <= Poles.UpperCol(); j++) {
1451 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1452 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1456 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1458 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1459 Poles(i,j) = temp (i,j-Poles.LowerCol());
1465 //=======================================================================
1466 //function : Reverse
1468 //=======================================================================
1470 void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1471 const Standard_Integer Last,
1472 const Standard_Boolean UDirection)
1474 Standard_Integer i,j, l = Last;
1476 l = Weights.LowerRow() +
1477 (l - Weights.LowerRow())%(Weights.ColLength());
1478 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1479 Weights.LowerCol(), Weights.UpperCol());
1481 for (i = Weights.LowerRow(); i <= l; i++) {
1483 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1484 temp(l-i,j) = Weights(i,j);
1488 for (i = l+1; i <= Weights.UpperRow(); i++) {
1490 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1491 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1495 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1497 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1498 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1503 l = Weights.LowerCol() +
1504 (l - Weights.LowerCol())%(Weights.RowLength());
1505 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1506 0, Weights.RowLength()-1);
1508 for (j = Weights.LowerCol(); j <= l; j++) {
1510 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1511 temp(i,l-j) = Weights(i,j);
1515 for (j = l+1; j <= Weights.UpperCol(); j++) {
1517 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1518 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1522 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1524 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1525 Weights(i,j) = temp (i,j-Weights.LowerCol());
1531 //=======================================================================
1532 //function : IsRational
1534 //=======================================================================
1536 Standard_Boolean BSplSLib::IsRational
1537 (const TColStd_Array2OfReal& Weights,
1538 const Standard_Integer I1,
1539 const Standard_Integer I2,
1540 const Standard_Integer J1,
1541 const Standard_Integer J2,
1542 const Standard_Real Epsi)
1544 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1545 Standard_Integer i,j;
1546 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1547 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1549 for (i = I1 - fi; i < I2 - fi; i++) {
1551 for (j = J1 - fj; j < J2 - fj; j++) {
1552 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1553 return Standard_True;
1556 return Standard_False;
1559 //=======================================================================
1560 //function : SetPoles
1562 //=======================================================================
1564 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1565 TColStd_Array1OfReal& FP,
1566 const Standard_Boolean UDirection)
1568 Standard_Integer i, j, l = FP.Lower();
1569 Standard_Integer PLowerRow = Poles.LowerRow();
1570 Standard_Integer PUpperRow = Poles.UpperRow();
1571 Standard_Integer PLowerCol = Poles.LowerCol();
1572 Standard_Integer PUpperCol = Poles.UpperCol();
1575 for ( i = PLowerRow; i <= PUpperRow; i++) {
1577 for ( j = PLowerCol; j <= PUpperCol; j++) {
1578 const gp_Pnt& P = Poles.Value(i,j);
1587 for ( j = PLowerCol; j <= PUpperCol; j++) {
1589 for ( i = PLowerRow; i <= PUpperRow; i++) {
1590 const gp_Pnt& P = Poles.Value(i,j);
1599 //=======================================================================
1600 //function : SetPoles
1602 //=======================================================================
1604 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1605 const TColStd_Array2OfReal& Weights,
1606 TColStd_Array1OfReal& FP,
1607 const Standard_Boolean UDirection)
1609 Standard_Integer i, j, l = FP.Lower();
1610 Standard_Integer PLowerRow = Poles.LowerRow();
1611 Standard_Integer PUpperRow = Poles.UpperRow();
1612 Standard_Integer PLowerCol = Poles.LowerCol();
1613 Standard_Integer PUpperCol = Poles.UpperCol();
1616 for ( i = PLowerRow; i <= PUpperRow; i++) {
1618 for ( j = PLowerCol; j <= PUpperCol; j++) {
1619 const gp_Pnt& P = Poles .Value(i,j);
1620 Standard_Real w = Weights.Value(i,j);
1621 FP(l) = P.X() * w; l++;
1622 FP(l) = P.Y() * w; l++;
1623 FP(l) = P.Z() * w; l++;
1630 for ( j = PLowerCol; j <= PUpperCol; j++) {
1632 for ( i = PLowerRow; i <= PUpperRow; i++) {
1633 const gp_Pnt& P = Poles .Value(i,j);
1634 Standard_Real w = Weights.Value(i,j);
1635 FP(l) = P.X() * w; l++;
1636 FP(l) = P.Y() * w; l++;
1637 FP(l) = P.Z() * w; l++;
1644 //=======================================================================
1645 //function : GetPoles
1647 //=======================================================================
1649 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1650 TColgp_Array2OfPnt& Poles,
1651 const Standard_Boolean UDirection)
1653 Standard_Integer i, j, l = FP.Lower();
1654 Standard_Integer PLowerRow = Poles.LowerRow();
1655 Standard_Integer PUpperRow = Poles.UpperRow();
1656 Standard_Integer PLowerCol = Poles.LowerCol();
1657 Standard_Integer PUpperCol = Poles.UpperCol();
1660 for ( i = PLowerRow; i <= PUpperRow; i++) {
1662 for ( j = PLowerCol; j <= PUpperCol; j++) {
1663 gp_Pnt& P = Poles.ChangeValue(i,j);
1672 for ( j = PLowerCol; j <= PUpperCol; j++) {
1674 for ( i = PLowerRow; i <= PUpperRow; i++) {
1675 gp_Pnt& P = Poles.ChangeValue(i,j);
1684 //=======================================================================
1685 //function : GetPoles
1687 //=======================================================================
1689 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1690 TColgp_Array2OfPnt& Poles,
1691 TColStd_Array2OfReal& Weights,
1692 const Standard_Boolean UDirection)
1694 Standard_Integer i, j, l = FP.Lower();
1695 Standard_Integer PLowerRow = Poles.LowerRow();
1696 Standard_Integer PUpperRow = Poles.UpperRow();
1697 Standard_Integer PLowerCol = Poles.LowerCol();
1698 Standard_Integer PUpperCol = Poles.UpperCol();
1701 for ( i = PLowerRow; i <= PUpperRow; i++) {
1703 for ( j = PLowerCol; j <= PUpperCol; j++) {
1704 Standard_Real w = FP( l + 3);
1706 gp_Pnt& P = Poles.ChangeValue(i,j);
1707 P.SetX(FP(l) / w); l++;
1708 P.SetY(FP(l) / w); l++;
1709 P.SetZ(FP(l) / w); l++;
1716 for ( j = PLowerCol; j <= PUpperCol; j++) {
1718 for ( i = PLowerRow; i <= PUpperRow; i++) {
1719 Standard_Real w = FP( l + 3);
1721 gp_Pnt& P = Poles.ChangeValue(i,j);
1722 P.SetX(FP(l) / w); l++;
1723 P.SetY(FP(l) / w); l++;
1724 P.SetZ(FP(l) / w); l++;
1731 //=======================================================================
1732 //function : InsertKnots
1734 //=======================================================================
1736 void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1737 const Standard_Integer Degree,
1738 const Standard_Boolean Periodic,
1739 const TColgp_Array2OfPnt& Poles,
1740 const TColStd_Array2OfReal& Weights,
1741 const TColStd_Array1OfReal& Knots,
1742 const TColStd_Array1OfInteger& Mults,
1743 const TColStd_Array1OfReal& AddKnots,
1744 const TColStd_Array1OfInteger& AddMults,
1745 TColgp_Array2OfPnt& NewPoles,
1746 TColStd_Array2OfReal& NewWeights,
1747 TColStd_Array1OfReal& NewKnots,
1748 TColStd_Array1OfInteger& NewMults,
1749 const Standard_Real Epsilon,
1750 const Standard_Boolean Add )
1752 Standard_Boolean rational = &Weights != NULL;
1753 Standard_Integer dim = 3;
1754 if (rational) dim++;
1756 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1757 TColStd_Array1OfReal
1758 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1760 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1761 else SetPoles(Poles,poles,UDirection);
1764 dim *= Poles.RowLength();
1767 dim *= Poles.ColLength();
1769 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1770 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1773 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1774 else GetPoles(newpoles,NewPoles,UDirection);
1777 //=======================================================================
1778 //function : RemoveKnot
1780 //=======================================================================
1782 Standard_Boolean BSplSLib::RemoveKnot
1783 (const Standard_Boolean UDirection,
1784 const Standard_Integer Index,
1785 const Standard_Integer Mult,
1786 const Standard_Integer Degree,
1787 const Standard_Boolean Periodic,
1788 const TColgp_Array2OfPnt& Poles,
1789 const TColStd_Array2OfReal& Weights,
1790 const TColStd_Array1OfReal& Knots,
1791 const TColStd_Array1OfInteger& Mults,
1792 TColgp_Array2OfPnt& NewPoles,
1793 TColStd_Array2OfReal& NewWeights,
1794 TColStd_Array1OfReal& NewKnots,
1795 TColStd_Array1OfInteger& NewMults,
1796 const Standard_Real Tolerance)
1798 Standard_Boolean rational = &Weights != NULL;
1799 Standard_Integer dim = 3;
1800 if (rational) dim++;
1802 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1803 TColStd_Array1OfReal
1804 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1806 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1807 else SetPoles(Poles,poles,UDirection);
1810 dim *= Poles.RowLength();
1813 dim *= Poles.ColLength();
1816 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1817 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1819 return Standard_False;
1821 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1822 else GetPoles(newpoles,NewPoles,UDirection);
1823 return Standard_True;
1826 //=======================================================================
1827 //function : IncreaseDegree
1829 //=======================================================================
1831 void BSplSLib::IncreaseDegree
1832 (const Standard_Boolean UDirection,
1833 const Standard_Integer Degree,
1834 const Standard_Integer NewDegree,
1835 const Standard_Boolean Periodic,
1836 const TColgp_Array2OfPnt& Poles,
1837 const TColStd_Array2OfReal& Weights,
1838 const TColStd_Array1OfReal& Knots,
1839 const TColStd_Array1OfInteger& Mults,
1840 TColgp_Array2OfPnt& NewPoles,
1841 TColStd_Array2OfReal& NewWeights,
1842 TColStd_Array1OfReal& NewKnots,
1843 TColStd_Array1OfInteger& NewMults)
1845 Standard_Boolean rational = &Weights != NULL;
1846 Standard_Integer dim = 3;
1847 if (rational) dim++;
1849 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1850 TColStd_Array1OfReal
1851 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1853 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1854 else SetPoles(Poles,poles,UDirection);
1857 dim *= Poles.RowLength();
1860 dim *= Poles.ColLength();
1863 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1864 newpoles,NewKnots,NewMults);
1866 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1867 else GetPoles(newpoles,NewPoles,UDirection);
1870 //=======================================================================
1871 //function : Unperiodize
1873 //=======================================================================
1875 void BSplSLib::Unperiodize
1876 (const Standard_Boolean UDirection,
1877 const Standard_Integer Degree,
1878 const TColStd_Array1OfInteger& Mults,
1879 const TColStd_Array1OfReal& Knots,
1880 const TColgp_Array2OfPnt& Poles,
1881 const TColStd_Array2OfReal& Weights,
1882 TColStd_Array1OfInteger& NewMults,
1883 TColStd_Array1OfReal& NewKnots,
1884 TColgp_Array2OfPnt& NewPoles,
1885 TColStd_Array2OfReal& NewWeights)
1887 Standard_Boolean rational = &Weights != NULL;
1888 Standard_Integer dim = 3;
1889 if (rational) dim++;
1891 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1892 TColStd_Array1OfReal
1893 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1895 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1896 else SetPoles(Poles,poles,UDirection);
1899 dim *= Poles.RowLength();
1902 dim *= Poles.ColLength();
1905 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1906 NewMults, NewKnots, newpoles);
1908 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1909 else GetPoles(newpoles,NewPoles,UDirection);
1912 //=======================================================================
1913 //function : BuildCache
1914 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1915 // Cache : in case of a rational function the Poles are
1916 // stored in homogeneous form
1917 //=======================================================================
1919 void BSplSLib::BuildCache
1920 (const Standard_Real U,
1921 const Standard_Real V,
1922 const Standard_Real USpanDomain,
1923 const Standard_Real VSpanDomain,
1924 const Standard_Boolean UPeriodic,
1925 const Standard_Boolean VPeriodic,
1926 const Standard_Integer UDegree,
1927 const Standard_Integer VDegree,
1928 const Standard_Integer UIndex,
1929 const Standard_Integer VIndex,
1930 const TColStd_Array1OfReal& UFlatKnots,
1931 const TColStd_Array1OfReal& VFlatKnots,
1932 const TColgp_Array2OfPnt& Poles,
1933 const TColStd_Array2OfReal& Weights,
1934 TColgp_Array2OfPnt& CachePoles,
1935 TColStd_Array2OfReal& CacheWeights)
1937 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1938 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1939 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1940 if (&Weights != NULL)
1941 rational_u = rational_v = Standard_True;
1943 rational_u = rational_v = Standard_False;
1944 BSplSLib_DataContainer dc (UDegree, VDegree);
1960 (BSplCLib::NoMults()),
1961 (BSplCLib::NoMults()),
1971 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1973 for (kk = 0; kk <= d1 ; kk++)
1974 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1976 min_degree_domain = USpanDomain ;
1977 max_degree_domain = VSpanDomain ;
1980 min_degree_domain = VSpanDomain ;
1981 max_degree_domain = USpanDomain ;
1985 for (ii = 0 ; ii <= d2 ; ii++) {
1989 for (jj = 0 ; jj <= d1 ; jj++) {
1991 Index = jj * d2p1 + ii ;
1993 gp_Pnt& P = CachePoles(iii,jjj);
1994 f = factor[0] * factor[1];
1995 P.SetX( f * dc.poles[Index]); Index++;
1996 P.SetY( f * dc.poles[Index]); Index++;
1997 P.SetZ( f * dc.poles[Index]); Index++;
1998 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1999 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
2001 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
2005 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
2007 for (kk = 0; kk <= d1 ; kk++)
2008 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
2010 min_degree_domain = USpanDomain ;
2011 max_degree_domain = VSpanDomain ;
2014 min_degree_domain = VSpanDomain ;
2015 max_degree_domain = USpanDomain ;
2019 for (ii = 0 ; ii <= d2 ; ii++) {
2023 for (jj = 0 ; jj <= d1 ; jj++) {
2025 Index = jj * d2p1 + ii ;
2026 Index = (Index << 1) + Index;
2027 gp_Pnt& P = CachePoles(iii,jjj);
2028 f = factor[0] * factor[1];
2029 P.SetX( f * dc.poles[Index]); Index++;
2030 P.SetY( f * dc.poles[Index]); Index++;
2031 P.SetZ( f * dc.poles[Index]);
2032 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
2034 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
2036 if (&Weights != NULL) {
2038 // means that PrepareEval did found out that the surface was
2039 // locally polynomial but since the surface is constructed
2040 // with some weights we need to set the weight polynomial to constant
2043 for (ii = 1 ; ii <= d2p1 ; ii++) {
2045 for (jj = 1 ; jj <= d1p1 ; jj++) {
2046 CacheWeights(ii,jj) = 0.0e0 ;
2049 CacheWeights(1,1) = 1.0e0 ;
2054 //=======================================================================
2055 //function : CacheD0
2056 //purpose : Evaluates the polynomial cache of the Bspline Curve
2058 //=======================================================================
2060 void BSplSLib::CacheD0(const Standard_Real UParameter,
2061 const Standard_Real VParameter,
2062 const Standard_Integer UDegree,
2063 const Standard_Integer VDegree,
2064 const Standard_Real UCacheParameter,
2065 const Standard_Real VCacheParameter,
2066 const Standard_Real USpanLenght,
2067 const Standard_Real VSpanLenght,
2068 const TColgp_Array2OfPnt& PolesArray,
2069 const TColStd_Array2OfReal& WeightsArray,
2073 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2075 // the SpanLenght is the normalizing factor so that the polynomial is between
2088 PArray = (Standard_Real *)
2089 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2091 myPoint = (Standard_Real *) &aPoint ;
2092 if (UDegree <= VDegree) {
2093 min_degree = UDegree ;
2094 max_degree = VDegree ;
2095 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
2096 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
2097 dimension = 3 * (UDegree + 1) ;
2100 min_degree = VDegree ;
2101 max_degree = UDegree ;
2102 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
2103 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
2104 dimension = 3 * (VDegree + 1) ;
2106 NCollection_LocalArray<Standard_Real> locpoles(dimension);
2108 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2111 max_degree*dimension,
2115 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2118 (min_degree << 1) + min_degree,
2121 if (&WeightsArray != NULL) {
2122 dimension = min_degree + 1 ;
2124 WArray = (Standard_Real *)
2125 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2126 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2129 max_degree*dimension,
2133 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2139 inverse = 1.0e0/ inverse ;
2141 myPoint[0] *= inverse ;
2142 myPoint[1] *= inverse ;
2143 myPoint[2] *= inverse ;
2147 //=======================================================================
2148 //function : CacheD1
2149 //purpose : Evaluates the polynomial cache of the Bspline Curve
2151 //=======================================================================
2153 void BSplSLib::CacheD1(const Standard_Real UParameter,
2154 const Standard_Real VParameter,
2155 const Standard_Integer UDegree,
2156 const Standard_Integer VDegree,
2157 const Standard_Real UCacheParameter,
2158 const Standard_Real VCacheParameter,
2159 const Standard_Real USpanLenght,
2160 const Standard_Real VSpanLenght,
2161 const TColgp_Array2OfPnt& PolesArray,
2162 const TColStd_Array2OfReal& WeightsArray,
2168 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2170 // the SpanLenght is the normalizing factor so that the polynomial is between
2186 PArray = (Standard_Real *)
2187 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2188 Standard_Real local_poles_array[2][2][3],
2189 local_poles_and_weights_array[2][2][4],
2190 local_weights_array[2][2] ;
2191 Standard_Real * my_vec_min,
2194 my_point = (Standard_Real *) &aPoint ;
2196 // initialize in case of rational evaluation
2197 // because RationalDerivative will use all
2201 if (&WeightsArray != NULL) {
2203 local_poles_array [0][0][0] = 0.0e0 ;
2204 local_poles_array [0][0][1] = 0.0e0 ;
2205 local_poles_array [0][0][2] = 0.0e0 ;
2206 local_weights_array [0][0] = 0.0e0 ;
2207 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2208 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2209 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2210 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2212 local_poles_array [0][1][0] = 0.0e0 ;
2213 local_poles_array [0][1][1] = 0.0e0 ;
2214 local_poles_array [0][1][2] = 0.0e0 ;
2215 local_weights_array [0][1] = 0.0e0 ;
2216 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2217 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2218 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2219 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2221 local_poles_array [1][0][0] = 0.0e0 ;
2222 local_poles_array [1][0][1] = 0.0e0 ;
2223 local_poles_array [1][0][2] = 0.0e0 ;
2224 local_weights_array [1][0] = 0.0e0 ;
2225 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2226 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2227 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2228 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2230 local_poles_array [1][1][0] = 0.0e0 ;
2231 local_poles_array [1][1][1] = 0.0e0 ;
2232 local_poles_array [1][1][2] = 0.0e0 ;
2233 local_weights_array [1][1] = 0.0e0 ;
2234 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2235 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2236 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2237 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2240 if (UDegree <= VDegree) {
2241 min_degree = UDegree ;
2242 max_degree = VDegree ;
2243 inverse_min = 1.0e0/USpanLenght ;
2244 inverse_max = 1.0e0/VSpanLenght ;
2245 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2246 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2248 dimension = 3 * (UDegree + 1) ;
2249 my_vec_min = (Standard_Real *) &aVecU ;
2250 my_vec_max = (Standard_Real *) &aVecV ;
2253 min_degree = VDegree ;
2254 max_degree = UDegree ;
2255 inverse_min = 1.0e0/VSpanLenght ;
2256 inverse_max = 1.0e0/USpanLenght ;
2257 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2258 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2259 dimension = 3 * (VDegree + 1) ;
2260 my_vec_min = (Standard_Real *) &aVecV ;
2261 my_vec_max = (Standard_Real *) &aVecU ;
2264 NCollection_LocalArray<Standard_Real> locpoles (2 * dimension);
2266 PLib::EvalPolynomial(new_parameter[0],
2273 PLib::EvalPolynomial(new_parameter[1],
2278 local_poles_array[0][0][0]) ;
2279 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2282 (min_degree << 1) + min_degree,
2283 locpoles[dimension],
2284 local_poles_array[1][0][0]) ;
2286 if (&WeightsArray != NULL) {
2287 dimension = min_degree + 1 ;
2289 WArray = (Standard_Real *)
2290 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2291 PLib::EvalPolynomial(new_parameter[0],
2298 PLib::EvalPolynomial(new_parameter[1],
2303 local_weights_array[0][0]) ;
2304 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2308 locpoles[dimension],
2309 local_weights_array[1][0]) ;
2311 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2312 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2313 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2314 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2316 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2317 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2318 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2319 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2321 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2322 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2323 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2324 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2326 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2327 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2328 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2329 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2331 BSplSLib::RationalDerivative(1,
2335 local_poles_and_weights_array[0][0][0],
2336 local_poles_array[0][0][0]) ;
2339 my_point [0] = local_poles_array [0][0][0] ;
2340 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2341 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2343 my_point [1] = local_poles_array [0][0][1] ;
2344 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2345 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2347 my_point [2] = local_poles_array [0][0][2] ;
2348 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2349 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2352 //=======================================================================
2353 //function : CacheD2
2354 //purpose : Evaluates the polynomial cache of the Bspline Curve
2356 //=======================================================================
2358 void BSplSLib::CacheD2(const Standard_Real UParameter,
2359 const Standard_Real VParameter,
2360 const Standard_Integer UDegree,
2361 const Standard_Integer VDegree,
2362 const Standard_Real UCacheParameter,
2363 const Standard_Real VCacheParameter,
2364 const Standard_Real USpanLenght,
2365 const Standard_Real VSpanLenght,
2366 const TColgp_Array2OfPnt& PolesArray,
2367 const TColStd_Array2OfReal& WeightsArray,
2376 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2378 // the SpanLenght is the normalizing factor so that the polynomial is between
2395 PArray = (Standard_Real *)
2396 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2397 Standard_Real local_poles_array[3][3][3],
2398 local_poles_and_weights_array[3][3][4],
2399 local_weights_array[3][3] ;
2400 Standard_Real * my_vec_min,
2406 my_point = (Standard_Real *) &aPoint ;
2409 // initialize in case the min and max degree are less than 2
2411 local_poles_array[0][0][0] = 0.0e0 ;
2412 local_poles_array[0][0][1] = 0.0e0 ;
2413 local_poles_array[0][0][2] = 0.0e0 ;
2414 local_poles_array[0][1][0] = 0.0e0 ;
2415 local_poles_array[0][1][1] = 0.0e0 ;
2416 local_poles_array[0][1][2] = 0.0e0 ;
2417 local_poles_array[0][2][0] = 0.0e0 ;
2418 local_poles_array[0][2][1] = 0.0e0 ;
2419 local_poles_array[0][2][2] = 0.0e0 ;
2421 local_poles_array[1][0][0] = 0.0e0 ;
2422 local_poles_array[1][0][1] = 0.0e0 ;
2423 local_poles_array[1][0][2] = 0.0e0 ;
2424 local_poles_array[1][1][0] = 0.0e0 ;
2425 local_poles_array[1][1][1] = 0.0e0 ;
2426 local_poles_array[1][1][2] = 0.0e0 ;
2427 local_poles_array[1][2][0] = 0.0e0 ;
2428 local_poles_array[1][2][1] = 0.0e0 ;
2429 local_poles_array[1][2][2] = 0.0e0 ;
2431 local_poles_array[2][0][0] = 0.0e0 ;
2432 local_poles_array[2][0][1] = 0.0e0 ;
2433 local_poles_array[2][0][2] = 0.0e0 ;
2434 local_poles_array[2][1][0] = 0.0e0 ;
2435 local_poles_array[2][1][1] = 0.0e0 ;
2436 local_poles_array[2][1][2] = 0.0e0 ;
2437 local_poles_array[2][2][0] = 0.0e0 ;
2438 local_poles_array[2][2][1] = 0.0e0 ;
2439 local_poles_array[2][2][2] = 0.0e0 ;
2441 // initialize in case of rational evaluation
2442 // because RationalDerivative will use all
2446 if (&WeightsArray != NULL) {
2448 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2449 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2450 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2451 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2452 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2453 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2454 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2455 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2456 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2458 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2459 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2460 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2461 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2462 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2463 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2464 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2465 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2466 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2468 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2469 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2470 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2471 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2472 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2473 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2474 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2475 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2476 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2478 local_poles_and_weights_array[0][0][3] =
2479 local_weights_array[0][0] = 0.0e0 ;
2480 local_poles_and_weights_array[0][1][3] =
2481 local_weights_array[0][1] = 0.0e0 ;
2482 local_poles_and_weights_array[0][2][3] =
2483 local_weights_array[0][2] = 0.0e0 ;
2484 local_poles_and_weights_array[1][0][3] =
2485 local_weights_array[1][0] = 0.0e0 ;
2486 local_poles_and_weights_array[1][1][3] =
2487 local_weights_array[1][1] = 0.0e0 ;
2488 local_poles_and_weights_array[1][2][3] =
2489 local_weights_array[1][2] = 0.0e0 ;
2490 local_poles_and_weights_array[2][0][3] =
2491 local_weights_array[2][0] = 0.0e0 ;
2492 local_poles_and_weights_array[2][1][3] =
2493 local_weights_array[2][1] = 0.0e0 ;
2494 local_poles_and_weights_array[2][2][3] =
2495 local_weights_array[2][2] = 0.0e0 ;
2498 if (UDegree <= VDegree) {
2499 min_degree = UDegree ;
2500 max_degree = VDegree ;
2501 inverse_min = 1.0e0/USpanLenght ;
2502 inverse_max = 1.0e0/VSpanLenght ;
2503 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2504 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2506 dimension = 3 * (UDegree + 1) ;
2507 my_vec_min = (Standard_Real *) &aVecU ;
2508 my_vec_max = (Standard_Real *) &aVecV ;
2509 my_vec_min_min = (Standard_Real *) &aVecUU ;
2510 my_vec_min_max = (Standard_Real *) &aVecUV ;
2511 my_vec_max_max = (Standard_Real *) &aVecVV ;
2514 min_degree = VDegree ;
2515 max_degree = UDegree ;
2516 inverse_min = 1.0e0/VSpanLenght ;
2517 inverse_max = 1.0e0/USpanLenght ;
2518 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2519 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2520 dimension = 3 * (VDegree + 1) ;
2521 my_vec_min = (Standard_Real *) &aVecV ;
2522 my_vec_max = (Standard_Real *) &aVecU ;
2523 my_vec_min_min = (Standard_Real *) &aVecVV ;
2524 my_vec_min_max = (Standard_Real *) &aVecUV ;
2525 my_vec_max_max = (Standard_Real *) &aVecUU ;
2528 NCollection_LocalArray<Standard_Real> locpoles (3 * dimension);
2531 // initialize in case min or max degree are less than 2
2533 Standard_Integer MinIndMax = 2;
2534 if ( max_degree < 2) MinIndMax = max_degree;
2535 Standard_Integer MinIndMin = 2;
2536 if ( min_degree < 2) MinIndMin = min_degree;
2538 index = MinIndMax * dimension ;
2540 for (ii = MinIndMax ; ii < 3 ; ii++) {
2542 for (kk = 0 ; kk < dimension ; kk++) {
2543 locpoles[index] = 0.0e0 ;
2548 PLib::EvalPolynomial(new_parameter[0],
2555 PLib::EvalPolynomial(new_parameter[1],
2560 local_poles_array[0][0][0]) ;
2561 PLib::EvalPolynomial(new_parameter[1],
2565 locpoles[dimension],
2566 local_poles_array[1][0][0]) ;
2568 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2571 (min_degree << 1) + min_degree,
2572 locpoles[dimension + dimension],
2573 local_poles_array[2][0][0]) ;
2575 if (&WeightsArray != NULL) {
2576 dimension = min_degree + 1 ;
2578 WArray = (Standard_Real *)
2579 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2580 PLib::EvalPolynomial(new_parameter[0],
2587 PLib::EvalPolynomial(new_parameter[1],
2592 local_weights_array[0][0]) ;
2593 PLib::EvalPolynomial(new_parameter[1],
2597 locpoles[dimension],
2598 local_weights_array[1][0]) ;
2599 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2603 locpoles[dimension + dimension],
2604 local_weights_array[2][0]) ;
2607 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2608 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2609 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2610 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2611 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2612 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2613 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2614 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2615 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2617 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2618 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2619 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2620 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2621 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2622 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2623 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2624 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2625 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2627 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2628 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2629 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2630 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2631 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2632 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2633 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2634 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2635 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2638 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2639 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2640 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2641 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2642 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2643 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2644 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2645 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2646 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2648 BSplSLib::RationalDerivative(2,
2652 local_poles_and_weights_array[0][0][0],
2653 local_poles_array[0][0][0]) ;
2657 Standard_Real minmin = inverse_min * inverse_min;
2658 Standard_Real minmax = inverse_min * inverse_max;
2659 Standard_Real maxmax = inverse_max * inverse_max;
2661 my_point [0] = local_poles_array [0][0][0] ;
2662 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2663 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2664 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2665 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2666 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2668 my_point [1] = local_poles_array [0][0][1] ;
2669 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2670 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2671 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2672 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2673 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2675 my_point [2] = local_poles_array [0][0][2] ;
2676 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2677 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2678 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2679 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2680 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2683 //=======================================================================
2684 //function : MovePoint
2685 //purpose : Find the new poles which allows an old point (with a
2686 // given u and v as parameters) to reach a new position
2687 //=======================================================================
2689 void BSplSLib::MovePoint (const Standard_Real U,
2690 const Standard_Real V,
2691 const gp_Vec& Displ,
2692 const Standard_Integer UIndex1,
2693 const Standard_Integer UIndex2,
2694 const Standard_Integer VIndex1,
2695 const Standard_Integer VIndex2,
2696 const Standard_Integer UDegree,
2697 const Standard_Integer VDegree,
2698 const Standard_Boolean Rational,
2699 const TColgp_Array2OfPnt& Poles,
2700 const TColStd_Array2OfReal& Weights,
2701 const TColStd_Array1OfReal& UFlatKnots,
2702 const TColStd_Array1OfReal& VFlatKnots,
2703 Standard_Integer& UFirstIndex,
2704 Standard_Integer& ULastIndex,
2705 Standard_Integer& VFirstIndex,
2706 Standard_Integer& VLastIndex,
2707 TColgp_Array2OfPnt& NewPoles)
2709 // calculate the UBSplineBasis in the parameter U
2710 Standard_Integer UFirstNonZeroBsplineIndex;
2711 math_Matrix UBSplineBasis(1, 1,
2713 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2718 UFirstNonZeroBsplineIndex,
2720 // calculate the VBSplineBasis in the parameter V
2721 Standard_Integer VFirstNonZeroBsplineIndex;
2722 math_Matrix VBSplineBasis(1, 1,
2724 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2729 VFirstNonZeroBsplineIndex,
2731 if (ErrorCod1 || ErrorCod2) {
2739 // find the span which is predominant for parameter U
2740 UFirstIndex = UFirstNonZeroBsplineIndex;
2741 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2742 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2743 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2745 Standard_Real maxValue = 0.0;
2746 Standard_Integer i, ukk1=0, ukk2;
2748 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2749 if (UBSplineBasis(1,i) > maxValue) {
2750 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2751 maxValue = UBSplineBasis(1,i);
2755 // find a ukk2 if symetriy
2757 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2758 if ((ukk1+1) <= ULastIndex) {
2759 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2764 // find the span which is predominant for parameter V
2765 VFirstIndex = VFirstNonZeroBsplineIndex;
2766 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2768 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2769 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2772 Standard_Integer j, vkk1=0, vkk2;
2774 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2775 if (VBSplineBasis(1,j) > maxValue) {
2776 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2777 maxValue = VBSplineBasis(1,j);
2781 // find a vkk2 if symetriy
2783 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2784 if ((vkk1+1) <= VLastIndex) {
2785 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2790 // compute the vector of displacement
2791 Standard_Real D1 = 0.0;
2792 Standard_Real D2 = 0.0;
2793 Standard_Real hN, Coef, DvalU, DvalV;
2795 Standard_Integer ii, jj;
2797 for (i = 1; i <= UDegree+1; i++) {
2798 ii = i + UFirstNonZeroBsplineIndex - 1;
2802 else if (ii > ukk2) {
2809 for (j = 1; j <= VDegree+1; j++) {
2810 jj = j + VFirstNonZeroBsplineIndex - 1;
2812 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2816 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2818 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2822 else if (jj > vkk2) {
2828 D1 += 1./(DvalU + DvalV + 1.) * hN;
2840 // compute the new poles
2842 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2846 else if (i > ukk2) {
2853 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2854 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2858 else if (j > vkk2) {
2864 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2867 NewPoles(i,j) = Poles(i,j);
2873 //=======================================================================
2874 // function : Resolution
2875 // purpose : this computes an estimate for the maximum of the
2876 // partial derivatives both in U and in V
2879 // The calculation resembles at the calculation of curves with
2880 // additional index for the control point. Let Si,j be the
2881 // control points for ls surface and Di,j the weights.
2882 // The checking of upper bounds for the partial derivatives
2883 // will be omitted and Su is the next upper bound in the polynomial case :
2887 // | Si,j - Si-1,j |
2888 // d * Max | ------------- |
2889 // i = 2,n | ti+d - ti |
2893 // and in the rational case :
2897 // Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2898 // Max Max d * -----------------------------------------------
2899 // k=1,n i dans Rj ti+d - ti
2901 // ----------------------------------------------------------------------
2909 // with Rj = {j-d, ...., j+d+d+1}.
2912 //=======================================================================
2914 void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
2915 const TColStd_Array2OfReal& Weights,
2916 const TColStd_Array1OfReal& UKnots,
2917 const TColStd_Array1OfReal& VKnots,
2918 const TColStd_Array1OfInteger& UMults,
2919 const TColStd_Array1OfInteger& VMults,
2920 const Standard_Integer UDegree,
2921 const Standard_Integer VDegree,
2922 const Standard_Boolean URational,
2923 const Standard_Boolean VRational,
2924 const Standard_Boolean UPeriodic,
2925 const Standard_Boolean VPeriodic,
2926 const Standard_Real Tolerance3D,
2927 Standard_Real& UTolerance,
2928 Standard_Real& VTolerance)
2930 Standard_Real Wij,Wmj,Wji,Wjm;
2931 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2932 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2933 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2934 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2936 max_derivative[0] = max_derivative[1] = 0.0e0 ;
2938 Standard_Integer PRowLength, PColLength;
2939 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2940 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2941 Standard_Integer num_poles[2],num_flat_knots[2];
2944 BSplCLib::KnotSequenceLength(UMults,
2948 BSplCLib::KnotSequenceLength(VMults,
2951 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
2952 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
2953 BSplCLib::KnotSequence(UKnots,
2958 BSplCLib::KnotSequence(VKnots,
2963 PRowLength = Poles.RowLength();
2964 PColLength = Poles.ColLength();
2965 if (URational || VRational) {
2966 Standard_Integer Wsize = PRowLength * PColLength;
2967 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2968 min_weights = WG[0];
2970 for (ii = 1 ; ii < Wsize ; ii++) {
2972 if (min_weights > min) min_weights = min;
2975 Standard_Integer UD1 = UDegree + 1;
2976 Standard_Integer VD1 = VDegree + 1;
2977 num_poles[0] = num_flat_knots[0] - UD1;
2978 num_poles[1] = num_flat_knots[1] - VD1;
2979 poles_length[0] = PColLength;
2980 poles_length[1] = PRowLength;
2982 Standard_Integer UD2 = UDegree << 1;
2983 Standard_Integer VD2 = VDegree << 1;
2985 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2986 ii_index = (ii - 1) % poles_length[0] + 1 ;
2987 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2988 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2989 inverse = 1.0e0 / inverse ;
2990 lower[0] = ii - UD1;
2991 if (lower[0] < 1) lower[0] = 1;
2992 upper[0] = ii + UD2 + 1;
2993 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2995 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2996 jj_index = (jj - 1) % poles_length[1] + 1 ;
2997 lower[1] = jj - VD1;
2998 if (lower[1] < 1) lower[1] = 1;
2999 upper[1] = jj + VD2 + 1;
3000 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
3001 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
3002 Wij = Weights.Value(ii_index,jj_index);
3003 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
3004 Wmj = Weights.Value(ii_minus,jj_index);
3012 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
3013 pp_index = (pp - 1) % poles_length[0] + 1 ;
3015 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
3017 qq_index = (qq - 1) % poles_length[1] + 1 ;
3018 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
3022 factor = (Xpq - Xij) * Wij;
3023 factor -= (Xpq - Xmj) * Wmj;
3024 if (factor < 0) factor = - factor;
3026 factor = (Ypq - Yij) * Wij;
3027 factor -= (Ypq - Ymj) * Wmj;
3028 if (factor < 0) factor = - factor;
3030 factor = (Zpq - Zij) * Wij;
3031 factor -= (Zpq - Zmj) * Wmj;
3032 if (factor < 0) factor = - factor;
3035 if (max_derivative[0] < value) max_derivative[0] = value ;
3040 max_derivative[0] /= min_weights ;
3044 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3045 ii_index = (ii - 1) % poles_length[0] + 1 ;
3046 ii_minus = (ii - 2) % poles_length[0] + 1 ;
3047 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3048 inverse = 1.0e0 / inverse ;
3050 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3051 jj_index = (jj - 1) % poles_length[1] + 1 ;
3053 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
3054 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
3055 factor = Pij.X() - Pmj.X();
3056 if (factor < 0) factor = - factor;
3058 factor = Pij.Y() - Pmj.Y();
3059 if (factor < 0) factor = - factor;
3061 factor = Pij.Z() - Pmj.Z();
3062 if (factor < 0) factor = - factor;
3065 if (max_derivative[0] < value) max_derivative[0] = value ;
3069 max_derivative[0] *= UDegree ;
3071 Standard_Integer UD2 = UDegree << 1;
3072 Standard_Integer VD2 = VDegree << 1;
3074 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3075 ii_index = (ii - 1) % poles_length[1] + 1 ;
3076 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3077 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3078 inverse = 1.0e0 / inverse ;
3079 lower[0] = ii - VD1;
3080 if (lower[0] < 1) lower[0] = 1;
3081 upper[0] = ii + VD2 + 1;
3082 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
3084 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3085 jj_index = (jj - 1) % poles_length[0] + 1 ;
3086 lower[1] = jj - UD1;
3087 if (lower[1] < 1) lower[1] = 1;
3088 upper[1] = jj + UD2 + 1;
3089 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
3090 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
3091 Wji = Weights.Value(jj_index,ii_index);
3092 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
3093 Wjm = Weights.Value(jj_index,ii_minus);
3101 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
3102 pp_index = (pp - 1) % poles_length[1] + 1 ;
3104 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
3106 qq_index = (qq - 1) % poles_length[0] + 1 ;
3107 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
3111 factor = (Xqp - Xji) * Wji;
3112 factor -= (Xqp - Xjm) * Wjm;
3113 if (factor < 0) factor = - factor;
3115 factor = (Yqp - Yji) * Wji;
3116 factor -= (Yqp - Yjm) * Wjm;
3117 if (factor < 0) factor = - factor;
3119 factor = (Zqp - Zji) * Wji;
3120 factor -= (Zqp - Zjm) * Wjm;
3121 if (factor < 0) factor = - factor;
3124 if (max_derivative[1] < value) max_derivative[1] = value ;
3129 max_derivative[1] /= min_weights ;
3133 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3134 ii_index = (ii - 1) % poles_length[1] + 1 ;
3135 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3136 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3137 inverse = 1.0e0 / inverse ;
3139 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3140 jj_index = (jj - 1) % poles_length[0] + 1 ;
3142 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3143 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3144 factor = Pji.X() - Pjm.X() ;
3145 if (factor < 0) factor = - factor;
3147 factor = Pji.Y() - Pjm.Y() ;
3148 if (factor < 0) factor = - factor;
3150 factor = Pji.Z() - Pjm.Z() ;
3151 if (factor < 0) factor = - factor;
3154 if (max_derivative[1] < value) max_derivative[1] = value ;
3158 max_derivative[1] *= VDegree ;
3159 max_derivative[0] *= M_SQRT2 ;
3160 max_derivative[1] *= M_SQRT2 ;
3161 if(max_derivative[0] && max_derivative[1]) {
3162 UTolerance = Tolerance3D / max_derivative[0] ;
3163 VTolerance = Tolerance3D / max_derivative[1] ;
3166 UTolerance=VTolerance=0.0;
3168 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3173 //=======================================================================
3174 //function : Interpolate
3176 //=======================================================================
3178 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3179 const Standard_Integer VDegree,
3180 const TColStd_Array1OfReal& UFlatKnots,
3181 const TColStd_Array1OfReal& VFlatKnots,
3182 const TColStd_Array1OfReal& UParameters,
3183 const TColStd_Array1OfReal& VParameters,
3184 TColgp_Array2OfPnt& Poles,
3185 TColStd_Array2OfReal& Weights,
3186 Standard_Integer& InversionProblem)
3188 Standard_Integer ii, jj, ll, kk, dimension;
3189 Standard_Integer ULength = UParameters.Length();
3190 Standard_Integer VLength = VParameters.Length();
3191 Standard_Real * poles_array;
3193 // extraction of iso u
3194 dimension = 4*ULength;
3195 TColStd_Array2OfReal Points(1, VLength,
3198 Handle(TColStd_HArray1OfInteger) ContactOrder =
3199 new (TColStd_HArray1OfInteger)(1, VLength);
3200 ContactOrder->Init(0);
3202 for (ii=1; ii <= VLength; ii++) {
3204 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3205 Points(ii,ll) = Poles(jj, ii).X();
3206 Points(ii,ll+1) = Poles(jj, ii).Y();
3207 Points(ii,ll+2) = Poles(jj, ii).Z();
3208 Points(ii,ll+3) = Weights(jj,ii) ;
3212 // interpolation of iso u
3213 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3214 BSplCLib::Interpolate(VDegree,
3217 ContactOrder->Array1(),
3221 if (InversionProblem != 0) return;
3223 // extraction of iso v
3225 dimension = VLength*4;
3226 TColStd_Array2OfReal IsoPoles(1, ULength,
3229 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3230 ContactOrder->Init(0);
3231 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3233 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3235 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3236 IsoPoles (ii,ll) = Points(jj, kk);
3237 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3238 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3239 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3242 // interpolation of iso v
3243 BSplCLib::Interpolate(UDegree,
3246 ContactOrder->Array1(),
3253 for (ii=1; ii <= ULength; ii++) {
3255 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3256 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3257 Poles.SetValue(ii, jj, Pnt);
3258 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3263 //=======================================================================
3264 //function : Interpolate
3266 //=======================================================================
3268 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3269 const Standard_Integer VDegree,
3270 const TColStd_Array1OfReal& UFlatKnots,
3271 const TColStd_Array1OfReal& VFlatKnots,
3272 const TColStd_Array1OfReal& UParameters,
3273 const TColStd_Array1OfReal& VParameters,
3274 TColgp_Array2OfPnt& Poles,
3275 Standard_Integer& InversionProblem)
3277 Standard_Integer ii, jj, ll, kk, dimension;
3278 Standard_Integer ULength = UParameters.Length();
3279 Standard_Integer VLength = VParameters.Length();
3280 Standard_Real * poles_array;
3282 // extraction of iso u
3283 dimension = 3*ULength;
3284 TColStd_Array2OfReal Points(1, VLength,
3287 Handle(TColStd_HArray1OfInteger) ContactOrder =
3288 new (TColStd_HArray1OfInteger)(1, VLength);
3289 ContactOrder->Init(0);
3291 for (ii=1; ii <= VLength; ii++) {
3293 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3294 Points(ii,ll) = Poles(jj, ii).X();
3295 Points(ii,ll+1) = Poles(jj, ii).Y();
3296 Points(ii,ll+2) = Poles(jj, ii).Z();
3300 // interpolation of iso u
3301 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3302 BSplCLib::Interpolate(VDegree,
3305 ContactOrder->Array1(),
3309 if (InversionProblem != 0) return;
3311 // extraction of iso v
3313 dimension = VLength*3;
3314 TColStd_Array2OfReal IsoPoles(1, ULength,
3317 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3318 ContactOrder->Init(0);
3319 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3321 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3323 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3324 IsoPoles (ii,ll) = Points(jj, kk);
3325 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3326 IsoPoles (ii,ll+2) = Points(jj, kk+2);
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+=3) {
3343 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3344 Poles.SetValue(ii, jj, Pnt);
3349 //=======================================================================
3350 //function : FunctionMultiply
3352 //=======================================================================
3354 void BSplSLib::FunctionMultiply
3355 (const BSplSLib_EvaluatorFunction& Function,
3356 const Standard_Integer UBSplineDegree,
3357 const Standard_Integer VBSplineDegree,
3358 const TColStd_Array1OfReal& UBSplineKnots,
3359 const TColStd_Array1OfReal& VBSplineKnots,
3360 const TColStd_Array1OfInteger & UMults,
3361 const TColStd_Array1OfInteger & VMults,
3362 const TColgp_Array2OfPnt& Poles,
3363 const TColStd_Array2OfReal& Weights,
3364 const TColStd_Array1OfReal& UFlatKnots,
3365 const TColStd_Array1OfReal& VFlatKnots,
3366 const Standard_Integer UNewDegree,
3367 const Standard_Integer VNewDegree,
3368 TColgp_Array2OfPnt& NewNumerator,
3369 TColStd_Array2OfReal& NewDenominator,
3370 Standard_Integer& Status)
3372 Standard_Integer num_uparameters,
3377 Standard_Real result ;
3379 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3380 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3381 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3382 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3384 if ((NewNumerator.ColLength() == num_uparameters) &&
3385 (NewNumerator.RowLength() == num_vparameters) &&
3386 (NewDenominator.ColLength() == num_uparameters) &&
3387 (NewDenominator.RowLength() == num_vparameters)) {
3390 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3394 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3398 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3400 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3401 HomogeneousD0(UParameters(ii),
3417 NewDenominator(ii,jj),
3418 NewNumerator(ii,jj)) ;
3420 Function.Evaluate (0,
3426 Standard_ConstructionError::Raise();
3428 gp_Pnt& P = NewNumerator(ii,jj);
3429 P.SetX(P.X() * result);
3430 P.SetY(P.Y() * result);
3431 P.SetZ(P.Z() * result);
3432 NewDenominator(ii,jj) *= result ;
3435 Interpolate(UNewDegree,
3446 Standard_ConstructionError::Raise();