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 <PLib_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 Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
57 VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
58 "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 typedef PLib_LocalArray BSplSLib_LocalArray;
69 //**************************************************************************
71 //**************************************************************************
73 //=======================================================================
74 //function : RationalDerivative
75 //purpose : computes the rational derivatives when whe have the
76 // the derivatives of the homogeneous numerator and the
77 // the derivatives of the denominator
78 //=======================================================================
80 void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
81 const Standard_Integer VDeg,
82 const Standard_Integer N,
83 const Standard_Integer M,
84 Standard_Real& HDerivatives,
85 Standard_Real& RDerivatives,
86 const Standard_Boolean All)
89 // if All is True all derivatives are computed. if Not only
90 // the requested N, M is computed
93 // let f(u,v) be a rational function = ------------------
97 // Let (N,M) the order of the derivatives we want : then since
100 // Numerator = f * Denominator
104 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
105 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
106 // (0,0) ( p<N q<M p q )
116 // HDerivatives is an array where derivatives are stored in the following form
117 // Numerator is assumee to have 3 functions that is a vector of dimension
120 // (0,0) (0,0) (0, DegV) (0, DegV)
121 // Numerator Denominator ... Numerator Denominator
123 // (1,0) (1,0) (1, DegV) (1, DegV)
124 // Numerator Denominator ... Numerator Denominator
126 // ...........................................................
129 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
130 // Numerator Denominator ... Numerator Denominator
133 Standard_Integer ii,jj,pp,qq,index,index1,index2;
134 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
135 Standard_Integer MinN,MinN1,MinM,MinM1;
136 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
142 M4 = (VDeg + 1) << 2;
144 BSplSLib_LocalArray StoreDerivatives (All ? 0 : ii * 3);
145 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
146 BSplSLib_LocalArray StoreW (ii);
147 Standard_Real *HomogeneousArray = &HDerivatives;
148 Standard_Real denominator,Pii,Pip,Pjq;
150 denominator = 1.0e0 / HomogeneousArray[3];
153 if (UDeg < N) MinN = UDeg;
155 if (VDeg < M) MinM = VDeg;
161 for (ii = 0 ; ii < MinN1 ; ii++) {
167 for (jj = 0 ; jj < MinM1 ; jj++) {
168 RArray[index_v++] = HomogeneousArray[index_v1++];
169 RArray[index_v++] = HomogeneousArray[index_v1++];
170 RArray[index_v++] = HomogeneousArray[index_v1++];
171 StoreW[index_w++] = HomogeneousArray[index_v1++];
174 for (jj = MinM1 ; jj < M1 ; jj++) {
175 RArray[index_v++] = 0.;
176 RArray[index_v++] = 0.;
177 RArray[index_v++] = 0.;
178 StoreW[index_w++] = 0.;
183 index_v = MinN1 * M3;
184 index_w = MinN1 * M1;
186 for (ii = MinN1 ; ii < N1 ; ii++) {
188 for (jj = 0 ; jj < M1 ; jj++) {
189 RArray[index_v++] = 0.0e0;
190 RArray[index_v++] = 0.0e0;
191 RArray[index_v++] = 0.0e0;
192 StoreW[index_w++] = 0.0e0;
196 // --------------- Calculation ----------------
201 for (ii = 0 ; ii <= N ; ii++) {
207 for (jj = 0 ; jj <= M ; jj++) {
213 for (pp = 0 ; pp < ii ; pp++) {
217 index2 = jjM1 - ppM1;
218 Pip = PLib::Bin(ii,pp);
220 for (qq = 0 ; qq <= jj ; qq++) {
222 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
223 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
224 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
225 RArray[index1] -= Pjq * RArray[index]; index++;
231 Pii = PLib::Bin(ii,ii);
233 for (qq = 0 ; qq < jj ; qq++) {
235 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
236 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
237 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
238 RArray[index1] -= Pjq * RArray[index]; index++;
241 RArray[index1] *= denominator; index1++;
242 RArray[index1] *= denominator; index1++;
243 RArray[index1] *= denominator;
248 RArray = &RDerivatives;
250 index = (index << 1) + index;
251 RArray[0] = StoreDerivatives[index]; index++;
252 RArray[1] = StoreDerivatives[index]; index++;
253 RArray[2] = StoreDerivatives[index];
257 //=======================================================================
258 //function : PrepareEval
260 //=======================================================================
265 // Prepare all data for computing points :
266 // local arrays of knots
267 // local array of poles (multiplied by the weights if rational)
269 // The first direction to compute (smaller degree) is returned
270 // and the poles are stored according to this direction.
272 static Standard_Boolean PrepareEval
273 (const Standard_Real U,
274 const Standard_Real V,
275 const Standard_Integer Uindex,
276 const Standard_Integer Vindex,
277 const Standard_Integer UDegree,
278 const Standard_Integer VDegree,
279 const Standard_Boolean URat,
280 const Standard_Boolean VRat,
281 const Standard_Boolean UPer,
282 const Standard_Boolean VPer,
283 const TColgp_Array2OfPnt& Poles,
284 const TColStd_Array2OfReal& Weights,
285 const TColStd_Array1OfReal& UKnots,
286 const TColStd_Array1OfReal& VKnots,
287 const TColStd_Array1OfInteger& UMults,
288 const TColStd_Array1OfInteger& VMults,
289 Standard_Real& u1, // first parameter to use
290 Standard_Real& u2, // second parameter to use
291 Standard_Integer& d1, // first degree
292 Standard_Integer& d2, // second degree
293 Standard_Boolean& rational,
294 BSplSLib_DataContainer& dc)
296 rational = URat || VRat;
297 Standard_Integer uindex = Uindex;
298 Standard_Integer vindex = Vindex;
299 Standard_Integer UKLower = UKnots.Lower();
300 Standard_Integer UKUpper = UKnots.Upper();
301 Standard_Integer VKLower = VKnots.Lower();
302 Standard_Integer VKUpper = VKnots.Upper();
303 if (UDegree <= VDegree) {
304 // compute the indices
305 if (uindex < UKLower || uindex > UKUpper)
306 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
308 if (vindex < VKLower || vindex > VKUpper)
309 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
314 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
315 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
316 if (&UMults == NULL) uindex -= UKLower + UDegree;
317 else uindex = BSplCLib::PoleIndex
318 (UDegree,uindex,UPer,UMults);
319 if (&VMults == NULL) vindex -= VKLower + VDegree;
320 else vindex = BSplCLib::PoleIndex
321 (VDegree,vindex,VPer,VMults);
323 // Standard_Integer i,j,k,ip,jp;
324 Standard_Integer i,j,ip,jp;
325 Standard_Real w, *pole = dc.poles;
328 Standard_Integer PLowerRow = Poles.LowerRow();
329 Standard_Integer PUpperRow = Poles.UpperRow();
330 Standard_Integer PLowerCol = Poles.LowerCol();
331 Standard_Integer PUpperCol = Poles.UpperCol();
332 if (rational) { // verify if locally non rational
333 rational = Standard_False;
334 ip = PLowerRow + uindex;
335 jp = PLowerCol + vindex;
336 w = Weights.Value(ip,jp);
337 Standard_Real eps = Epsilon(w);
340 for (i = 0; i <= UDegree && !rational; i++) {
341 jp = PLowerCol + vindex;
343 for (j = 0; j <= VDegree && !rational; j++) {
344 dw = Weights.Value(ip,jp) - w;
345 if (dw < 0) dw = - dw;
348 if (jp > PUpperCol) jp = PLowerCol;
351 if (ip > PUpperRow) ip = PLowerRow;
355 ip = PLowerRow + uindex;
358 for (i = 0; i <= d1; i++) {
359 jp = PLowerCol + vindex;
361 for (j = 0; j <= d2; j++) {
362 const gp_Pnt& P = Poles .Value(ip,jp);
363 pole[3] = w = Weights.Value(ip,jp);
369 if (jp > PUpperCol) jp = PLowerCol;
372 if (ip > PUpperRow) ip = PLowerRow;
377 for (i = 0; i <= d1; i++) {
378 jp = PLowerCol + vindex;
380 for (j = 0; j <= d2; j++) {
381 const gp_Pnt& P = Poles.Value(ip,jp);
387 if (jp > PUpperCol) jp = PLowerCol;
390 if (ip > PUpperRow) ip = PLowerRow;
393 return Standard_True;
396 // compute the indices
397 if (uindex < UKLower || uindex > UKUpper)
398 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
400 if (vindex < VKLower || vindex > VKUpper)
401 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
406 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
407 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
408 if (&UMults == NULL) uindex -= UKLower + UDegree;
409 else uindex = BSplCLib::PoleIndex
410 (UDegree,uindex,UPer,UMults);
411 if (&VMults == NULL) vindex -= VKLower + VDegree;
412 else vindex = BSplCLib::PoleIndex
413 (VDegree,vindex,VPer,VMults);
415 // Standard_Integer i,j,k,ip,jp;
416 Standard_Integer i,j,ip,jp;
417 Standard_Real w, *pole = dc.poles;
420 Standard_Integer PLowerRow = Poles.LowerRow();
421 Standard_Integer PUpperRow = Poles.UpperRow();
422 Standard_Integer PLowerCol = Poles.LowerCol();
423 Standard_Integer PUpperCol = Poles.UpperCol();
424 if (rational) { // verify if locally non rational
425 rational = Standard_False;
426 ip = PLowerRow + uindex;
427 jp = PLowerCol + vindex;
428 w = Weights.Value(ip,jp);
429 Standard_Real eps = Epsilon(w);
432 for (i = 0; i <= UDegree && !rational; i++) {
433 jp = PLowerCol + vindex;
435 for (j = 0; j <= VDegree && !rational; j++) {
436 dw = Weights.Value(ip,jp) - w;
437 if (dw < 0) dw = - dw;
440 if (jp > PUpperCol) jp = PLowerCol;
443 if (ip > PUpperRow) ip = PLowerRow;
447 jp = PLowerCol + vindex;
450 for (i = 0; i <= d1; i++) {
451 ip = PLowerRow + uindex;
453 for (j = 0; j <= d2; j++) {
454 const gp_Pnt& P = Poles .Value(ip,jp);
455 pole[3] = w = Weights.Value(ip,jp);
461 if (ip > PUpperRow) ip = PLowerRow;
464 if (jp > PUpperCol) jp = PLowerCol;
469 for (i = 0; i <= d1; i++) {
470 ip = PLowerRow + uindex;
472 for (j = 0; j <= d2; j++) {
473 const gp_Pnt& P = Poles.Value(ip,jp);
479 if (ip > PUpperRow) ip = PLowerRow;
482 if (jp > PUpperCol) jp = PLowerCol;
485 return Standard_False;
489 //=======================================================================
492 //=======================================================================
495 (const Standard_Real U,
496 const Standard_Real V,
497 const Standard_Integer UIndex,
498 const Standard_Integer VIndex,
499 const TColgp_Array2OfPnt& Poles,
500 const TColStd_Array2OfReal& Weights,
501 const TColStd_Array1OfReal& UKnots,
502 const TColStd_Array1OfReal& VKnots,
503 const TColStd_Array1OfInteger& UMults,
504 const TColStd_Array1OfInteger& VMults,
505 const Standard_Integer UDegree,
506 const Standard_Integer VDegree,
507 const Standard_Boolean URat,
508 const Standard_Boolean VRat,
509 const Standard_Boolean UPer,
510 const Standard_Boolean VPer,
513 // Standard_Integer k ;
538 //=======================================================================
541 //=======================================================================
543 void BSplSLib::HomogeneousD0
544 (const Standard_Real U,
545 const Standard_Real V,
546 const Standard_Integer UIndex,
547 const Standard_Integer VIndex,
548 const TColgp_Array2OfPnt& Poles,
549 const TColStd_Array2OfReal& Weights,
550 const TColStd_Array1OfReal& UKnots,
551 const TColStd_Array1OfReal& VKnots,
552 const TColStd_Array1OfInteger& UMults,
553 const TColStd_Array1OfInteger& VMults,
554 const Standard_Integer UDegree,
555 const Standard_Integer VDegree,
556 const Standard_Boolean URat,
557 const Standard_Boolean VRat,
558 const Standard_Boolean UPer,
559 const Standard_Boolean VPer,
563 Standard_Boolean rational;
564 // Standard_Integer k,dim;
565 Standard_Integer dim;
567 Standard_Integer d1,d2;
570 BSplSLib_DataContainer dc (UDegree, VDegree);
571 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
572 Poles,Weights,UKnots,VKnots,UMults,VMults,
573 u1,u2,d1,d2,rational,dc);
576 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
577 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
585 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
586 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
593 //=======================================================================
596 //=======================================================================
599 (const Standard_Real U,
600 const Standard_Real V,
601 const Standard_Integer UIndex,
602 const Standard_Integer VIndex,
603 const TColgp_Array2OfPnt& Poles,
604 const TColStd_Array2OfReal& Weights,
605 const TColStd_Array1OfReal& UKnots,
606 const TColStd_Array1OfReal& VKnots,
607 const TColStd_Array1OfInteger& UMults,
608 const TColStd_Array1OfInteger& VMults,
609 const Standard_Integer UDegree,
610 const Standard_Integer VDegree,
611 const Standard_Boolean URat,
612 const Standard_Boolean VRat,
613 const Standard_Boolean UPer,
614 const Standard_Boolean VPer,
619 Standard_Boolean rational;
620 // Standard_Integer k,dim,dim2;
621 Standard_Integer dim,dim2;
623 Standard_Integer d1,d2;
624 Standard_Real *result, *resVu, *resVv;
625 BSplSLib_DataContainer dc (UDegree, VDegree);
627 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
628 Poles,Weights,UKnots,VKnots,UMults,VMults,
629 u1,u2,d1,d2,rational,dc)) {
632 dim2 = (d2 + 1) << 2;
633 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
634 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
635 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
636 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
644 dim2 = (dim2 << 1) + dim2;
645 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
646 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
647 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
649 resVu = result + dim2;
656 dim2 = (d2 + 1) << 2;
657 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
658 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
659 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
660 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
668 dim2 = (dim2 << 1) + dim2;
669 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
670 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
671 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
674 resVv = result + dim2;
691 //=======================================================================
694 //=======================================================================
696 void BSplSLib::HomogeneousD1
697 (const Standard_Real U,
698 const Standard_Real V,
699 const Standard_Integer UIndex,
700 const Standard_Integer VIndex,
701 const TColgp_Array2OfPnt& Poles,
702 const TColStd_Array2OfReal& Weights,
703 const TColStd_Array1OfReal& UKnots,
704 const TColStd_Array1OfReal& VKnots,
705 const TColStd_Array1OfInteger& UMults,
706 const TColStd_Array1OfInteger& VMults,
707 const Standard_Integer UDegree,
708 const Standard_Integer VDegree,
709 const Standard_Boolean URat,
710 const Standard_Boolean VRat,
711 const Standard_Boolean UPer,
712 const Standard_Boolean VPer,
720 Standard_Boolean rational;
721 // Standard_Integer k,dim;
722 Standard_Integer dim;
724 Standard_Integer d1,d2;
729 BSplSLib_DataContainer dc (UDegree, VDegree);
730 Standard_Boolean ufirst = PrepareEval
731 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
732 Poles,Weights,UKnots,VKnots,UMults,VMults,
733 u1,u2,d1,d2,rational,dc);
734 dim = rational ? 4 : 3;
736 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
737 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
738 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
740 Standard_Real *result, *resVu, *resVv;
742 resVu = result + (ufirst ? dim*(d2+1) : dim);
743 resVv = result + (ufirst ? dim : dim*(d2+1));
764 //=======================================================================
767 //=======================================================================
770 (const Standard_Real U,
771 const Standard_Real V,
772 const Standard_Integer UIndex,
773 const Standard_Integer VIndex,
774 const TColgp_Array2OfPnt& Poles,
775 const TColStd_Array2OfReal& Weights,
776 const TColStd_Array1OfReal& UKnots,
777 const TColStd_Array1OfReal& VKnots,
778 const TColStd_Array1OfInteger& UMults,
779 const TColStd_Array1OfInteger& VMults,
780 const Standard_Integer UDegree,
781 const Standard_Integer VDegree,
782 const Standard_Boolean URat,
783 const Standard_Boolean VRat,
784 const Standard_Boolean UPer,
785 const Standard_Boolean VPer,
793 Standard_Boolean rational;
794 // Standard_Integer k,dim,dim2;
795 Standard_Integer dim,dim2;
797 Standard_Integer d1,d2;
798 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
799 BSplSLib_DataContainer dc (UDegree, VDegree);
801 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
802 Poles,Weights,UKnots,VKnots,UMults,VMults,
803 u1,u2,d1,d2,rational,dc)) {
806 dim2 = (d2 + 1) << 2;
807 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
808 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
809 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
811 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
812 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
816 resVuu = result + 18;
818 resVuv = result + 12;
823 dim2 = (dim2 << 1) + dim2;
824 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
825 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
826 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
828 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
830 resVu = result + dim2;
832 if (UDegree <= 1) resVuu = BSplSLib_zero;
833 else resVuu = result + (dim2 << 1);
834 if (VDegree <= 1) resVvv = BSplSLib_zero;
835 else resVvv = result + 6;
836 resVuv = result + (d2 << 1) + d2 + 6;
842 dim2 = (d2 + 1) << 2;
843 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
844 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
845 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
847 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
848 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
853 resVvv = result + 18;
854 resVuv = result + 12;
859 dim2 = (dim2 << 1) + dim2;
860 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
861 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
862 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
864 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
867 resVv = result + dim2;
868 if (UDegree <= 1) resVuu = BSplSLib_zero;
869 else resVuu = result + 6;
870 if (VDegree <= 1) resVvv = BSplSLib_zero;
871 else resVvv = result + (dim2 << 1);
872 resVuv = result + (d2 << 1) + d2 + 6;
898 //=======================================================================
901 //=======================================================================
904 (const Standard_Real U,
905 const Standard_Real V,
906 const Standard_Integer UIndex,
907 const Standard_Integer VIndex,
908 const TColgp_Array2OfPnt& Poles,
909 const TColStd_Array2OfReal& Weights,
910 const TColStd_Array1OfReal& UKnots,
911 const TColStd_Array1OfReal& VKnots,
912 const TColStd_Array1OfInteger& UMults,
913 const TColStd_Array1OfInteger& VMults,
914 const Standard_Integer UDegree,
915 const Standard_Integer VDegree,
916 const Standard_Boolean URat,
917 const Standard_Boolean VRat,
918 const Standard_Boolean UPer,
919 const Standard_Boolean VPer,
931 Standard_Boolean rational;
932 // Standard_Integer k,dim,dim2;
933 Standard_Integer dim,dim2;
935 Standard_Integer d1,d2;
936 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
937 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
938 BSplSLib_DataContainer dc (UDegree, VDegree);
940 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
941 Poles,Weights,UKnots,VKnots,UMults,VMults,
942 u1,u2,d1,d2,rational,dc)) {
945 dim2 = (d2 + 1) << 2;
946 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
947 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
948 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
950 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
952 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
953 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
957 resVuu = result + 24;
959 resVuv = result + 15;
960 resVuuu = result + 36;
961 resVvvv = result + 9;
962 resVuuv = result + 27;
963 resVuvv = result + 18;
968 dim2 = (dim2 << 1) + dim2;
969 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
970 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
971 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
973 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
975 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
977 resVu = result + dim2;
980 resVuu = BSplSLib_zero;
981 resVuuv = BSplSLib_zero;
984 resVuu = result + (dim2 << 1);
985 resVuuv = result + (dim2 << 1) + 3;
988 resVvv = BSplSLib_zero;
989 resVuvv = BSplSLib_zero;
993 resVuvv = result + dim2 + 6;
995 resVuv = result + (d2 << 1) + d2 + 6;
996 if (UDegree <= 2) resVuuu = BSplSLib_zero;
997 else resVuuu = result + (dim2 << 1) + dim2;
998 if (VDegree <= 2) resVvvv = BSplSLib_zero;
999 else resVvvv = result + 9;
1005 dim2 = (d2 + 1) << 2;
1006 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1007 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1008 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1010 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1012 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1013 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1016 resVv = result + 12;
1017 resVuu = result + 6;
1018 resVvv = result + 24;
1019 resVuv = result + 15;
1020 resVuuu = result + 9;
1021 resVvvv = result + 36;
1022 resVuuv = result + 18;
1023 resVuvv = result + 27;
1028 dim2 = (dim2 << 1) + dim2;
1029 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1030 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1031 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1033 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1035 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1038 resVv = result + dim2;
1040 resVuu = BSplSLib_zero;
1041 resVuuv = BSplSLib_zero;
1044 resVuu = result + 6;
1045 resVuuv = result + dim2 + 6;
1048 resVvv = BSplSLib_zero;
1049 resVuvv = BSplSLib_zero;
1052 resVvv = result + (dim2 << 1);
1053 resVuvv = result + (dim2 << 1) + 3;
1055 resVuv = result + (d2 << 1) + d2 + 6;
1056 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1057 else resVuuu = result + 9;
1058 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1059 else resVvvv = result + (dim2 << 1) + dim2;
1063 P .SetX(result [0]);
1064 Vu .SetX(resVu [0]);
1065 Vv .SetX(resVv [0]);
1066 Vuu .SetX(resVuu [0]);
1067 Vvv .SetX(resVvv [0]);
1068 Vuv .SetX(resVuv [0]);
1069 Vuuu.SetX(resVuuu[0]);
1070 Vvvv.SetX(resVvvv[0]);
1071 Vuuv.SetX(resVuuv[0]);
1072 Vuvv.SetX(resVuvv[0]);
1074 P .SetY(result [1]);
1075 Vu .SetY(resVu [1]);
1076 Vv .SetY(resVv [1]);
1077 Vuu .SetY(resVuu [1]);
1078 Vvv .SetY(resVvv [1]);
1079 Vuv .SetY(resVuv [1]);
1080 Vuuu.SetY(resVuuu[1]);
1081 Vvvv.SetY(resVvvv[1]);
1082 Vuuv.SetY(resVuuv[1]);
1083 Vuvv.SetY(resVuvv[1]);
1085 P .SetZ(result [2]);
1086 Vu .SetZ(resVu [2]);
1087 Vv .SetZ(resVv [2]);
1088 Vuu .SetZ(resVuu [2]);
1089 Vvv .SetZ(resVvv [2]);
1090 Vuv .SetZ(resVuv [2]);
1091 Vuuu.SetZ(resVuuu[2]);
1092 Vvvv.SetZ(resVvvv[2]);
1093 Vuuv.SetZ(resVuuv[2]);
1094 Vuvv.SetZ(resVuvv[2]);
1097 //=======================================================================
1100 //=======================================================================
1103 (const Standard_Real U,
1104 const Standard_Real V,
1105 const Standard_Integer Nu,
1106 const Standard_Integer Nv,
1107 const Standard_Integer UIndex,
1108 const Standard_Integer VIndex,
1109 const TColgp_Array2OfPnt& Poles,
1110 const TColStd_Array2OfReal& Weights,
1111 const TColStd_Array1OfReal& UKnots,
1112 const TColStd_Array1OfReal& VKnots,
1113 const TColStd_Array1OfInteger& UMults,
1114 const TColStd_Array1OfInteger& VMults,
1115 const Standard_Integer UDegree,
1116 const Standard_Integer VDegree,
1117 const Standard_Boolean URat,
1118 const Standard_Boolean VRat,
1119 const Standard_Boolean UPer,
1120 const Standard_Boolean VPer,
1123 Standard_Boolean rational;
1124 Standard_Integer k,dim;
1125 Standard_Real u1,u2;
1126 Standard_Integer d1,d2;
1128 BSplSLib_DataContainer dc (UDegree, VDegree);
1129 Standard_Boolean ufirst = PrepareEval
1130 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1131 Poles,Weights,UKnots,VKnots,UMults,VMults,
1132 u1,u2,d1,d2,rational,dc);
1133 dim = rational ? 4 : 3;
1136 if ((Nu > UDegree) || (Nv > VDegree)) {
1144 Standard_Integer n1 = ufirst ? Nu : Nv;
1145 Standard_Integer n2 = ufirst ? Nv : Nu;
1147 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1149 for (k = 0; k <= Min(n1,d1); k++)
1150 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1152 Standard_Real *result;
1154 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1155 result = dc.ders; // because Standard_False ci-dessus.
1159 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1168 // Surface modifications
1170 // a surface is processed as a curve of curves.
1171 // i.e : a curve of parameter u where the current point is the set of poles
1175 //=======================================================================
1178 //=======================================================================
1180 void BSplSLib::Iso(const Standard_Real Param,
1181 const Standard_Boolean IsU,
1182 const TColgp_Array2OfPnt& Poles,
1183 const TColStd_Array2OfReal& Weights,
1184 const TColStd_Array1OfReal& Knots,
1185 const TColStd_Array1OfInteger& Mults,
1186 const Standard_Integer Degree,
1187 const Standard_Boolean Periodic,
1188 TColgp_Array1OfPnt& CPoles,
1189 TColStd_Array1OfReal& CWeights)
1191 Standard_Integer index = 0;
1192 Standard_Real u = Param;
1193 Standard_Boolean rational = &Weights != NULL;
1194 Standard_Integer dim = rational ? 4 : 3;
1196 // compute local knots
1198 BSplSLib_LocalArray locknots1 (2*Degree);
1199 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1200 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1202 index -= Knots.Lower() + Degree;
1204 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1207 // copy the local poles
1209 // Standard_Integer f1,l1,f2,l2,i,j,k;
1210 Standard_Integer f1,l1,f2,l2,i,j;
1213 f1 = Poles.LowerRow();
1214 l1 = Poles.UpperRow();
1215 f2 = Poles.LowerCol();
1216 l2 = Poles.UpperCol();
1219 f1 = Poles.LowerCol();
1220 l1 = Poles.UpperCol();
1221 f2 = Poles.LowerRow();
1222 l2 = Poles.UpperRow();
1225 BSplSLib_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim);
1227 Standard_Real w, *pole = locpoles;
1230 for (i = 0; i <= Degree; i++) {
1232 for (j = f2; j <= l2; j++) {
1234 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1236 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1237 pole[0] = P.X() * w;
1238 pole[1] = P.Y() * w;
1239 pole[2] = P.Z() * w;
1249 if (index > l1) index = f1;
1253 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1258 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1259 gp_Pnt& P = CPoles(i);
1261 CWeights(i) = w = pole[3];
1262 P.SetX( pole[0] / w);
1263 P.SetY( pole[1] / w);
1264 P.SetZ( pole[2] / w);
1274 // if the input is not rational but weights are wanted
1275 if (!rational && (&CWeights != NULL)) {
1277 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1282 //=======================================================================
1283 //function : Reverse
1285 //=======================================================================
1287 void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1288 const Standard_Integer Last,
1289 const Standard_Boolean UDirection)
1291 Standard_Integer i,j, l = Last;
1293 l = Poles.LowerRow() +
1294 (l - Poles.LowerRow())%(Poles.ColLength());
1295 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1296 Poles.LowerCol(), Poles.UpperCol());
1298 for (i = Poles.LowerRow(); i <= l; i++) {
1300 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1301 temp(l-i,j) = Poles(i,j);
1305 for (i = l+1; i <= Poles.UpperRow(); i++) {
1307 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1308 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1312 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1314 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1315 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1320 l = Poles.LowerCol() +
1321 (l - Poles.LowerCol())%(Poles.RowLength());
1322 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1323 0, Poles.RowLength()-1);
1325 for (j = Poles.LowerCol(); j <= l; j++) {
1327 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1328 temp(i,l-j) = Poles(i,j);
1332 for (j = l+1; j <= Poles.UpperCol(); j++) {
1334 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1335 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1339 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1341 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1342 Poles(i,j) = temp (i,j-Poles.LowerCol());
1348 //=======================================================================
1349 //function : Reverse
1351 //=======================================================================
1353 void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1354 const Standard_Integer Last,
1355 const Standard_Boolean UDirection)
1357 Standard_Integer i,j, l = Last;
1359 l = Weights.LowerRow() +
1360 (l - Weights.LowerRow())%(Weights.ColLength());
1361 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1362 Weights.LowerCol(), Weights.UpperCol());
1364 for (i = Weights.LowerRow(); i <= l; i++) {
1366 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1367 temp(l-i,j) = Weights(i,j);
1371 for (i = l+1; i <= Weights.UpperRow(); i++) {
1373 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1374 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1378 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1380 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1381 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1386 l = Weights.LowerCol() +
1387 (l - Weights.LowerCol())%(Weights.RowLength());
1388 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1389 0, Weights.RowLength()-1);
1391 for (j = Weights.LowerCol(); j <= l; j++) {
1393 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1394 temp(i,l-j) = Weights(i,j);
1398 for (j = l+1; j <= Weights.UpperCol(); j++) {
1400 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1401 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1405 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1407 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1408 Weights(i,j) = temp (i,j-Weights.LowerCol());
1414 //=======================================================================
1415 //function : IsRational
1417 //=======================================================================
1419 Standard_Boolean BSplSLib::IsRational
1420 (const TColStd_Array2OfReal& Weights,
1421 const Standard_Integer I1,
1422 const Standard_Integer I2,
1423 const Standard_Integer J1,
1424 const Standard_Integer J2,
1425 const Standard_Real Epsi)
1427 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1428 Standard_Integer i,j;
1429 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1430 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1432 for (i = I1 - fi; i < I2 - fi; i++) {
1434 for (j = J1 - fj; j < J2 - fj; j++) {
1435 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1436 return Standard_True;
1439 return Standard_False;
1442 //=======================================================================
1443 //function : SetPoles
1445 //=======================================================================
1447 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1448 TColStd_Array1OfReal& FP,
1449 const Standard_Boolean UDirection)
1451 Standard_Integer i, j, l = FP.Lower();
1452 Standard_Integer PLowerRow = Poles.LowerRow();
1453 Standard_Integer PUpperRow = Poles.UpperRow();
1454 Standard_Integer PLowerCol = Poles.LowerCol();
1455 Standard_Integer PUpperCol = Poles.UpperCol();
1458 for ( i = PLowerRow; i <= PUpperRow; i++) {
1460 for ( j = PLowerCol; j <= PUpperCol; j++) {
1461 const gp_Pnt& P = Poles.Value(i,j);
1470 for ( j = PLowerCol; j <= PUpperCol; j++) {
1472 for ( i = PLowerRow; i <= PUpperRow; i++) {
1473 const gp_Pnt& P = Poles.Value(i,j);
1482 //=======================================================================
1483 //function : SetPoles
1485 //=======================================================================
1487 void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1488 const TColStd_Array2OfReal& Weights,
1489 TColStd_Array1OfReal& FP,
1490 const Standard_Boolean UDirection)
1492 Standard_Integer i, j, l = FP.Lower();
1493 Standard_Integer PLowerRow = Poles.LowerRow();
1494 Standard_Integer PUpperRow = Poles.UpperRow();
1495 Standard_Integer PLowerCol = Poles.LowerCol();
1496 Standard_Integer PUpperCol = Poles.UpperCol();
1499 for ( i = PLowerRow; i <= PUpperRow; i++) {
1501 for ( j = PLowerCol; j <= PUpperCol; j++) {
1502 const gp_Pnt& P = Poles .Value(i,j);
1503 Standard_Real w = Weights.Value(i,j);
1504 FP(l) = P.X() * w; l++;
1505 FP(l) = P.Y() * w; l++;
1506 FP(l) = P.Z() * w; l++;
1513 for ( j = PLowerCol; j <= PUpperCol; j++) {
1515 for ( i = PLowerRow; i <= PUpperRow; i++) {
1516 const gp_Pnt& P = Poles .Value(i,j);
1517 Standard_Real w = Weights.Value(i,j);
1518 FP(l) = P.X() * w; l++;
1519 FP(l) = P.Y() * w; l++;
1520 FP(l) = P.Z() * w; l++;
1527 //=======================================================================
1528 //function : GetPoles
1530 //=======================================================================
1532 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1533 TColgp_Array2OfPnt& Poles,
1534 const Standard_Boolean UDirection)
1536 Standard_Integer i, j, l = FP.Lower();
1537 Standard_Integer PLowerRow = Poles.LowerRow();
1538 Standard_Integer PUpperRow = Poles.UpperRow();
1539 Standard_Integer PLowerCol = Poles.LowerCol();
1540 Standard_Integer PUpperCol = Poles.UpperCol();
1543 for ( i = PLowerRow; i <= PUpperRow; i++) {
1545 for ( j = PLowerCol; j <= PUpperCol; j++) {
1546 gp_Pnt& P = Poles.ChangeValue(i,j);
1555 for ( j = PLowerCol; j <= PUpperCol; j++) {
1557 for ( i = PLowerRow; i <= PUpperRow; i++) {
1558 gp_Pnt& P = Poles.ChangeValue(i,j);
1567 //=======================================================================
1568 //function : GetPoles
1570 //=======================================================================
1572 void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1573 TColgp_Array2OfPnt& Poles,
1574 TColStd_Array2OfReal& Weights,
1575 const Standard_Boolean UDirection)
1577 Standard_Integer i, j, l = FP.Lower();
1578 Standard_Integer PLowerRow = Poles.LowerRow();
1579 Standard_Integer PUpperRow = Poles.UpperRow();
1580 Standard_Integer PLowerCol = Poles.LowerCol();
1581 Standard_Integer PUpperCol = Poles.UpperCol();
1584 for ( i = PLowerRow; i <= PUpperRow; i++) {
1586 for ( j = PLowerCol; j <= PUpperCol; j++) {
1587 Standard_Real w = FP( l + 3);
1589 gp_Pnt& P = Poles.ChangeValue(i,j);
1590 P.SetX(FP(l) / w); l++;
1591 P.SetY(FP(l) / w); l++;
1592 P.SetZ(FP(l) / w); l++;
1599 for ( j = PLowerCol; j <= PUpperCol; j++) {
1601 for ( i = PLowerRow; i <= PUpperRow; i++) {
1602 Standard_Real w = FP( l + 3);
1604 gp_Pnt& P = Poles.ChangeValue(i,j);
1605 P.SetX(FP(l) / w); l++;
1606 P.SetY(FP(l) / w); l++;
1607 P.SetZ(FP(l) / w); l++;
1614 //=======================================================================
1615 //function : InsertKnots
1617 //=======================================================================
1619 void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1620 const Standard_Integer Degree,
1621 const Standard_Boolean Periodic,
1622 const TColgp_Array2OfPnt& Poles,
1623 const TColStd_Array2OfReal& Weights,
1624 const TColStd_Array1OfReal& Knots,
1625 const TColStd_Array1OfInteger& Mults,
1626 const TColStd_Array1OfReal& AddKnots,
1627 const TColStd_Array1OfInteger& AddMults,
1628 TColgp_Array2OfPnt& NewPoles,
1629 TColStd_Array2OfReal& NewWeights,
1630 TColStd_Array1OfReal& NewKnots,
1631 TColStd_Array1OfInteger& NewMults,
1632 const Standard_Real Epsilon,
1633 const Standard_Boolean Add )
1635 Standard_Boolean rational = &Weights != NULL;
1636 Standard_Integer dim = 3;
1637 if (rational) dim++;
1639 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1640 TColStd_Array1OfReal
1641 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1643 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1644 else SetPoles(Poles,poles,UDirection);
1647 dim *= Poles.RowLength();
1650 dim *= Poles.ColLength();
1652 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1653 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1656 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1657 else GetPoles(newpoles,NewPoles,UDirection);
1660 //=======================================================================
1661 //function : RemoveKnot
1663 //=======================================================================
1665 Standard_Boolean BSplSLib::RemoveKnot
1666 (const Standard_Boolean UDirection,
1667 const Standard_Integer Index,
1668 const Standard_Integer Mult,
1669 const Standard_Integer Degree,
1670 const Standard_Boolean Periodic,
1671 const TColgp_Array2OfPnt& Poles,
1672 const TColStd_Array2OfReal& Weights,
1673 const TColStd_Array1OfReal& Knots,
1674 const TColStd_Array1OfInteger& Mults,
1675 TColgp_Array2OfPnt& NewPoles,
1676 TColStd_Array2OfReal& NewWeights,
1677 TColStd_Array1OfReal& NewKnots,
1678 TColStd_Array1OfInteger& NewMults,
1679 const Standard_Real Tolerance)
1681 Standard_Boolean rational = &Weights != NULL;
1682 Standard_Integer dim = 3;
1683 if (rational) dim++;
1685 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1686 TColStd_Array1OfReal
1687 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1689 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1690 else SetPoles(Poles,poles,UDirection);
1693 dim *= Poles.RowLength();
1696 dim *= Poles.ColLength();
1699 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1700 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1702 return Standard_False;
1704 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1705 else GetPoles(newpoles,NewPoles,UDirection);
1706 return Standard_True;
1709 //=======================================================================
1710 //function : IncreaseDegree
1712 //=======================================================================
1714 void BSplSLib::IncreaseDegree
1715 (const Standard_Boolean UDirection,
1716 const Standard_Integer Degree,
1717 const Standard_Integer NewDegree,
1718 const Standard_Boolean Periodic,
1719 const TColgp_Array2OfPnt& Poles,
1720 const TColStd_Array2OfReal& Weights,
1721 const TColStd_Array1OfReal& Knots,
1722 const TColStd_Array1OfInteger& Mults,
1723 TColgp_Array2OfPnt& NewPoles,
1724 TColStd_Array2OfReal& NewWeights,
1725 TColStd_Array1OfReal& NewKnots,
1726 TColStd_Array1OfInteger& NewMults)
1728 Standard_Boolean rational = &Weights != NULL;
1729 Standard_Integer dim = 3;
1730 if (rational) dim++;
1732 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1733 TColStd_Array1OfReal
1734 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1736 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1737 else SetPoles(Poles,poles,UDirection);
1740 dim *= Poles.RowLength();
1743 dim *= Poles.ColLength();
1746 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1747 newpoles,NewKnots,NewMults);
1749 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1750 else GetPoles(newpoles,NewPoles,UDirection);
1753 //=======================================================================
1754 //function : Unperiodize
1756 //=======================================================================
1758 void BSplSLib::Unperiodize
1759 (const Standard_Boolean UDirection,
1760 const Standard_Integer Degree,
1761 const TColStd_Array1OfInteger& Mults,
1762 const TColStd_Array1OfReal& Knots,
1763 const TColgp_Array2OfPnt& Poles,
1764 const TColStd_Array2OfReal& Weights,
1765 TColStd_Array1OfInteger& NewMults,
1766 TColStd_Array1OfReal& NewKnots,
1767 TColgp_Array2OfPnt& NewPoles,
1768 TColStd_Array2OfReal& NewWeights)
1770 Standard_Boolean rational = &Weights != NULL;
1771 Standard_Integer dim = 3;
1772 if (rational) dim++;
1774 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1775 TColStd_Array1OfReal
1776 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1778 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1779 else SetPoles(Poles,poles,UDirection);
1782 dim *= Poles.RowLength();
1785 dim *= Poles.ColLength();
1788 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1789 NewMults, NewKnots, newpoles);
1791 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1792 else GetPoles(newpoles,NewPoles,UDirection);
1795 //=======================================================================
1796 //function : BuildCache
1797 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1798 // Cache : in case of a rational function the Poles are
1799 // stored in homogeneous form
1800 //=======================================================================
1802 void BSplSLib::BuildCache
1803 (const Standard_Real U,
1804 const Standard_Real V,
1805 const Standard_Real USpanDomain,
1806 const Standard_Real VSpanDomain,
1807 const Standard_Boolean UPeriodic,
1808 const Standard_Boolean VPeriodic,
1809 const Standard_Integer UDegree,
1810 const Standard_Integer VDegree,
1811 const Standard_Integer UIndex,
1812 const Standard_Integer VIndex,
1813 const TColStd_Array1OfReal& UFlatKnots,
1814 const TColStd_Array1OfReal& VFlatKnots,
1815 const TColgp_Array2OfPnt& Poles,
1816 const TColStd_Array2OfReal& Weights,
1817 TColgp_Array2OfPnt& CachePoles,
1818 TColStd_Array2OfReal& CacheWeights)
1820 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1821 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1822 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1823 if (&Weights != NULL)
1824 rational_u = rational_v = Standard_True;
1826 rational_u = rational_v = Standard_False;
1827 BSplSLib_DataContainer dc (UDegree, VDegree);
1843 (BSplCLib::NoMults()),
1844 (BSplCLib::NoMults()),
1854 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1856 for (kk = 0; kk <= d1 ; kk++)
1857 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1859 min_degree_domain = USpanDomain ;
1860 max_degree_domain = VSpanDomain ;
1863 min_degree_domain = VSpanDomain ;
1864 max_degree_domain = USpanDomain ;
1868 for (ii = 0 ; ii <= d2 ; ii++) {
1872 for (jj = 0 ; jj <= d1 ; jj++) {
1874 Index = jj * d2p1 + ii ;
1876 gp_Pnt& P = CachePoles(iii,jjj);
1877 f = factor[0] * factor[1];
1878 P.SetX( f * dc.poles[Index]); Index++;
1879 P.SetY( f * dc.poles[Index]); Index++;
1880 P.SetZ( f * dc.poles[Index]); Index++;
1881 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1882 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1884 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1888 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1890 for (kk = 0; kk <= d1 ; kk++)
1891 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
1893 min_degree_domain = USpanDomain ;
1894 max_degree_domain = VSpanDomain ;
1897 min_degree_domain = VSpanDomain ;
1898 max_degree_domain = USpanDomain ;
1902 for (ii = 0 ; ii <= d2 ; ii++) {
1906 for (jj = 0 ; jj <= d1 ; jj++) {
1908 Index = jj * d2p1 + ii ;
1909 Index = (Index << 1) + Index;
1910 gp_Pnt& P = CachePoles(iii,jjj);
1911 f = factor[0] * factor[1];
1912 P.SetX( f * dc.poles[Index]); Index++;
1913 P.SetY( f * dc.poles[Index]); Index++;
1914 P.SetZ( f * dc.poles[Index]);
1915 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1917 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1919 if (&Weights != NULL) {
1921 // means that PrepareEval did found out that the surface was
1922 // locally polynomial but since the surface is constructed
1923 // with some weights we need to set the weight polynomial to constant
1926 for (ii = 1 ; ii <= d2p1 ; ii++) {
1928 for (jj = 1 ; jj <= d1p1 ; jj++) {
1929 CacheWeights(ii,jj) = 0.0e0 ;
1932 CacheWeights(1,1) = 1.0e0 ;
1937 //=======================================================================
1938 //function : CacheD0
1939 //purpose : Evaluates the polynomial cache of the Bspline Curve
1941 //=======================================================================
1943 void BSplSLib::CacheD0(const Standard_Real UParameter,
1944 const Standard_Real VParameter,
1945 const Standard_Integer UDegree,
1946 const Standard_Integer VDegree,
1947 const Standard_Real UCacheParameter,
1948 const Standard_Real VCacheParameter,
1949 const Standard_Real USpanLenght,
1950 const Standard_Real VSpanLenght,
1951 const TColgp_Array2OfPnt& PolesArray,
1952 const TColStd_Array2OfReal& WeightsArray,
1956 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
1958 // the SpanLenght is the normalizing factor so that the polynomial is between
1971 PArray = (Standard_Real *)
1972 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
1974 myPoint = (Standard_Real *) &aPoint ;
1975 if (UDegree <= VDegree) {
1976 min_degree = UDegree ;
1977 max_degree = VDegree ;
1978 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
1979 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
1980 dimension = 3 * (UDegree + 1) ;
1983 min_degree = VDegree ;
1984 max_degree = UDegree ;
1985 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
1986 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
1987 dimension = 3 * (VDegree + 1) ;
1989 BSplSLib_LocalArray locpoles(dimension);
1991 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
1994 max_degree*dimension,
1998 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2001 (min_degree << 1) + min_degree,
2004 if (&WeightsArray != NULL) {
2005 dimension = min_degree + 1 ;
2007 WArray = (Standard_Real *)
2008 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2009 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2012 max_degree*dimension,
2016 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2022 inverse = 1.0e0/ inverse ;
2024 myPoint[0] *= inverse ;
2025 myPoint[1] *= inverse ;
2026 myPoint[2] *= inverse ;
2030 //=======================================================================
2031 //function : CacheD1
2032 //purpose : Evaluates the polynomial cache of the Bspline Curve
2034 //=======================================================================
2036 void BSplSLib::CacheD1(const Standard_Real UParameter,
2037 const Standard_Real VParameter,
2038 const Standard_Integer UDegree,
2039 const Standard_Integer VDegree,
2040 const Standard_Real UCacheParameter,
2041 const Standard_Real VCacheParameter,
2042 const Standard_Real USpanLenght,
2043 const Standard_Real VSpanLenght,
2044 const TColgp_Array2OfPnt& PolesArray,
2045 const TColStd_Array2OfReal& WeightsArray,
2051 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2053 // the SpanLenght is the normalizing factor so that the polynomial is between
2069 PArray = (Standard_Real *)
2070 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2071 Standard_Real local_poles_array[2][2][3],
2072 local_poles_and_weights_array[2][2][4],
2073 local_weights_array[2][2] ;
2074 Standard_Real * my_vec_min,
2077 my_point = (Standard_Real *) &aPoint ;
2079 // initialize in case of rational evaluation
2080 // because RationalDerivative will use all
2084 if (&WeightsArray != NULL) {
2086 local_poles_array [0][0][0] = 0.0e0 ;
2087 local_poles_array [0][0][1] = 0.0e0 ;
2088 local_poles_array [0][0][2] = 0.0e0 ;
2089 local_weights_array [0][0] = 0.0e0 ;
2090 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2091 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2092 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2093 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2095 local_poles_array [0][1][0] = 0.0e0 ;
2096 local_poles_array [0][1][1] = 0.0e0 ;
2097 local_poles_array [0][1][2] = 0.0e0 ;
2098 local_weights_array [0][1] = 0.0e0 ;
2099 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2100 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2101 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2102 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2104 local_poles_array [1][0][0] = 0.0e0 ;
2105 local_poles_array [1][0][1] = 0.0e0 ;
2106 local_poles_array [1][0][2] = 0.0e0 ;
2107 local_weights_array [1][0] = 0.0e0 ;
2108 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2109 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2110 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2111 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2113 local_poles_array [1][1][0] = 0.0e0 ;
2114 local_poles_array [1][1][1] = 0.0e0 ;
2115 local_poles_array [1][1][2] = 0.0e0 ;
2116 local_weights_array [1][1] = 0.0e0 ;
2117 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2118 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2119 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2120 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2123 if (UDegree <= VDegree) {
2124 min_degree = UDegree ;
2125 max_degree = VDegree ;
2126 inverse_min = 1.0e0/USpanLenght ;
2127 inverse_max = 1.0e0/VSpanLenght ;
2128 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2129 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2131 dimension = 3 * (UDegree + 1) ;
2132 my_vec_min = (Standard_Real *) &aVecU ;
2133 my_vec_max = (Standard_Real *) &aVecV ;
2136 min_degree = VDegree ;
2137 max_degree = UDegree ;
2138 inverse_min = 1.0e0/VSpanLenght ;
2139 inverse_max = 1.0e0/USpanLenght ;
2140 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2141 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2142 dimension = 3 * (VDegree + 1) ;
2143 my_vec_min = (Standard_Real *) &aVecV ;
2144 my_vec_max = (Standard_Real *) &aVecU ;
2147 BSplSLib_LocalArray locpoles (2 * dimension);
2149 PLib::EvalPolynomial(new_parameter[0],
2156 PLib::EvalPolynomial(new_parameter[1],
2161 local_poles_array[0][0][0]) ;
2162 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2165 (min_degree << 1) + min_degree,
2166 locpoles[dimension],
2167 local_poles_array[1][0][0]) ;
2169 if (&WeightsArray != NULL) {
2170 dimension = min_degree + 1 ;
2172 WArray = (Standard_Real *)
2173 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2174 PLib::EvalPolynomial(new_parameter[0],
2181 PLib::EvalPolynomial(new_parameter[1],
2186 local_weights_array[0][0]) ;
2187 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2191 locpoles[dimension],
2192 local_weights_array[1][0]) ;
2194 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2195 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2196 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2197 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2199 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2200 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2201 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2202 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2204 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2205 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2206 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2207 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2209 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2210 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2211 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2212 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2214 BSplSLib::RationalDerivative(1,
2218 local_poles_and_weights_array[0][0][0],
2219 local_poles_array[0][0][0]) ;
2222 my_point [0] = local_poles_array [0][0][0] ;
2223 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2224 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2226 my_point [1] = local_poles_array [0][0][1] ;
2227 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2228 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2230 my_point [2] = local_poles_array [0][0][2] ;
2231 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2232 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2235 //=======================================================================
2236 //function : CacheD2
2237 //purpose : Evaluates the polynomial cache of the Bspline Curve
2239 //=======================================================================
2241 void BSplSLib::CacheD2(const Standard_Real UParameter,
2242 const Standard_Real VParameter,
2243 const Standard_Integer UDegree,
2244 const Standard_Integer VDegree,
2245 const Standard_Real UCacheParameter,
2246 const Standard_Real VCacheParameter,
2247 const Standard_Real USpanLenght,
2248 const Standard_Real VSpanLenght,
2249 const TColgp_Array2OfPnt& PolesArray,
2250 const TColStd_Array2OfReal& WeightsArray,
2259 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2261 // the SpanLenght is the normalizing factor so that the polynomial is between
2278 PArray = (Standard_Real *)
2279 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2280 Standard_Real local_poles_array[3][3][3],
2281 local_poles_and_weights_array[3][3][4],
2282 local_weights_array[3][3] ;
2283 Standard_Real * my_vec_min,
2289 my_point = (Standard_Real *) &aPoint ;
2292 // initialize in case the min and max degree are less than 2
2294 local_poles_array[0][0][0] = 0.0e0 ;
2295 local_poles_array[0][0][1] = 0.0e0 ;
2296 local_poles_array[0][0][2] = 0.0e0 ;
2297 local_poles_array[0][1][0] = 0.0e0 ;
2298 local_poles_array[0][1][1] = 0.0e0 ;
2299 local_poles_array[0][1][2] = 0.0e0 ;
2300 local_poles_array[0][2][0] = 0.0e0 ;
2301 local_poles_array[0][2][1] = 0.0e0 ;
2302 local_poles_array[0][2][2] = 0.0e0 ;
2304 local_poles_array[1][0][0] = 0.0e0 ;
2305 local_poles_array[1][0][1] = 0.0e0 ;
2306 local_poles_array[1][0][2] = 0.0e0 ;
2307 local_poles_array[1][1][0] = 0.0e0 ;
2308 local_poles_array[1][1][1] = 0.0e0 ;
2309 local_poles_array[1][1][2] = 0.0e0 ;
2310 local_poles_array[1][2][0] = 0.0e0 ;
2311 local_poles_array[1][2][1] = 0.0e0 ;
2312 local_poles_array[1][2][2] = 0.0e0 ;
2314 local_poles_array[2][0][0] = 0.0e0 ;
2315 local_poles_array[2][0][1] = 0.0e0 ;
2316 local_poles_array[2][0][2] = 0.0e0 ;
2317 local_poles_array[2][1][0] = 0.0e0 ;
2318 local_poles_array[2][1][1] = 0.0e0 ;
2319 local_poles_array[2][1][2] = 0.0e0 ;
2320 local_poles_array[2][2][0] = 0.0e0 ;
2321 local_poles_array[2][2][1] = 0.0e0 ;
2322 local_poles_array[2][2][2] = 0.0e0 ;
2324 // initialize in case of rational evaluation
2325 // because RationalDerivative will use all
2329 if (&WeightsArray != NULL) {
2331 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2332 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2333 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2334 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2335 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2336 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2337 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2338 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2339 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2341 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2342 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2343 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2344 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2345 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2346 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2347 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2348 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2349 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2351 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2352 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2353 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2354 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2355 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2356 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2357 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2358 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2359 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2361 local_poles_and_weights_array[0][0][3] =
2362 local_weights_array[0][0] = 0.0e0 ;
2363 local_poles_and_weights_array[0][1][3] =
2364 local_weights_array[0][1] = 0.0e0 ;
2365 local_poles_and_weights_array[0][2][3] =
2366 local_weights_array[0][2] = 0.0e0 ;
2367 local_poles_and_weights_array[1][0][3] =
2368 local_weights_array[1][0] = 0.0e0 ;
2369 local_poles_and_weights_array[1][1][3] =
2370 local_weights_array[1][1] = 0.0e0 ;
2371 local_poles_and_weights_array[1][2][3] =
2372 local_weights_array[1][2] = 0.0e0 ;
2373 local_poles_and_weights_array[2][0][3] =
2374 local_weights_array[2][0] = 0.0e0 ;
2375 local_poles_and_weights_array[2][1][3] =
2376 local_weights_array[2][1] = 0.0e0 ;
2377 local_poles_and_weights_array[2][2][3] =
2378 local_weights_array[2][2] = 0.0e0 ;
2381 if (UDegree <= VDegree) {
2382 min_degree = UDegree ;
2383 max_degree = VDegree ;
2384 inverse_min = 1.0e0/USpanLenght ;
2385 inverse_max = 1.0e0/VSpanLenght ;
2386 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2387 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2389 dimension = 3 * (UDegree + 1) ;
2390 my_vec_min = (Standard_Real *) &aVecU ;
2391 my_vec_max = (Standard_Real *) &aVecV ;
2392 my_vec_min_min = (Standard_Real *) &aVecUU ;
2393 my_vec_min_max = (Standard_Real *) &aVecUV ;
2394 my_vec_max_max = (Standard_Real *) &aVecVV ;
2397 min_degree = VDegree ;
2398 max_degree = UDegree ;
2399 inverse_min = 1.0e0/VSpanLenght ;
2400 inverse_max = 1.0e0/USpanLenght ;
2401 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2402 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2403 dimension = 3 * (VDegree + 1) ;
2404 my_vec_min = (Standard_Real *) &aVecV ;
2405 my_vec_max = (Standard_Real *) &aVecU ;
2406 my_vec_min_min = (Standard_Real *) &aVecVV ;
2407 my_vec_min_max = (Standard_Real *) &aVecUV ;
2408 my_vec_max_max = (Standard_Real *) &aVecUU ;
2411 BSplSLib_LocalArray locpoles (3 * dimension);
2414 // initialize in case min or max degree are less than 2
2416 Standard_Integer MinIndMax = 2;
2417 if ( max_degree < 2) MinIndMax = max_degree;
2418 Standard_Integer MinIndMin = 2;
2419 if ( min_degree < 2) MinIndMin = min_degree;
2421 index = MinIndMax * dimension ;
2423 for (ii = MinIndMax ; ii < 3 ; ii++) {
2425 for (kk = 0 ; kk < dimension ; kk++) {
2426 locpoles[index] = 0.0e0 ;
2431 PLib::EvalPolynomial(new_parameter[0],
2438 PLib::EvalPolynomial(new_parameter[1],
2443 local_poles_array[0][0][0]) ;
2444 PLib::EvalPolynomial(new_parameter[1],
2448 locpoles[dimension],
2449 local_poles_array[1][0][0]) ;
2451 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2454 (min_degree << 1) + min_degree,
2455 locpoles[dimension + dimension],
2456 local_poles_array[2][0][0]) ;
2458 if (&WeightsArray != NULL) {
2459 dimension = min_degree + 1 ;
2461 WArray = (Standard_Real *)
2462 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2463 PLib::EvalPolynomial(new_parameter[0],
2470 PLib::EvalPolynomial(new_parameter[1],
2475 local_weights_array[0][0]) ;
2476 PLib::EvalPolynomial(new_parameter[1],
2480 locpoles[dimension],
2481 local_weights_array[1][0]) ;
2482 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2486 locpoles[dimension + dimension],
2487 local_weights_array[2][0]) ;
2490 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2491 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2492 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2493 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2494 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2495 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2496 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2497 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2498 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2500 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2501 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2502 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2503 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2504 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2505 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2506 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2507 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2508 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2510 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2511 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2512 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2513 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2514 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2515 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2516 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2517 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2518 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2521 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2522 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2523 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2524 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2525 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2526 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2527 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2528 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2529 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2531 BSplSLib::RationalDerivative(2,
2535 local_poles_and_weights_array[0][0][0],
2536 local_poles_array[0][0][0]) ;
2540 Standard_Real minmin = inverse_min * inverse_min;
2541 Standard_Real minmax = inverse_min * inverse_max;
2542 Standard_Real maxmax = inverse_max * inverse_max;
2544 my_point [0] = local_poles_array [0][0][0] ;
2545 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2546 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2547 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2548 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2549 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2551 my_point [1] = local_poles_array [0][0][1] ;
2552 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2553 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2554 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2555 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2556 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2558 my_point [2] = local_poles_array [0][0][2] ;
2559 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2560 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2561 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2562 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2563 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2566 //=======================================================================
2567 //function : MovePoint
2568 //purpose : Find the new poles which allows an old point (with a
2569 // given u and v as parameters) to reach a new position
2570 //=======================================================================
2572 void BSplSLib::MovePoint (const Standard_Real U,
2573 const Standard_Real V,
2574 const gp_Vec& Displ,
2575 const Standard_Integer UIndex1,
2576 const Standard_Integer UIndex2,
2577 const Standard_Integer VIndex1,
2578 const Standard_Integer VIndex2,
2579 const Standard_Integer UDegree,
2580 const Standard_Integer VDegree,
2581 const Standard_Boolean Rational,
2582 const TColgp_Array2OfPnt& Poles,
2583 const TColStd_Array2OfReal& Weights,
2584 const TColStd_Array1OfReal& UFlatKnots,
2585 const TColStd_Array1OfReal& VFlatKnots,
2586 Standard_Integer& UFirstIndex,
2587 Standard_Integer& ULastIndex,
2588 Standard_Integer& VFirstIndex,
2589 Standard_Integer& VLastIndex,
2590 TColgp_Array2OfPnt& NewPoles)
2592 // calculate the UBSplineBasis in the parameter U
2593 Standard_Integer UFirstNonZeroBsplineIndex;
2594 math_Matrix UBSplineBasis(1, 1,
2596 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2601 UFirstNonZeroBsplineIndex,
2603 // calculate the VBSplineBasis in the parameter V
2604 Standard_Integer VFirstNonZeroBsplineIndex;
2605 math_Matrix VBSplineBasis(1, 1,
2607 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2612 VFirstNonZeroBsplineIndex,
2614 if (ErrorCod1 || ErrorCod2) {
2622 // find the span which is predominant for parameter U
2623 UFirstIndex = UFirstNonZeroBsplineIndex;
2624 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2625 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2626 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2628 Standard_Real maxValue = 0.0;
2629 Standard_Integer i, ukk1=0, ukk2;
2631 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2632 if (UBSplineBasis(1,i) > maxValue) {
2633 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2634 maxValue = UBSplineBasis(1,i);
2638 // find a ukk2 if symetriy
2640 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2641 if ((ukk1+1) <= ULastIndex) {
2642 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2647 // find the span which is predominant for parameter V
2648 VFirstIndex = VFirstNonZeroBsplineIndex;
2649 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2651 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2652 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2655 Standard_Integer j, vkk1=0, vkk2;
2657 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2658 if (VBSplineBasis(1,j) > maxValue) {
2659 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2660 maxValue = VBSplineBasis(1,j);
2664 // find a vkk2 if symetriy
2666 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2667 if ((vkk1+1) <= VLastIndex) {
2668 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2673 // compute the vector of displacement
2674 Standard_Real D1 = 0.0;
2675 Standard_Real D2 = 0.0;
2676 Standard_Real hN, Coef, DvalU, DvalV;
2678 Standard_Integer ii, jj;
2680 for (i = 1; i <= UDegree+1; i++) {
2681 ii = i + UFirstNonZeroBsplineIndex - 1;
2685 else if (ii > ukk2) {
2692 for (j = 1; j <= VDegree+1; j++) {
2693 jj = j + VFirstNonZeroBsplineIndex - 1;
2695 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2699 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2701 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2705 else if (jj > vkk2) {
2711 D1 += 1./(DvalU + DvalV + 1.) * hN;
2723 // compute the new poles
2725 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2729 else if (i > ukk2) {
2736 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2737 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2741 else if (j > vkk2) {
2747 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2750 NewPoles(i,j) = Poles(i,j);
2756 //=======================================================================
2757 // function : Resolution
2758 // purpose : this computes an estimate for the maximum of the
2759 // partial derivatives both in U and in V
2762 // The calculation resembles at the calculation of curves with
2763 // additional index for the control point. Let Si,j be the
2764 // control points for ls surface and Di,j the weights.
2765 // The checking of upper bounds for the partial derivatives
2766 // will be omitted and Su is the next upper bound in the polynomial case :
2770 // | Si,j - Si-1,j |
2771 // d * Max | ------------- |
2772 // i = 2,n | ti+d - ti |
2776 // and in the rational case :
2780 // Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2781 // Max Max d * -----------------------------------------------
2782 // k=1,n i dans Rj ti+d - ti
2784 // ----------------------------------------------------------------------
2792 // with Rj = {j-d, ...., j+d+d+1}.
2795 //=======================================================================
2797 void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
2798 const TColStd_Array2OfReal& Weights,
2799 const TColStd_Array1OfReal& UKnots,
2800 const TColStd_Array1OfReal& VKnots,
2801 const TColStd_Array1OfInteger& UMults,
2802 const TColStd_Array1OfInteger& VMults,
2803 const Standard_Integer UDegree,
2804 const Standard_Integer VDegree,
2805 const Standard_Boolean URational,
2806 const Standard_Boolean VRational,
2807 const Standard_Boolean UPeriodic,
2808 const Standard_Boolean VPeriodic,
2809 const Standard_Real Tolerance3D,
2810 Standard_Real& UTolerance,
2811 Standard_Real& VTolerance)
2813 Standard_Real Wij,Wmj,Wji,Wjm;
2814 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2815 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2816 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2817 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2819 max_derivative[0] = max_derivative[1] = 0.0e0 ;
2821 Standard_Integer PRowLength, PColLength;
2822 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2823 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2824 Standard_Integer num_poles[2],num_flat_knots[2];
2827 BSplCLib::KnotSequenceLength(UMults,
2831 BSplCLib::KnotSequenceLength(VMults,
2834 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
2835 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
2836 BSplCLib::KnotSequence(UKnots,
2841 BSplCLib::KnotSequence(VKnots,
2846 PRowLength = Poles.RowLength();
2847 PColLength = Poles.ColLength();
2848 if (URational || VRational) {
2849 Standard_Integer Wsize = PRowLength * PColLength;
2850 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2851 min_weights = WG[0];
2853 for (ii = 1 ; ii < Wsize ; ii++) {
2855 if (min_weights > min) min_weights = min;
2858 Standard_Integer UD1 = UDegree + 1;
2859 Standard_Integer VD1 = VDegree + 1;
2860 num_poles[0] = num_flat_knots[0] - UD1;
2861 num_poles[1] = num_flat_knots[1] - VD1;
2862 poles_length[0] = PColLength;
2863 poles_length[1] = PRowLength;
2865 Standard_Integer UD2 = UDegree << 1;
2866 Standard_Integer VD2 = VDegree << 1;
2868 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2869 ii_index = (ii - 1) % poles_length[0] + 1 ;
2870 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2871 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2872 inverse = 1.0e0 / inverse ;
2873 lower[0] = ii - UD1;
2874 if (lower[0] < 1) lower[0] = 1;
2875 upper[0] = ii + UD2 + 1;
2876 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2878 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2879 jj_index = (jj - 1) % poles_length[1] + 1 ;
2880 lower[1] = jj - VD1;
2881 if (lower[1] < 1) lower[1] = 1;
2882 upper[1] = jj + VD2 + 1;
2883 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
2884 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
2885 Wij = Weights.Value(ii_index,jj_index);
2886 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
2887 Wmj = Weights.Value(ii_minus,jj_index);
2895 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
2896 pp_index = (pp - 1) % poles_length[0] + 1 ;
2898 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
2900 qq_index = (qq - 1) % poles_length[1] + 1 ;
2901 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
2905 factor = (Xpq - Xij) * Wij;
2906 factor -= (Xpq - Xmj) * Wmj;
2907 if (factor < 0) factor = - factor;
2909 factor = (Ypq - Yij) * Wij;
2910 factor -= (Ypq - Ymj) * Wmj;
2911 if (factor < 0) factor = - factor;
2913 factor = (Zpq - Zij) * Wij;
2914 factor -= (Zpq - Zmj) * Wmj;
2915 if (factor < 0) factor = - factor;
2918 if (max_derivative[0] < value) max_derivative[0] = value ;
2923 max_derivative[0] /= min_weights ;
2927 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2928 ii_index = (ii - 1) % poles_length[0] + 1 ;
2929 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2930 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2931 inverse = 1.0e0 / inverse ;
2933 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2934 jj_index = (jj - 1) % poles_length[1] + 1 ;
2936 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
2937 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
2938 factor = Pij.X() - Pmj.X();
2939 if (factor < 0) factor = - factor;
2941 factor = Pij.Y() - Pmj.Y();
2942 if (factor < 0) factor = - factor;
2944 factor = Pij.Z() - Pmj.Z();
2945 if (factor < 0) factor = - factor;
2948 if (max_derivative[0] < value) max_derivative[0] = value ;
2952 max_derivative[0] *= UDegree ;
2954 Standard_Integer UD2 = UDegree << 1;
2955 Standard_Integer VD2 = VDegree << 1;
2957 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
2958 ii_index = (ii - 1) % poles_length[1] + 1 ;
2959 ii_minus = (ii - 2) % poles_length[1] + 1 ;
2960 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
2961 inverse = 1.0e0 / inverse ;
2962 lower[0] = ii - VD1;
2963 if (lower[0] < 1) lower[0] = 1;
2964 upper[0] = ii + VD2 + 1;
2965 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
2967 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
2968 jj_index = (jj - 1) % poles_length[0] + 1 ;
2969 lower[1] = jj - UD1;
2970 if (lower[1] < 1) lower[1] = 1;
2971 upper[1] = jj + UD2 + 1;
2972 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
2973 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
2974 Wji = Weights.Value(jj_index,ii_index);
2975 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
2976 Wjm = Weights.Value(jj_index,ii_minus);
2984 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
2985 pp_index = (pp - 1) % poles_length[1] + 1 ;
2987 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
2989 qq_index = (qq - 1) % poles_length[0] + 1 ;
2990 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
2994 factor = (Xqp - Xji) * Wji;
2995 factor -= (Xqp - Xjm) * Wjm;
2996 if (factor < 0) factor = - factor;
2998 factor = (Yqp - Yji) * Wji;
2999 factor -= (Yqp - Yjm) * Wjm;
3000 if (factor < 0) factor = - factor;
3002 factor = (Zqp - Zji) * Wji;
3003 factor -= (Zqp - Zjm) * Wjm;
3004 if (factor < 0) factor = - factor;
3007 if (max_derivative[1] < value) max_derivative[1] = value ;
3012 max_derivative[1] /= min_weights ;
3016 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3017 ii_index = (ii - 1) % poles_length[1] + 1 ;
3018 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3019 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3020 inverse = 1.0e0 / inverse ;
3022 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3023 jj_index = (jj - 1) % poles_length[0] + 1 ;
3025 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3026 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3027 factor = Pji.X() - Pjm.X() ;
3028 if (factor < 0) factor = - factor;
3030 factor = Pji.Y() - Pjm.Y() ;
3031 if (factor < 0) factor = - factor;
3033 factor = Pji.Z() - Pjm.Z() ;
3034 if (factor < 0) factor = - factor;
3037 if (max_derivative[1] < value) max_derivative[1] = value ;
3041 max_derivative[1] *= VDegree ;
3042 max_derivative[0] *= M_SQRT2 ;
3043 max_derivative[1] *= M_SQRT2 ;
3044 if(max_derivative[0] && max_derivative[1]) {
3045 UTolerance = Tolerance3D / max_derivative[0] ;
3046 VTolerance = Tolerance3D / max_derivative[1] ;
3049 UTolerance=VTolerance=0.0;
3051 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3056 //=======================================================================
3057 //function : Interpolate
3059 //=======================================================================
3061 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3062 const Standard_Integer VDegree,
3063 const TColStd_Array1OfReal& UFlatKnots,
3064 const TColStd_Array1OfReal& VFlatKnots,
3065 const TColStd_Array1OfReal& UParameters,
3066 const TColStd_Array1OfReal& VParameters,
3067 TColgp_Array2OfPnt& Poles,
3068 TColStd_Array2OfReal& Weights,
3069 Standard_Integer& InversionProblem)
3071 Standard_Integer ii, jj, ll, kk, dimension;
3072 Standard_Integer ULength = UParameters.Length();
3073 Standard_Integer VLength = VParameters.Length();
3074 Standard_Real * poles_array;
3076 // extraction of iso u
3077 dimension = 4*ULength;
3078 TColStd_Array2OfReal Points(1, VLength,
3081 Handle(TColStd_HArray1OfInteger) ContactOrder =
3082 new (TColStd_HArray1OfInteger)(1, VLength);
3083 ContactOrder->Init(0);
3085 for (ii=1; ii <= VLength; ii++) {
3087 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3088 Points(ii,ll) = Poles(jj, ii).X();
3089 Points(ii,ll+1) = Poles(jj, ii).Y();
3090 Points(ii,ll+2) = Poles(jj, ii).Z();
3091 Points(ii,ll+3) = Weights(jj,ii) ;
3095 // interpolation of iso u
3096 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3097 BSplCLib::Interpolate(VDegree,
3100 ContactOrder->Array1(),
3104 if (InversionProblem != 0) return;
3106 // extraction of iso v
3108 dimension = VLength*4;
3109 TColStd_Array2OfReal IsoPoles(1, ULength,
3112 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3113 ContactOrder->Init(0);
3114 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3116 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3118 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3119 IsoPoles (ii,ll) = Points(jj, kk);
3120 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3121 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3122 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3125 // interpolation of iso v
3126 BSplCLib::Interpolate(UDegree,
3129 ContactOrder->Array1(),
3136 for (ii=1; ii <= ULength; ii++) {
3138 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3139 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3140 Poles.SetValue(ii, jj, Pnt);
3141 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3146 //=======================================================================
3147 //function : Interpolate
3149 //=======================================================================
3151 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3152 const Standard_Integer VDegree,
3153 const TColStd_Array1OfReal& UFlatKnots,
3154 const TColStd_Array1OfReal& VFlatKnots,
3155 const TColStd_Array1OfReal& UParameters,
3156 const TColStd_Array1OfReal& VParameters,
3157 TColgp_Array2OfPnt& Poles,
3158 Standard_Integer& InversionProblem)
3160 Standard_Integer ii, jj, ll, kk, dimension;
3161 Standard_Integer ULength = UParameters.Length();
3162 Standard_Integer VLength = VParameters.Length();
3163 Standard_Real * poles_array;
3165 // extraction of iso u
3166 dimension = 3*ULength;
3167 TColStd_Array2OfReal Points(1, VLength,
3170 Handle(TColStd_HArray1OfInteger) ContactOrder =
3171 new (TColStd_HArray1OfInteger)(1, VLength);
3172 ContactOrder->Init(0);
3174 for (ii=1; ii <= VLength; ii++) {
3176 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3177 Points(ii,ll) = Poles(jj, ii).X();
3178 Points(ii,ll+1) = Poles(jj, ii).Y();
3179 Points(ii,ll+2) = Poles(jj, ii).Z();
3183 // interpolation of iso u
3184 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3185 BSplCLib::Interpolate(VDegree,
3188 ContactOrder->Array1(),
3192 if (InversionProblem != 0) return;
3194 // extraction of iso v
3196 dimension = VLength*3;
3197 TColStd_Array2OfReal IsoPoles(1, ULength,
3200 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3201 ContactOrder->Init(0);
3202 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3204 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3206 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3207 IsoPoles (ii,ll) = Points(jj, kk);
3208 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3209 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3212 // interpolation of iso v
3213 BSplCLib::Interpolate(UDegree,
3216 ContactOrder->Array1(),
3223 for (ii=1; ii <= ULength; ii++) {
3225 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3226 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3227 Poles.SetValue(ii, jj, Pnt);
3232 //=======================================================================
3233 //function : FunctionMultiply
3235 //=======================================================================
3237 void BSplSLib::FunctionMultiply
3238 (const BSplSLib_EvaluatorFunction& Function,
3239 const Standard_Integer UBSplineDegree,
3240 const Standard_Integer VBSplineDegree,
3241 const TColStd_Array1OfReal& UBSplineKnots,
3242 const TColStd_Array1OfReal& VBSplineKnots,
3243 const TColStd_Array1OfInteger & UMults,
3244 const TColStd_Array1OfInteger & VMults,
3245 const TColgp_Array2OfPnt& Poles,
3246 const TColStd_Array2OfReal& Weights,
3247 const TColStd_Array1OfReal& UFlatKnots,
3248 const TColStd_Array1OfReal& VFlatKnots,
3249 const Standard_Integer UNewDegree,
3250 const Standard_Integer VNewDegree,
3251 TColgp_Array2OfPnt& NewNumerator,
3252 TColStd_Array2OfReal& NewDenominator,
3253 Standard_Integer& Status)
3255 Standard_Integer num_uparameters,
3260 Standard_Real result ;
3262 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3263 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3264 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3265 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3267 if ((NewNumerator.ColLength() == num_uparameters) &&
3268 (NewNumerator.RowLength() == num_vparameters) &&
3269 (NewDenominator.ColLength() == num_uparameters) &&
3270 (NewDenominator.RowLength() == num_vparameters)) {
3273 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3277 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3281 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3283 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3284 HomogeneousD0(UParameters(ii),
3300 NewDenominator(ii,jj),
3301 NewNumerator(ii,jj)) ;
3303 Function.Evaluate (0,
3309 Standard_ConstructionError::Raise();
3311 gp_Pnt& P = NewNumerator(ii,jj);
3312 P.SetX(P.X() * result);
3313 P.SetY(P.Y() * result);
3314 P.SetZ(P.Z() * result);
3315 NewDenominator(ii,jj) *= result ;
3318 Interpolate(UNewDegree,
3329 Standard_ConstructionError::Raise();