projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022312: Translation of french commentaries in OCCT files
[occt.git] / src / BSplSLib / BSplSLib.cxx
CommitLineData
7fd59977 1// File: BSplSLib.cxx
2// Created : Mon Aug 26 07:39:13 1991
3// Author: JCV
4
5// Modifed RLE Aug 93 - Complete rewrite
6// xab 21-Mar-95 implemented cache mecanism
7// pmn 25-09-96 Interpolation
8// jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
9// pmn 07-10-96 : Correction de DN dans le cas rationnal.
10// pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
11
12#define No_Standard_RangeError
13#define No_Standard_OutOfRange
14
15#include <BSplSLib.ixx>
16#include <PLib.hxx>
17#include <BSplCLib.hxx>
18#include <TColgp_Array2OfXYZ.hxx>
19#include <TColgp_Array1OfXYZ.hxx>
20#include <TColStd_HArray1OfInteger.hxx>
21#include <Standard_NotImplemented.hxx>
22#include <Standard_ConstructionError.hxx>
23#include <math_Matrix.hxx>
24
25// for null derivatives
26static Standard_Real BSplSLib_zero[3] = {0.,0.,0.};
27#ifdef WNT
28#define M_SQRT2 1.4142135623730950488016887
29#endif
30
31//=======================================================================
32//struct : BSplCLib_DataContainer
33//purpose: Auxiliary structure providing buffers for poles and knots used in
34// evaluation of bspline (allocated in the stack)
35//=======================================================================
36
37struct BSplSLib_DataContainer
38{
39 BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
40 {
41 if ( UDegree > BSplCLib::MaxDegree() ||
42 VDegree > BSplCLib::MaxDegree() ||
43 BSplCLib::MaxDegree() != 25 )
44 Standard_OutOfRange::Raise ("BSplCLib: bspline degree is greater than maximum supported");
45 }
46
47 Standard_Real poles[4*(25+1)*(25+1)];
48 Standard_Real knots1[2*25];
49 Standard_Real knots2[2*25];
50 Standard_Real ders[48];
51};
52
53//=======================================================================
54//class : BSplSLib_LocalArray
55//purpose: Auxiliary class optimizing creation of array buffer for
56// evaluation of bspline (using stack allocation for small arrays)
57//=======================================================================
58
59#define LOCARRAY_BUFFER 1024
60class BSplSLib_LocalArray
61{
62 public:
63 BSplSLib_LocalArray (Standard_Integer Size)
64 : myPtr(myBuffer)
65 {
66 if ( Size > LOCARRAY_BUFFER )
67 myPtr = (Standard_Real*)Standard::Allocate (Size * sizeof(Standard_Real));
68 }
69
70 ~BSplSLib_LocalArray ()
71 {
72 if ( myPtr != myBuffer )
73 Standard::Free (*(Standard_Address*)&myPtr);
74 }
75
76 operator Standard_Real* () { return myPtr; }
77
78 private:
79 Standard_Real myBuffer[LOCARRAY_BUFFER];
80 Standard_Real* myPtr;
81};
82
83//**************************************************************************
84// Evaluation methods
85//**************************************************************************
86
87//=======================================================================
88//function : RationalDerivative
89//purpose : computes the rational derivatives when whe have the
90// the derivatives of the homogeneous numerator and the
91// the derivatives of the denominator
92//=======================================================================
93
94void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
95 const Standard_Integer VDeg,
96 const Standard_Integer N,
97 const Standard_Integer M,
98 Standard_Real& HDerivatives,
99 Standard_Real& RDerivatives,
100 const Standard_Boolean All)
101{
102 //
103 // if All is True all derivatives are computed. if Not only
104 // the requested N, M is computed
105 //
106 // Numerator(u,v)
107 // let f(u,v) be a rational function = ------------------
108 // Denominator(u,v)
109 //
110 //
111 // Let (N,M) the order of the derivatives we want : then since
112 // we have :
113 //
114 // Numerator = f * Denominator
115 //
116 // we derive :
117 //
118 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
119 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
120 // (0,0) ( p<N q<M p q )
121 // Denominator
122 //
123 // with :
124 //
125 // ( N ) ( M )
126 // a = ( ) ( )
127 // p q ( p ) ( q )
128 //
129 //
130 // HDerivatives is an array where derivatives are stored in the following form
131 // Numerator is assumee to have 3 functions that is a vector of dimension
132 // 3
133 //
134 // (0,0) (0,0) (0, DegV) (0, DegV)
135 // Numerator Denominator ... Numerator Denominator
136 //
137 // (1,0) (1,0) (1, DegV) (1, DegV)
138 // Numerator Denominator ... Numerator Denominator
139 //
140 // ...........................................................
141 //
142 //
143 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
144 // Numerator Denominator ... Numerator Denominator
145 //
146 //
147 Standard_Integer ii,jj,pp,qq,index,index1,index2;
148 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
149 Standard_Integer MinN,MinN1,MinM,MinM1;
150 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
151
152 M1 = M + 1;
153 N1 = N + 1;
154 ii = N1 * M1;
155 M3 = (M1 << 1) + M1;
156 M4 = (VDeg + 1) << 2;
157
158 BSplSLib_LocalArray StoreDerivatives (All ? 0 : ii * 3);
159 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
160 BSplSLib_LocalArray StoreW (ii);
161 Standard_Real *HomogeneousArray = &HDerivatives;
162 Standard_Real denominator,Pii,Pip,Pjq;
163
164 denominator = 1.0e0 / HomogeneousArray[3];
165 index_u = 0;
166 index_u1 = 0;
167 if (UDeg < N) MinN = UDeg;
168 else MinN = N;
169 if (VDeg < M) MinM = VDeg;
170 else MinM = M;
171 MinN1 = MinN + 1;
172 MinM1 = MinM + 1;
173 iiM1 = - M1;
174
175 for (ii = 0 ; ii < MinN1 ; ii++) {
176 iiM1 += M1;
177 index_v = index_u;
178 index_v1 = index_u1;
179 index_w = iiM1;
180
181 for (jj = 0 ; jj < MinM1 ; jj++) {
182 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
183 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
184 RArray[index_v] = HomogeneousArray[index_v1]; index_v++; index_v1++;
185 StoreW[index_w] = HomogeneousArray[index_v1]; index_w++; index_v1++;
186 }
187
188 for (jj = MinM1 ; jj < M1 ; jj++) {
189 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
190 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
191 RArray[index_v] = 0.0e0 ; index_v++; index_v1++;
192 StoreW[index_w] = HomogeneousArray[index_v1]; index_w++; index_v1++;
193 }
194 index_u1 += M4;
195 index_u += M3;
196 }
197 index_v = MinN1 * M3;
198 index_w = MinN1 * M1;
199
200 for (ii = MinN1 ; ii < N1 ; ii++) {
201
202 for (jj = 0 ; jj < M1 ; jj++) {
203 RArray[index_v] = 0.0e0; index_v++;
204 RArray[index_v] = 0.0e0; index_v++;
205 RArray[index_v] = 0.0e0; index_v++;
206 StoreW[index_w] = 0.0e0; index_w++;
207 }
208 }
209
0d969553 210 // --------------- Calculation ----------------
7fd59977 211
212 iiM1 = - M1;
213 iiM3 = - M3;
214 PLib::Binomial(N);
215 PLib::Binomial(M);
216
217 for (ii = 0 ; ii <= N ; ii++) {
218 iiM1 += M1;
219 iiM3 += M3;
220 index1 = iiM3 - 3;
221 jjM1 = iiM1;
222
223 for (jj = 0 ; jj <= M ; jj++) {
224 jjM1 ++;
225 ppM1 = - M1;
226 ppM3 = - M3;
227 index1 += 3;
228
229 for (pp = 0 ; pp < ii ; pp++) {
230 ppM1 += M1;
231 ppM3 += M3;
232 index = ppM3;
233 index2 = jjM1 - ppM1;
234 Pip = PLib::Bin(ii,pp);
235
236 for (qq = 0 ; qq <= jj ; qq++) {
237 index2--;
238 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
239 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
240 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
241 RArray[index1] -= Pjq * RArray[index]; index++;
242 index1 -= 2;
243 }
244 }
245 index = iiM3;
246 index2 = jj + 1;
247 Pii = PLib::Bin(ii,ii);
248
249 for (qq = 0 ; qq < jj ; qq++) {
250 index2--;
251 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
252 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
253 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
254 RArray[index1] -= Pjq * RArray[index]; index++;
255 index1 -= 2;
256 }
257 RArray[index1] *= denominator; index1++;
258 RArray[index1] *= denominator; index1++;
259 RArray[index1] *= denominator;
260 index1 -= 2;
261 }
262 }
263 if (!All) {
264 RArray = &RDerivatives;
265 index = N * M1 + M;
266 index = (index << 1) + index;
267 RArray[0] = StoreDerivatives[index]; index++;
268 RArray[1] = StoreDerivatives[index]; index++;
269 RArray[2] = StoreDerivatives[index];
270 }
271}
272
273//=======================================================================
274//function : PrepareEval
275//purpose :
276//=======================================================================
277
278//
279// PrepareEval :
280//
0d969553 281// Prepare all data for computing points :
7fd59977 282// local arrays of knots
283// local array of poles (multiplied by the weights if rational)
284//
285// The first direction to compute (smaller degree) is returned
286// and the poles are stored according to this direction.
287
288static Standard_Boolean PrepareEval
289(const Standard_Real U,
290 const Standard_Real V,
291 const Standard_Integer Uindex,
292 const Standard_Integer Vindex,
293 const Standard_Integer UDegree,
294 const Standard_Integer VDegree,
295 const Standard_Boolean URat,
296 const Standard_Boolean VRat,
297 const Standard_Boolean UPer,
298 const Standard_Boolean VPer,
299 const TColgp_Array2OfPnt& Poles,
300 const TColStd_Array2OfReal& Weights,
301 const TColStd_Array1OfReal& UKnots,
302 const TColStd_Array1OfReal& VKnots,
303 const TColStd_Array1OfInteger& UMults,
304 const TColStd_Array1OfInteger& VMults,
305 Standard_Real& u1, // first parameter to use
306 Standard_Real& u2, // second parameter to use
307 Standard_Integer& d1, // first degree
308 Standard_Integer& d2, // second degree
309 Standard_Boolean& rational,
310 BSplSLib_DataContainer& dc)
311{
312 rational = URat || VRat;
313 Standard_Integer uindex = Uindex;
314 Standard_Integer vindex = Vindex;
315 Standard_Integer UKLower = UKnots.Lower();
316 Standard_Integer UKUpper = UKnots.Upper();
317 Standard_Integer VKLower = VKnots.Lower();
318 Standard_Integer VKUpper = VKnots.Upper();
319 if (UDegree <= VDegree) {
320 // compute the indices
321 if (uindex < UKLower || uindex > UKUpper)
322 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
323 else u1 = U;
324 if (vindex < VKLower || vindex > VKUpper)
325 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
326 else u2 = V;
327 // get the knots
328 d1 = UDegree;
329 d2 = VDegree;
330 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
331 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
332 if (&UMults == NULL) uindex -= UKLower + UDegree;
333 else uindex = BSplCLib::PoleIndex
334 (UDegree,uindex,UPer,UMults);
335 if (&VMults == NULL) vindex -= VKLower + VDegree;
336 else vindex = BSplCLib::PoleIndex
337 (VDegree,vindex,VPer,VMults);
338 // get the poles
339// Standard_Integer i,j,k,ip,jp;
340 Standard_Integer i,j,ip,jp;
341 Standard_Real w, *pole = dc.poles;
342 d1 = UDegree;
343 d2 = VDegree;
344 Standard_Integer PLowerRow = Poles.LowerRow();
345 Standard_Integer PUpperRow = Poles.UpperRow();
346 Standard_Integer PLowerCol = Poles.LowerCol();
347 Standard_Integer PUpperCol = Poles.UpperCol();
348 if (rational) { // verify if locally non rational
349 rational = Standard_False;
350 ip = PLowerRow + uindex;
351 jp = PLowerCol + vindex;
352 w = Weights.Value(ip,jp);
353 Standard_Real eps = Epsilon(w);
354 Standard_Real dw;
355
356 for (i = 0; i <= UDegree && !rational; i++) {
357 jp = PLowerCol + vindex;
358
359 for (j = 0; j <= VDegree && !rational; j++) {
360 dw = Weights.Value(ip,jp) - w;
361 if (dw < 0) dw = - dw;
362 rational = dw > eps;
363 jp++;
364 if (jp > PUpperCol) jp = PLowerCol;
365 }
366 ip++;
367 if (ip > PUpperRow) ip = PLowerRow;
368 }
369 }
370 // copy the poles
371 ip = PLowerRow + uindex;
372 if (rational) {
373
374 for (i = 0; i <= d1; i++) {
375 jp = PLowerCol + vindex;
376
377 for (j = 0; j <= d2; j++) {
378 const gp_Pnt& P = Poles .Value(ip,jp);
379 pole[3] = w = Weights.Value(ip,jp);
380 pole[0] = P.X() * w;
381 pole[1] = P.Y() * w;
382 pole[2] = P.Z() * w;
383 pole += 4;
384 jp++;
385 if (jp > PUpperCol) jp = PLowerCol;
386 }
387 ip++;
388 if (ip > PUpperRow) ip = PLowerRow;
389 }
390 }
391 else {
392
393 for (i = 0; i <= d1; i++) {
394 jp = PLowerCol + vindex;
395
396 for (j = 0; j <= d2; j++) {
397 const gp_Pnt& P = Poles.Value(ip,jp);
398 pole[0] = P.X();
399 pole[1] = P.Y();
400 pole[2] = P.Z();
401 pole += 3;
402 jp++;
403 if (jp > PUpperCol) jp = PLowerCol;
404 }
405 ip++;
406 if (ip > PUpperRow) ip = PLowerRow;
407 }
408 }
409 return Standard_True;
410 }
411 else {
412 // compute the indices
413 if (uindex < UKLower || uindex > UKUpper)
414 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
415 else u2 = U;
416 if (vindex < VKLower || vindex > VKUpper)
417 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
418 else u1 = V;
419 // get the knots
420 d2 = UDegree;
421 d1 = VDegree;
422 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
423 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
424 if (&UMults == NULL) uindex -= UKLower + UDegree;
425 else uindex = BSplCLib::PoleIndex
426 (UDegree,uindex,UPer,UMults);
427 if (&VMults == NULL) vindex -= VKLower + VDegree;
428 else vindex = BSplCLib::PoleIndex
429 (VDegree,vindex,VPer,VMults);
430 // get the poles
431// Standard_Integer i,j,k,ip,jp;
432 Standard_Integer i,j,ip,jp;
433 Standard_Real w, *pole = dc.poles;
434 d1 = VDegree;
435 d2 = UDegree;
436 Standard_Integer PLowerRow = Poles.LowerRow();
437 Standard_Integer PUpperRow = Poles.UpperRow();
438 Standard_Integer PLowerCol = Poles.LowerCol();
439 Standard_Integer PUpperCol = Poles.UpperCol();
440 if (rational) { // verify if locally non rational
441 rational = Standard_False;
442 ip = PLowerRow + uindex;
443 jp = PLowerCol + vindex;
444 w = Weights.Value(ip,jp);
445 Standard_Real eps = Epsilon(w);
446 Standard_Real dw;
447
448 for (i = 0; i <= UDegree && !rational; i++) {
449 jp = PLowerCol + vindex;
450
451 for (j = 0; j <= VDegree && !rational; j++) {
452 dw = Weights.Value(ip,jp) - w;
453 if (dw < 0) dw = - dw;
454 rational = dw > eps;
455 jp++;
456 if (jp > PUpperCol) jp = PLowerCol;
457 }
458 ip++;
459 if (ip > PUpperRow) ip = PLowerRow;
460 }
461 }
462 // copy the poles
463 jp = PLowerCol + vindex;
464 if (rational) {
465
466 for (i = 0; i <= d1; i++) {
467 ip = PLowerRow + uindex;
468
469 for (j = 0; j <= d2; j++) {
470 const gp_Pnt& P = Poles .Value(ip,jp);
471 pole[3] = w = Weights.Value(ip,jp);
472 pole[0] = P.X() * w;
473 pole[1] = P.Y() * w;
474 pole[2] = P.Z() * w;
475 pole += 4;
476 ip++;
477 if (ip > PUpperRow) ip = PLowerRow;
478 }
479 jp++;
480 if (jp > PUpperCol) jp = PLowerCol;
481 }
482 }
483 else {
484
485 for (i = 0; i <= d1; i++) {
486 ip = PLowerRow + uindex;
487
488 for (j = 0; j <= d2; j++) {
489 const gp_Pnt& P = Poles.Value(ip,jp);
490 pole[0] = P.X();
491 pole[1] = P.Y();
492 pole[2] = P.Z();
493 pole += 3;
494 ip++;
495 if (ip > PUpperRow) ip = PLowerRow;
496 }
497 jp++;
498 if (jp > PUpperCol) jp = PLowerCol;
499 }
500 }
501 return Standard_False;
502 }
503}
504
505//=======================================================================
506//function : D0
507//purpose :
508//=======================================================================
509
510void BSplSLib::D0
511(const Standard_Real U,
512 const Standard_Real V,
513 const Standard_Integer UIndex,
514 const Standard_Integer VIndex,
515 const TColgp_Array2OfPnt& Poles,
516 const TColStd_Array2OfReal& Weights,
517 const TColStd_Array1OfReal& UKnots,
518 const TColStd_Array1OfReal& VKnots,
519 const TColStd_Array1OfInteger& UMults,
520 const TColStd_Array1OfInteger& VMults,
521 const Standard_Integer UDegree,
522 const Standard_Integer VDegree,
523 const Standard_Boolean URat,
524 const Standard_Boolean VRat,
525 const Standard_Boolean UPer,
526 const Standard_Boolean VPer,
527 gp_Pnt& P)
528{
529// Standard_Integer k ;
530 Standard_Real W ;
531 HomogeneousD0(U,
532 V,
533 UIndex,
534 VIndex,
535 Poles,
536 Weights,
537 UKnots,
538 VKnots,
539 UMults,
540 VMults,
541 UDegree,
542 VDegree,
543 URat,
544 VRat,
545 UPer,
546 VPer,
547 W,
548 P) ;
549 P.SetX(P.X() / W);
550 P.SetY(P.Y() / W);
551 P.SetZ(P.Z() / W);
552}
553
554//=======================================================================
555//function : D0
556//purpose :
557//=======================================================================
558
559void BSplSLib::HomogeneousD0
560(const Standard_Real U,
561 const Standard_Real V,
562 const Standard_Integer UIndex,
563 const Standard_Integer VIndex,
564 const TColgp_Array2OfPnt& Poles,
565 const TColStd_Array2OfReal& Weights,
566 const TColStd_Array1OfReal& UKnots,
567 const TColStd_Array1OfReal& VKnots,
568 const TColStd_Array1OfInteger& UMults,
569 const TColStd_Array1OfInteger& VMults,
570 const Standard_Integer UDegree,
571 const Standard_Integer VDegree,
572 const Standard_Boolean URat,
573 const Standard_Boolean VRat,
574 const Standard_Boolean UPer,
575 const Standard_Boolean VPer,
576 Standard_Real & W,
577 gp_Pnt& P)
578{
579 Standard_Boolean rational;
580// Standard_Integer k,dim;
581 Standard_Integer dim;
582 Standard_Real u1,u2;
583 Standard_Integer d1,d2;
584 W = 1.0e0 ;
585
586 BSplSLib_DataContainer dc (UDegree, VDegree);
587 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
588 Poles,Weights,UKnots,VKnots,UMults,VMults,
589 u1,u2,d1,d2,rational,dc);
590 if (rational) {
591 dim = 4;
592 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
593 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
594 W = dc.poles[3];
595 P.SetX(dc.poles[0]);
596 P.SetY(dc.poles[1]);
597 P.SetZ(dc.poles[2]);
598 }
599 else {
600 dim = 3;
601 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
602 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
603 P.SetX(dc.poles[0]);
604 P.SetY(dc.poles[1]);
605 P.SetZ(dc.poles[2]);
606 }
607}
608
609//=======================================================================
610//function : D1
611//purpose :
612//=======================================================================
613
614void BSplSLib::D1
615(const Standard_Real U,
616 const Standard_Real V,
617 const Standard_Integer UIndex,
618 const Standard_Integer VIndex,
619 const TColgp_Array2OfPnt& Poles,
620 const TColStd_Array2OfReal& Weights,
621 const TColStd_Array1OfReal& UKnots,
622 const TColStd_Array1OfReal& VKnots,
623 const TColStd_Array1OfInteger& UMults,
624 const TColStd_Array1OfInteger& VMults,
625 const Standard_Integer UDegree,
626 const Standard_Integer VDegree,
627 const Standard_Boolean URat,
628 const Standard_Boolean VRat,
629 const Standard_Boolean UPer,
630 const Standard_Boolean VPer,
631 gp_Pnt& P,
632 gp_Vec& Vu,
633 gp_Vec& Vv)
634{
635 Standard_Boolean rational;
636// Standard_Integer k,dim,dim2;
637 Standard_Integer dim,dim2;
638 Standard_Real u1,u2;
639 Standard_Integer d1,d2;
640 Standard_Real *result, *resVu, *resVv;
641 BSplSLib_DataContainer dc (UDegree, VDegree);
642 if (PrepareEval
643 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
644 Poles,Weights,UKnots,VKnots,UMults,VMults,
645 u1,u2,d1,d2,rational,dc)) {
646 if (rational) {
647 dim = 4;
648 dim2 = (d2 + 1) << 2;
649 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
650 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
651 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
652 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
653 result = dc.ders;
654 resVu = result + 6;
655 resVv = result + 3;
656 }
657 else {
658 dim = 3;
659 dim2 = d2 + 1;
660 dim2 = (dim2 << 1) + dim2;
661 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
662 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
663 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
664 result = dc.poles;
665 resVu = result + dim2;
666 resVv = result + 3;
667 }
668 }
669 else {
670 if (rational) {
671 dim = 4;
672 dim2 = (d2 + 1) << 2;
673 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
674 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
675 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
676 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
677 result = dc.ders;
678 resVu = result + 3;
679 resVv = result + 6;
680 }
681 else {
682 dim = 3;
683 dim2 = d2 + 1;
684 dim2 = (dim2 << 1) + dim2;
685 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
686 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
687 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
688 result = dc.poles;
689 resVu = result + 3;
690 resVv = result + dim2;
691 }
692 }
693
694 P .SetX(result[0]);
695 Vu.SetX(resVu [0]);
696 Vv.SetX(resVv [0]);
697
698 P .SetY(result[1]);
699 Vu.SetY(resVu [1]);
700 Vv.SetY(resVv [1]);
701
702 P .SetZ(result[2]);
703 Vu.SetZ(resVu [2]);
704 Vv.SetZ(resVv [2]);
705}
706
707//=======================================================================
708//function : D1
709//purpose :
710//=======================================================================
711
712void BSplSLib::HomogeneousD1
713(const Standard_Real U,
714 const Standard_Real V,
715 const Standard_Integer UIndex,
716 const Standard_Integer VIndex,
717 const TColgp_Array2OfPnt& Poles,
718 const TColStd_Array2OfReal& Weights,
719 const TColStd_Array1OfReal& UKnots,
720 const TColStd_Array1OfReal& VKnots,
721 const TColStd_Array1OfInteger& UMults,
722 const TColStd_Array1OfInteger& VMults,
723 const Standard_Integer UDegree,
724 const Standard_Integer VDegree,
725 const Standard_Boolean URat,
726 const Standard_Boolean VRat,
727 const Standard_Boolean UPer,
728 const Standard_Boolean VPer,
729 gp_Pnt& N,
730 gp_Vec& Nu,
731 gp_Vec& Nv,
732 Standard_Real& D,
733 Standard_Real& Du,
734 Standard_Real& Dv)
735{
736 Standard_Boolean rational;
737// Standard_Integer k,dim;
738 Standard_Integer dim;
739 Standard_Real u1,u2;
740 Standard_Integer d1,d2;
741
742 D = 1.0e0 ;
743 Du = 0.0e0 ;
744 Dv = 0.0e0 ;
745 BSplSLib_DataContainer dc (UDegree, VDegree);
746 Standard_Boolean ufirst = PrepareEval
747 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
748 Poles,Weights,UKnots,VKnots,UMults,VMults,
749 u1,u2,d1,d2,rational,dc);
750 dim = rational ? 4 : 3;
751
752 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
753 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
754 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
755
756 Standard_Real *result, *resVu, *resVv;
757 result = dc.poles;
758 resVu = result + (ufirst ? dim*(d2+1) : dim);
759 resVv = result + (ufirst ? dim : dim*(d2+1));
760
761 N .SetX(result[0]);
762 Nu.SetX(resVu [0]);
763 Nv.SetX(resVv [0]);
764
765 N .SetY(result[1]);
766 Nu.SetY(resVu [1]);
767 Nv.SetY(resVv [1]);
768
769 N .SetZ(result[2]);
770 Nu.SetZ(resVu [2]);
771 Nv.SetZ(resVv [2]);
772
773 if (rational) {
774 D = result[3];
775 Du = resVu [3];
776 Dv = resVv [3];
777 }
778}
779
780//=======================================================================
781//function : D2
782//purpose :
783//=======================================================================
784
785void BSplSLib::D2
786(const Standard_Real U,
787 const Standard_Real V,
788 const Standard_Integer UIndex,
789 const Standard_Integer VIndex,
790 const TColgp_Array2OfPnt& Poles,
791 const TColStd_Array2OfReal& Weights,
792 const TColStd_Array1OfReal& UKnots,
793 const TColStd_Array1OfReal& VKnots,
794 const TColStd_Array1OfInteger& UMults,
795 const TColStd_Array1OfInteger& VMults,
796 const Standard_Integer UDegree,
797 const Standard_Integer VDegree,
798 const Standard_Boolean URat,
799 const Standard_Boolean VRat,
800 const Standard_Boolean UPer,
801 const Standard_Boolean VPer,
802 gp_Pnt& P,
803 gp_Vec& Vu,
804 gp_Vec& Vv,
805 gp_Vec& Vuu,
806 gp_Vec& Vvv,
807 gp_Vec& Vuv)
808{
809 Standard_Boolean rational;
810// Standard_Integer k,dim,dim2;
811 Standard_Integer dim,dim2;
812 Standard_Real u1,u2;
813 Standard_Integer d1,d2;
814 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
815 BSplSLib_DataContainer dc (UDegree, VDegree);
816 if (PrepareEval
817 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
818 Poles,Weights,UKnots,VKnots,UMults,VMults,
819 u1,u2,d1,d2,rational,dc)) {
820 if (rational) {
821 dim = 4;
822 dim2 = (d2 + 1) << 2;
823 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
824 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
825 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
826 if (d1 > 1)
827 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
828 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
829 result = dc.ders;
830 resVu = result + 9;
831 resVv = result + 3;
832 resVuu = result + 18;
833 resVvv = result + 6;
834 resVuv = result + 12;
835 }
836 else {
837 dim = 3;
838 dim2 = d2 + 1;
839 dim2 = (dim2 << 1) + dim2;
840 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
841 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
842 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
843 if (d1 > 1)
844 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
845 result = dc.poles;
846 resVu = result + dim2;
847 resVv = result + 3;
848 if (UDegree <= 1) resVuu = BSplSLib_zero;
849 else resVuu = result + (dim2 << 1);
850 if (VDegree <= 1) resVvv = BSplSLib_zero;
851 else resVvv = result + 6;
852 resVuv = result + (d2 << 1) + d2 + 6;
853 }
854 }
855 else {
856 if (rational) {
857 dim = 4;
858 dim2 = (d2 + 1) << 2;
859 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
860 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
861 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
862 if (d1 > 1)
863 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
864 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
865 result = dc.ders;
866 resVu = result + 3;
867 resVv = result + 9;
868 resVuu = result + 6;
869 resVvv = result + 18;
870 resVuv = result + 12;
871 }
872 else {
873 dim = 3;
874 dim2 = d2 + 1;
875 dim2 = (dim2 << 1) + dim2;
876 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
877 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
878 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
879 if (d1 > 1)
880 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
881 result = dc.poles;
882 resVu = result + 3;
883 resVv = result + dim2;
884 if (UDegree <= 1) resVuu = BSplSLib_zero;
885 else resVuu = result + 6;
886 if (VDegree <= 1) resVvv = BSplSLib_zero;
887 else resVvv = result + (dim2 << 1);
888 resVuv = result + (d2 << 1) + d2 + 6;
889 }
890 }
891
892 P .SetX(result[0]);
893 Vu .SetX(resVu [0]);
894 Vv .SetX(resVv [0]);
895 Vuu.SetX(resVuu[0]);
896 Vvv.SetX(resVvv[0]);
897 Vuv.SetX(resVuv[0]);
898
899 P .SetY(result[1]);
900 Vu .SetY(resVu [1]);
901 Vv .SetY(resVv [1]);
902 Vuu.SetY(resVuu[1]);
903 Vvv.SetY(resVvv[1]);
904 Vuv.SetY(resVuv[1]);
905
906 P .SetZ(result[2]);
907 Vu .SetZ(resVu [2]);
908 Vv .SetZ(resVv [2]);
909 Vuu.SetZ(resVuu[2]);
910 Vvv.SetZ(resVvv[2]);
911 Vuv.SetZ(resVuv[2]);
912}
913
914//=======================================================================
915//function : D3
916//purpose :
917//=======================================================================
918
919void BSplSLib::D3
920(const Standard_Real U,
921 const Standard_Real V,
922 const Standard_Integer UIndex,
923 const Standard_Integer VIndex,
924 const TColgp_Array2OfPnt& Poles,
925 const TColStd_Array2OfReal& Weights,
926 const TColStd_Array1OfReal& UKnots,
927 const TColStd_Array1OfReal& VKnots,
928 const TColStd_Array1OfInteger& UMults,
929 const TColStd_Array1OfInteger& VMults,
930 const Standard_Integer UDegree,
931 const Standard_Integer VDegree,
932 const Standard_Boolean URat,
933 const Standard_Boolean VRat,
934 const Standard_Boolean UPer,
935 const Standard_Boolean VPer,
936 gp_Pnt& P,
937 gp_Vec& Vu,
938 gp_Vec& Vv,
939 gp_Vec& Vuu,
940 gp_Vec& Vvv,
941 gp_Vec& Vuv,
942 gp_Vec& Vuuu,
943 gp_Vec& Vvvv,
944 gp_Vec& Vuuv,
945 gp_Vec& Vuvv)
946{
947 Standard_Boolean rational;
948// Standard_Integer k,dim,dim2;
949 Standard_Integer dim,dim2;
950 Standard_Real u1,u2;
951 Standard_Integer d1,d2;
952 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
953 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
954 BSplSLib_DataContainer dc (UDegree, VDegree);
955 if (PrepareEval
956 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
957 Poles,Weights,UKnots,VKnots,UMults,VMults,
958 u1,u2,d1,d2,rational,dc)) {
959 if (rational) {
960 dim = 4;
961 dim2 = (d2 + 1) << 2;
962 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
963 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
964 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
965 if (d1 > 1)
966 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
967 if (d1 > 2)
968 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
969 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
970 result = dc.ders;
971 resVu = result + 12;
972 resVv = result + 3;
973 resVuu = result + 24;
974 resVvv = result + 6;
975 resVuv = result + 15;
976 resVuuu = result + 36;
977 resVvvv = result + 9;
978 resVuuv = result + 27;
979 resVuvv = result + 18;
980 }
981 else {
982 dim = 3;
983 dim2 = (d2 + 1);
984 dim2 = (dim2 << 1) + dim2;
985 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
986 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
987 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
988 if (d1 > 1)
989 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
990 if (d1 > 2)
991 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
992 result = dc.poles;
993 resVu = result + dim2;
994 resVv = result + 3;
995 if (UDegree <= 1) {
996 resVuu = BSplSLib_zero;
997 resVuuv = BSplSLib_zero;
998 }
999 else {
1000 resVuu = result + (dim2 << 1);
1001 resVuuv = result + (dim2 << 1) + 3;
1002 }
1003 if (VDegree <= 1) {
1004 resVvv = BSplSLib_zero;
1005 resVuvv = BSplSLib_zero;
1006 }
1007 else {
1008 resVvv = result + 6;
1009 resVuvv = result + dim2 + 6;
1010 }
1011 resVuv = result + (d2 << 1) + d2 + 6;
1012 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1013 else resVuuu = result + (dim2 << 1) + dim2;
1014 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1015 else resVvvv = result + 9;
1016 }
1017 }
1018 else {
1019 if (rational) {
1020 dim = 4;
1021 dim2 = (d2 + 1) << 2;
1022 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1023 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1024 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1025 if (d1 > 1)
1026 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1027 if (d1 > 2)
1028 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1029 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1030 result = dc.ders;
1031 resVu = result + 3;
1032 resVv = result + 12;
1033 resVuu = result + 6;
1034 resVvv = result + 24;
1035 resVuv = result + 15;
1036 resVuuu = result + 9;
1037 resVvvv = result + 36;
1038 resVuuv = result + 18;
1039 resVuvv = result + 27;
1040 }
1041 else {
1042 dim = 3;
1043 dim2 = (d2 + 1);
1044 dim2 = (dim2 << 1) + dim2;
1045 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1046 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1047 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1048 if (d1 > 1)
1049 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1050 if (d1 > 2)
1051 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1052 result = dc.poles;
1053 resVu = result + 3;
1054 resVv = result + dim2;
1055 if (UDegree <= 1) {
1056 resVuu = BSplSLib_zero;
1057 resVuuv = BSplSLib_zero;
1058 }
1059 else {
1060 resVuu = result + 6;
1061 resVuuv = result + dim2 + 6;
1062 }
1063 if (VDegree <= 1) {
1064 resVvv = BSplSLib_zero;
1065 resVuvv = BSplSLib_zero;
1066 }
1067 else {
1068 resVvv = result + (dim2 << 1);
1069 resVuvv = result + (dim2 << 1) + 3;
1070 }
1071 resVuv = result + (d2 << 1) + d2 + 6;
1072 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1073 else resVuuu = result + 9;
1074 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1075 else resVvvv = result + (dim2 << 1) + dim2;
1076 }
1077 }
1078
1079 P .SetX(result [0]);
1080 Vu .SetX(resVu [0]);
1081 Vv .SetX(resVv [0]);
1082 Vuu .SetX(resVuu [0]);
1083 Vvv .SetX(resVvv [0]);
1084 Vuv .SetX(resVuv [0]);
1085 Vuuu.SetX(resVuuu[0]);
1086 Vvvv.SetX(resVvvv[0]);
1087 Vuuv.SetX(resVuuv[0]);
1088 Vuvv.SetX(resVuvv[0]);
1089
1090 P .SetY(result [1]);
1091 Vu .SetY(resVu [1]);
1092 Vv .SetY(resVv [1]);
1093 Vuu .SetY(resVuu [1]);
1094 Vvv .SetY(resVvv [1]);
1095 Vuv .SetY(resVuv [1]);
1096 Vuuu.SetY(resVuuu[1]);
1097 Vvvv.SetY(resVvvv[1]);
1098 Vuuv.SetY(resVuuv[1]);
1099 Vuvv.SetY(resVuvv[1]);
1100
1101 P .SetZ(result [2]);
1102 Vu .SetZ(resVu [2]);
1103 Vv .SetZ(resVv [2]);
1104 Vuu .SetZ(resVuu [2]);
1105 Vvv .SetZ(resVvv [2]);
1106 Vuv .SetZ(resVuv [2]);
1107 Vuuu.SetZ(resVuuu[2]);
1108 Vvvv.SetZ(resVvvv[2]);
1109 Vuuv.SetZ(resVuuv[2]);
1110 Vuvv.SetZ(resVuvv[2]);
1111}
1112
1113//=======================================================================
1114//function : DN
1115//purpose :
1116//=======================================================================
1117
1118void BSplSLib::DN
1119(const Standard_Real U,
1120 const Standard_Real V,
1121 const Standard_Integer Nu,
1122 const Standard_Integer Nv,
1123 const Standard_Integer UIndex,
1124 const Standard_Integer VIndex,
1125 const TColgp_Array2OfPnt& Poles,
1126 const TColStd_Array2OfReal& Weights,
1127 const TColStd_Array1OfReal& UKnots,
1128 const TColStd_Array1OfReal& VKnots,
1129 const TColStd_Array1OfInteger& UMults,
1130 const TColStd_Array1OfInteger& VMults,
1131 const Standard_Integer UDegree,
1132 const Standard_Integer VDegree,
1133 const Standard_Boolean URat,
1134 const Standard_Boolean VRat,
1135 const Standard_Boolean UPer,
1136 const Standard_Boolean VPer,
1137 gp_Vec& Vn)
1138{
1139 Standard_Boolean rational;
1140 Standard_Integer k,dim;
1141 Standard_Real u1,u2;
1142 Standard_Integer d1,d2;
1143
1144 BSplSLib_DataContainer dc (UDegree, VDegree);
1145 Standard_Boolean ufirst = PrepareEval
1146 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1147 Poles,Weights,UKnots,VKnots,UMults,VMults,
1148 u1,u2,d1,d2,rational,dc);
1149 dim = rational ? 4 : 3;
1150
1151 if (!rational) {
1152 if ((Nu > UDegree) || (Nv > VDegree)) {
1153 Vn.SetX(0.);
1154 Vn.SetY(0.);
1155 Vn.SetZ(0.);
1156 return;
1157 }
1158 }
1159
1160 Standard_Integer n1 = ufirst ? Nu : Nv;
1161 Standard_Integer n2 = ufirst ? Nv : Nu;
1162
1163 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1164
1165 for (k = 0; k <= Min(n1,d1); k++)
1166 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1167
1168 Standard_Real *result;
1169 if (rational) {
1170 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1171 result = dc.ders; // because Standard_False ci-dessus.
1172
1173 }
1174 else {
1175 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1176 }
1177
1178 Vn.SetX(result[0]);
1179 Vn.SetY(result[1]);
1180 Vn.SetZ(result[2]);
1181}
1182
1183//
1184// Surface modifications
1185//
1186// a surface is processed as a curve of curves.
1187// i.e : a curve of parameter u where the current point is the set of poles
1188// of the iso.
1189//
1190
1191//=======================================================================
1192//function : Iso
1193//purpose :
1194//=======================================================================
1195
1196void BSplSLib::Iso(const Standard_Real Param,
1197 const Standard_Boolean IsU,
1198 const TColgp_Array2OfPnt& Poles,
1199 const TColStd_Array2OfReal& Weights,
1200 const TColStd_Array1OfReal& Knots,
1201 const TColStd_Array1OfInteger& Mults,
1202 const Standard_Integer Degree,
1203 const Standard_Boolean Periodic,
1204 TColgp_Array1OfPnt& CPoles,
1205 TColStd_Array1OfReal& CWeights)
1206{
1207 Standard_Integer index = 0;
1208 Standard_Real u = Param;
1209 Standard_Boolean rational = &Weights != NULL;
1210 Standard_Integer dim = rational ? 4 : 3;
1211
1212 // compute local knots
1213
1214 BSplSLib_LocalArray locknots1 (2*Degree);
1215 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1216 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1217 if (&Mults == NULL)
1218 index -= Knots.Lower() + Degree;
1219 else
1220 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1221
1222
1223 // copy the local poles
1224
1225// Standard_Integer f1,l1,f2,l2,i,j,k;
1226 Standard_Integer f1,l1,f2,l2,i,j;
1227
1228 if (IsU) {
1229 f1 = Poles.LowerRow();
1230 l1 = Poles.UpperRow();
1231 f2 = Poles.LowerCol();
1232 l2 = Poles.UpperCol();
1233 }
1234 else {
1235 f1 = Poles.LowerCol();
1236 l1 = Poles.UpperCol();
1237 f2 = Poles.LowerRow();
1238 l2 = Poles.UpperRow();
1239 }
1240
1241 BSplSLib_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim);
1242
1243 Standard_Real w, *pole = locpoles;
1244 index += f1;
1245
1246 for (i = 0; i <= Degree; i++) {
1247
1248 for (j = f2; j <= l2; j++) {
1249
1250 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1251 if (rational) {
1252 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1253 pole[0] = P.X() * w;
1254 pole[1] = P.Y() * w;
1255 pole[2] = P.Z() * w;
1256 }
1257 else {
1258 pole[0] = P.X();
1259 pole[1] = P.Y();
1260 pole[2] = P.Z();
1261 }
1262 pole += dim;
1263 }
1264 index++;
1265 if (index > l1) index = f1;
1266 }
1267
1268 // compute the iso
1269 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1270
1271 // get the result
1272 pole = locpoles;
1273
1274 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1275 gp_Pnt& P = CPoles(i);
1276 if (rational) {
1277 CWeights(i) = w = pole[3];
1278 P.SetX( pole[0] / w);
1279 P.SetY( pole[1] / w);
1280 P.SetZ( pole[2] / w);
1281 }
1282 else {
1283 P.SetX( pole[0]);
1284 P.SetY( pole[1]);
1285 P.SetZ( pole[2]);
1286 }
1287 pole += dim;
1288 }
1289
1290 // if the input is not rational but weights are wanted
1291 if (!rational && (&CWeights != NULL)) {
1292
1293 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1294 CWeights(i) = 1.;
1295 }
1296}
1297
1298//=======================================================================
1299//function : Reverse
1300//purpose :
1301//=======================================================================
1302
1303void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1304 const Standard_Integer Last,
1305 const Standard_Boolean UDirection)
1306{
1307 Standard_Integer i,j, l = Last;
1308 if ( UDirection) {
1309 l = Poles.LowerRow() +
1310 (l - Poles.LowerRow())%(Poles.ColLength());
1311 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1312 Poles.LowerCol(), Poles.UpperCol());
1313
1314 for (i = Poles.LowerRow(); i <= l; i++) {
1315
1316 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1317 temp(l-i,j) = Poles(i,j);
1318 }
1319 }
1320
1321 for (i = l+1; i <= Poles.UpperRow(); i++) {
1322
1323 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1324 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1325 }
1326 }
1327
1328 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1329
1330 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1331 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1332 }
1333 }
1334 }
1335 else {
1336 l = Poles.LowerCol() +
1337 (l - Poles.LowerCol())%(Poles.RowLength());
1338 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1339 0, Poles.RowLength()-1);
1340
1341 for (j = Poles.LowerCol(); j <= l; j++) {
1342
1343 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1344 temp(i,l-j) = Poles(i,j);
1345 }
1346 }
1347
1348 for (j = l+1; j <= Poles.UpperCol(); j++) {
1349
1350 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1351 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1352 }
1353 }
1354
1355 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1356
1357 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1358 Poles(i,j) = temp (i,j-Poles.LowerCol());
1359 }
1360 }
1361 }
1362}
1363
1364//=======================================================================
1365//function : Reverse
1366//purpose :
1367//=======================================================================
1368
1369void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1370 const Standard_Integer Last,
1371 const Standard_Boolean UDirection)
1372{
1373 Standard_Integer i,j, l = Last;
1374 if ( UDirection) {
1375 l = Weights.LowerRow() +
1376 (l - Weights.LowerRow())%(Weights.ColLength());
1377 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1378 Weights.LowerCol(), Weights.UpperCol());
1379
1380 for (i = Weights.LowerRow(); i <= l; i++) {
1381
1382 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1383 temp(l-i,j) = Weights(i,j);
1384 }
1385 }
1386
1387 for (i = l+1; i <= Weights.UpperRow(); i++) {
1388
1389 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1390 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1391 }
1392 }
1393
1394 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1395
1396 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1397 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1398 }
1399 }
1400 }
1401 else {
1402 l = Weights.LowerCol() +
1403 (l - Weights.LowerCol())%(Weights.RowLength());
1404 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1405 0, Weights.RowLength()-1);
1406
1407 for (j = Weights.LowerCol(); j <= l; j++) {
1408
1409 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1410 temp(i,l-j) = Weights(i,j);
1411 }
1412 }
1413
1414 for (j = l+1; j <= Weights.UpperCol(); j++) {
1415
1416 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1417 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1418 }
1419 }
1420
1421 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1422
1423 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1424 Weights(i,j) = temp (i,j-Weights.LowerCol());
1425 }
1426 }
1427 }
1428}
1429
1430//=======================================================================
1431//function : IsRational
1432//purpose :
1433//=======================================================================
1434
1435Standard_Boolean BSplSLib::IsRational
1436(const TColStd_Array2OfReal& Weights,
1437 const Standard_Integer I1,
1438 const Standard_Integer I2,
1439 const Standard_Integer J1,
1440 const Standard_Integer J2,
1441 const Standard_Real Epsi)
1442{
1443 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1444 Standard_Integer i,j;
1445 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1446 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1447
1448 for (i = I1 - fi; i < I2 - fi; i++) {
1449
1450 for (j = J1 - fj; j < J2 - fj; j++) {
1451 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1452 return Standard_True;
1453 }
1454 }
1455 return Standard_False;
1456}
1457
1458//=======================================================================
1459//function : SetPoles
1460//purpose :
1461//=======================================================================
1462
1463void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1464 TColStd_Array1OfReal& FP,
1465 const Standard_Boolean UDirection)
1466{
1467 Standard_Integer i, j, l = FP.Lower();
1468 Standard_Integer PLowerRow = Poles.LowerRow();
1469 Standard_Integer PUpperRow = Poles.UpperRow();
1470 Standard_Integer PLowerCol = Poles.LowerCol();
1471 Standard_Integer PUpperCol = Poles.UpperCol();
1472 if (UDirection) {
1473
1474 for ( i = PLowerRow; i <= PUpperRow; i++) {
1475
1476 for ( j = PLowerCol; j <= PUpperCol; j++) {
1477 const gp_Pnt& P = Poles.Value(i,j);
1478 FP(l) = P.X(); l++;
1479 FP(l) = P.Y(); l++;
1480 FP(l) = P.Z(); l++;
1481 }
1482 }
1483 }
1484 else {
1485
1486 for ( j = PLowerCol; j <= PUpperCol; j++) {
1487
1488 for ( i = PLowerRow; i <= PUpperRow; i++) {
1489 const gp_Pnt& P = Poles.Value(i,j);
1490 FP(l) = P.X(); l++;
1491 FP(l) = P.Y(); l++;
1492 FP(l) = P.Z(); l++;
1493 }
1494 }
1495 }
1496}
1497
1498//=======================================================================
1499//function : SetPoles
1500//purpose :
1501//=======================================================================
1502
1503void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1504 const TColStd_Array2OfReal& Weights,
1505 TColStd_Array1OfReal& FP,
1506 const Standard_Boolean UDirection)
1507{
1508 Standard_Integer i, j, l = FP.Lower();
1509 Standard_Integer PLowerRow = Poles.LowerRow();
1510 Standard_Integer PUpperRow = Poles.UpperRow();
1511 Standard_Integer PLowerCol = Poles.LowerCol();
1512 Standard_Integer PUpperCol = Poles.UpperCol();
1513 if (UDirection) {
1514
1515 for ( i = PLowerRow; i <= PUpperRow; i++) {
1516
1517 for ( j = PLowerCol; j <= PUpperCol; j++) {
1518 const gp_Pnt& P = Poles .Value(i,j);
1519 Standard_Real w = Weights.Value(i,j);
1520 FP(l) = P.X() * w; l++;
1521 FP(l) = P.Y() * w; l++;
1522 FP(l) = P.Z() * w; l++;
1523 FP(l) = w; l++;
1524 }
1525 }
1526 }
1527 else {
1528
1529 for ( j = PLowerCol; j <= PUpperCol; j++) {
1530
1531 for ( i = PLowerRow; i <= PUpperRow; i++) {
1532 const gp_Pnt& P = Poles .Value(i,j);
1533 Standard_Real w = Weights.Value(i,j);
1534 FP(l) = P.X() * w; l++;
1535 FP(l) = P.Y() * w; l++;
1536 FP(l) = P.Z() * w; l++;
1537 FP(l) = w; l++;
1538 }
1539 }
1540 }
1541}
1542
1543//=======================================================================
1544//function : GetPoles
1545//purpose :
1546//=======================================================================
1547
1548void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1549 TColgp_Array2OfPnt& Poles,
1550 const Standard_Boolean UDirection)
1551{
1552 Standard_Integer i, j, l = FP.Lower();
1553 Standard_Integer PLowerRow = Poles.LowerRow();
1554 Standard_Integer PUpperRow = Poles.UpperRow();
1555 Standard_Integer PLowerCol = Poles.LowerCol();
1556 Standard_Integer PUpperCol = Poles.UpperCol();
1557 if (UDirection) {
1558
1559 for ( i = PLowerRow; i <= PUpperRow; i++) {
1560
1561 for ( j = PLowerCol; j <= PUpperCol; j++) {
1562 gp_Pnt& P = Poles.ChangeValue(i,j);
1563 P.SetX(FP(l)); l++;
1564 P.SetY(FP(l)); l++;
1565 P.SetZ(FP(l)); l++;
1566 }
1567 }
1568 }
1569 else {
1570
1571 for ( j = PLowerCol; j <= PUpperCol; j++) {
1572
1573 for ( i = PLowerRow; i <= PUpperRow; i++) {
1574 gp_Pnt& P = Poles.ChangeValue(i,j);
1575 P.SetX(FP(l)); l++;
1576 P.SetY(FP(l)); l++;
1577 P.SetZ(FP(l)); l++;
1578 }
1579 }
1580 }
1581}
1582
1583//=======================================================================
1584//function : GetPoles
1585//purpose :
1586//=======================================================================
1587
1588void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1589 TColgp_Array2OfPnt& Poles,
1590 TColStd_Array2OfReal& Weights,
1591 const Standard_Boolean UDirection)
1592{
1593 Standard_Integer i, j, l = FP.Lower();
1594 Standard_Integer PLowerRow = Poles.LowerRow();
1595 Standard_Integer PUpperRow = Poles.UpperRow();
1596 Standard_Integer PLowerCol = Poles.LowerCol();
1597 Standard_Integer PUpperCol = Poles.UpperCol();
1598 if (UDirection) {
1599
1600 for ( i = PLowerRow; i <= PUpperRow; i++) {
1601
1602 for ( j = PLowerCol; j <= PUpperCol; j++) {
1603 Standard_Real w = FP( l + 3);
1604 Weights(i,j) = w;
1605 gp_Pnt& P = Poles.ChangeValue(i,j);
1606 P.SetX(FP(l) / w); l++;
1607 P.SetY(FP(l) / w); l++;
1608 P.SetZ(FP(l) / w); l++;
1609 l++;
1610 }
1611 }
1612 }
1613 else {
1614
1615 for ( j = PLowerCol; j <= PUpperCol; j++) {
1616
1617 for ( i = PLowerRow; i <= PUpperRow; i++) {
1618 Standard_Real w = FP( l + 3);
1619 Weights(i,j) = w;
1620 gp_Pnt& P = Poles.ChangeValue(i,j);
1621 P.SetX(FP(l) / w); l++;
1622 P.SetY(FP(l) / w); l++;
1623 P.SetZ(FP(l) / w); l++;
1624 l++;
1625 }
1626 }
1627 }
1628}
1629
1630//=======================================================================
1631//function : InsertKnots
1632//purpose :
1633//=======================================================================
1634
1635void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1636 const Standard_Integer Degree,
1637 const Standard_Boolean Periodic,
1638 const TColgp_Array2OfPnt& Poles,
1639 const TColStd_Array2OfReal& Weights,
1640 const TColStd_Array1OfReal& Knots,
1641 const TColStd_Array1OfInteger& Mults,
1642 const TColStd_Array1OfReal& AddKnots,
1643 const TColStd_Array1OfInteger& AddMults,
1644 TColgp_Array2OfPnt& NewPoles,
1645 TColStd_Array2OfReal& NewWeights,
1646 TColStd_Array1OfReal& NewKnots,
1647 TColStd_Array1OfInteger& NewMults,
1648 const Standard_Real Epsilon,
1649 const Standard_Boolean Add )
1650{
1651 Standard_Boolean rational = &Weights != NULL;
1652 Standard_Integer dim = 3;
1653 if (rational) dim++;
1654
1655 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1656 TColStd_Array1OfReal
1657 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1658
1659 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1660 else SetPoles(Poles,poles,UDirection);
1661
1662 if (UDirection) {
1663 dim *= Poles.RowLength();
1664 }
1665 else {
1666 dim *= Poles.ColLength();
1667 }
1668 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1669 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1670 Epsilon,Add);
1671
1672 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1673 else GetPoles(newpoles,NewPoles,UDirection);
1674}
1675
1676//=======================================================================
1677//function : RemoveKnot
1678//purpose :
1679//=======================================================================
1680
1681Standard_Boolean BSplSLib::RemoveKnot
1682(const Standard_Boolean UDirection,
1683 const Standard_Integer Index,
1684 const Standard_Integer Mult,
1685 const Standard_Integer Degree,
1686 const Standard_Boolean Periodic,
1687 const TColgp_Array2OfPnt& Poles,
1688 const TColStd_Array2OfReal& Weights,
1689 const TColStd_Array1OfReal& Knots,
1690 const TColStd_Array1OfInteger& Mults,
1691 TColgp_Array2OfPnt& NewPoles,
1692 TColStd_Array2OfReal& NewWeights,
1693 TColStd_Array1OfReal& NewKnots,
1694 TColStd_Array1OfInteger& NewMults,
1695 const Standard_Real Tolerance)
1696{
1697 Standard_Boolean rational = &Weights != NULL;
1698 Standard_Integer dim = 3;
1699 if (rational) dim++;
1700
1701 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1702 TColStd_Array1OfReal
1703 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1704
1705 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1706 else SetPoles(Poles,poles,UDirection);
1707
1708 if (UDirection) {
1709 dim *= Poles.RowLength();
1710 }
1711 else {
1712 dim *= Poles.ColLength();
1713 }
1714
1715 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1716 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1717 Tolerance))
1718 return Standard_False;
1719
1720 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1721 else GetPoles(newpoles,NewPoles,UDirection);
1722 return Standard_True;
1723}
1724
1725//=======================================================================
1726//function : IncreaseDegree
1727//purpose :
1728//=======================================================================
1729
1730void BSplSLib::IncreaseDegree
1731(const Standard_Boolean UDirection,
1732 const Standard_Integer Degree,
1733 const Standard_Integer NewDegree,
1734 const Standard_Boolean Periodic,
1735 const TColgp_Array2OfPnt& Poles,
1736 const TColStd_Array2OfReal& Weights,
1737 const TColStd_Array1OfReal& Knots,
1738 const TColStd_Array1OfInteger& Mults,
1739 TColgp_Array2OfPnt& NewPoles,
1740 TColStd_Array2OfReal& NewWeights,
1741 TColStd_Array1OfReal& NewKnots,
1742 TColStd_Array1OfInteger& NewMults)
1743{
1744 Standard_Boolean rational = &Weights != NULL;
1745 Standard_Integer dim = 3;
1746 if (rational) dim++;
1747
1748 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1749 TColStd_Array1OfReal
1750 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1751
1752 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1753 else SetPoles(Poles,poles,UDirection);
1754
1755 if (UDirection) {
1756 dim *= Poles.RowLength();
1757 }
1758 else {
1759 dim *= Poles.ColLength();
1760 }
1761
1762 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1763 newpoles,NewKnots,NewMults);
1764
1765 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1766 else GetPoles(newpoles,NewPoles,UDirection);
1767}
1768
1769//=======================================================================
1770//function : Unperiodize
1771//purpose :
1772//=======================================================================
1773
1774void BSplSLib::Unperiodize
1775(const Standard_Boolean UDirection,
1776 const Standard_Integer Degree,
1777 const TColStd_Array1OfInteger& Mults,
1778 const TColStd_Array1OfReal& Knots,
1779 const TColgp_Array2OfPnt& Poles,
1780 const TColStd_Array2OfReal& Weights,
1781 TColStd_Array1OfInteger& NewMults,
1782 TColStd_Array1OfReal& NewKnots,
1783 TColgp_Array2OfPnt& NewPoles,
1784 TColStd_Array2OfReal& NewWeights)
1785{
1786 Standard_Boolean rational = &Weights != NULL;
1787 Standard_Integer dim = 3;
1788 if (rational) dim++;
1789
1790 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1791 TColStd_Array1OfReal
1792 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1793
1794 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1795 else SetPoles(Poles,poles,UDirection);
1796
1797 if (UDirection) {
1798 dim *= Poles.RowLength();
1799 }
1800 else {
1801 dim *= Poles.ColLength();
1802 }
1803
1804 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1805 NewMults, NewKnots, newpoles);
1806
1807 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1808 else GetPoles(newpoles,NewPoles,UDirection);
1809}
1810
1811//=======================================================================
1812//function : BuildCache
1813//purpose : Stores theTaylor expansion normalized between 0,1 in the
1814// Cache : in case of a rational function the Poles are
1815// stored in homogeneous form
1816//=======================================================================
1817
1818void BSplSLib::BuildCache
1819(const Standard_Real U,
1820 const Standard_Real V,
1821 const Standard_Real USpanDomain,
1822 const Standard_Real VSpanDomain,
1823 const Standard_Boolean UPeriodic,
1824 const Standard_Boolean VPeriodic,
1825 const Standard_Integer UDegree,
1826 const Standard_Integer VDegree,
1827 const Standard_Integer UIndex,
1828 const Standard_Integer VIndex,
1829 const TColStd_Array1OfReal& UFlatKnots,
1830 const TColStd_Array1OfReal& VFlatKnots,
1831 const TColgp_Array2OfPnt& Poles,
1832 const TColStd_Array2OfReal& Weights,
1833 TColgp_Array2OfPnt& CachePoles,
1834 TColStd_Array2OfReal& CacheWeights)
1835{
1836 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1837 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1838 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1839 if (&Weights != NULL)
1840 rational_u = rational_v = Standard_True;
1841 else
1842 rational_u = rational_v = Standard_False;
1843 BSplSLib_DataContainer dc (UDegree, VDegree);
1844 flag_u_or_v =
1845 PrepareEval (U,
1846 V,
1847 UIndex,
1848 VIndex,
1849 UDegree,
1850 VDegree,
1851 rational_u,
1852 rational_v,
1853 UPeriodic,
1854 VPeriodic,
1855 Poles,
1856 Weights,
1857 UFlatKnots,
1858 VFlatKnots,
1859 (BSplCLib::NoMults()),
1860 (BSplCLib::NoMults()),
1861 u1,
1862 u2,
1863 d1,
1864 d2,
1865 rational,
1866 dc);
1867 d1p1 = d1 + 1;
1868 d2p1 = d2 + 1;
1869 if (rational) {
1870 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1871
1872 for (kk = 0; kk <= d1 ; kk++)
1873 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1874 if (flag_u_or_v) {
1875 min_degree_domain = USpanDomain ;
1876 max_degree_domain = VSpanDomain ;
1877 }
1878 else {
1879 min_degree_domain = VSpanDomain ;
1880 max_degree_domain = USpanDomain ;
1881 }
1882 factor[0] = 1.0e0 ;
1883
1884 for (ii = 0 ; ii <= d2 ; ii++) {
1885 iii = ii + 1;
1886 factor[1] = 1.0e0 ;
1887
1888 for (jj = 0 ; jj <= d1 ; jj++) {
1889 jjj = jj + 1;
1890 Index = jj * d2p1 + ii ;
1891 Index = Index << 2;
1892 gp_Pnt& P = CachePoles(iii,jjj);
1893 f = factor[0] * factor[1];
1894 P.SetX( f * dc.poles[Index]); Index++;
1895 P.SetY( f * dc.poles[Index]); Index++;
1896 P.SetZ( f * dc.poles[Index]); Index++;
1897 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1898 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1899 }
1900 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1901 }
1902 }
1903 else {
1904 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1905
1906 for (kk = 0; kk <= d1 ; kk++)
1907 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
1908 if (flag_u_or_v) {
1909 min_degree_domain = USpanDomain ;
1910 max_degree_domain = VSpanDomain ;
1911 }
1912 else {
1913 min_degree_domain = VSpanDomain ;
1914 max_degree_domain = USpanDomain ;
1915 }
1916 factor[0] = 1.0e0 ;
1917
1918 for (ii = 0 ; ii <= d2 ; ii++) {
1919 iii = ii + 1;
1920 factor[1] = 1.0e0 ;
1921
1922 for (jj = 0 ; jj <= d1 ; jj++) {
1923 jjj = jj + 1;
1924 Index = jj * d2p1 + ii ;
1925 Index = (Index << 1) + Index;
1926 gp_Pnt& P = CachePoles(iii,jjj);
1927 f = factor[0] * factor[1];
1928 P.SetX( f * dc.poles[Index]); Index++;
1929 P.SetY( f * dc.poles[Index]); Index++;
1930 P.SetZ( f * dc.poles[Index]);
1931 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1932 }
1933 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1934 }
1935 if (&Weights != NULL) {
1936 //
1937 // means that PrepareEval did found out that the surface was
1938 // locally polynomial but since the surface is constructed
1939 // with some weights we need to set the weight polynomial to constant
1940 //
1941
1942 for (ii = 1 ; ii <= d2p1 ; ii++) {
1943
1944 for (jj = 1 ; jj <= d1p1 ; jj++) {
1945 CacheWeights(ii,jj) = 0.0e0 ;
1946 }
1947 }
1948 CacheWeights(1,1) = 1.0e0 ;
1949 }
1950 }
1951}
1952
1953//=======================================================================
1954//function : CacheD0
1955//purpose : Evaluates the polynomial cache of the Bspline Curve
1956//
1957//=======================================================================
1958
1959void BSplSLib::CacheD0(const Standard_Real UParameter,
1960 const Standard_Real VParameter,
1961 const Standard_Integer UDegree,
1962 const Standard_Integer VDegree,
1963 const Standard_Real UCacheParameter,
1964 const Standard_Real VCacheParameter,
1965 const Standard_Real USpanLenght,
1966 const Standard_Real VSpanLenght,
1967 const TColgp_Array2OfPnt& PolesArray,
1968 const TColStd_Array2OfReal& WeightsArray,
1969 gp_Pnt& aPoint)
1970{
1971 //
1972 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
1973 // form
1974 // the SpanLenght is the normalizing factor so that the polynomial is between
1975 // 0 and 1
1976 Standard_Integer
1977// ii,
1978 dimension,
1979 min_degree,
1980 max_degree ;
1981
1982 Standard_Real
1983 new_parameter[2] ,
1984 inverse ;
1985
1986 Standard_Real *
1987 PArray = (Standard_Real *)
1988 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
1989 Standard_Real *
1990 myPoint = (Standard_Real *) &aPoint ;
1991 if (UDegree <= VDegree) {
1992 min_degree = UDegree ;
1993 max_degree = VDegree ;
1994 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
1995 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
1996 dimension = 3 * (UDegree + 1) ;
1997 }
1998 else {
1999 min_degree = VDegree ;
2000 max_degree = UDegree ;
2001 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
2002 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
2003 dimension = 3 * (VDegree + 1) ;
2004 }
2005 BSplSLib_LocalArray locpoles(dimension);
2006
2007 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2008 max_degree,
2009 dimension,
2010 max_degree*dimension,
2011 PArray[0],
2012 locpoles[0]) ;
2013
2014 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2015 min_degree,
2016 3,
2017 (min_degree << 1) + min_degree,
2018 locpoles[0],
2019 myPoint[0]) ;
2020 if (&WeightsArray != NULL) {
2021 dimension = min_degree + 1 ;
2022 Standard_Real *
2023 WArray = (Standard_Real *)
2024 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2025 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2026 max_degree,
2027 dimension,
2028 max_degree*dimension,
2029 WArray[0],
2030 locpoles[0]) ;
2031
2032 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2033 min_degree,
2034 1,
2035 min_degree,
2036 locpoles[0],
2037 inverse) ;
2038 inverse = 1.0e0/ inverse ;
2039
2040 myPoint[0] *= inverse ;
2041 myPoint[1] *= inverse ;
2042 myPoint[2] *= inverse ;
2043 }
2044}
2045
2046//=======================================================================
2047//function : CacheD1
2048//purpose : Evaluates the polynomial cache of the Bspline Curve
2049//
2050//=======================================================================
2051
2052void BSplSLib::CacheD1(const Standard_Real UParameter,
2053 const Standard_Real VParameter,
2054 const Standard_Integer UDegree,
2055 const Standard_Integer VDegree,
2056 const Standard_Real UCacheParameter,
2057 const Standard_Real VCacheParameter,
2058 const Standard_Real USpanLenght,
2059 const Standard_Real VSpanLenght,
2060 const TColgp_Array2OfPnt& PolesArray,
2061 const TColStd_Array2OfReal& WeightsArray,
2062 gp_Pnt& aPoint,
2063 gp_Vec& aVecU,
2064 gp_Vec& aVecV)
2065{
2066 //
2067 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2068 // form
2069 // the SpanLenght is the normalizing factor so that the polynomial is between
2070 // 0 and 1
2071 Standard_Integer
2072// ii,
2073// jj,
2074// kk,
2075 dimension,
2076 min_degree,
2077 max_degree ;
2078
2079 Standard_Real
2080 inverse_min,
2081 inverse_max,
2082 new_parameter[2] ;
2083
2084 Standard_Real *
2085 PArray = (Standard_Real *)
2086 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2087 Standard_Real local_poles_array[2][2][3],
2088 local_poles_and_weights_array[2][2][4],
2089 local_weights_array[2][2] ;
2090 Standard_Real * my_vec_min,
2091 * my_vec_max,
2092 * my_point ;
2093 my_point = (Standard_Real *) &aPoint ;
2094 //
2095 // initialize in case of rational evaluation
2096 // because RationalDerivative will use all
2097 // the coefficients
2098 //
2099 //
2100 if (&WeightsArray != NULL) {
2101
2102 local_poles_array [0][0][0] = 0.0e0 ;
2103 local_poles_array [0][0][1] = 0.0e0 ;
2104 local_poles_array [0][0][2] = 0.0e0 ;
2105 local_weights_array [0][0] = 0.0e0 ;
2106 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2107 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2108 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2109 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2110
2111 local_poles_array [0][1][0] = 0.0e0 ;
2112 local_poles_array [0][1][1] = 0.0e0 ;
2113 local_poles_array [0][1][2] = 0.0e0 ;
2114 local_weights_array [0][1] = 0.0e0 ;
2115 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2116 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2117 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2118 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2119
2120 local_poles_array [1][0][0] = 0.0e0 ;
2121 local_poles_array [1][0][1] = 0.0e0 ;
2122 local_poles_array [1][0][2] = 0.0e0 ;
2123 local_weights_array [1][0] = 0.0e0 ;
2124 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2125 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2126 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2127 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2128
2129 local_poles_array [1][1][0] = 0.0e0 ;
2130 local_poles_array [1][1][1] = 0.0e0 ;
2131 local_poles_array [1][1][2] = 0.0e0 ;
2132 local_weights_array [1][1] = 0.0e0 ;
2133 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2134 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2135 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2136 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2137 }
2138
2139 if (UDegree <= VDegree) {
2140 min_degree = UDegree ;
2141 max_degree = VDegree ;
2142 inverse_min = 1.0e0/USpanLenght ;
2143 inverse_max = 1.0e0/VSpanLenght ;
2144 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2145 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2146
2147 dimension = 3 * (UDegree + 1) ;
2148 my_vec_min = (Standard_Real *) &aVecU ;
2149 my_vec_max = (Standard_Real *) &aVecV ;
2150 }
2151 else {
2152 min_degree = VDegree ;
2153 max_degree = UDegree ;
2154 inverse_min = 1.0e0/VSpanLenght ;
2155 inverse_max = 1.0e0/USpanLenght ;
2156 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2157 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2158 dimension = 3 * (VDegree + 1) ;
2159 my_vec_min = (Standard_Real *) &aVecV ;
2160 my_vec_max = (Standard_Real *) &aVecU ;
2161 }
2162
2163 BSplSLib_LocalArray locpoles (2 * dimension);
2164
2165 PLib::EvalPolynomial(new_parameter[0],
2166 1,
2167 max_degree,
2168 dimension,
2169 PArray[0],
2170 locpoles[0]) ;
2171
2172 PLib::EvalPolynomial(new_parameter[1],
2173 1,
2174 min_degree,
2175 3,
2176 locpoles[0],
2177 local_poles_array[0][0][0]) ;
2178 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2179 min_degree,
2180 3,
2181 (min_degree << 1) + min_degree,
2182 locpoles[dimension],
2183 local_poles_array[1][0][0]) ;
2184
2185 if (&WeightsArray != NULL) {
2186 dimension = min_degree + 1 ;
2187 Standard_Real *
2188 WArray = (Standard_Real *)
2189 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2190 PLib::EvalPolynomial(new_parameter[0],
2191 1,
2192 max_degree,
2193 dimension,
2194 WArray[0],
2195 locpoles[0]) ;
2196
2197 PLib::EvalPolynomial(new_parameter[1],
2198 1,
2199 min_degree,
2200 1,
2201 locpoles[0],
2202 local_weights_array[0][0]) ;
2203 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2204 min_degree,
2205 1,
2206 min_degree,
2207 locpoles[dimension],
2208 local_weights_array[1][0]) ;
2209
2210 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2211 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2212 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2213 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2214
2215 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2216 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2217 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2218 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2219
2220 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2221 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2222 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2223 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2224
2225 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2226 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2227 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2228 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2229
2230 BSplSLib::RationalDerivative(1,
2231 1,
2232 1,
2233 1,
2234 local_poles_and_weights_array[0][0][0],
2235 local_poles_array[0][0][0]) ;
2236 }
2237
2238 my_point [0] = local_poles_array [0][0][0] ;
2239 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2240 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2241
2242 my_point [1] = local_poles_array [0][0][1] ;
2243 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2244 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2245
2246 my_point [2] = local_poles_array [0][0][2] ;
2247 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2248 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2249}
2250
2251//=======================================================================
2252//function : CacheD2
2253//purpose : Evaluates the polynomial cache of the Bspline Curve
2254//
2255//=======================================================================
2256
2257void BSplSLib::CacheD2(const Standard_Real UParameter,
2258 const Standard_Real VParameter,
2259 const Standard_Integer UDegree,
2260 const Standard_Integer VDegree,
2261 const Standard_Real UCacheParameter,
2262 const Standard_Real VCacheParameter,
2263 const Standard_Real USpanLenght,
2264 const Standard_Real VSpanLenght,
2265 const TColgp_Array2OfPnt& PolesArray,
2266 const TColStd_Array2OfReal& WeightsArray,
2267 gp_Pnt& aPoint,
2268 gp_Vec& aVecU,
2269 gp_Vec& aVecV,
2270 gp_Vec& aVecUU,
2271 gp_Vec& aVecUV,
2272 gp_Vec& aVecVV)
2273{
2274 //
2275 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2276 // form
2277 // the SpanLenght is the normalizing factor so that the polynomial is between
2278 // 0 and 1
2279 Standard_Integer
2280 ii,
2281// jj,
2282 kk,
2283 index,
2284 dimension,
2285 min_degree,
2286 max_degree ;
2287
2288 Standard_Real
2289 inverse_min,
2290 inverse_max,
2291 new_parameter[2] ;
2292
2293 Standard_Real *
2294 PArray = (Standard_Real *)
2295 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2296 Standard_Real local_poles_array[3][3][3],
2297 local_poles_and_weights_array[3][3][4],
2298 local_weights_array[3][3] ;
2299 Standard_Real * my_vec_min,
2300 * my_vec_max,
2301 * my_vec_min_min,
2302 * my_vec_max_max,
2303 * my_vec_min_max,
2304 * my_point ;
2305 my_point = (Standard_Real *) &aPoint ;
2306
2307 //
2308 // initialize in case the min and max degree are less than 2
2309 //
2310 local_poles_array[0][0][0] = 0.0e0 ;
2311 local_poles_array[0][0][1] = 0.0e0 ;
2312 local_poles_array[0][0][2] = 0.0e0 ;
2313 local_poles_array[0][1][0] = 0.0e0 ;
2314 local_poles_array[0][1][1] = 0.0e0 ;
2315 local_poles_array[0][1][2] = 0.0e0 ;
2316 local_poles_array[0][2][0] = 0.0e0 ;
2317 local_poles_array[0][2][1] = 0.0e0 ;
2318 local_poles_array[0][2][2] = 0.0e0 ;
2319
2320 local_poles_array[1][0][0] = 0.0e0 ;
2321 local_poles_array[1][0][1] = 0.0e0 ;
2322 local_poles_array[1][0][2] = 0.0e0 ;
2323 local_poles_array[1][1][0] = 0.0e0 ;
2324 local_poles_array[1][1][1] = 0.0e0 ;
2325 local_poles_array[1][1][2] = 0.0e0 ;
2326 local_poles_array[1][2][0] = 0.0e0 ;
2327 local_poles_array[1][2][1] = 0.0e0 ;
2328 local_poles_array[1][2][2] = 0.0e0 ;
2329
2330 local_poles_array[2][0][0] = 0.0e0 ;
2331 local_poles_array[2][0][1] = 0.0e0 ;
2332 local_poles_array[2][0][2] = 0.0e0 ;
2333 local_poles_array[2][1][0] = 0.0e0 ;
2334 local_poles_array[2][1][1] = 0.0e0 ;
2335 local_poles_array[2][1][2] = 0.0e0 ;
2336 local_poles_array[2][2][0] = 0.0e0 ;
2337 local_poles_array[2][2][1] = 0.0e0 ;
2338 local_poles_array[2][2][2] = 0.0e0 ;
2339 //
2340 // initialize in case of rational evaluation
2341 // because RationalDerivative will use all
2342 // the coefficients
2343 //
2344 //
2345 if (&WeightsArray != NULL) {
2346
2347 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2348 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2349 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2350 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2351 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2352 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2353 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2354 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2355 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2356
2357 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2358 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2359 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2360 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2361 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2362 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2363 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2364 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2365 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2366
2367 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2368 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2369 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2370 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2371 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2372 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2373 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2374 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2375 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2376
2377 local_poles_and_weights_array[0][0][3] =
2378 local_weights_array[0][0] = 0.0e0 ;
2379 local_poles_and_weights_array[0][1][3] =
2380 local_weights_array[0][1] = 0.0e0 ;
2381 local_poles_and_weights_array[0][2][3] =
2382 local_weights_array[0][2] = 0.0e0 ;
2383 local_poles_and_weights_array[1][0][3] =
2384 local_weights_array[1][0] = 0.0e0 ;
2385 local_poles_and_weights_array[1][1][3] =
2386 local_weights_array[1][1] = 0.0e0 ;
2387 local_poles_and_weights_array[1][2][3] =
2388 local_weights_array[1][2] = 0.0e0 ;
2389 local_poles_and_weights_array[2][0][3] =
2390 local_weights_array[2][0] = 0.0e0 ;
2391 local_poles_and_weights_array[2][1][3] =
2392 local_weights_array[2][1] = 0.0e0 ;
2393 local_poles_and_weights_array[2][2][3] =
2394 local_weights_array[2][2] = 0.0e0 ;
2395 }
2396
2397 if (UDegree <= VDegree) {
2398 min_degree = UDegree ;
2399 max_degree = VDegree ;
2400 inverse_min = 1.0e0/USpanLenght ;
2401 inverse_max = 1.0e0/VSpanLenght ;
2402 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2403 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2404
2405 dimension = 3 * (UDegree + 1) ;
2406 my_vec_min = (Standard_Real *) &aVecU ;
2407 my_vec_max = (Standard_Real *) &aVecV ;
2408 my_vec_min_min = (Standard_Real *) &aVecUU ;
2409 my_vec_min_max = (Standard_Real *) &aVecUV ;
2410 my_vec_max_max = (Standard_Real *) &aVecVV ;
2411 }
2412 else {
2413 min_degree = VDegree ;
2414 max_degree = UDegree ;
2415 inverse_min = 1.0e0/VSpanLenght ;
2416 inverse_max = 1.0e0/USpanLenght ;
2417 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2418 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2419 dimension = 3 * (VDegree + 1) ;
2420 my_vec_min = (Standard_Real *) &aVecV ;
2421 my_vec_max = (Standard_Real *) &aVecU ;
2422 my_vec_min_min = (Standard_Real *) &aVecVV ;
2423 my_vec_min_max = (Standard_Real *) &aVecUV ;
2424 my_vec_max_max = (Standard_Real *) &aVecUU ;
2425 }
2426
2427 BSplSLib_LocalArray locpoles (3 * dimension);
2428
2429 //
2430 // initialize in case min or max degree are less than 2
2431 //
2432 Standard_Integer MinIndMax = 2;
2433 if ( max_degree < 2) MinIndMax = max_degree;
2434 Standard_Integer MinIndMin = 2;
2435 if ( min_degree < 2) MinIndMin = min_degree;
2436
2437 index = MinIndMax * dimension ;
2438
2439 for (ii = MinIndMax ; ii < 3 ; ii++) {
2440
2441 for (kk = 0 ; kk < dimension ; kk++) {
2442 locpoles[index] = 0.0e0 ;
2443 index += 1 ;
2444 }
2445 }
2446
2447 PLib::EvalPolynomial(new_parameter[0],
2448 MinIndMax,
2449 max_degree,
2450 dimension,
2451 PArray[0],
2452 locpoles[0]) ;
2453
2454 PLib::EvalPolynomial(new_parameter[1],
2455 MinIndMin,
2456 min_degree,
2457 3,
2458 locpoles[0],
2459 local_poles_array[0][0][0]) ;
2460 PLib::EvalPolynomial(new_parameter[1],
2461 1,
2462 min_degree,
2463 3,
2464 locpoles[dimension],
2465 local_poles_array[1][0][0]) ;
2466
2467 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2468 min_degree,
2469 3,
2470 (min_degree << 1) + min_degree,
2471 locpoles[dimension + dimension],
2472 local_poles_array[2][0][0]) ;
2473
2474 if (&WeightsArray != NULL) {
2475 dimension = min_degree + 1 ;
2476 Standard_Real *
2477 WArray = (Standard_Real *)
2478 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2479 PLib::EvalPolynomial(new_parameter[0],
2480 MinIndMax,
2481 max_degree,
2482 dimension,
2483 WArray[0],
2484 locpoles[0]) ;
2485
2486 PLib::EvalPolynomial(new_parameter[1],
2487 MinIndMin,
2488 min_degree,
2489 1,
2490 locpoles[0],
2491 local_weights_array[0][0]) ;
2492 PLib::EvalPolynomial(new_parameter[1],
2493 1,
2494 min_degree,
2495 1,
2496 locpoles[dimension],
2497 local_weights_array[1][0]) ;
2498 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2499 min_degree,
2500 1,
2501 min_degree,
2502 locpoles[dimension + dimension],
2503 local_weights_array[2][0]) ;
2504
2505
2506 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2507 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2508 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2509 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2510 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2511 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2512 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2513 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2514 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2515
2516 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2517 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2518 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2519 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2520 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2521 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2522 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2523 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2524 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2525
2526 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2527 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2528 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2529 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2530 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2531 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2532 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2533 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2534 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2535
2536
2537 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2538 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2539 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2540 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2541 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2542 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2543 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2544 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2545 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2546
2547 BSplSLib::RationalDerivative(2,
2548 2,
2549 2,
2550 2,
2551 local_poles_and_weights_array[0][0][0],
2552 local_poles_array[0][0][0]) ;
2553 }
2554
2555
2556 Standard_Real minmin = inverse_min * inverse_min;
2557 Standard_Real minmax = inverse_min * inverse_max;
2558 Standard_Real maxmax = inverse_max * inverse_max;
2559
2560 my_point [0] = local_poles_array [0][0][0] ;
2561 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2562 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2563 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2564 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2565 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2566
2567 my_point [1] = local_poles_array [0][0][1] ;
2568 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2569 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2570 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2571 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2572 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2573
2574 my_point [2] = local_poles_array [0][0][2] ;
2575 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2576 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2577 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2578 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2579 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2580}
2581
2582//=======================================================================
2583//function : MovePoint
2584//purpose : Find the new poles which allows an old point (with a
2585// given u and v as parameters) to reach a new position
2586//=======================================================================
2587
2588void BSplSLib::MovePoint (const Standard_Real U,
2589 const Standard_Real V,
2590 const gp_Vec& Displ,
2591 const Standard_Integer UIndex1,
2592 const Standard_Integer UIndex2,
2593 const Standard_Integer VIndex1,
2594 const Standard_Integer VIndex2,
2595 const Standard_Integer UDegree,
2596 const Standard_Integer VDegree,
2597 const Standard_Boolean Rational,
2598 const TColgp_Array2OfPnt& Poles,
2599 const TColStd_Array2OfReal& Weights,
2600 const TColStd_Array1OfReal& UFlatKnots,
2601 const TColStd_Array1OfReal& VFlatKnots,
2602 Standard_Integer& UFirstIndex,
2603 Standard_Integer& ULastIndex,
2604 Standard_Integer& VFirstIndex,
2605 Standard_Integer& VLastIndex,
2606 TColgp_Array2OfPnt& NewPoles)
2607{
2608 // calculate the UBSplineBasis in the parameter U
2609 Standard_Integer UFirstNonZeroBsplineIndex;
2610 math_Matrix UBSplineBasis(1, 1,
2611 1, UDegree+1);
2612 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2613 0,
2614 UDegree+1,
2615 UFlatKnots,
2616 U,
2617 UFirstNonZeroBsplineIndex,
2618 UBSplineBasis);
2619 // calculate the VBSplineBasis in the parameter V
2620 Standard_Integer VFirstNonZeroBsplineIndex;
2621 math_Matrix VBSplineBasis(1, 1,
2622 1, VDegree+1);
2623 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2624 0,
2625 VDegree+1,
2626 VFlatKnots,
2627 V,
2628 VFirstNonZeroBsplineIndex,
2629 VBSplineBasis);
2630 if (ErrorCod1 || ErrorCod2) {
2631 UFirstIndex = 0;
2632 ULastIndex = 0;
2633 VFirstIndex = 0;
2634 VLastIndex = 0;
2635 return;
2636 }
2637
2638 // find the span which is predominant for parameter U
2639 UFirstIndex = UFirstNonZeroBsplineIndex;
2640 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2641 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2642 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2643
2644 Standard_Real maxValue = 0.0;
2645 Standard_Integer i, ukk1=0, ukk2;
2646
2647 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2648 if (UBSplineBasis(1,i) > maxValue) {
2649 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2650 maxValue = UBSplineBasis(1,i);
2651 }
2652 }
2653
2654 // find a ukk2 if symetriy
2655 ukk2 = ukk1;
2656 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2657 if ((ukk1+1) <= ULastIndex) {
2658 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2659 ukk2 = ukk1+1;
2660 }
2661 }
2662
2663 // find the span which is predominant for parameter V
2664 VFirstIndex = VFirstNonZeroBsplineIndex;
2665 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2666
2667 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2668 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2669
2670 maxValue = 0.0;
2671 Standard_Integer j, vkk1=0, vkk2;
2672
2673 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2674 if (VBSplineBasis(1,j) > maxValue) {
2675 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2676 maxValue = VBSplineBasis(1,j);
2677 }
2678 }
2679
2680 // find a vkk2 if symetriy
2681 vkk2 = vkk1;
2682 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2683 if ((vkk1+1) <= VLastIndex) {
2684 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2685 vkk2 = vkk1+1;
2686 }
2687 }
2688
2689 // compute the vector of displacement
2690 Standard_Real D1 = 0.0;
2691 Standard_Real D2 = 0.0;
2692 Standard_Real hN, Coef, DvalU, DvalV;
2693
2694 Standard_Integer ii, jj;
2695
2696 for (i = 1; i <= UDegree+1; i++) {
2697 ii = i + UFirstNonZeroBsplineIndex - 1;
2698 if (ii < ukk1) {
2699 DvalU = ukk1-ii;
2700 }
2701 else if (ii > ukk2) {
2702 DvalU = ii - ukk2;
2703 }
2704 else {
2705 DvalU = 0.0;
2706 }
2707
2708 for (j = 1; j <= VDegree+1; j++) {
2709 jj = j + VFirstNonZeroBsplineIndex - 1;
2710 if (Rational) {
2711 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2712 D2 += hN;
2713 }
2714 else {
2715 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2716 }
2717 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2718 if (jj < vkk1) {
2719 DvalV = vkk1-jj;
2720 }
2721 else if (jj > vkk2) {
2722 DvalV = jj - vkk2;
2723 }
2724 else {
2725 DvalV = 0.0;
2726 }
2727 D1 += 1./(DvalU + DvalV + 1.) * hN;
2728 }
2729 }
2730 }
2731
2732 if (Rational) {
2733 Coef = D2/D1;
2734 }
2735 else {
2736 Coef = 1./D1;
2737 }
2738
2739 // compute the new poles
2740
2741 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2742 if (i < ukk1) {
2743 DvalU = ukk1-i;
2744 }
2745 else if (i > ukk2) {
2746 DvalU = i - ukk2;
2747 }
2748 else {
2749 DvalU = 0.0;
2750 }
2751
2752 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2753 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2754 if (j < vkk1) {
2755 DvalV = vkk1-j;
2756 }
2757 else if (j > vkk2) {
2758 DvalV = j - vkk2;
2759 }
2760 else {
2761 DvalV = 0.0;
2762 }
2763 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2764 }
2765 else {
2766 NewPoles(i,j) = Poles(i,j);
2767 }
2768 }
2769 }
2770}
2771
2772//=======================================================================
0d969553
Y		 2773// function : Resolution
2774// purpose : this computes an estimate for the maximum of the
7fd59977 2775// partial derivatives both in U and in V
2776//
2777//
0d969553
Y		 2778// The calculation resembles at the calculation of curves with
2779// additional index for the control point. Let Si,j be the
2780// control points for ls surface and Di,j the weights.
2781// The checking of upper bounds for the partial derivatives
2782// will be omitted and Su is the next upper bound in the polynomial case :
7fd59977 2783//
2784//
2785//
2786// | Si,j - Si-1,j |
2787// d * Max | ------------- |
2788// i = 2,n | ti+d - ti |
2789// i=1.m
2790//
2791//
0d969553 2792// and in the rational case :
7fd59977 2793//
2794//
2795//
2796// Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2797// Max Max d * -----------------------------------------------
2798// k=1,n i dans Rj ti+d - ti
2799// j=1,m
2800// ----------------------------------------------------------------------
2801//
2802// Min Di,j
2803// i=1,n
2804// j=1,m
2805//
2806//
2807//
0d969553 2808// with Rj = {j-d, ...., j+d+d+1}.
7fd59977 2809//
2810//
2811//=======================================================================
2812
2813void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
2814 const TColStd_Array2OfReal& Weights,
2815 const TColStd_Array1OfReal& UKnots,
2816 const TColStd_Array1OfReal& VKnots,
2817 const TColStd_Array1OfInteger& UMults,
2818 const TColStd_Array1OfInteger& VMults,
2819 const Standard_Integer UDegree,
2820 const Standard_Integer VDegree,
2821 const Standard_Boolean URational,
2822 const Standard_Boolean VRational,
2823 const Standard_Boolean UPeriodic,
2824 const Standard_Boolean VPeriodic,
2825 const Standard_Real Tolerance3D,
2826 Standard_Real& UTolerance,
2827 Standard_Real& VTolerance)
2828{
2829 Standard_Real Wij,Wmj,Wji,Wjm;
2830 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2831 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2832 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2833 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2834
2835 max_derivative[0] = max_derivative[1] = 0.0e0 ;
2836
2837 Standard_Integer PRowLength, PColLength;
2838 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2839 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2840 Standard_Integer num_poles[2],num_flat_knots[2];
2841
2842 num_flat_knots[0] =
2843 BSplCLib::KnotSequenceLength(UMults,
2844 UDegree,
2845 UPeriodic) ;
2846 num_flat_knots[1] =
2847 BSplCLib::KnotSequenceLength(VMults,
2848 VDegree,
2849 VPeriodic) ;
2850 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
2851 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
2852 BSplCLib::KnotSequence(UKnots,
2853 UMults,
2854 UDegree,
2855 UPeriodic,
2856 flat_knots_in_u) ;
2857 BSplCLib::KnotSequence(VKnots,
2858 VMults,
2859 VDegree,
2860 VPeriodic,
2861 flat_knots_in_v) ;
2862 PRowLength = Poles.RowLength();
2863 PColLength = Poles.ColLength();
2864 if (URational || VRational) {
2865 Standard_Integer Wsize = PRowLength * PColLength;
2866 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2867 min_weights = WG[0];
2868
2869 for (ii = 1 ; ii < Wsize ; ii++) {
2870 min = WG[ii];
2871 if (min_weights > min) min_weights = min;
2872 }
2873 }
2874 Standard_Integer UD1 = UDegree + 1;
2875 Standard_Integer VD1 = VDegree + 1;
2876 num_poles[0] = num_flat_knots[0] - UD1;
2877 num_poles[1] = num_flat_knots[1] - VD1;
2878 poles_length[0] = PColLength;
2879 poles_length[1] = PRowLength;
2880 if(URational) {
2881 Standard_Integer UD2 = UDegree << 1;
2882 Standard_Integer VD2 = VDegree << 1;
2883
2884 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2885 ii_index = (ii - 1) % poles_length[0] + 1 ;
2886 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2887 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2888 inverse = 1.0e0 / inverse ;
2889 lower[0] = ii - UD1;
2890 if (lower[0] < 1) lower[0] = 1;
2891 upper[0] = ii + UD2 + 1;
2892 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2893
2894 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2895 jj_index = (jj - 1) % poles_length[1] + 1 ;
2896 lower[1] = jj - VD1;
2897 if (lower[1] < 1) lower[1] = 1;
2898 upper[1] = jj + VD2 + 1;
2899 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
2900 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
2901 Wij = Weights.Value(ii_index,jj_index);
2902 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
2903 Wmj = Weights.Value(ii_minus,jj_index);
2904 Xij = Pij.X();
2905 Yij = Pij.Y();
2906 Zij = Pij.Z();
2907 Xmj = Pmj.X();
2908 Ymj = Pmj.Y();
2909 Zmj = Pmj.Z();
2910
2911 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
2912 pp_index = (pp - 1) % poles_length[0] + 1 ;
2913
2914 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
2915 value = 0.0e0 ;
2916 qq_index = (qq - 1) % poles_length[1] + 1 ;
2917 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
2918 Xpq = Ppq.X();
2919 Ypq = Ppq.Y();
2920 Zpq = Ppq.Z();
2921 factor = (Xpq - Xij) * Wij;
2922 factor -= (Xpq - Xmj) * Wmj;
2923 if (factor < 0) factor = - factor;
2924 value += factor ;
2925 factor = (Ypq - Yij) * Wij;
2926 factor -= (Ypq - Ymj) * Wmj;
2927 if (factor < 0) factor = - factor;
2928 value += factor ;
2929 factor = (Zpq - Zij) * Wij;
2930 factor -= (Zpq - Zmj) * Wmj;
2931 if (factor < 0) factor = - factor;
2932 value += factor ;
2933 value *= inverse ;
2934 if (max_derivative[0] < value) max_derivative[0] = value ;
2935 }
2936 }
2937 }
2938 }
2939 max_derivative[0] /= min_weights ;
2940 }
2941 else {
2942
2943 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2944 ii_index = (ii - 1) % poles_length[0] + 1 ;
2945 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2946 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2947 inverse = 1.0e0 / inverse ;
2948
2949 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2950 jj_index = (jj - 1) % poles_length[1] + 1 ;
2951 value = 0.0e0 ;
2952 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
2953 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
2954 factor = Pij.X() - Pmj.X();
2955 if (factor < 0) factor = - factor;
2956 value += factor;
2957 factor = Pij.Y() - Pmj.Y();
2958 if (factor < 0) factor = - factor;
2959 value += factor;
2960 factor = Pij.Z() - Pmj.Z();
2961 if (factor < 0) factor = - factor;
2962 value += factor;
2963 value *= inverse ;
2964 if (max_derivative[0] < value) max_derivative[0] = value ;
2965 }
2966 }
2967 }
2968 max_derivative[0] *= UDegree ;
2969 if(VRational) {
2970 Standard_Integer UD2 = UDegree << 1;
2971 Standard_Integer VD2 = VDegree << 1;
2972
2973 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
2974 ii_index = (ii - 1) % poles_length[1] + 1 ;
2975 ii_minus = (ii - 2) % poles_length[1] + 1 ;
2976 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
2977 inverse = 1.0e0 / inverse ;
2978 lower[0] = ii - VD1;
2979 if (lower[0] < 1) lower[0] = 1;
2980 upper[0] = ii + VD2 + 1;
2981 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
2982
2983 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
2984 jj_index = (jj - 1) % poles_length[0] + 1 ;
2985 lower[1] = jj - UD1;
2986 if (lower[1] < 1) lower[1] = 1;
2987 upper[1] = jj + UD2 + 1;
2988 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
2989 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
2990 Wji = Weights.Value(jj_index,ii_index);
2991 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
2992 Wjm = Weights.Value(jj_index,ii_minus);
2993 Xji = Pji.X();
2994 Yji = Pji.Y();
2995 Zji = Pji.Z();
2996 Xjm = Pjm.X();
2997 Yjm = Pjm.Y();
2998 Zjm = Pjm.Z();
2999
3000 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
3001 pp_index = (pp - 1) % poles_length[1] + 1 ;
3002
3003 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
3004 value = 0.0e0 ;
3005 qq_index = (qq - 1) % poles_length[0] + 1 ;
3006 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
3007 Xqp = Pqp.X();
3008 Yqp = Pqp.Y();
3009 Zqp = Pqp.Z();
3010 factor = (Xqp - Xji) * Wji;
3011 factor -= (Xqp - Xjm) * Wjm;
3012 if (factor < 0) factor = - factor;
3013 value += factor ;
3014 factor = (Yqp - Yji) * Wji;
3015 factor -= (Yqp - Yjm) * Wjm;
3016 if (factor < 0) factor = - factor;
3017 value += factor ;
3018 factor = (Zqp - Zji) * Wji;
3019 factor -= (Zqp - Zjm) * Wjm;
3020 if (factor < 0) factor = - factor;
3021 value += factor ;
3022 value *= inverse ;
3023 if (max_derivative[1] < value) max_derivative[1] = value ;
3024 }
3025 }
3026 }
3027 }
3028 max_derivative[1] /= min_weights ;
3029 }
3030 else {
3031
3032 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3033 ii_index = (ii - 1) % poles_length[1] + 1 ;
3034 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3035 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3036 inverse = 1.0e0 / inverse ;
3037
3038 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3039 jj_index = (jj - 1) % poles_length[0] + 1 ;
3040 value = 0.0e0 ;
3041 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3042 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3043 factor = Pji.X() - Pjm.X() ;
3044 if (factor < 0) factor = - factor;
3045 value += factor;
3046 factor = Pji.Y() - Pjm.Y() ;
3047 if (factor < 0) factor = - factor;
3048 value += factor;
3049 factor = Pji.Z() - Pjm.Z() ;
3050 if (factor < 0) factor = - factor;
3051 value += factor;
3052 value *= inverse ;
3053 if (max_derivative[1] < value) max_derivative[1] = value ;
3054 }
3055 }
3056 }
3057 max_derivative[1] *= VDegree ;
3058 max_derivative[0] *= M_SQRT2 ;
3059 max_derivative[1] *= M_SQRT2 ;
3060 if(max_derivative[0] && max_derivative[1]) {
3061 UTolerance = Tolerance3D / max_derivative[0] ;
3062 VTolerance = Tolerance3D / max_derivative[1] ;
3063 }
3064 else {
3065 UTolerance=VTolerance=0.0;
3066#ifdef DEB
3067 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3068#endif
3069 }
3070}
3071
3072//=======================================================================
3073//function : Interpolate
3074//purpose :
3075//=======================================================================
3076
3077void BSplSLib::Interpolate(const Standard_Integer UDegree,
3078 const Standard_Integer VDegree,
3079 const TColStd_Array1OfReal& UFlatKnots,
3080 const TColStd_Array1OfReal& VFlatKnots,
3081 const TColStd_Array1OfReal& UParameters,
3082 const TColStd_Array1OfReal& VParameters,
3083 TColgp_Array2OfPnt& Poles,
3084 TColStd_Array2OfReal& Weights,
3085 Standard_Integer& InversionProblem)
3086{
3087 Standard_Integer ii, jj, ll, kk, dimension;
3088 Standard_Integer ULength = UParameters.Length();
3089 Standard_Integer VLength = VParameters.Length();
3090 Standard_Real * poles_array;
3091
0d969553 3092 // extraction of iso u
7fd59977 3093 dimension = 4*ULength;
3094 TColStd_Array2OfReal Points(1, VLength,
3095 1, dimension);
3096
3097 Handle(TColStd_HArray1OfInteger) ContactOrder =
3098 new (TColStd_HArray1OfInteger)(1, VLength);
3099 ContactOrder->Init(0);
3100
3101 for (ii=1; ii <= VLength; ii++) {
3102
3103 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3104 Points(ii,ll) = Poles(jj, ii).X();
3105 Points(ii,ll+1) = Poles(jj, ii).Y();
3106 Points(ii,ll+2) = Poles(jj, ii).Z();
3107 Points(ii,ll+3) = Weights(jj,ii) ;
3108 }
3109 }
3110
0d969553 3111 // interpolation of iso u
7fd59977 3112 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3113 BSplCLib::Interpolate(VDegree,
3114 VFlatKnots,
3115 VParameters,
3116 ContactOrder->Array1(),
3117 dimension,
3118 poles_array[0],
3119 InversionProblem) ;
3120 if (InversionProblem != 0) return;
3121
0d969553 3122 // extraction of iso v
7fd59977 3123
3124 dimension = VLength*4;
3125 TColStd_Array2OfReal IsoPoles(1, ULength,
3126 1, dimension);
3127
3128 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3129 ContactOrder->Init(0);
3130 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3131
3132 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3133
3134 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3135 IsoPoles (ii,ll) = Points(jj, kk);
3136 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3137 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3138 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3139 }
3140 }
0d969553 3141 // interpolation of iso v
7fd59977 3142 BSplCLib::Interpolate(UDegree,
3143 UFlatKnots,
3144 UParameters,
3145 ContactOrder->Array1(),
3146 dimension,
3147 poles_array[0],
3148 InversionProblem);
3149
0d969553 3150 // return results
7fd59977 3151
3152 for (ii=1; ii <= ULength; ii++) {
3153
3154 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3155 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3156 Poles.SetValue(ii, jj, Pnt);
3157 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3158 }
3159 }
3160}
3161
3162//=======================================================================
3163//function : Interpolate
3164//purpose :
3165//=======================================================================
3166
3167void BSplSLib::Interpolate(const Standard_Integer UDegree,
3168 const Standard_Integer VDegree,
3169 const TColStd_Array1OfReal& UFlatKnots,
3170 const TColStd_Array1OfReal& VFlatKnots,
3171 const TColStd_Array1OfReal& UParameters,
3172 const TColStd_Array1OfReal& VParameters,
3173 TColgp_Array2OfPnt& Poles,
3174 Standard_Integer& InversionProblem)
3175{
3176 Standard_Integer ii, jj, ll, kk, dimension;
3177 Standard_Integer ULength = UParameters.Length();
3178 Standard_Integer VLength = VParameters.Length();
3179 Standard_Real * poles_array;
3180
0d969553 3181 // extraction of iso u
7fd59977 3182 dimension = 3*ULength;
3183 TColStd_Array2OfReal Points(1, VLength,
3184 1, dimension);
3185
3186 Handle(TColStd_HArray1OfInteger) ContactOrder =
3187 new (TColStd_HArray1OfInteger)(1, VLength);
3188 ContactOrder->Init(0);
3189
3190 for (ii=1; ii <= VLength; ii++) {
3191
3192 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3193 Points(ii,ll) = Poles(jj, ii).X();
3194 Points(ii,ll+1) = Poles(jj, ii).Y();
3195 Points(ii,ll+2) = Poles(jj, ii).Z();
3196 }
3197 }
3198
0d969553 3199 // interpolation of iso u
7fd59977 3200 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3201 BSplCLib::Interpolate(VDegree,
3202 VFlatKnots,
3203 VParameters,
3204 ContactOrder->Array1(),
3205 dimension,
3206 poles_array[0],
3207 InversionProblem) ;
3208 if (InversionProblem != 0) return;
3209
0d969553 3210 // extraction of iso v
7fd59977 3211
3212 dimension = VLength*3;
3213 TColStd_Array2OfReal IsoPoles(1, ULength,
3214 1, dimension);
3215
3216 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3217 ContactOrder->Init(0);
3218 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3219
3220 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3221
3222 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3223 IsoPoles (ii,ll) = Points(jj, kk);
3224 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3225 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3226 }
3227 }
0d969553 3228 // interpolation of iso v
7fd59977 3229 BSplCLib::Interpolate(UDegree,
3230 UFlatKnots,
3231 UParameters,
3232 ContactOrder->Array1(),
3233 dimension,
3234 poles_array[0],
3235 InversionProblem);
3236
0d969553 3237 // return results
7fd59977 3238
3239 for (ii=1; ii <= ULength; ii++) {
3240
3241 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3242 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3243 Poles.SetValue(ii, jj, Pnt);
3244 }
3245 }
3246}
3247
3248//=======================================================================
3249//function : FunctionMultiply
3250//purpose :
3251//=======================================================================
3252
3253void BSplSLib::FunctionMultiply
3254(const BSplSLib_EvaluatorFunction& Function,
3255 const Standard_Integer UBSplineDegree,
3256 const Standard_Integer VBSplineDegree,
3257 const TColStd_Array1OfReal& UBSplineKnots,
3258 const TColStd_Array1OfReal& VBSplineKnots,
3259 const TColStd_Array1OfInteger & UMults,
3260 const TColStd_Array1OfInteger & VMults,
3261 const TColgp_Array2OfPnt& Poles,
3262 const TColStd_Array2OfReal& Weights,
3263 const TColStd_Array1OfReal& UFlatKnots,
3264 const TColStd_Array1OfReal& VFlatKnots,
3265 const Standard_Integer UNewDegree,
3266 const Standard_Integer VNewDegree,
3267 TColgp_Array2OfPnt& NewNumerator,
3268 TColStd_Array2OfReal& NewDenominator,
3269 Standard_Integer& Status)
3270{
3271 Standard_Integer num_uparameters,
3272// ii,jj,kk,
3273 ii,jj,
3274 error_code,
3275 num_vparameters ;
3276 Standard_Real result ;
3277
3278 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3279 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3280 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3281 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3282
3283 if ((NewNumerator.ColLength() == num_uparameters) &&
3284 (NewNumerator.RowLength() == num_vparameters) &&
3285 (NewDenominator.ColLength() == num_uparameters) &&
3286 (NewDenominator.RowLength() == num_vparameters)) {
3287
3288
3289 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3290 UFlatKnots,
3291 UParameters) ;
3292
3293 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3294 VFlatKnots,
3295 VParameters) ;
3296
3297 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3298
3299 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3300 HomogeneousD0(UParameters(ii),
3301 VParameters(jj),
3302 0,
3303 0,
3304 Poles,
3305 Weights,
3306 UBSplineKnots,
3307 VBSplineKnots,
3308 UMults,
3309 VMults,
3310 UBSplineDegree,
3311 VBSplineDegree,
3312 Standard_True,
3313 Standard_True,
3314 Standard_False,
3315 Standard_False,
3316 NewDenominator(ii,jj),
3317 NewNumerator(ii,jj)) ;
3318
3319 Function(0,
3320 UParameters(ii),
3321 VParameters(jj),
3322 result,
3323 error_code) ;
3324 if (error_code) {
3325 Standard_ConstructionError::Raise();
3326 }
3327 gp_Pnt& P = NewNumerator(ii,jj);
3328 P.SetX(P.X() * result);
3329 P.SetY(P.Y() * result);
3330 P.SetZ(P.Z() * result);
3331 NewDenominator(ii,jj) *= result ;
3332 }
3333 }
3334 Interpolate(UNewDegree,
3335 VNewDegree,
3336 UFlatKnots,
3337 VFlatKnots,
3338 UParameters,
3339 VParameters,
3340 NewNumerator,
3341 NewDenominator,
3342 Status) ;
3343 }
3344 else {
3345 Standard_ConstructionError::Raise();
3346 }
3347}
3348
Open CASCADE Technology repository