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