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