0022550: Fixing data races
[occt.git] / src / BSplCLib / BSplCLib.cxx
CommitLineData
7fd59977 1// File: BSplCLib.cxx
2// Created: Fri Aug 9 16:32:46 1991
3// Author: JCV
4
5//
6// Modified RLE 9 Sep 1993
7//
8// pmn : modified 28-01-97 : fixed a mistake in LocateParameter (PRO6973)
9// pmn : modified 4-11-96 : fixed a mistake in BuildKnots (PRO6124)
10// pmn : modified 28-Jun-96 : fixed a mistake in AntiBoorScheme
11// xab : modified 15-Jun-95 : fixed a mistake in IsRational
12// xab : modified 15-Mar-95 : removed Epsilon comparison in IsRational
13// added RationalDerivatives.
14// xab : 30-Mar-95 : fixed coupling with lti in RationalDerivatives
15// xab : 15-Mar-96 : fixed a typo in Eval with extrapolation
16// jct : 15-Apr-97 : added TangExtendToConstraint
17// jct : 24-Apr-97 : correction on computation of Tbord and NewFlatKnots
18// in TangExtendToConstraint; Continuity can be equal to 0
19//
20
7fd59977 21#include <BSplCLib.ixx>
22#include <PLib.hxx>
41194117 23#include <PLib_LocalArray.hxx>
7fd59977 24#include <Precision.hxx>
25#include <Standard_NotImplemented.hxx>
26
27typedef gp_Pnt Pnt;
28typedef gp_Vec Vec;
29typedef TColgp_Array1OfPnt Array1OfPnt;
30typedef TColStd_Array1OfReal Array1OfReal;
31typedef TColStd_Array1OfInteger Array1OfInteger;
32
33//=======================================================================
34//class : BSplCLib_LocalMatrix
35//purpose: Auxiliary class optimizing creation of matrix buffer for
36// evaluation of bspline (using stack allocation for main matrix)
37//=======================================================================
38
39class BSplCLib_LocalMatrix : public math_Matrix
40{
41194117 41public:
7fd59977 42 BSplCLib_LocalMatrix (Standard_Integer DerivativeRequest, Standard_Integer Order)
43 : math_Matrix (myBuffer, 1, DerivativeRequest + 1, 1, Order)
44 {
41194117
K
45 Standard_OutOfRange_Raise_if (DerivativeRequest > BSplCLib::MaxDegree() ||
46 Order > BSplCLib::MaxDegree()+1 || BSplCLib::MaxDegree() > 25,
47 "BSplCLib: bspline degree is greater than maximum supported");
7fd59977 48 }
41194117 49
7fd59977 50 private:
51 // local buffer, to be sufficient for addressing by index [Degree+1][Degree+1]
52 // (see math_Matrix implementation)
53 Standard_Real myBuffer[27*27];
54};
55
41194117 56typedef PLib_LocalArray BSplCLib_LocalArray;
7fd59977 57
58//=======================================================================
59//function : Hunt
60//purpose :
61//=======================================================================
62
63void BSplCLib::Hunt (const Array1OfReal& XX,
64 const Standard_Real X,
65 Standard_Integer& Ilc)
66{
67 // replaced by simple dichotomy (RLE)
68 Ilc = XX.Lower();
69 const Standard_Real *px = &XX(Ilc);
70 px -= Ilc;
71
72 if (X < px[Ilc]) {
73 Ilc--;
74 return;
75 }
76 Standard_Integer Ihi = XX.Upper();
77 if (X > px[Ihi]) {
78 Ilc = Ihi + 1;
79 return;
80 }
81 Standard_Integer Im;
82
83 while (Ihi - Ilc != 1) {
84 Im = (Ihi + Ilc) >> 1;
85 if (X > px[Im]) Ilc = Im;
86 else Ihi = Im;
87 }
88}
89
90//=======================================================================
91//function : FirstUKnotIndex
92//purpose :
93//=======================================================================
94
95Standard_Integer BSplCLib::FirstUKnotIndex (const Standard_Integer Degree,
96 const TColStd_Array1OfInteger& Mults)
97{
98 Standard_Integer Index = Mults.Lower();
99 Standard_Integer SigmaMult = Mults(Index);
100
101 while (SigmaMult <= Degree) {
102 Index++;
103 SigmaMult += Mults (Index);
104 }
105 return Index;
106}
107
108//=======================================================================
109//function : LastUKnotIndex
110//purpose :
111//=======================================================================
112
113Standard_Integer BSplCLib::LastUKnotIndex (const Standard_Integer Degree,
114 const Array1OfInteger& Mults)
115{
116 Standard_Integer Index = Mults.Upper();
117 Standard_Integer SigmaMult = Mults(Index);
118
119 while (SigmaMult <= Degree) {
120 Index--;
121 SigmaMult += Mults.Value (Index);
122 }
123 return Index;
124}
125
126//=======================================================================
127//function : FlatIndex
128//purpose :
129//=======================================================================
130
131Standard_Integer BSplCLib::FlatIndex
132 (const Standard_Integer Degree,
133 const Standard_Integer Index,
134 const TColStd_Array1OfInteger& Mults,
135 const Standard_Boolean Periodic)
136{
137 Standard_Integer i, index = Index;
138 const Standard_Integer MLower = Mults.Lower();
139 const Standard_Integer *pmu = &Mults(MLower);
140 pmu -= MLower;
141
142 for (i = MLower + 1; i <= Index; i++)
143 index += pmu[i] - 1;
144 if ( Periodic)
145 index += Degree;
146 else
147 index += pmu[MLower] - 1;
148 return index;
149}
150
151//=======================================================================
152//function : LocateParameter
0d969553
Y
153//purpose : Processing of nodes with multiplicities
154//pmn 28-01-97 -> compute eventual of the period.
7fd59977 155//=======================================================================
156
157void BSplCLib::LocateParameter
158(const Standard_Integer , //Degree,
159 const Array1OfReal& Knots,
160 const Array1OfInteger& , //Mults,
161 const Standard_Real U,
162 const Standard_Boolean IsPeriodic,
163 const Standard_Integer FromK1,
164 const Standard_Integer ToK2,
165 Standard_Integer& KnotIndex,
166 Standard_Real& NewU)
167{
168 Standard_Real uf = 0, ul=1;
169 if (IsPeriodic) {
170 uf = Knots(Knots.Lower());
171 ul = Knots(Knots.Upper());
172 }
173 BSplCLib::LocateParameter(Knots,U,IsPeriodic,FromK1,ToK2,
174 KnotIndex,NewU, uf, ul);
175}
176
177//=======================================================================
178//function : LocateParameter
0d969553
Y
179//purpose : For plane nodes
180// pmn 28-01-97 -> There is a need of the degre to calculate
181// the eventual period
7fd59977 182//=======================================================================
183
184void BSplCLib::LocateParameter
185(const Standard_Integer Degree,
186 const Array1OfReal& Knots,
187 const Standard_Real U,
188 const Standard_Boolean IsPeriodic,
189 const Standard_Integer FromK1,
190 const Standard_Integer ToK2,
191 Standard_Integer& KnotIndex,
192 Standard_Real& NewU)
193{
194 if (IsPeriodic)
195 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
196 KnotIndex, NewU,
197 Knots(Knots.Lower() + Degree),
198 Knots(Knots.Upper() - Degree));
199 else
200 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
201 KnotIndex, NewU,
202 0.,
203 1.);
204}
205
206//=======================================================================
207//function : LocateParameter
0d969553
Y
208//purpose : Effective computation
209// pmn 28-01-97 : Add limits of the period as input argument,
210// as it is imposible to produce them at this level.
7fd59977 211//=======================================================================
212
213void BSplCLib::LocateParameter
214(const TColStd_Array1OfReal& Knots,
215 const Standard_Real U,
216 const Standard_Boolean IsPeriodic,
217 const Standard_Integer FromK1,
218 const Standard_Integer ToK2,
219 Standard_Integer& KnotIndex,
220 Standard_Real& NewU,
221 const Standard_Real UFirst,
222 const Standard_Real ULast)
223{
224 Standard_Integer First,Last;
225 if (FromK1 < ToK2) {
226 First = FromK1;
227 Last = ToK2;
228 }
229 else {
230 First = ToK2;
231 Last = FromK1;
232 }
233 Standard_Integer Last1 = Last - 1;
234 NewU = U;
235 if (IsPeriodic) {
236 Standard_Real Period = ULast - UFirst;
237
238 while (NewU > ULast )
239 NewU -= Period;
240
241 while (NewU < UFirst)
242 NewU += Period;
243 }
244
245 BSplCLib::Hunt (Knots, NewU, KnotIndex);
246
247 Standard_Real Eps = Epsilon(U);
248 Standard_Real val;
249 if (Eps < 0) Eps = - Eps;
250 Standard_Integer KLower = Knots.Lower();
251 const Standard_Real *knots = &Knots(KLower);
252 knots -= KLower;
253 if ( KnotIndex < Knots.Upper()) {
254 val = NewU - knots[KnotIndex + 1];
255 if (val < 0) val = - val;
0d969553 256 // <= to be coherent with Segment where Eps corresponds to a bit of error.
7fd59977 257 if (val <= Eps) KnotIndex++;
258 }
259 if (KnotIndex < First) KnotIndex = First;
260 if (KnotIndex > Last1) KnotIndex = Last1;
261
262 if (KnotIndex != Last1) {
263 Standard_Real K1 = knots[KnotIndex];
264 Standard_Real K2 = knots[KnotIndex + 1];
265 val = K2 - K1;
266 if (val < 0) val = - val;
267
268 while (val <= Eps) {
269 KnotIndex++;
270 K1 = K2;
271 K2 = knots[KnotIndex + 1];
272 val = K2 - K1;
273 if (val < 0) val = - val;
274 }
275 }
276}
277
278//=======================================================================
279//function : LocateParameter
280//purpose : the index is recomputed only if out of range
0d969553 281//pmn 28-01-97 -> eventual computation of the period.
7fd59977 282//=======================================================================
283
284void BSplCLib::LocateParameter
285(const Standard_Integer Degree,
286 const TColStd_Array1OfReal& Knots,
287 const TColStd_Array1OfInteger& Mults,
288 const Standard_Real U,
289 const Standard_Boolean Periodic,
290 Standard_Integer& KnotIndex,
291 Standard_Real& NewU)
292{
293 Standard_Integer first,last;
294 if (&Mults) {
295 if (Periodic) {
296 first = Knots.Lower();
297 last = Knots.Upper();
298 }
299 else {
300 first = FirstUKnotIndex(Degree,Mults);
301 last = LastUKnotIndex (Degree,Mults);
302 }
303 }
304 else {
305 first = Knots.Lower() + Degree;
306 last = Knots.Upper() - Degree;
307 }
308 if ( KnotIndex < first || KnotIndex > last)
309 BSplCLib::LocateParameter(Knots, U, Periodic, first, last,
310 KnotIndex, NewU, Knots(first), Knots(last));
311 else
312 NewU = U;
313}
314
315//=======================================================================
316//function : MaxKnotMult
317//purpose :
318//=======================================================================
319
320Standard_Integer BSplCLib::MaxKnotMult
321(const Array1OfInteger& Mults,
322 const Standard_Integer FromK1,
323 const Standard_Integer ToK2)
324{
325 Standard_Integer MLower = Mults.Lower();
326 const Standard_Integer *pmu = &Mults(MLower);
327 pmu -= MLower;
328 Standard_Integer MaxMult = pmu[FromK1];
329
330 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
331 if (MaxMult < pmu[i]) MaxMult = pmu[i];
332 }
333 return MaxMult;
334}
335
336//=======================================================================
337//function : MinKnotMult
338//purpose :
339//=======================================================================
340
341Standard_Integer BSplCLib::MinKnotMult
342(const Array1OfInteger& Mults,
343 const Standard_Integer FromK1,
344 const Standard_Integer ToK2)
345{
346 Standard_Integer MLower = Mults.Lower();
347 const Standard_Integer *pmu = &Mults(MLower);
348 pmu -= MLower;
349 Standard_Integer MinMult = pmu[FromK1];
350
351 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
352 if (MinMult > pmu[i]) MinMult = pmu[i];
353 }
354 return MinMult;
355}
356
357//=======================================================================
358//function : NbPoles
359//purpose :
360//=======================================================================
361
362Standard_Integer BSplCLib::NbPoles(const Standard_Integer Degree,
363 const Standard_Boolean Periodic,
364 const TColStd_Array1OfInteger& Mults)
365{
366 Standard_Integer i,sigma = 0;
367 Standard_Integer f = Mults.Lower();
368 Standard_Integer l = Mults.Upper();
369 const Standard_Integer * pmu = &Mults(f);
370 pmu -= f;
371 Standard_Integer Mf = pmu[f];
372 Standard_Integer Ml = pmu[l];
373 if (Mf <= 0) return 0;
374 if (Ml <= 0) return 0;
375 if (Periodic) {
376 if (Mf > Degree) return 0;
377 if (Ml > Degree) return 0;
378 if (Mf != Ml ) return 0;
379 sigma = Mf;
380 }
381 else {
382 Standard_Integer Deg1 = Degree + 1;
383 if (Mf > Deg1) return 0;
384 if (Ml > Deg1) return 0;
385 sigma = Mf + Ml - Deg1;
386 }
387
388 for (i = f + 1; i < l; i++) {
389 if (pmu[i] <= 0 ) return 0;
390 if (pmu[i] > Degree) return 0;
391 sigma += pmu[i];
392 }
393 return sigma;
394}
395
396//=======================================================================
397//function : KnotSequenceLength
398//purpose :
399//=======================================================================
400
401Standard_Integer BSplCLib::KnotSequenceLength
402(const TColStd_Array1OfInteger& Mults,
403 const Standard_Integer Degree,
404 const Standard_Boolean Periodic)
405{
406 Standard_Integer i,l = 0;
407 Standard_Integer MLower = Mults.Lower();
408 Standard_Integer MUpper = Mults.Upper();
409 const Standard_Integer * pmu = &Mults(MLower);
410 pmu -= MLower;
411
412 for (i = MLower; i <= MUpper; i++)
413 l += pmu[i];
414 if (Periodic) l += 2 * (Degree + 1 - pmu[MLower]);
415 return l;
416}
417
418//=======================================================================
419//function : KnotSequence
420//purpose :
421//=======================================================================
422
423void BSplCLib::KnotSequence
424(const TColStd_Array1OfReal& Knots,
425 const TColStd_Array1OfInteger& Mults,
426 TColStd_Array1OfReal& KnotSeq)
427{
428 BSplCLib::KnotSequence(Knots,Mults,0,Standard_False,KnotSeq);
429}
430
431//=======================================================================
432//function : KnotSequence
433//purpose :
434//=======================================================================
435
436void BSplCLib::KnotSequence
437(const TColStd_Array1OfReal& Knots,
438 const TColStd_Array1OfInteger& Mults,
439 const Standard_Integer Degree,
440 const Standard_Boolean Periodic,
441 TColStd_Array1OfReal& KnotSeq)
442{
443 Standard_Real K;
444 Standard_Integer Mult;
445 Standard_Integer MLower = Mults.Lower();
446 const Standard_Integer * pmu = &Mults(MLower);
447 pmu -= MLower;
448 Standard_Integer KLower = Knots.Lower();
449 Standard_Integer KUpper = Knots.Upper();
450 const Standard_Real * pkn = &Knots(KLower);
451 pkn -= KLower;
452 Standard_Integer M1 = Degree + 1 - pmu[MLower]; // for periodic
453 Standard_Integer i,j,index = Periodic ? M1 + 1 : 1;
454
455 for (i = KLower; i <= KUpper; i++) {
456 Mult = pmu[i];
457 K = pkn[i];
458
459 for (j = 1; j <= Mult; j++) {
460 KnotSeq (index) = K;
461 index++;
462 }
463 }
464 if (Periodic) {
465 Standard_Real period = pkn[KUpper] - pkn[KLower];
466 Standard_Integer m;
467 m = 1;
468 j = KUpper - 1;
469
470 for (i = M1; i >= 1; i--) {
471 KnotSeq(i) = pkn[j] - period;
472 m++;
473 if (m > pmu[j]) {
474 j--;
475 m = 1;
476 }
477 }
478 m = 1;
479 j = KLower + 1;
480
481 for (i = index; i <= KnotSeq.Upper(); i++) {
482 KnotSeq(i) = pkn[j] + period;
483 m++;
484 if (m > pmu[j]) {
485 j++;
486 m = 1;
487 }
488 }
489 }
490}
491
492//=======================================================================
493//function : KnotsLength
494//purpose :
495//=======================================================================
496 Standard_Integer BSplCLib::KnotsLength(const TColStd_Array1OfReal& SeqKnots,
497// const Standard_Boolean Periodic)
498 const Standard_Boolean )
499{
500 Standard_Integer sizeMult = 1;
501 Standard_Real val = SeqKnots(1);
502 for (Standard_Integer jj=2;
503 jj<=SeqKnots.Length();jj++)
504 {
505 // test on strict equality on nodes
506 if (SeqKnots(jj)!=val)
507 {
508 val = SeqKnots(jj);
509 sizeMult++;
510 }
511 }
512 return sizeMult;
513}
514
515//=======================================================================
516//function : Knots
517//purpose :
518//=======================================================================
519void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots,
520 TColStd_Array1OfReal &knots,
521 TColStd_Array1OfInteger &mult,
522// const Standard_Boolean Periodic)
523 const Standard_Boolean )
524{
525 Standard_Real val = SeqKnots(1);
526 Standard_Integer kk=1;
527 knots(kk) = val;
528 mult(kk) = 1;
529
530 for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++)
531 {
532 // test on strict equality on nodes
533 if (SeqKnots(jj)!=val)
534 {
535 val = SeqKnots(jj);
536 kk++;
537 knots(kk) = val;
538 mult(kk) = 1;
539 }
540 else
541 {
542 mult(kk)++;
543 }
544 }
545}
546
547//=======================================================================
548//function : KnotForm
549//purpose :
550//=======================================================================
551
552BSplCLib_KnotDistribution BSplCLib::KnotForm
553(const Array1OfReal& Knots,
554 const Standard_Integer FromK1,
555 const Standard_Integer ToK2)
556{
557 Standard_Real DU0,DU1,Ui,Uj,Eps0,val;
558 BSplCLib_KnotDistribution KForm = BSplCLib_Uniform;
559
560 Standard_Integer KLower = Knots.Lower();
561 const Standard_Real * pkn = &Knots(KLower);
562 pkn -= KLower;
563 Ui = pkn[FromK1];
564 if (Ui < 0) Ui = - Ui;
565 Uj = pkn[FromK1 + 1];
566 if (Uj < 0) Uj = - Uj;
567 DU0 = Uj - Ui;
568 if (DU0 < 0) DU0 = - DU0;
569 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
570 Standard_Integer i = FromK1 + 1;
571
572 while (KForm != BSplCLib_NonUniform && i < ToK2) {
573 Ui = pkn[i];
574 if (Ui < 0) Ui = - Ui;
575 i++;
576 Uj = pkn[i];
577 if (Uj < 0) Uj = - Uj;
578 DU1 = Uj - Ui;
579 if (DU1 < 0) DU1 = - DU1;
580 val = DU1 - DU0;
581 if (val < 0) val = -val;
582 if (val > Eps0) KForm = BSplCLib_NonUniform;
583 DU0 = DU1;
584 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
585 }
586 return KForm;
587}
588
589//=======================================================================
590//function : MultForm
591//purpose :
592//=======================================================================
593
594BSplCLib_MultDistribution BSplCLib::MultForm
595(const Array1OfInteger& Mults,
596 const Standard_Integer FromK1,
597 const Standard_Integer ToK2)
598{
599 Standard_Integer First,Last;
600 if (FromK1 < ToK2) {
601 First = FromK1;
602 Last = ToK2;
603 }
604 else {
605 First = ToK2;
606 Last = FromK1;
607 }
608 Standard_Integer MLower = Mults.Lower();
609 const Standard_Integer *pmu = &Mults(MLower);
610 pmu -= MLower;
611 Standard_Integer FirstMult = pmu[First];
612 BSplCLib_MultDistribution MForm = BSplCLib_Constant;
613 Standard_Integer i = First + 1;
614 Standard_Integer Mult = pmu[i];
615
616// while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR????????
617 while (MForm != BSplCLib_NonConstant && i <= Last) {
618 if (i == First + 1) {
619 if (Mult != FirstMult) MForm = BSplCLib_QuasiConstant;
620 }
621 else if (i == Last) {
622 if (MForm == BSplCLib_QuasiConstant) {
623 if (FirstMult != pmu[i]) MForm = BSplCLib_NonConstant;
624 }
625 else {
626 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
627 }
628 }
629 else {
630 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
631 Mult = pmu[i];
632 }
633 i++;
634 }
635 return MForm;
636}
637
638//=======================================================================
639//function : Reparametrize
640//purpose :
641//=======================================================================
642
643void BSplCLib::Reparametrize
644(const Standard_Real U1,
645 const Standard_Real U2,
646 Array1OfReal& Knots)
647{
648 Standard_Integer Lower = Knots.Lower();
649 Standard_Integer Upper = Knots.Upper();
650 Standard_Real UFirst = Min (U1, U2);
651 Standard_Real ULast = Max (U1, U2);
652 Standard_Real NewLength = ULast - UFirst;
653 BSplCLib_KnotDistribution KSet = BSplCLib::KnotForm (Knots, Lower, Upper);
654 if (KSet == BSplCLib_Uniform) {
655 Standard_Real DU = NewLength / (Upper - Lower);
656 Knots (Lower) = UFirst;
657
658 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
659 Knots (i) = Knots (i-1) + DU;
660 }
661 }
662 else {
663 Standard_Real K2;
664 Standard_Real Ratio;
665 Standard_Real K1 = Knots (Lower);
666 Standard_Real Length = Knots (Upper) - Knots (Lower);
667 Knots (Lower) = UFirst;
668
669 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
670 K2 = Knots (i);
671 Ratio = (K2 - K1) / Length;
672 Knots (i) = Knots (i-1) + (NewLength * Ratio);
673
674 //for CheckCurveData
675 Standard_Real Eps = Epsilon( Abs(Knots(i-1)) );
676 if (Knots(i) - Knots(i-1) <= Eps)
677 Knots(i) += 1.1*Eps;
678
679 K1 = K2;
680 }
681 }
682}
683
684//=======================================================================
685//function : Reverse
686//purpose :
687//=======================================================================
688
689void BSplCLib::Reverse(TColStd_Array1OfReal& Knots)
690{
691 Standard_Integer first = Knots.Lower();
692 Standard_Integer last = Knots.Upper();
693 Standard_Real kfirst = Knots(first);
694 Standard_Real klast = Knots(last);
695 Standard_Real tfirst = kfirst;
696 Standard_Real tlast = klast;
697 first++;
698 last--;
699
700 while (first <= last) {
701 tfirst += klast - Knots(last);
702 tlast -= Knots(first) - kfirst;
703 kfirst = Knots(first);
704 klast = Knots(last);
705 Knots(first) = tfirst;
706 Knots(last) = tlast;
707 first++;
708 last--;
709 }
710}
711
712//=======================================================================
713//function : Reverse
714//purpose :
715//=======================================================================
716
717void BSplCLib::Reverse(TColStd_Array1OfInteger& Mults)
718{
719 Standard_Integer first = Mults.Lower();
720 Standard_Integer last = Mults.Upper();
721 Standard_Integer temp;
722
723 while (first < last) {
724 temp = Mults(first);
725 Mults(first) = Mults(last);
726 Mults(last) = temp;
727 first++;
728 last--;
729 }
730}
731
732//=======================================================================
733//function : Reverse
734//purpose :
735//=======================================================================
736
737void BSplCLib::Reverse(TColStd_Array1OfReal& Weights,
738 const Standard_Integer L)
739{
740 Standard_Integer i, l = L;
741 l = Weights.Lower()+(l-Weights.Lower())%(Weights.Upper()-Weights.Lower()+1);
742
743 TColStd_Array1OfReal temp(0,Weights.Length()-1);
744
745 for (i = Weights.Lower(); i <= l; i++)
746 temp(l-i) = Weights(i);
747
748 for (i = l+1; i <= Weights.Upper(); i++)
749 temp(l-Weights.Lower()+Weights.Upper()-i+1) = Weights(i);
750
751 for (i = Weights.Lower(); i <= Weights.Upper(); i++)
752 Weights(i) = temp(i-Weights.Lower());
753}
754
755//=======================================================================
756//function : IsRational
757//purpose :
758//=======================================================================
759
760Standard_Boolean BSplCLib::IsRational(const TColStd_Array1OfReal& Weights,
761 const Standard_Integer I1,
762 const Standard_Integer I2,
763// const Standard_Real Epsi)
764 const Standard_Real )
765{
766 Standard_Integer i, f = Weights.Lower(), l = Weights.Length();
767 Standard_Integer I3 = I2 - f;
768 const Standard_Real * WG = &Weights(f);
769 WG -= f;
770
771 for (i = I1 - f; i < I3; i++) {
772 if (WG[f + (i % l)] != WG[f + ((i + 1) % l)]) return Standard_True;
773 }
774 return Standard_False ;
775}
776
777//=======================================================================
778//function : Eval
779//purpose : evaluate point and derivatives
780//=======================================================================
781
782void BSplCLib::Eval(const Standard_Real U,
783 const Standard_Integer Degree,
784 Standard_Real& Knots,
785 const Standard_Integer Dimension,
786 Standard_Real& Poles)
787{
788 Standard_Integer step,i,Dms,Dm1,Dpi,Sti;
789 Standard_Real X, Y, *poles, *knots = &Knots;
790 Dm1 = Dms = Degree;
791 Dm1--;
792 Dms++;
793 switch (Dimension) {
794
795 case 1 : {
796
797 for (step = - 1; step < Dm1; step++) {
798 Dms--;
799 poles = &Poles;
800 Dpi = Dm1;
801 Sti = step;
802
803 for (i = 0; i < Dms; i++) {
804 Dpi++;
805 Sti++;
806 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
807 Y = 1 - X;
808 poles[0] *= X; poles[0] += Y * poles[1];
809 poles += 1;
810 }
811 }
812 break;
813 }
814 case 2 : {
815
816 for (step = - 1; step < Dm1; step++) {
817 Dms--;
818 poles = &Poles;
819 Dpi = Dm1;
820 Sti = step;
821
822 for (i = 0; i < Dms; i++) {
823 Dpi++;
824 Sti++;
825 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
826 Y = 1 - X;
827 poles[0] *= X; poles[0] += Y * poles[2];
828 poles[1] *= X; poles[1] += Y * poles[3];
829 poles += 2;
830 }
831 }
832 break;
833 }
834 case 3 : {
835
836 for (step = - 1; step < Dm1; step++) {
837 Dms--;
838 poles = &Poles;
839 Dpi = Dm1;
840 Sti = step;
841
842 for (i = 0; i < Dms; i++) {
843 Dpi++;
844 Sti++;
845 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
846 Y = 1 - X;
847 poles[0] *= X; poles[0] += Y * poles[3];
848 poles[1] *= X; poles[1] += Y * poles[4];
849 poles[2] *= X; poles[2] += Y * poles[5];
850 poles += 3;
851 }
852 }
853 break;
854 }
855 case 4 : {
856
857 for (step = - 1; step < Dm1; step++) {
858 Dms--;
859 poles = &Poles;
860 Dpi = Dm1;
861 Sti = step;
862
863 for (i = 0; i < Dms; i++) {
864 Dpi++;
865 Sti++;
866 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
867 Y = 1 - X;
868 poles[0] *= X; poles[0] += Y * poles[4];
869 poles[1] *= X; poles[1] += Y * poles[5];
870 poles[2] *= X; poles[2] += Y * poles[6];
871 poles[3] *= X; poles[3] += Y * poles[7];
872 poles += 4;
873 }
874 }
875 break;
876 }
877 default : {
878 Standard_Integer k;
879
880 for (step = - 1; step < Dm1; step++) {
881 Dms--;
882 poles = &Poles;
883 Dpi = Dm1;
884 Sti = step;
885
886 for (i = 0; i < Dms; i++) {
887 Dpi++;
888 Sti++;
889 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
890 Y = 1 - X;
891
892 for (k = 0; k < Dimension; k++) {
893 poles[k] *= X;
894 poles[k] += Y * poles[k + Dimension];
895 }
896 poles += Dimension;
897 }
898 }
899 }
900 }
901}
902
903//=======================================================================
904//function : BoorScheme
905//purpose :
906//=======================================================================
907
908void BSplCLib::BoorScheme(const Standard_Real U,
909 const Standard_Integer Degree,
910 Standard_Real& Knots,
911 const Standard_Integer Dimension,
912 Standard_Real& Poles,
913 const Standard_Integer Depth,
914 const Standard_Integer Length)
915{
916 //
917 // Compute the values
918 //
919 // P(i,j) (U).
920 //
921 // for i = 0 to Depth,
922 // j = 0 to Length - i
923 //
924 // The Boor scheme is :
925 //
926 // P(0,j) = Pole(j)
927 // P(i,j) = x * P(i-1,j) + (1-x) * P(i-1,j+1)
928 //
929 // where x = (knot(i+j+Degree) - U) / (knot(i+j+Degree) - knot(i+j))
930 //
931 //
932 // The values are stored in the array Poles
933 // They are alternatively written if the odd and even positions.
934 //
935 // The successives contents of the array are
936 // ***** means unitialised, l = Degree + Length
937 //
938 // P(0,0) ****** P(0,1) ...... P(0,l-1) ******** P(0,l)
939 // P(0,0) P(1,0) P(0,1) ...... P(0,l-1) P(1,l-1) P(0,l)
940 // P(0,0) P(1,0) P(2,0) ...... P(2,l-1) P(1,l-1) P(0,l)
941 //
942
943 Standard_Integer i,k,step;
944 Standard_Real *knots = &Knots;
945 Standard_Real *pole, *firstpole = &Poles - 2 * Dimension;
946 // the steps of recursion
947
948 for (step = 0; step < Depth; step++) {
949 firstpole += Dimension;
950 pole = firstpole;
951 // compute the new row of poles
952
953 for (i = step; i < Length; i++) {
954 pole += 2 * Dimension;
955 // coefficient
956 Standard_Real X = (knots[i+Degree-step] - U)
957 / (knots[i+Degree-step] - knots[i]);
958 Standard_Real Y = 1. - X;
959 // Value
960 // P(i,j) = X * P(i-1,j) + (1-X) * P(i-1,j+1)
961
962 for (k = 0; k < Dimension; k++)
963 pole[k] = X * pole[k - Dimension] + Y * pole[k + Dimension];
964 }
965 }
966}
967
968//=======================================================================
969//function : AntiBoorScheme
970//purpose :
971//=======================================================================
972
973Standard_Boolean BSplCLib::AntiBoorScheme(const Standard_Real U,
974 const Standard_Integer Degree,
975 Standard_Real& Knots,
976 const Standard_Integer Dimension,
977 Standard_Real& Poles,
978 const Standard_Integer Depth,
979 const Standard_Integer Length,
980 const Standard_Real Tolerance)
981{
982 // do the Boor scheme reverted.
983
984 Standard_Integer i,k,step, half_length;
985 Standard_Real *knots = &Knots;
986 Standard_Real z,X,Y,*pole, *firstpole = &Poles + (Depth-1) * Dimension;
987
988 // Test the special case length = 1
989 // only verification of the central point
990
991 if (Length == 1) {
992 X = (knots[Degree] - U) / (knots[Degree] - knots[0]);
993 Y = 1. - X;
994
995 for (k = 0; k < Dimension; k++) {
996 z = X * firstpole[k] + Y * firstpole[k+2*Dimension];
997 if (Abs(z - firstpole[k+Dimension]) > Tolerance)
998 return Standard_False;
999 }
1000 return Standard_True;
1001 }
1002
1003 // General case
1004 // the steps of recursion
1005
1006 for (step = Depth-1; step >= 0; step--) {
1007 firstpole -= Dimension;
1008 pole = firstpole;
1009
1010 // first step from left to right
1011
1012 for (i = step; i < Length-1; i++) {
1013 pole += 2 * Dimension;
1014
1015 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1016 Y = 1. - X;
1017
1018 for (k = 0; k < Dimension; k++)
1019 pole[k+Dimension] = (pole[k] - X*pole[k-Dimension]) / Y;
1020
1021 }
1022
1023 // second step from right to left
1024 pole += 4* Dimension;
1025 half_length = (Length - 1 + step) / 2 ;
1026 //
1027 // only do half of the way from right to left
1028 // otherwise it start degenerating because of
1029 // overflows
1030 //
1031
1032 for (i = Length-1; i > half_length ; i--) {
1033 pole -= 2 * Dimension;
1034
1035 // coefficient
1036 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1037 Y = 1. - X;
1038
1039 for (k = 0; k < Dimension; k++) {
1040 z = (pole[k] - Y * pole[k+Dimension]) / X;
1041 if (Abs(z-pole[k-Dimension]) > Tolerance)
1042 return Standard_False;
1043 pole[k-Dimension] += z;
1044 pole[k-Dimension] /= 2.;
1045 }
1046 }
1047 }
1048 return Standard_True;
1049}
1050
1051//=======================================================================
1052//function : Derivative
1053//purpose :
1054//=======================================================================
1055
1056void BSplCLib::Derivative(const Standard_Integer Degree,
1057 Standard_Real& Knots,
1058 const Standard_Integer Dimension,
1059 const Standard_Integer Length,
1060 const Standard_Integer Order,
1061 Standard_Real& Poles)
1062{
1063 Standard_Integer i,k,step,span = Degree;
1064 Standard_Real *knot = &Knots;
1065
1066 for (step = 1; step <= Order; step++) {
1067 Standard_Real* pole = &Poles;
1068
1069 for (i = step; i < Length; i++) {
1070 Standard_Real coef = - span / (knot[i+span] - knot[i]);
1071
1072 for (k = 0; k < Dimension; k++) {
1073 pole[k] -= pole[k+Dimension];
1074 pole[k] *= coef;
1075 }
1076 pole += Dimension;
1077 }
1078 span--;
1079 }
1080}
1081
1082//=======================================================================
1083//function : Bohm
1084//purpose :
1085//=======================================================================
1086
1087void BSplCLib::Bohm(const Standard_Real U,
1088 const Standard_Integer Degree,
1089 const Standard_Integer N,
1090 Standard_Real& Knots,
1091 const Standard_Integer Dimension,
1092 Standard_Real& Poles)
1093{
1094 // First phase independant of U, compute the poles of the derivatives
1095 Standard_Integer i,j,iDim,min,Dmi,DDmi,jDmi,Degm1;
1096 Standard_Real *knot = &Knots, *pole, coef, *tbis, *psav, *psDD, *psDDmDim;
1097 psav = &Poles;
1098 if (N < Degree) min = N;
1099 else min = Degree;
1100 Degm1 = Degree - 1;
1101 DDmi = (Degree << 1) + 1;
1102 switch (Dimension) {
1103 case 1 : {
1104 psDD = psav + Degree;
1105 psDDmDim = psDD - 1;
1106
1107 for (i = 0; i < Degree; i++) {
1108 DDmi--;
1109 pole = psDD;
1110 tbis = psDDmDim;
1111 jDmi = DDmi;
1112
1113 for (j = Degm1; j >= i; j--) {
1114 jDmi--;
1115 *pole -= *tbis; *pole /= (knot[jDmi] - knot[j]);
1116 pole--;
1117 tbis--;
1118 }
1119 }
1120 // Second phase, dependant of U
1121 iDim = - 1;
1122
1123 for (i = 0; i < Degree; i++) {
1124 iDim += 1;
1125 pole = psav + iDim;
1126 tbis = pole + 1;
1127 coef = U - knot[i];
1128
1129 for (j = i; j >= 0; j--) {
1130 *pole += coef * (*tbis);
1131 pole--;
1132 tbis--;
1133 }
1134 }
1135 // multiply by the degrees
1136 coef = Degree;
1137 Dmi = Degree;
1138 pole = psav + 1;
1139
1140 for (i = 1; i <= min; i++) {
1141 *pole *= coef; pole++;
1142 Dmi--;
1143 coef *= Dmi;
1144 }
1145 break;
1146 }
1147 case 2 : {
1148 psDD = psav + (Degree << 1);
1149 psDDmDim = psDD - 2;
1150
1151 for (i = 0; i < Degree; i++) {
1152 DDmi--;
1153 pole = psDD;
1154 tbis = psDDmDim;
1155 jDmi = DDmi;
1156
1157 for (j = Degm1; j >= i; j--) {
1158 jDmi--;
1159 coef = 1. / (knot[jDmi] - knot[j]);
1160 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1161 *pole -= *tbis; *pole *= coef;
1162 pole -= 3;
1163 tbis -= 3;
1164 }
1165 }
1166 // Second phase, dependant of U
1167 iDim = - 2;
1168
1169 for (i = 0; i < Degree; i++) {
1170 iDim += 2;
1171 pole = psav + iDim;
1172 tbis = pole + 2;
1173 coef = U - knot[i];
1174
1175 for (j = i; j >= 0; j--) {
1176 *pole += coef * (*tbis); pole++; tbis++;
1177 *pole += coef * (*tbis);
1178 pole -= 3;
1179 tbis -= 3;
1180 }
1181 }
1182 // multiply by the degrees
1183 coef = Degree;
1184 Dmi = Degree;
1185 pole = psav + 2;
1186
1187 for (i = 1; i <= min; i++) {
1188 *pole *= coef; pole++;
1189 *pole *= coef; pole++;
1190 Dmi--;
1191 coef *= Dmi;
1192 }
1193 break;
1194 }
1195 case 3 : {
1196 psDD = psav + (Degree << 1) + Degree;
1197 psDDmDim = psDD - 3;
1198
1199 for (i = 0; i < Degree; i++) {
1200 DDmi--;
1201 pole = psDD;
1202 tbis = psDDmDim;
1203 jDmi = DDmi;
1204
1205 for (j = Degm1; j >= i; j--) {
1206 jDmi--;
1207 coef = 1. / (knot[jDmi] - knot[j]);
1208 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1209 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1210 *pole -= *tbis; *pole *= coef;
1211 pole -= 5;
1212 tbis -= 5;
1213 }
1214 }
1215 // Second phase, dependant of U
1216 iDim = - 3;
1217
1218 for (i = 0; i < Degree; i++) {
1219 iDim += 3;
1220 pole = psav + iDim;
1221 tbis = pole + 3;
1222 coef = U - knot[i];
1223
1224 for (j = i; j >= 0; j--) {
1225 *pole += coef * (*tbis); pole++; tbis++;
1226 *pole += coef * (*tbis); pole++; tbis++;
1227 *pole += coef * (*tbis);
1228 pole -= 5;
1229 tbis -= 5;
1230 }
1231 }
1232 // multiply by the degrees
1233 coef = Degree;
1234 Dmi = Degree;
1235 pole = psav + 3;
1236
1237 for (i = 1; i <= min; i++) {
1238 *pole *= coef; pole++;
1239 *pole *= coef; pole++;
1240 *pole *= coef; pole++;
1241 Dmi--;
1242 coef *= Dmi;
1243 }
1244 break;
1245 }
1246 case 4 : {
1247 psDD = psav + (Degree << 2);
1248 psDDmDim = psDD - 4;
1249
1250 for (i = 0; i < Degree; i++) {
1251 DDmi--;
1252 pole = psDD;
1253 tbis = psDDmDim;
1254 jDmi = DDmi;
1255
1256 for (j = Degm1; j >= i; j--) {
1257 jDmi--;
1258 coef = 1. / (knot[jDmi] - knot[j]);
1259 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1260 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1261 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1262 *pole -= *tbis; *pole *= coef;
1263 pole -= 7;
1264 tbis -= 7;
1265 }
1266 }
1267 // Second phase, dependant of U
1268 iDim = - 4;
1269
1270 for (i = 0; i < Degree; i++) {
1271 iDim += 4;
1272 pole = psav + iDim;
1273 tbis = pole + 4;
1274 coef = U - knot[i];
1275
1276 for (j = i; j >= 0; j--) {
1277 *pole += coef * (*tbis); pole++; tbis++;
1278 *pole += coef * (*tbis); pole++; tbis++;
1279 *pole += coef * (*tbis); pole++; tbis++;
1280 *pole += coef * (*tbis);
1281 pole -= 7;
1282 tbis -= 7;
1283 }
1284 }
1285 // multiply by the degrees
1286 coef = Degree;
1287 Dmi = Degree;
1288 pole = psav + 4;
1289
1290 for (i = 1; i <= min; i++) {
1291 *pole *= coef; pole++;
1292 *pole *= coef; pole++;
1293 *pole *= coef; pole++;
1294 *pole *= coef; pole++;
1295 Dmi--;
1296 coef *= Dmi;
1297 }
1298 break;
1299 }
1300 default : {
1301 Standard_Integer k;
1302 Standard_Integer Dim2 = Dimension << 1;
1303 psDD = psav + Degree * Dimension;
1304 psDDmDim = psDD - Dimension;
1305
1306 for (i = 0; i < Degree; i++) {
1307 DDmi--;
1308 pole = psDD;
1309 tbis = psDDmDim;
1310 jDmi = DDmi;
1311
1312 for (j = Degm1; j >= i; j--) {
1313 jDmi--;
1314 coef = 1. / (knot[jDmi] - knot[j]);
1315
1316 for (k = 0; k < Dimension; k++) {
1317 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1318 }
1319 pole -= Dim2;
1320 tbis -= Dim2;
1321 }
1322 }
1323 // Second phase, dependant of U
1324 iDim = - Dimension;
1325
1326 for (i = 0; i < Degree; i++) {
1327 iDim += Dimension;
1328 pole = psav + iDim;
1329 tbis = pole + Dimension;
1330 coef = U - knot[i];
1331
1332 for (j = i; j >= 0; j--) {
1333
1334 for (k = 0; k < Dimension; k++) {
1335 *pole += coef * (*tbis); pole++; tbis++;
1336 }
1337 pole -= Dim2;
1338 tbis -= Dim2;
1339 }
1340 }
1341 // multiply by the degrees
1342 coef = Degree;
1343 Dmi = Degree;
1344 pole = psav + Dimension;
1345
1346 for (i = 1; i <= min; i++) {
1347
1348 for (k = 0; k < Dimension; k++) {
1349 *pole *= coef; pole++;
1350 }
1351 Dmi--;
1352 coef *= Dmi;
1353 }
1354 }
1355 }
1356}
1357
1358//=======================================================================
1359//function : BuildKnots
1360//purpose :
1361//=======================================================================
1362
1363void BSplCLib::BuildKnots(const Standard_Integer Degree,
1364 const Standard_Integer Index,
1365 const Standard_Boolean Periodic,
1366 const TColStd_Array1OfReal& Knots,
1367 const TColStd_Array1OfInteger& Mults,
1368 Standard_Real& LK)
1369{
1370 Standard_Integer KLower = Knots.Lower();
1371 const Standard_Real * pkn = &Knots(KLower);
1372 pkn -= KLower;
1373 Standard_Real *knot = &LK;
1374 if (&Mults == NULL) {
1375 switch (Degree) {
1376 case 1 : {
1377 Standard_Integer j = Index ;
1378 knot[0] = pkn[j]; j++;
1379 knot[1] = pkn[j];
1380 break;
1381 }
1382 case 2 : {
1383 Standard_Integer j = Index - 1;
1384 knot[0] = pkn[j]; j++;
1385 knot[1] = pkn[j]; j++;
1386 knot[2] = pkn[j]; j++;
1387 knot[3] = pkn[j];
1388 break;
1389 }
1390 case 3 : {
1391 Standard_Integer j = Index - 2;
1392 knot[0] = pkn[j]; j++;
1393 knot[1] = pkn[j]; j++;
1394 knot[2] = pkn[j]; j++;
1395 knot[3] = pkn[j]; j++;
1396 knot[4] = pkn[j]; j++;
1397 knot[5] = pkn[j];
1398 break;
1399 }
1400 case 4 : {
1401 Standard_Integer j = Index - 3;
1402 knot[0] = pkn[j]; j++;
1403 knot[1] = pkn[j]; j++;
1404 knot[2] = pkn[j]; j++;
1405 knot[3] = pkn[j]; j++;
1406 knot[4] = pkn[j]; j++;
1407 knot[5] = pkn[j]; j++;
1408 knot[6] = pkn[j]; j++;
1409 knot[7] = pkn[j];
1410 break;
1411 }
1412 case 5 : {
1413 Standard_Integer j = Index - 4;
1414 knot[0] = pkn[j]; j++;
1415 knot[1] = pkn[j]; j++;
1416 knot[2] = pkn[j]; j++;
1417 knot[3] = pkn[j]; j++;
1418 knot[4] = pkn[j]; j++;
1419 knot[5] = pkn[j]; j++;
1420 knot[6] = pkn[j]; j++;
1421 knot[7] = pkn[j]; j++;
1422 knot[8] = pkn[j]; j++;
1423 knot[9] = pkn[j];
1424 break;
1425 }
1426 case 6 : {
1427 Standard_Integer j = Index - 5;
1428 knot[ 0] = pkn[j]; j++;
1429 knot[ 1] = pkn[j]; j++;
1430 knot[ 2] = pkn[j]; j++;
1431 knot[ 3] = pkn[j]; j++;
1432 knot[ 4] = pkn[j]; j++;
1433 knot[ 5] = pkn[j]; j++;
1434 knot[ 6] = pkn[j]; j++;
1435 knot[ 7] = pkn[j]; j++;
1436 knot[ 8] = pkn[j]; j++;
1437 knot[ 9] = pkn[j]; j++;
1438 knot[10] = pkn[j]; j++;
1439 knot[11] = pkn[j];
1440 break;
1441 }
1442 default : {
1443 Standard_Integer i,j;
1444 Standard_Integer Deg2 = Degree << 1;
1445 j = Index - Degree;
1446
1447 for (i = 0; i < Deg2; i++) {
1448 j++;
1449 knot[i] = pkn[j];
1450 }
1451 }
1452 }
1453 }
1454 else {
1455 Standard_Integer i;
1456 Standard_Integer Deg1 = Degree - 1;
1457 Standard_Integer KUpper = Knots.Upper();
1458 Standard_Integer MLower = Mults.Lower();
1459 Standard_Integer MUpper = Mults.Upper();
1460 const Standard_Integer * pmu = &Mults(MLower);
1461 pmu -= MLower;
1462 Standard_Real dknot = 0;
1463 Standard_Integer ilow = Index , mlow = 0;
1464 Standard_Integer iupp = Index + 1, mupp = 0;
1465 Standard_Real loffset = 0., uoffset = 0.;
1466 Standard_Boolean getlow = Standard_True, getupp = Standard_True;
1467 if (Periodic) {
1468 dknot = pkn[KUpper] - pkn[KLower];
1469 if (iupp > MUpper) {
1470 iupp = MLower + 1;
1471 uoffset = dknot;
1472 }
1473 }
1474 // Find the knots around Index
1475
1476 for (i = 0; i < Degree; i++) {
1477 if (getlow) {
1478 mlow++;
1479 if (mlow > pmu[ilow]) {
1480 mlow = 1;
1481 ilow--;
1482 getlow = (ilow >= MLower);
1483 if (Periodic && !getlow) {
1484 ilow = MUpper - 1;
1485 loffset = dknot;
1486 getlow = Standard_True;
1487 }
1488 }
1489 if (getlow)
1490 knot[Deg1 - i] = pkn[ilow] - loffset;
1491 }
1492 if (getupp) {
1493 mupp++;
1494 if (mupp > pmu[iupp]) {
1495 mupp = 1;
1496 iupp++;
1497 getupp = (iupp <= MUpper);
1498 if (Periodic && !getupp) {
1499 iupp = MLower + 1;
1500 uoffset = dknot;
1501 getupp = Standard_True;
1502 }
1503 }
1504 if (getupp)
1505 knot[Degree + i] = pkn[iupp] + uoffset;
1506 }
1507 }
1508 }
1509}
1510
1511//=======================================================================
1512//function : PoleIndex
1513//purpose :
1514//=======================================================================
1515
1516Standard_Integer BSplCLib::PoleIndex(const Standard_Integer Degree,
1517 const Standard_Integer Index,
1518 const Standard_Boolean Periodic,
1519 const TColStd_Array1OfInteger& Mults)
1520{
1521 Standard_Integer i, pindex = 0;
1522
1523 for (i = Mults.Lower(); i <= Index; i++)
1524 pindex += Mults(i);
1525 if (Periodic)
1526 pindex -= Mults(Mults.Lower());
1527 else
1528 pindex -= Degree + 1;
1529
1530 return pindex;
1531}
1532
1533//=======================================================================
1534//function : BuildBoor
1535//purpose : builds the local array for boor
1536//=======================================================================
1537
1538void BSplCLib::BuildBoor(const Standard_Integer Index,
1539 const Standard_Integer Length,
1540 const Standard_Integer Dimension,
1541 const TColStd_Array1OfReal& Poles,
1542 Standard_Real& LP)
1543{
1544 Standard_Real *poles = &LP;
1545 Standard_Integer i,k, ip = Poles.Lower() + Index * Dimension;
1546
1547 for (i = 0; i < Length+1; i++) {
1548
1549 for (k = 0; k < Dimension; k++) {
1550 poles[k] = Poles(ip);
1551 ip++;
1552 if (ip > Poles.Upper()) ip = Poles.Lower();
1553 }
1554 poles += 2 * Dimension;
1555 }
1556}
1557
1558//=======================================================================
1559//function : BoorIndex
1560//purpose :
1561//=======================================================================
1562
1563Standard_Integer BSplCLib::BoorIndex(const Standard_Integer Index,
1564 const Standard_Integer Length,
1565 const Standard_Integer Depth)
1566{
1567 if (Index <= Depth) return Index;
1568 if (Index <= Length) return 2 * Index - Depth;
1569 return Length + Index - Depth;
1570}
1571
1572//=======================================================================
1573//function : GetPole
1574//purpose :
1575//=======================================================================
1576
1577void BSplCLib::GetPole(const Standard_Integer Index,
1578 const Standard_Integer Length,
1579 const Standard_Integer Depth,
1580 const Standard_Integer Dimension,
1581 Standard_Real& LP,
1582 Standard_Integer& Position,
1583 TColStd_Array1OfReal& Pole)
1584{
1585 Standard_Integer k;
1586 Standard_Real* pole = &LP + BoorIndex(Index,Length,Depth) * Dimension;
1587
1588 for (k = 0; k < Dimension; k++) {
1589 Pole(Position) = pole[k];
1590 Position++;
1591 }
1592 if (Position > Pole.Upper()) Position = Pole.Lower();
1593}
1594
1595//=======================================================================
1596//function : PrepareInsertKnots
1597//purpose :
1598//=======================================================================
1599
1600Standard_Boolean BSplCLib::PrepareInsertKnots
1601(const Standard_Integer Degree,
1602 const Standard_Boolean Periodic,
1603 const TColStd_Array1OfReal& Knots,
1604 const TColStd_Array1OfInteger& Mults,
1605 const TColStd_Array1OfReal& AddKnots,
1606 const TColStd_Array1OfInteger& AddMults,
1607 Standard_Integer& NbPoles,
1608 Standard_Integer& NbKnots,
1609 const Standard_Real Tolerance,
1610 const Standard_Boolean Add)
1611{
1612 Standard_Boolean addflat = &AddMults == NULL;
1613
1614 Standard_Integer first,last;
1615 if (Periodic) {
1616 first = Knots.Lower();
1617 last = Knots.Upper();
1618 }
1619 else {
1620 first = FirstUKnotIndex(Degree,Mults);
1621 last = LastUKnotIndex(Degree,Mults);
1622 }
1623 Standard_Real adeltaK1 = Knots(first)-AddKnots(AddKnots.Lower());
1624 Standard_Real adeltaK2 = AddKnots(AddKnots.Upper())-Knots(last);
1625 if (adeltaK1 > Tolerance) return Standard_False;
1626 if (adeltaK2 > Tolerance) return Standard_False;
1627
1628 Standard_Integer sigma = 0, mult, amult, lastmult = 0;
1629 NbKnots = 0;
1630 Standard_Integer k = Knots.Lower() - 1;
1631 Standard_Integer ak = AddKnots.Lower();
1632
1633 if(Periodic && AddKnots.Length() > 1)
1634 {
1635 //gka for case when segments was produced on full period only one knot
1636 //was added in the end of curve
1637 if(fabs(adeltaK1) <= Precision::PConfusion() &&
1638 fabs(adeltaK2) <= Precision::PConfusion())
1639 ak++;
1640 }
1641
1642 Standard_Real au,oldau = AddKnots(ak),Eps;
1643
1644 while (ak <= AddKnots.Upper()) {
1645 au = AddKnots(ak);
1646 if (au < oldau) return Standard_False;
1647 oldau = au;
1648
1649 Eps = Max(Tolerance,Epsilon(au));
1650
1651 while ((k < Knots.Upper()) && (Knots(k+1) - au <= Eps)) {
1652 k++;
1653 NbKnots++;
1654 sigma += Mults(k);
1655 }
1656
1657 if (addflat) amult = 1;
1658 else amult = Max(0,AddMults(ak));
1659
1660 while ((ak < AddKnots.Upper()) &&
1661 (Abs(au - AddKnots(ak+1)) <= Eps)) {
1662 ak++;
1663 if (Add) {
1664 if (addflat) amult++;
1665 else amult += Max(0,AddMults(ak));
1666 }
1667 }
1668
1669
1670 if (Abs(au - Knots(k)) <= Eps) {
1671 // identic to existing knot
1672 mult = Mults(k);
1673 lastmult = mult;//gka
1674 if (Add) {
1675 if (mult + amult > Degree)
1676 amult = Max(0,Degree - mult);
1677 sigma += amult;
1678 //lastmult = mult + amult;
1679
1680 }
1681 else if (amult > mult) {
1682 if (amult > Degree) amult = Degree;
1683 sigma += amult - mult;
1684 //lastmult = amult;//gka modified
1685 }
1686 /*
1687 // on periodic curves if this is the last knot
1688 // the multiplicity is added twice to account for the first knot
1689 if (k == Knots.Upper() && Periodic) {
1690 if (Add)
1691 sigma += amult;
1692 else
1693 sigma += amult - mult;
1694 }
1695 */
1696 }
1697 else {
1698 // not identic to existing knot
1699 if (amult > 0) {
1700 if (amult > Degree) amult = Degree;
1701 NbKnots++;
1702 //lastmult = amult;
1703 sigma += amult;
1704 }
1705 }
1706
1707 ak++;
1708 }
1709
1710 // count the last knots
1711 if (lastmult == 0)// || k < Knots.Upper())
1712 lastmult = Mults(Knots.Upper());
1713
1714 while (k < Knots.Upper()) {
1715 k++;
1716 NbKnots++;
1717 sigma += Mults(k);
1718 }
1719
1720 if (Periodic) {
1721 NbPoles = sigma - lastmult;
1722 }
1723 else {
1724 NbPoles = sigma - Degree - 1;
1725 }
1726
1727 return Standard_True;
1728}
1729
1730//=======================================================================
1731//function : Copy
1732//purpose : copy reals from an array to an other
1733//
1734// NbValues are copied from OldPoles(OldFirst)
1735// to NewPoles(NewFirst)
1736//
1737// Periodicity is handled.
1738// OldFirst and NewFirst are updated
1739// to the position after the last copied pole.
1740//
1741//=======================================================================
1742
1743static void Copy(const Standard_Integer NbPoles,
1744 Standard_Integer& OldFirst,
1745 const TColStd_Array1OfReal& OldPoles,
1746 Standard_Integer& NewFirst,
1747 TColStd_Array1OfReal& NewPoles)
1748{
1749 // reset the index in the range for periodicity
1750
1751 OldFirst = OldPoles.Lower() +
1752 (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1);
1753
1754 NewFirst = NewPoles.Lower() +
1755 (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1);
1756
1757 // copy
1758 Standard_Integer i;
1759
1760 for (i = 1; i <= NbPoles; i++) {
1761 NewPoles(NewFirst) = OldPoles(OldFirst);
1762 OldFirst++;
1763 if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower();
1764 NewFirst++;
1765 if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower();
1766 }
1767}
1768
1769//=======================================================================
1770//function : InsertKnots
1771//purpose : insert an array of knots and multiplicities
1772//=======================================================================
1773
1774void BSplCLib::InsertKnots
1775(const Standard_Integer Degree,
1776 const Standard_Boolean Periodic,
1777 const Standard_Integer Dimension,
1778 const TColStd_Array1OfReal& Poles,
1779 const TColStd_Array1OfReal& Knots,
1780 const TColStd_Array1OfInteger& Mults,
1781 const TColStd_Array1OfReal& AddKnots,
1782 const TColStd_Array1OfInteger& AddMults,
1783 TColStd_Array1OfReal& NewPoles,
1784 TColStd_Array1OfReal& NewKnots,
1785 TColStd_Array1OfInteger& NewMults,
1786 const Standard_Real Tolerance,
1787 const Standard_Boolean Add)
1788{
1789 Standard_Boolean addflat = &AddMults == NULL;
1790
1791 Standard_Integer i,k,mult,firstmult;
1792 Standard_Integer index,kn,curnk,curk;
1793 Standard_Integer p,np, curp, curnp, length, depth;
1794 Standard_Real u;
1795 Standard_Integer need;
1796 Standard_Real Eps;
1797
1798 // -------------------
1799 // create local arrays
1800 // -------------------
1801
1802 Standard_Real *knots = new Standard_Real[2*Degree];
1803 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
1804
1805 //----------------------------
1806 // loop on the knots to insert
1807 //----------------------------
1808
1809 curk = Knots.Lower()-1; // current position in Knots
1810 curnk = NewKnots.Lower()-1; // current position in NewKnots
1811 curp = Poles.Lower(); // current position in Poles
1812 curnp = NewPoles.Lower(); // current position in NewPoles
1813
1814 // NewKnots, NewMults, NewPoles contains the current state of the curve
1815
1816 // index is the first pole of the current curve for insertion schema
1817
1818 if (Periodic) index = -Mults(Mults.Lower());
1819 else index = -Degree-1;
1820
1821 // on Periodic curves the first knot and the last knot are inserted later
1822 // (they are the same knot)
1823 firstmult = 0; // multiplicity of the first-last knot for periodic
1824
1825
1826 // kn current knot to insert in AddKnots
1827
1828 for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) {
1829
1830 u = AddKnots(kn);
1831 Eps = Max(Tolerance,Epsilon(u));
1832
1833 //-----------------------------------
1834 // find the position in the old knots
1835 // and copy to the new knots
1836 //-----------------------------------
1837
1838 while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) {
1839 curk++; curnk++;
1840 NewKnots(curnk) = Knots(curk);
1841 index += NewMults(curnk) = Mults(curk);
1842 }
1843
1844 //-----------------------------------
1845 // Slice the knots and the mults
1846 // to the current size of the new curve
1847 //-----------------------------------
1848
1849 i = curnk + Knots.Upper() - curk;
1850 TColStd_Array1OfReal nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i);
1851 TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i);
1852
1853 //-----------------------------------
1854 // copy enough knots
1855 // to compute the insertion schema
1856 //-----------------------------------
1857
1858 k = curk;
1859 i = curnk;
1860 mult = 0;
1861
1862 while (mult < Degree && k < Knots.Upper()) {
1863 k++; i++;
1864 nknots(i) = Knots(k);
1865 mult += nmults(i) = Mults(k);
1866 }
1867
1868 // copy knots at the end for periodic curve
1869 if (Periodic) {
1870 mult = 0;
1871 k = Knots.Upper();
1872 i = nknots.Upper();
1873
1874 while (mult < Degree && i > curnk) {
1875 nknots(i) = Knots(k);
1876 mult += nmults(i) = Mults(k);
1877 k--;
1878 i--;
1879 }
1880 nmults(nmults.Upper()) = nmults(nmults.Lower());
1881 }
1882
1883
1884
1885 //------------------------------------
1886 // do the boor scheme on the new curve
1887 // to insert the new knot
1888 //------------------------------------
1889
1890 Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps);
1891
1892 if (sameknot) length = Max(0,Degree - NewMults(curnk));
1893 else length = Degree;
1894
1895 if (addflat) depth = 1;
1896 else depth = Min(Degree,AddMults(kn));
1897
1898 if (sameknot) {
1899 if (Add) {
1900 if ((NewMults(curnk) + depth) > Degree)
1901 depth = Degree - NewMults(curnk);
1902 }
1903 else {
1904 depth = Max(0,depth-NewMults(curnk));
1905 }
1906
1907 if (Periodic) {
1908 // on periodic curve the first and last knot are delayed to the end
1909 if (curk == Knots.Lower() || (curk == Knots.Upper())) {
1910 firstmult += depth;
1911 depth = 0;
1912 }
1913 }
1914 }
1915 if (depth <= 0) continue;
1916
1917 BuildKnots(Degree,curnk,Periodic,nknots,nmults,*knots);
1918
1919 // copy the poles
1920
1921 need = NewPoles.Lower()+(index+length+1)*Dimension - curnp;
1922 need = Min(need,Poles.Upper() - curp + 1);
1923
1924 p = curp;
1925 np = curnp;
1926 Copy(need,p,Poles,np,NewPoles);
1927 curp += need;
1928 curnp += need;
1929
1930 // slice the poles to the current number of poles in case of periodic
1931 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1932
1933 BuildBoor(index,length,Dimension,npoles,*poles);
1934 BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length);
1935
1936 //-------------------
1937 // copy the new poles
1938 //-------------------
1939
1940 curnp += depth * Dimension; // number of poles is increased by depth
1941 TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1942 np = NewKnots.Lower()+(index+1)*Dimension;
1943
1944 for (i = 1; i <= length + depth; i++)
1945 GetPole(i,length,depth,Dimension,*poles,np,ThePoles);
1946
1947 //-------------------
1948 // insert the knot
1949 //-------------------
1950
1951 index += depth;
1952 if (sameknot) {
1953 NewMults(curnk) += depth;
1954 }
1955 else {
1956 curnk++;
1957 NewKnots(curnk) = u;
1958 NewMults(curnk) = depth;
1959 }
1960 }
1961
1962 //------------------------------
1963 // copy the last poles and knots
1964 //------------------------------
1965
1966 Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles);
1967
1968 while (curk < Knots.Upper()) {
1969 curk++; curnk++;
1970 NewKnots(curnk) = Knots(curk);
1971 NewMults(curnk) = Mults(curk);
1972 }
1973
1974 //------------------------------
1975 // process the first-last knot
1976 // on periodic curves
1977 //------------------------------
1978
1979 if (firstmult > 0) {
1980 curnk = NewKnots.Lower();
1981 if (NewMults(curnk) + firstmult > Degree) {
1982 firstmult = Degree - NewMults(curnk);
1983 }
1984 if (firstmult > 0) {
1985
1986 length = Degree - NewMults(curnk);
1987 depth = firstmult;
1988
1989 BuildKnots(Degree,curnk,Periodic,NewKnots,NewMults,*knots);
1990 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),
1991 NewPoles.Lower(),
1992 NewPoles.Upper()-depth*Dimension);
1993 BuildBoor(0,length,Dimension,npoles,*poles);
1994 BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length);
1995
1996 //---------------------------
1997 // copy the new poles
1998 // but rotate them with depth
1999 //---------------------------
2000
2001 np = NewPoles.Lower();
2002
2003 for (i = depth; i < length + depth; i++)
2004 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2005
2006 np = NewPoles.Upper() - depth*Dimension + 1;
2007
2008 for (i = 0; i < depth; i++)
2009 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2010
2011 NewMults(NewMults.Lower()) += depth;
2012 NewMults(NewMults.Upper()) += depth;
2013 }
2014 }
2015 // free local arrays
2016 delete [] knots;
2017 delete [] poles;
2018}
2019
2020//=======================================================================
2021//function : RemoveKnot
2022//purpose :
2023//=======================================================================
2024
2025Standard_Boolean BSplCLib::RemoveKnot
2026(const Standard_Integer Index,
2027 const Standard_Integer Mult,
2028 const Standard_Integer Degree,
2029 const Standard_Boolean Periodic,
2030 const Standard_Integer Dimension,
2031 const TColStd_Array1OfReal& Poles,
2032 const TColStd_Array1OfReal& Knots,
2033 const TColStd_Array1OfInteger& Mults,
2034 TColStd_Array1OfReal& NewPoles,
2035 TColStd_Array1OfReal& NewKnots,
2036 TColStd_Array1OfInteger& NewMults,
2037 const Standard_Real Tolerance)
2038{
2039 Standard_Integer index,i,j,k,p,np;
2040
2041 Standard_Integer TheIndex = Index;
2042
2043 // protection
2044 Standard_Integer first,last;
2045 if (Periodic) {
2046 first = Knots.Lower();
2047 last = Knots.Upper();
2048 }
2049 else {
2050 first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1;
2051 last = BSplCLib::LastUKnotIndex(Degree,Mults) - 1;
2052 }
2053 if (Index < first) return Standard_False;
2054 if (Index > last) return Standard_False;
2055
2056 if ( Periodic && (Index == first)) TheIndex = last;
2057
2058 Standard_Integer depth = Mults(TheIndex) - Mult;
2059 Standard_Integer length = Degree - Mult;
2060
2061 // -------------------
2062 // create local arrays
2063 // -------------------
2064
2065 Standard_Real *knots = new Standard_Real[4*Degree];
2066 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
2067
2068
2069 // ------------------------------------
2070 // build the knots for anti Boor Scheme
2071 // ------------------------------------
2072
2073 // the new sequence of knots
2074 // is obtained from the knots at Index-1 and Index
2075
2076 BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,Mults,*knots);
2077 index = PoleIndex(Degree,TheIndex-1,Periodic,Mults);
2078 BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,Mults,knots[2*Degree]);
2079
2080 index += Mult;
2081
2082 for (i = 0; i < Degree - Mult; i++)
2083 knots[i] = knots[i+Mult];
2084
2085 for (i = Degree-Mult; i < 2*Degree; i++)
2086 knots[i] = knots[2*Degree+i];
2087
2088
2089 // ------------------------------------
2090 // build the poles for anti Boor Scheme
2091 // ------------------------------------
2092
2093 p = Poles.Lower()+index * Dimension;
2094
2095 for (i = 0; i <= length + depth; i++) {
2096 j = Dimension * BoorIndex(i,length,depth);
2097
2098 for (k = 0; k < Dimension; k++) {
2099 poles[j+k] = Poles(p+k);
2100 }
2101 p += Dimension;
2102 if (p > Poles.Upper()) p = Poles.Lower();
2103 }
2104
2105
2106 // ----------------
2107 // Anti Boor Scheme
2108 // ----------------
2109
2110 Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots,
2111 Dimension,*poles,
2112 depth,length,Tolerance);
2113
2114 // ----------------
2115 // copy the results
2116 // ----------------
2117
2118 if (result) {
2119
2120 // poles
2121
2122 p = Poles.Lower();
2123 np = NewPoles.Lower();
2124
2125 // unmodified poles before
2126 Copy((index+1)*Dimension,p,Poles,np,NewPoles);
2127
2128
2129 // modified
2130
2131 for (i = 1; i <= length; i++)
2132 BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles);
2133 p += (length + depth) * Dimension ;
2134
2135 // unmodified poles after
2136 if (p != Poles.Lower()) {
2137 i = Poles.Upper() - p + 1;
2138 Copy(i,p,Poles,np,NewPoles);
2139 }
2140
2141 // knots and mults
2142
2143 if (Mult > 0) {
2144 NewKnots = Knots;
2145 NewMults = Mults;
2146 NewMults(TheIndex) = Mult;
2147 if (Periodic) {
2148 if (TheIndex == first) NewMults(last) = Mult;
2149 if (TheIndex == last) NewMults(first) = Mult;
2150 }
2151 }
2152 else {
2153 if (!Periodic || (TheIndex != first && TheIndex != last)) {
2154
2155 for (i = Knots.Lower(); i < TheIndex; i++) {
2156 NewKnots(i) = Knots(i);
2157 NewMults(i) = Mults(i);
2158 }
2159
2160 for (i = TheIndex+1; i <= Knots.Upper(); i++) {
2161 NewKnots(i-1) = Knots(i);
2162 NewMults(i-1) = Mults(i);
2163 }
2164 }
2165 else {
2166 // The interesting case of a Periodic curve
2167 // where the first and last knot is removed.
2168
2169 for (i = first; i < last-1; i++) {
2170 NewKnots(i) = Knots(i+1);
2171 NewMults(i) = Mults(i+1);
2172 }
2173 NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first);
2174 NewMults(last-1) = NewMults(first);
2175 }
2176 }
2177 }
2178
2179
2180 // free local arrays
2181 delete [] knots;
2182 delete [] poles;
2183
2184 return result;
2185}
2186
2187//=======================================================================
2188//function : IncreaseDegreeCountKnots
2189//purpose :
2190//=======================================================================
2191
2192Standard_Integer BSplCLib::IncreaseDegreeCountKnots
2193(const Standard_Integer Degree,
2194 const Standard_Integer NewDegree,
2195 const Standard_Boolean Periodic,
2196 const TColStd_Array1OfInteger& Mults)
2197{
2198 if (Periodic) return Mults.Length();
2199 Standard_Integer f = FirstUKnotIndex(Degree,Mults);
2200 Standard_Integer l = LastUKnotIndex(Degree,Mults);
2201 Standard_Integer m,i,removed = 0, step = NewDegree - Degree;
2202
2203 i = Mults.Lower();
2204 m = Degree + (f - i + 1) * step + 1;
2205
2206 while (m > NewDegree+1) {
2207 removed++;
2208 m -= Mults(i) + step;
2209 i++;
2210 }
2211 if (m < NewDegree+1) removed--;
2212
2213 i = Mults.Upper();
2214 m = Degree + (i - l + 1) * step + 1;
2215
2216 while (m > NewDegree+1) {
2217 removed++;
2218 m -= Mults(i) + step;
2219 i--;
2220 }
2221 if (m < NewDegree+1) removed--;
2222
2223 return Mults.Length() - removed;
2224}
2225
2226//=======================================================================
2227//function : IncreaseDegree
2228//purpose :
2229//=======================================================================
2230
2231void BSplCLib::IncreaseDegree
2232(const Standard_Integer Degree,
2233 const Standard_Integer NewDegree,
2234 const Standard_Boolean Periodic,
2235 const Standard_Integer Dimension,
2236 const TColStd_Array1OfReal& Poles,
2237 const TColStd_Array1OfReal& Knots,
2238 const TColStd_Array1OfInteger& Mults,
2239 TColStd_Array1OfReal& NewPoles,
2240 TColStd_Array1OfReal& NewKnots,
2241 TColStd_Array1OfInteger& NewMults)
2242{
2243 // Degree elevation of a BSpline Curve
2244
2245 // This algorithms loops on degree incrementation from Degree to NewDegree.
2246 // The variable curDeg is the current degree to increment.
2247
2248 // Before degree incrementations a "working curve" is created.
2249 // It has the same knot, poles and multiplicities.
2250
2251 // If the curve is periodic knots are added on the working curve before
2252 // and after the existing knots to make it a non-periodic curves.
2253 // The poles are also copied.
2254
2255 // The first and last multiplicity of the working curve are set to Degree+1,
2256 // null poles are added if necessary.
2257
2258 // Then the degree is incremented on the working curve.
2259 // The knots are unchanged but all multiplicities will be incremented.
2260
2261 // Each degree incrementation is achieved by averaging curDeg+1 curves.
2262
2263 // See : Degree elevation of B-spline curves
2264 // Hartmut PRAUTZSCH
2265 // CAGD 1 (1984)
2266
2267
2268 //-------------------------
2269 // create the working curve
2270 //-------------------------
2271
2272 Standard_Integer i,k,f,l,m,pf,pl,firstknot;
2273
2274 pf = 0; // number of null poles added at beginning
2275 pl = 0; // number of null poles added at end
2276
2277 Standard_Integer nbwknots = Knots.Length();
2278 f = FirstUKnotIndex(Degree,Mults);
2279 l = LastUKnotIndex (Degree,Mults);
2280
2281 if (Periodic) {
2282 // Periodic curves are transformed in non-periodic curves
2283
2284 nbwknots += f - Mults.Lower();
2285
2286 pf = -Degree - 1;
2287
2288 for (i = Mults.Lower(); i <= f; i++)
2289 pf += Mults(i);
2290
2291 nbwknots += Mults.Upper() - l;
2292
2293 pl = -Degree - 1;
2294
2295 for (i = l; i <= Mults.Upper(); i++)
2296 pl += Mults(i);
2297 }
2298
2299 // copy the knots and multiplicities
2300 TColStd_Array1OfReal wknots(1,nbwknots);
2301 TColStd_Array1OfInteger wmults(1,nbwknots);
2302 if (!Periodic) {
2303 wknots = Knots;
2304 wmults = Mults;
2305 }
2306 else {
2307 // copy the knots for a periodic curve
2308 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2309 i = 0;
2310
2311 for (k = l; k < Knots.Upper(); k++) {
2312 i++;
2313 wknots(i) = Knots(k) - period;
2314 wmults(i) = Mults(k);
2315 }
2316
2317 for (k = Knots.Lower(); k <= Knots.Upper(); k++) {
2318 i++;
2319 wknots(i) = Knots(k);
2320 wmults(i) = Mults(k);
2321 }
2322
2323 for (k = Knots.Lower()+1; k <= f; k++) {
2324 i++;
2325 wknots(i) = Knots(k) + period;
2326 wmults(i) = Mults(k);
2327 }
2328 }
2329
2330 // set the first and last mults to Degree+1
2331 // and add null poles
2332
2333 pf += Degree + 1 - wmults(1);
2334 wmults(1) = Degree + 1;
2335 pl += Degree + 1 - wmults(nbwknots);
2336 wmults(nbwknots) = Degree + 1;
2337
2338 //---------------------------
2339 // poles of the working curve
2340 //---------------------------
2341
2342 Standard_Integer nbwpoles = 0;
2343
2344 for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i);
2345 nbwpoles -= Degree + 1;
2346
2347 // we provide space for degree elevation
2348 TColStd_Array1OfReal
2349 wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension);
2350
2351 for (i = 1; i <= pf * Dimension; i++)
2352 wpoles(i) = 0;
2353
2354 k = Poles.Lower();
2355
2356 for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) {
2357 wpoles(i) = Poles(k);
2358 k++;
2359 if (k > Poles.Upper()) k = Poles.Lower();
2360 }
2361
2362 for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++)
2363 wpoles(i) = 0;
2364
2365
2366 //-----------------------------------------------------------
2367 // Declare the temporary arrays used in degree incrementation
2368 //-----------------------------------------------------------
2369
2370 Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree);
2371 // Arrays for storing the temporary curves
2372 TColStd_Array1OfReal tempc1(1,nbwp * Dimension);
2373 TColStd_Array1OfReal tempc2(1,nbwp * Dimension);
2374
2375 // Array for storing the knots to insert
2376 TColStd_Array1OfReal iknots(1,nbwknots);
2377
2378 // Arrays for receiving the knots after insertion
2379 TColStd_Array1OfReal nknots(1,nbwknots);
2380
2381
2382
2383 //------------------------------
2384 // Loop on degree incrementation
2385 //------------------------------
2386
2387 Standard_Integer step,curDeg;
2388 Standard_Integer nbp = nbwpoles;
2389 nbwp = nbp;
2390
2391 for (curDeg = Degree; curDeg < NewDegree; curDeg++) {
2392
2393 nbp = nbwp; // current number of poles
2394 nbwp = nbp + nbwknots - 1; // new number of poles
2395
2396 // For the averaging
2397 TColStd_Array1OfReal nwpoles(1,nbwp * Dimension);
2398 nwpoles.Init(0.0e0) ;
2399
2400
2401 for (step = 0; step <= curDeg; step++) {
2402
2403 // Compute the bspline of rank step.
2404
2405 // if not the first time, decrement the multiplicities back
2406 if (step != 0) {
2407 for (i = 1; i <= nbwknots; i++)
2408 wmults(i)--;
2409 }
2410
2411 // Poles are the current poles
2412 // but the poles congruent to step are duplicated.
2413
2414 Standard_Integer offset = 0;
2415
2416 for (i = 0; i < nbp; i++) {
2417 offset++;
2418
2419 for (k = 0; k < Dimension; k++) {
2420 tempc1((offset-1)*Dimension+k+1) =
2421 wpoles(NewPoles.Lower()+i*Dimension+k);
2422 }
2423 if (i % (curDeg+1) == step) {
2424 offset++;
2425
2426 for (k = 0; k < Dimension; k++) {
2427 tempc1((offset-1)*Dimension+k+1) =
2428 wpoles(NewPoles.Lower()+i*Dimension+k);
2429 }
2430 }
2431 }
2432
2433 // Knots multiplicities are increased
2434 // For each knot where the sum of multiplicities is congruent to step
2435
2436 Standard_Integer stepmult = step+1;
2437 Standard_Integer nbknots = 0;
2438 Standard_Integer smult = 0;
2439
2440 for (k = 1; k <= nbwknots; k++) {
2441 smult += wmults(k);
2442 if (smult >= stepmult) {
2443 // this knot is increased
2444 stepmult += curDeg+1;
2445 wmults(k)++;
2446 }
2447 else {
2448 // this knot is inserted
2449 nbknots++;
2450 iknots(nbknots) = wknots(k);
2451 }
2452 }
2453
2454 // the curve is obtained by inserting the knots
2455 // to raise the multiplicities
2456
2457 // we build "slices" of the arrays to set the correct size
2458 if (nbknots > 0) {
2459 TColStd_Array1OfReal aknots(iknots(1),1,nbknots);
2460 TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension);
2461 TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp * Dimension);
2462// InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2463// aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.));
2464
2465 InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2466 aknots,NoMults(),ncurve,nknots,wmults,0.0);
2467
2468 // add to the average
2469
2470 for (i = 1; i <= nbwp * Dimension; i++)
2471 nwpoles(i) += ncurve(i);
2472 }
2473 else {
2474 // add to the average
2475
2476 for (i = 1; i <= nbwp * Dimension; i++)
2477 nwpoles(i) += tempc1(i);
2478 }
2479 }
2480
2481 // The result is the average
2482
2483 for (i = 1; i <= nbwp * Dimension; i++) {
2484 wpoles(i) = nwpoles(i) / (curDeg+1);
2485 }
2486 }
2487
2488 //-----------------
2489 // Copy the results
2490 //-----------------
2491
2492 // index in new knots of the first knot of the curve
2493 if (Periodic)
2494 firstknot = Mults.Upper() - l + 1;
2495 else
2496 firstknot = f;
2497
2498 // the new curve starts at index firstknot
2499 // so we must remove knots until the sum of multiplicities
2500 // from the first to the start is NewDegree+1
2501
2502 // m is the current sum of multiplicities
2503 m = 0;
2504
2505 for (k = 1; k <= firstknot; k++)
2506 m += wmults(k);
2507
2508 // compute the new first knot (k), pf will be the index of the first pole
2509 k = 1;
2510 pf = 0;
2511
2512 while (m > NewDegree+1) {
2513 k++;
2514 m -= wmults(k);
2515 pf += wmults(k);
2516 }
2517 if (m < NewDegree+1) {
2518 k--;
2519 wmults(k) += m - NewDegree - 1;
2520 pf += m - NewDegree - 1;
2521 }
2522
2523 // on a periodic curve the knots start with firstknot
2524 if (Periodic)
2525 k = firstknot;
2526
2527 // copy knots
2528
2529 for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) {
2530 NewKnots(i) = wknots(k);
2531 NewMults(i) = wmults(k);
2532 k++;
2533 }
2534
2535 // copy poles
2536 pf *= Dimension;
2537
2538 for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) {
2539 pf++;
2540 NewPoles(i) = wpoles(pf);
2541 }
2542}
2543
2544//=======================================================================
2545//function : PrepareUnperiodize
2546//purpose :
2547//=======================================================================
2548
2549void BSplCLib::PrepareUnperiodize
2550(const Standard_Integer Degree,
2551 const TColStd_Array1OfInteger& Mults,
2552 Standard_Integer& NbKnots,
2553 Standard_Integer& NbPoles)
2554{
2555 Standard_Integer i;
2556 // initialize NbKnots and NbPoles
2557 NbKnots = Mults.Length();
2558 NbPoles = - Degree - 1;
2559
2560 for (i = Mults.Lower(); i <= Mults.Upper(); i++)
2561 NbPoles += Mults(i);
2562
2563 Standard_Integer sigma, k;
2564 // Add knots at the beginning of the curve to raise Multiplicities
2565 // to Degre + 1;
2566 sigma = Mults(Mults.Lower());
2567 k = Mults.Upper() - 1;
2568
2569 while ( sigma < Degree + 1) {
2570 sigma += Mults(k);
2571 NbPoles += Mults(k);
2572 k--;
2573 NbKnots++;
2574 }
2575 // We must add exactly until Degree + 1 ->
2576 // Supress the excedent.
2577 if ( sigma > Degree + 1)
2578 NbPoles -= sigma - Degree - 1;
2579
2580 // Add knots at the end of the curve to raise Multiplicities
2581 // to Degre + 1;
2582 sigma = Mults(Mults.Upper());
2583 k = Mults.Lower() + 1;
2584
2585 while ( sigma < Degree + 1) {
2586 sigma += Mults(k);
2587 NbPoles += Mults(k);
2588 k++;
2589 NbKnots++;
2590 }
2591 // We must add exactly until Degree + 1 ->
2592 // Supress the excedent.
2593 if ( sigma > Degree + 1)
2594 NbPoles -= sigma - Degree - 1;
2595}
2596
2597//=======================================================================
2598//function : Unperiodize
2599//purpose :
2600//=======================================================================
2601
2602void BSplCLib::Unperiodize
2603(const Standard_Integer Degree,
2604 const Standard_Integer , // Dimension,
2605 const TColStd_Array1OfInteger& Mults,
2606 const TColStd_Array1OfReal& Knots,
2607 const TColStd_Array1OfReal& Poles,
2608 TColStd_Array1OfInteger& NewMults,
2609 TColStd_Array1OfReal& NewKnots,
2610 TColStd_Array1OfReal& NewPoles)
2611{
2612 Standard_Integer sigma, k, index = 0;
2613 // evaluation of index : number of knots to insert before knot(1) to
2614 // raise sum of multiplicities to <Degree + 1>
2615 sigma = Mults(Mults.Lower());
2616 k = Mults.Upper() - 1;
2617
2618 while ( sigma < Degree + 1) {
2619 sigma += Mults(k);
2620 k--;
2621 index++;
2622 }
2623
2624 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2625
2626 // set the 'interior' knots;
2627
2628 for ( k = 1; k <= Knots.Length(); k++) {
2629 NewKnots ( k + index ) = Knots( k);
2630 NewMults ( k + index ) = Mults( k);
2631 }
2632
2633 // set the 'starting' knots;
2634
2635 for ( k = 1; k <= index; k++) {
2636 NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period;
2637 NewMults( k) = NewMults( k + Knots.Length() - 1);
2638 }
2639 NewMults( 1) -= sigma - Degree -1;
2640
2641 // set the 'ending' knots;
2642 sigma = NewMults( index + Knots.Length() );
2643
2644 for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) {
2645 NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period;
2646 NewMults( k) = NewMults( k - Knots.Length() + 1);
2647 sigma += NewMults( k - Knots.Length() + 1);
2648 }
2649 NewMults(NewMults.Length()) -= sigma - Degree - 1;
2650
2651 for ( k = 1 ; k <= NewPoles.Length(); k++) {
2652 NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1);
2653 }
2654}
2655
2656//=======================================================================
2657//function : PrepareTrimming
2658//purpose :
2659//=======================================================================
2660
2661void BSplCLib::PrepareTrimming(const Standard_Integer Degree,
2662 const Standard_Boolean Periodic,
2663 const TColStd_Array1OfReal& Knots,
2664 const TColStd_Array1OfInteger& Mults,
2665 const Standard_Real U1,
2666 const Standard_Real U2,
2667 Standard_Integer& NbKnots,
2668 Standard_Integer& NbPoles)
2669{
2670 Standard_Integer i;
2671 Standard_Real NewU1, NewU2;
2672 Standard_Integer index1 = 0, index2 = 0;
2673
2674 // Eval index1, index2 : position of U1 and U2 in the Array Knots
2675 // such as Knots(index1-1) <= U1 < Knots(index1)
2676 // Knots(index2-1) <= U2 < Knots(index2)
2677 LocateParameter( Degree, Knots, Mults, U1, Periodic,
2678 Knots.Lower(), Knots.Upper(), index1, NewU1);
2679 LocateParameter( Degree, Knots, Mults, U2, Periodic,
2680 Knots.Lower(), Knots.Upper(), index2, NewU2);
2681 index1++;
2682 if ( Abs(Knots(index2) - U2) <= Epsilon( U1))
2683 index2--;
2684
2685 // eval NbKnots:
2686 NbKnots = index2 - index1 + 3;
2687
2688 // eval NbPoles:
2689 NbPoles = Degree + 1;
2690
2691 for ( i = index1; i <= index2; i++)
2692 NbPoles += Mults(i);
2693}
2694
2695//=======================================================================
2696//function : Trimming
2697//purpose :
2698//=======================================================================
2699
2700void BSplCLib::Trimming(const Standard_Integer Degree,
2701 const Standard_Boolean Periodic,
2702 const Standard_Integer Dimension,
2703 const TColStd_Array1OfReal& Knots,
2704 const TColStd_Array1OfInteger& Mults,
2705 const TColStd_Array1OfReal& Poles,
2706 const Standard_Real U1,
2707 const Standard_Real U2,
2708 TColStd_Array1OfReal& NewKnots,
2709 TColStd_Array1OfInteger& NewMults,
2710 TColStd_Array1OfReal& NewPoles)
2711{
2712 Standard_Integer i, nbpoles, nbknots;
2713 Standard_Real kk[2];
2714 Standard_Integer mm[2];
2715 TColStd_Array1OfReal K( kk[0], 1, 2 );
2716 TColStd_Array1OfInteger M( mm[0], 1, 2 );
2717
2718 K(1) = U1; K(2) = U2;
2719 mm[0] = mm[1] = Degree;
2720 if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, M,
2721 nbpoles, nbknots, Epsilon( U1), 0))
2722 Standard_OutOfRange::Raise();
2723
2724 TColStd_Array1OfReal TempPoles(1, nbpoles*Dimension);
2725 TColStd_Array1OfReal TempKnots(1, nbknots);
2726 TColStd_Array1OfInteger TempMults(1, nbknots);
2727
2728//
2729// do not allow the multiplicities to Add : they must be less than Degree
2730//
2731 InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults,
2732 K, M, TempPoles, TempKnots, TempMults, Epsilon(U1),
2733 Standard_False);
2734
2735 // find in TempPoles the index of the pole corresponding to U1
2736 Standard_Integer Kindex = 0, Pindex;
2737 Standard_Real NewU1;
2738 LocateParameter( Degree, TempKnots, TempMults, U1, Periodic,
2739 TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1);
2740 Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults);
2741 Pindex *= Dimension;
2742
2743 for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i);
2744
2745 for ( i = 1; i <= NewKnots.Length(); i++) {
2746 NewKnots(i) = TempKnots( Kindex+i-1);
2747 NewMults(i) = TempMults( Kindex+i-1);
2748 }
2749 NewMults(1) = Min(Degree, NewMults(1)) + 1 ;
2750 NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ;
2751}
2752
2753//=======================================================================
2754//function : Solves a LU factored Matrix
2755//purpose :
2756//=======================================================================
2757
2758Standard_Integer
2759BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2760 const Standard_Integer UpperBandWidth,
2761 const Standard_Integer LowerBandWidth,
2762 const Standard_Integer ArrayDimension,
2763 Standard_Real& Array)
2764{
2765 Standard_Integer ii,
2766 jj,
2767 kk,
2768 MinIndex,
2769 MaxIndex,
2770 ReturnCode = 0 ;
2771
2772 Standard_Real *PolesArray = &Array ;
2773 Standard_Real Inverse ;
2774
2775
2776 if (Matrix.LowerCol() != 1 ||
2777 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2778 ReturnCode = 1 ;
2779 goto FINISH ;
2780 }
2781
2782 for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow() ; ii++) {
2783 MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ?
2784 ii - LowerBandWidth : Matrix.LowerRow()) ;
2785
2786 for ( jj = MinIndex ; jj < ii ; jj++) {
2787
2788 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2789 PolesArray[(ii-1) * ArrayDimension + kk] +=
2790 PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2791 }
2792 }
2793 }
2794
2795 for (ii = Matrix.UpperRow() ; ii >= Matrix.LowerRow() ; ii--) {
2796 MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ?
2797 ii + UpperBandWidth : Matrix.UpperRow()) ;
2798
2799 for (jj = MaxIndex ; jj > ii ; jj--) {
2800
2801 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2802 PolesArray[(ii-1) * ArrayDimension + kk] -=
2803 PolesArray[(jj - 1) * ArrayDimension + kk] *
2804 Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2805 }
2806 }
2807
2808 //fixing a bug PRO18577 to avoid divizion by zero
2809
2810 Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ;
2811 Standard_Real Toler = 1.0e-16;
2812 if ( Abs(divizor) > Toler )
2813 Inverse = 1.0e0 / divizor ;
2814 else {
2815 Inverse = 1.0e0;
2816// cout << " BSplCLib::SolveBandedSystem() : zero determinant " << endl;
2817 ReturnCode = 1;
2818 goto FINISH;
2819 }
2820
2821 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2822 PolesArray[(ii-1) * ArrayDimension + kk] *= Inverse ;
2823 }
2824 }
2825 FINISH :
2826 return (ReturnCode) ;
2827}
2828
2829//=======================================================================
2830//function : Solves a LU factored Matrix
2831//purpose :
2832//=======================================================================
2833
2834Standard_Integer
2835BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2836 const Standard_Integer UpperBandWidth,
2837 const Standard_Integer LowerBandWidth,
2838 const Standard_Boolean HomogeneousFlag,
2839 const Standard_Integer ArrayDimension,
2840 Standard_Real& Poles,
2841 Standard_Real& Weights)
2842{
2843 Standard_Integer ii,
2844 kk,
2845 ErrorCode = 0,
2846 ReturnCode = 0 ;
2847
2848 Standard_Real Inverse,
2849 *PolesArray = &Poles,
2850 *WeightsArray = &Weights ;
2851
2852 if (Matrix.LowerCol() != 1 ||
2853 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2854 ReturnCode = 1 ;
2855 goto FINISH ;
2856 }
2857 if (HomogeneousFlag == Standard_False) {
2858
2859 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) {
2860
2861 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2862 PolesArray[ii * ArrayDimension + kk] *=
2863 WeightsArray[ii] ;
2864 }
2865 }
2866 }
2867 ErrorCode =
2868 BSplCLib::SolveBandedSystem(Matrix,
2869 UpperBandWidth,
2870 LowerBandWidth,
2871 ArrayDimension,
2872 Poles) ;
2873 if (ErrorCode != 0) {
2874 ReturnCode = 2 ;
2875 goto FINISH ;
2876 }
2877 ErrorCode =
2878 BSplCLib::SolveBandedSystem(Matrix,
2879 UpperBandWidth,
2880 LowerBandWidth,
2881 1,
2882 Weights) ;
2883 if (ErrorCode != 0) {
2884 ReturnCode = 3 ;
2885 goto FINISH ;
2886 }
2887 if (HomogeneousFlag == Standard_False) {
2888
2889 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) {
2890 Inverse = 1.0e0 / WeightsArray[ii] ;
2891
2892 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2893 PolesArray[ii * ArrayDimension + kk] *= Inverse ;
2894 }
2895 }
2896 }
2897 FINISH : return (ReturnCode) ;
2898}
2899
2900//=======================================================================
2901//function : BuildSchoenbergPoints
2902//purpose :
2903//=======================================================================
2904
2905void BSplCLib::BuildSchoenbergPoints(const Standard_Integer Degree,
2906 const TColStd_Array1OfReal& FlatKnots,
2907 TColStd_Array1OfReal& Parameters)
2908{
2909 Standard_Integer ii,
2910 jj ;
2911 Standard_Real Inverse ;
2912 Inverse = 1.0e0 / (Standard_Real)Degree ;
2913
2914 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
2915 Parameters(ii) = 0.0e0 ;
2916
2917 for (jj = 1 ; jj <= Degree ; jj++) {
2918 Parameters(ii) += FlatKnots(jj + ii) ;
2919 }
2920 Parameters(ii) *= Inverse ;
2921 }
2922}
2923
2924//=======================================================================
2925//function : Interpolate
2926//purpose :
2927//=======================================================================
2928
2929void BSplCLib::Interpolate(const Standard_Integer Degree,
2930 const TColStd_Array1OfReal& FlatKnots,
2931 const TColStd_Array1OfReal& Parameters,
2932 const TColStd_Array1OfInteger& ContactOrderArray,
2933 const Standard_Integer ArrayDimension,
2934 Standard_Real& Poles,
2935 Standard_Integer& InversionProblem)
2936{
2937 Standard_Integer ErrorCode,
2938 UpperBandWidth,
2939 LowerBandWidth ;
2940// Standard_Real *PolesArray = &Poles ;
2941 math_Matrix InterpolationMatrix(1, Parameters.Length(),
2942 1, 2 * Degree + 1) ;
2943 ErrorCode =
2944 BSplCLib::BuildBSpMatrix(Parameters,
2945 ContactOrderArray,
2946 FlatKnots,
2947 Degree,
2948 InterpolationMatrix,
2949 UpperBandWidth,
2950 LowerBandWidth) ;
2951 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
2952
2953 ErrorCode =
2954 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
2955 UpperBandWidth,
2956 LowerBandWidth,
2957 InversionProblem) ;
2958 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
2959
2960 ErrorCode =
2961 BSplCLib::SolveBandedSystem(InterpolationMatrix,
2962 UpperBandWidth,
2963 LowerBandWidth,
2964 ArrayDimension,
2965 Poles) ;
2966
2967 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
2968}
2969
2970//=======================================================================
2971//function : Interpolate
2972//purpose :
2973//=======================================================================
2974
2975void BSplCLib::Interpolate(const Standard_Integer Degree,
2976 const TColStd_Array1OfReal& FlatKnots,
2977 const TColStd_Array1OfReal& Parameters,
2978 const TColStd_Array1OfInteger& ContactOrderArray,
2979 const Standard_Integer ArrayDimension,
2980 Standard_Real& Poles,
2981 Standard_Real& Weights,
2982 Standard_Integer& InversionProblem)
2983{
2984 Standard_Integer ErrorCode,
2985 UpperBandWidth,
2986 LowerBandWidth ;
2987
2988 math_Matrix InterpolationMatrix(1, Parameters.Length(),
2989 1, 2 * Degree + 1) ;
2990 ErrorCode =
2991 BSplCLib::BuildBSpMatrix(Parameters,
2992 ContactOrderArray,
2993 FlatKnots,
2994 Degree,
2995 InterpolationMatrix,
2996 UpperBandWidth,
2997 LowerBandWidth) ;
2998 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
2999
3000 ErrorCode =
3001 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3002 UpperBandWidth,
3003 LowerBandWidth,
3004 InversionProblem) ;
3005 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3006
3007 ErrorCode =
3008 BSplCLib::SolveBandedSystem(InterpolationMatrix,
3009 UpperBandWidth,
3010 LowerBandWidth,
3011 Standard_False,
3012 ArrayDimension,
3013 Poles,
3014 Weights) ;
3015
3016 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
3017}
3018
3019//=======================================================================
3020//function : Evaluates a Bspline function : uses the ExtrapMode
3021//purpose : the function is extrapolated using the Taylor expansion
3022// of degree ExtrapMode[0] to the left and the Taylor
3023// expansion of degree ExtrapMode[1] to the right
3024// this evaluates the numerator by multiplying by the weights
3025// and evaluating it but does not call RationalDerivatives after
3026//=======================================================================
3027
3028void BSplCLib::Eval
3029(const Standard_Real Parameter,
3030 const Standard_Boolean PeriodicFlag,
3031 const Standard_Integer DerivativeRequest,
3032 Standard_Integer& ExtrapMode,
3033 const Standard_Integer Degree,
3034 const TColStd_Array1OfReal& FlatKnots,
3035 const Standard_Integer ArrayDimension,
3036 Standard_Real& Poles,
3037 Standard_Real& Weights,
3038 Standard_Real& PolesResults,
3039 Standard_Real& WeightsResults)
3040{
3041 Standard_Integer ii,
3042 jj,
3043 kk=0,
3044 Index,
3045 Index1,
3046 Index2,
3047 *ExtrapModeArray,
3048 Modulus,
3049 NewRequest,
3050 ExtrapolatingFlag[2],
3051 ErrorCode,
3052 ReturnCode,
3053 Order = Degree + 1,
3054 FirstNonZeroBsplineIndex,
3055 LocalRequest = DerivativeRequest ;
3056 Standard_Real *PResultArray,
3057 *WResultArray,
3058 *PolesArray,
3059 *WeightsArray,
3060 LocalParameter,
3061 Period,
3062 Inverse,
3063 Delta ;
3064 PolesArray = &Poles ;
3065 WeightsArray = &Weights ;
3066 ExtrapModeArray = &ExtrapMode ;
3067 PResultArray = &PolesResults ;
3068 WResultArray = &WeightsResults ;
3069 LocalParameter = Parameter ;
3070 ExtrapolatingFlag[0] =
3071 ExtrapolatingFlag[1] = 0 ;
3072 //
3073 // check if we are extrapolating to a degree which is smaller than
3074 // the degree of the Bspline
3075 //
3076 if (PeriodicFlag) {
3077 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3078
3079 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3080 LocalParameter -= Period ;
3081 }
3082
3083 while (LocalParameter < FlatKnots(2)) {
3084 LocalParameter += Period ;
3085 }
3086 }
3087 if (Parameter < FlatKnots(2) &&
3088 LocalRequest < ExtrapModeArray[0] &&
3089 ExtrapModeArray[0] < Degree) {
3090 LocalRequest = ExtrapModeArray[0] ;
3091 LocalParameter = FlatKnots(2) ;
3092 ExtrapolatingFlag[0] = 1 ;
3093 }
3094 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3095 LocalRequest < ExtrapModeArray[1] &&
3096 ExtrapModeArray[1] < Degree) {
3097 LocalRequest = ExtrapModeArray[1] ;
3098 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3099 ExtrapolatingFlag[1] = 1 ;
3100 }
3101 Delta = Parameter - LocalParameter ;
3102 if (LocalRequest >= Order) {
3103 LocalRequest = Degree ;
3104 }
3105 if (PeriodicFlag) {
3106 Modulus = FlatKnots.Length() - Degree -1 ;
3107 }
3108 else {
3109 Modulus = FlatKnots.Length() - Degree ;
3110 }
3111
3112 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3113 ErrorCode =
3114 BSplCLib::EvalBsplineBasis(1,
3115 LocalRequest,
3116 Order,
3117 FlatKnots,
3118 LocalParameter,
3119 FirstNonZeroBsplineIndex,
3120 BsplineBasis) ;
3121 if (ErrorCode != 0) {
3122 ReturnCode = 1 ;
3123 goto FINISH ;
3124 }
3125 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3126 Index = 0 ;
3127 Index2 = 0 ;
3128
3129 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3130 Index1 = FirstNonZeroBsplineIndex ;
3131
3132 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3133 PResultArray[Index + kk] = 0.0e0 ;
3134 }
3135 WResultArray[Index] = 0.0e0 ;
3136
3137 for (jj = 1 ; jj <= Order ; jj++) {
3138
3139 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3140 PResultArray[Index + kk] +=
3141 PolesArray[(Index1-1) * ArrayDimension + kk]
3142 * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3143 }
3144 WResultArray[Index2] += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3145
3146 Index1 = Index1 % Modulus ;
3147 Index1 += 1 ;
3148 }
3149 Index += ArrayDimension ;
3150 Index2 += 1 ;
3151 }
3152 }
3153 else {
3154 //
3155 // store Taylor expansion in LocalRealArray
3156 //
3157 NewRequest = DerivativeRequest ;
3158 if (NewRequest > Degree) {
3159 NewRequest = Degree ;
3160 }
3161 BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension);
3162 Index = 0 ;
3163 Inverse = 1.0e0 ;
3164
3165 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3166 Index1 = FirstNonZeroBsplineIndex ;
3167
3168 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3169 LocalRealArray[Index + kk] = 0.0e0 ;
3170 }
3171
3172 for (jj = 1 ; jj <= Order ; jj++) {
3173
3174 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3175 LocalRealArray[Index + kk] +=
3176 PolesArray[(Index1-1)*ArrayDimension + kk] *
3177 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3178 }
3179 Index1 = Index1 % Modulus ;
3180 Index1 += 1 ;
3181 }
3182
3183 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3184 LocalRealArray[Index + kk] *= Inverse ;
3185 }
3186 Index += ArrayDimension ;
3187 Inverse /= (Standard_Real) ii ;
3188 }
3189 PLib::EvalPolynomial(Delta,
3190 NewRequest,
3191 Degree,
3192 ArrayDimension,
3193 LocalRealArray[0],
3194 PolesResults) ;
3195 Index = 0 ;
3196 Inverse = 1.0e0 ;
3197
3198 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3199 Index1 = FirstNonZeroBsplineIndex ;
3200 LocalRealArray[Index] = 0.0e0 ;
3201
3202 for (jj = 1 ; jj <= Order ; jj++) {
3203 LocalRealArray[Index] +=
3204 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3205 Index1 = Index1 % Modulus ;
3206 Index1 += 1 ;
3207 }
3208 LocalRealArray[Index + kk] *= Inverse ;
3209 Index += 1 ;
3210 Inverse /= (Standard_Real) ii ;
3211 }
3212 PLib::EvalPolynomial(Delta,
3213 NewRequest,
3214 Degree,
3215 1,
3216 LocalRealArray[0],
3217 WeightsResults) ;
3218 }
3219 FINISH : ;
3220}
3221
3222//=======================================================================
3223//function : Evaluates a Bspline function : uses the ExtrapMode
3224//purpose : the function is extrapolated using the Taylor expansion
3225// of degree ExtrapMode[0] to the left and the Taylor
3226// expansion of degree ExtrapMode[1] to the right
3227// WARNING : the array Results is supposed to have at least
3228// (DerivativeRequest + 1) * ArrayDimension slots and the
3229//
3230//=======================================================================
3231
3232void BSplCLib::Eval
3233(const Standard_Real Parameter,
3234 const Standard_Boolean PeriodicFlag,
3235 const Standard_Integer DerivativeRequest,
3236 Standard_Integer& ExtrapMode,
3237 const Standard_Integer Degree,
3238 const TColStd_Array1OfReal& FlatKnots,
3239 const Standard_Integer ArrayDimension,
3240 Standard_Real& Poles,
3241 Standard_Real& Results)
3242{
3243 Standard_Integer ii,
3244 jj,
3245 kk,
3246 Index,
3247 Index1,
3248 *ExtrapModeArray,
3249 Modulus,
3250 NewRequest,
3251 ExtrapolatingFlag[2],
3252 ErrorCode,
3253 ReturnCode,
3254 Order = Degree + 1,
3255 FirstNonZeroBsplineIndex,
3256 LocalRequest = DerivativeRequest ;
3257
3258 Standard_Real *ResultArray,
3259 *PolesArray,
3260 LocalParameter,
3261 Period,
3262 Inverse,
3263 Delta ;
3264
3265 PolesArray = &Poles ;
3266 ExtrapModeArray = &ExtrapMode ;
3267 ResultArray = &Results ;
3268 LocalParameter = Parameter ;
3269 ExtrapolatingFlag[0] =
3270 ExtrapolatingFlag[1] = 0 ;
3271 //
3272 // check if we are extrapolating to a degree which is smaller than
3273 // the degree of the Bspline
3274 //
3275 if (PeriodicFlag) {
3276 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3277
3278 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3279 LocalParameter -= Period ;
3280 }
3281
3282 while (LocalParameter < FlatKnots(2)) {
3283 LocalParameter += Period ;
3284 }
3285 }
3286 if (Parameter < FlatKnots(2) &&
3287 LocalRequest < ExtrapModeArray[0] &&
3288 ExtrapModeArray[0] < Degree) {
3289 LocalRequest = ExtrapModeArray[0] ;
3290 LocalParameter = FlatKnots(2) ;
3291 ExtrapolatingFlag[0] = 1 ;
3292 }
3293 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3294 LocalRequest < ExtrapModeArray[1] &&
3295 ExtrapModeArray[1] < Degree) {
3296 LocalRequest = ExtrapModeArray[1] ;
3297 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3298 ExtrapolatingFlag[1] = 1 ;
3299 }
3300 Delta = Parameter - LocalParameter ;
3301 if (LocalRequest >= Order) {
3302 LocalRequest = Degree ;
3303 }
3304
3305 if (PeriodicFlag) {
3306 Modulus = FlatKnots.Length() - Degree -1 ;
3307 }
3308 else {
3309 Modulus = FlatKnots.Length() - Degree ;
3310 }
3311
3312 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3313
3314 ErrorCode =
3315 BSplCLib::EvalBsplineBasis(1,
3316 LocalRequest,
3317 Order,
3318 FlatKnots,
3319 LocalParameter,
3320 FirstNonZeroBsplineIndex,
3321 BsplineBasis);
3322 if (ErrorCode != 0) {
3323 ReturnCode = 1 ;
3324 goto FINISH ;
3325 }
3326 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3327 Index = 0 ;
3328
3329 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3330 Index1 = FirstNonZeroBsplineIndex ;
3331
3332 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3333 ResultArray[Index + kk] = 0.0e0 ;
3334 }
3335
3336 for (jj = 1 ; jj <= Order ; jj++) {
3337
3338 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3339 ResultArray[Index + kk] +=
3340 PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3341 }
3342 Index1 = Index1 % Modulus ;
3343 Index1 += 1 ;
3344 }
3345 Index += ArrayDimension ;
3346 }
3347 }
3348 else {
3349 //
3350 // store Taylor expansion in LocalRealArray
3351 //
3352 NewRequest = DerivativeRequest ;
3353 if (NewRequest > Degree) {
3354 NewRequest = Degree ;
3355 }
3356 BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension);
3357
3358 Index = 0 ;
3359 Inverse = 1.0e0 ;
3360
3361 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3362 Index1 = FirstNonZeroBsplineIndex ;
3363
3364 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3365 LocalRealArray[Index + kk] = 0.0e0 ;
3366 }
3367
3368 for (jj = 1 ; jj <= Order ; jj++) {
3369
3370 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3371 LocalRealArray[Index + kk] +=
3372 PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3373 }
3374 Index1 = Index1 % Modulus ;
3375 Index1 += 1 ;
3376 }
3377
3378 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3379 LocalRealArray[Index + kk] *= Inverse ;
3380 }
3381 Index += ArrayDimension ;
3382 Inverse /= (Standard_Real) ii ;
3383 }
3384 PLib::EvalPolynomial(Delta,
3385 NewRequest,
3386 Degree,
3387 ArrayDimension,
3388 LocalRealArray[0],
3389 Results) ;
3390 }
3391 FINISH : ;
3392}
3393
3394//=======================================================================
3395//function : TangExtendToConstraint
3396//purpose : Extends a Bspline function using the tangency map
3397// WARNING :
3398//
3399//
3400//=======================================================================
3401
3402void BSplCLib::TangExtendToConstraint
3403(const TColStd_Array1OfReal& FlatKnots,
3404 const Standard_Real C1Coefficient,
3405 const Standard_Integer NumPoles,
3406 Standard_Real& Poles,
3407 const Standard_Integer CDimension,
3408 const Standard_Integer CDegree,
3409 const TColStd_Array1OfReal& ConstraintPoint,
3410 const Standard_Integer Continuity,
3411 const Standard_Boolean After,
3412 Standard_Integer& NbPolesResult,
3413 Standard_Integer& NbKnotsResult,
3414 Standard_Real& KnotsResult,
3415 Standard_Real& PolesResult)
3416{
3417#if DEB
3418 if (CDegree<Continuity+1) {
3419 cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3420 }
3421#endif
3422 Standard_Real * Padr = &Poles ;
3423 Standard_Real * KRadr = &KnotsResult ;
3424 Standard_Real * PRadr = &PolesResult ;
3425
3426////////////////////////////////////////////////////////////////////////
3427//
0d969553 3428// 1. calculation of extension nD
7fd59977 3429//
3430////////////////////////////////////////////////////////////////////////
3431
0d969553 3432// Hermite matrix
7fd59977 3433 Standard_Integer Csize = Continuity + 2;
3434 math_Matrix MatCoefs(1,Csize, 1,Csize);
3435 if (After) {
0d969553
Y
3436 PLib::HermiteCoefficients(0, 1, // Limits
3437 Continuity, 0, // Orders of constraints
7fd59977 3438 MatCoefs);
3439 }
3440 else {
0d969553
Y
3441 PLib::HermiteCoefficients(0, 1, // Limits
3442 0, Continuity, // Orders of constraints
7fd59977 3443 MatCoefs);
3444 }
3445
3446
0d969553 3447// position at the node of connection
7fd59977 3448 Standard_Real Tbord ;
3449 if (After) {
3450 Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3451 }
3452 else {
3453 Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3454 }
3455 Standard_Boolean periodic_flag = Standard_False ;
3456 Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3457 extrap_mode[0] = extrap_mode[1] = CDegree;
3458 TColStd_Array1OfReal EvalBS(1, CDimension * (derivative_request+1)) ;
3459 Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3460 BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3461 CDegree,FlatKnots,CDimension,Poles,*Eadr);
3462
0d969553 3463// norm of the tangent at the node of connection
7fd59977 3464 math_Vector Tgte(1,CDimension);
3465
3466 for (ipos=1;ipos<=CDimension;ipos++) {
3467 Tgte(ipos) = EvalBS(ipos+CDimension);
3468 }
3469 Standard_Real L1=Tgte.Norm();
3470
3471
0d969553 3472// matrix of constraints
7fd59977 3473 math_Matrix Contraintes(1,Csize,1,CDimension);
3474 if (After) {
3475
3476 for (ipos=1;ipos<=CDimension;ipos++) {
3477 Contraintes(1,ipos) = EvalBS(ipos);
3478 Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3479 if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3480 if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3481 Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3482 }
3483 }
3484 else {
3485
3486 for (ipos=1;ipos<=CDimension;ipos++) {
3487 Contraintes(1,ipos) = ConstraintPoint(ipos);
3488 Contraintes(2,ipos) = EvalBS(ipos);
3489 if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3490 if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3491 if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3492 }
3493 }
3494
0d969553 3495// calculate the coefficients of extension
7fd59977 3496 Standard_Integer ii, jj, kk;
3497 TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3498 ExtraCoeffs.Init(0.);
3499
3500 for (ii=1; ii<=Csize; ii++) {
3501
3502 for (jj=1; jj<=Csize; jj++) {
3503
3504 for (kk=1; kk<=CDimension; kk++) {
3505 ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3506 }
3507 }
3508 }
3509
0d969553 3510// calculate the poles of extension
7fd59977 3511 TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3512 Standard_Real * EPadr = &ExtrapPoles(1) ;
3513 PLib::CoefficientsPoles(CDimension,
3514 ExtraCoeffs, PLib::NoWeights(),
3515 ExtrapPoles, PLib::NoWeights());
3516
0d969553 3517// calculate the nodes of extension with multiplicities
7fd59977 3518 TColStd_Array1OfReal ExtrapNoeuds(1,2);
3519 ExtrapNoeuds(1) = 0.;
3520 ExtrapNoeuds(2) = 1.;
3521 TColStd_Array1OfInteger ExtrapMults(1,2);
3522 ExtrapMults(1) = Csize;
3523 ExtrapMults(2) = Csize;
3524
0d969553 3525// flat nodes of extension
7fd59977 3526 TColStd_Array1OfReal FK2(1, Csize*2);
3527 BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3528
0d969553 3529// norm of the tangent at the connection point
7fd59977 3530 if (After) {
3531 BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3532 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3533 }
3534 else {
3535 BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3536 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3537 }
3538
3539 for (ipos=1;ipos<=CDimension;ipos++) {
3540 Tgte(ipos) = EvalBS(ipos+CDimension);
3541 }
3542 Standard_Real L2 = Tgte.Norm();
3543
0d969553 3544// harmonisation of degrees
7fd59977 3545 TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3546 TColStd_Array1OfReal NewK2(1, 2);
3547 TColStd_Array1OfInteger NewM2(1, 2);
3548 if (Csize-1<CDegree) {
3549 BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3550 ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3551 NewP2,NewK2,NewM2);
3552 }
3553 else {
3554 NewP2 = ExtrapPoles;
3555 NewK2 = ExtrapNoeuds;
3556 NewM2 = ExtrapMults;
3557 }
3558
0d969553 3559// flat nodes of extension after harmonization of degrees
7fd59977 3560 TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3561 BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3562
3563
3564////////////////////////////////////////////////////////////////////////
3565//
3566// 2. concatenation C0
3567//
3568////////////////////////////////////////////////////////////////////////
3569
0d969553 3570// ratio of reparametrization
7fd59977 3571 Standard_Real Ratio=1, Delta;
3572 if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3573 Ratio = L2 / L1;
3574 }
3575 if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3576
3577 if (After) {
0d969553 3578// do not touch the first BSpline
7fd59977 3579 Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3580 }
3581 else {
0d969553 3582// do not touch the second BSpline
7fd59977 3583 Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower());
3584 }
3585
0d969553 3586// result of the concatenation
7fd59977 3587 Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1;
3588 Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1);
3589 TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension);
3590 TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2);
3591
0d969553 3592// poles
7fd59977 3593 Standard_Integer indNP, indP, indEP;
3594 if (After) {
3595
3596 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3597
3598 for (jj=1; jj<=CDimension; jj++) {
3599 indNP = (ii-1)*CDimension+jj;
3600 indP = (ii-1)*CDimension+jj-1;
3601 indEP = (ii-NbP1)*CDimension+jj;
3602 if (ii<NbP1) NewPoles(indNP) = Padr[indP];
3603 else NewPoles(indNP) = NewP2(indEP);
3604 }
3605 }
3606 }
3607 else {
3608
3609 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3610
3611 for (jj=1; jj<=CDimension; jj++) {
3612 indNP = (ii-1)*CDimension+jj;
3613 indEP = (ii-1)*CDimension+jj;
3614 indP = (ii-NbP2)*CDimension+jj-1;
3615 if (ii<NbP2) NewPoles(indNP) = NewP2(indEP);
3616 else NewPoles(indNP) = Padr[indP];
3617 }
3618 }
3619 }
3620
0d969553 3621// flat nodes
7fd59977 3622 if (After) {
0d969553 3623// start with the nodes of the initial surface
7fd59977 3624
3625 for (ii=1; ii<NbK1; ii++) {
3626 NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1);
3627 }
0d969553 3628// continue with the reparameterized nodes of the extension
7fd59977 3629
3630 for (ii=1; ii<=NbK2-CDegree-1; ii++) {
3631 NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta;
3632 }
3633 }
3634 else {
0d969553 3635// start with the reparameterized nodes of the extension
7fd59977 3636
3637 for (ii=1; ii<NbK2-CDegree; ii++) {
3638 NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta;
3639 }
0d969553 3640// continue with the nodes of the initial surface
7fd59977 3641
3642 for (ii=2; ii<=NbK1; ii++) {
3643 NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1);
3644 }
3645 }
3646
3647
3648////////////////////////////////////////////////////////////////////////
3649//
0d969553 3650// 3. reduction of multiplicite at the node of connection
7fd59977 3651//
3652////////////////////////////////////////////////////////////////////////
3653
0d969553 3654// number of separate nodes
7fd59977 3655 Standard_Integer KLength = 1;
3656
3657 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3658 if (NewFlats(ii) != NewFlats(ii-1)) KLength++;
3659 }
3660
0d969553 3661// flat nodes --> nodes + multiplicities
7fd59977 3662 TColStd_Array1OfReal NewKnots (1, KLength);
3663 TColStd_Array1OfInteger NewMults (1, KLength);
3664 NewMults.Init(1);
3665 jj = 1;
3666 NewKnots(jj) = NewFlats(1);
3667
3668 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3669 if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++;
3670 else {
3671 jj++;
3672 NewKnots(jj) = NewFlats(ii);
3673 }
3674 }
3675
0d969553 3676// reduction of multiplicity at the second or the last but one node
7fd59977 3677 Standard_Integer Index = 2, M = CDegree;
3678 if (After) Index = KLength-1;
3679 TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension);
3680 TColStd_Array1OfReal ResultKnots (1, KLength);
3681 TColStd_Array1OfInteger ResultMults (1, KLength);
3682 Standard_Real Tol = 1.e-6;
3683 Standard_Boolean Ok = Standard_True;
3684
3685 while ( (M>CDegree-Continuity) && Ok) {
3686 Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension,
3687 NewPoles, NewKnots, NewMults,
3688 ResultPoles, ResultKnots, ResultMults, Tol);
3689 if (Ok) M--;
3690 }
3691
3692 if (M == CDegree) {
0d969553 3693// number of poles of the concatenation
7fd59977 3694 NbPolesResult = NbP1 + NbP2 - 1;
0d969553 3695// the poles of the concatenation
7fd59977 3696 Standard_Integer PLength = NbPolesResult*CDimension;
3697
3698 for (jj=1; jj<=PLength; jj++) {
3699 PRadr[jj-1] = NewPoles(jj);
3700 }
3701
0d969553 3702// flat nodes of the concatenation
7fd59977 3703 Standard_Integer ideb = 0;
3704
3705 for (jj=0; jj<NewKnots.Length(); jj++) {
3706 for (ii=0; ii<NewMults(jj+1); ii++) {
3707 KRadr[ideb+ii] = NewKnots(jj+1);
3708 }
3709 ideb += NewMults(jj+1);
3710 }
3711 NbKnotsResult = ideb;
3712 }
3713
3714 else {
0d969553 3715// number of poles of the result
7fd59977 3716 NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M;
0d969553 3717// the poles of the result
7fd59977 3718 Standard_Integer PLength = NbPolesResult*CDimension;
3719
3720 for (jj=0; jj<PLength; jj++) {
3721 PRadr[jj] = ResultPoles(jj+1);
3722 }
3723
0d969553 3724// flat nodes of the result
7fd59977 3725 Standard_Integer ideb = 0;
3726
3727 for (jj=0; jj<ResultKnots.Length(); jj++) {
3728 for (ii=0; ii<ResultMults(jj+1); ii++) {
3729 KRadr[ideb+ii] = ResultKnots(jj+1);
3730 }
3731 ideb += ResultMults(jj+1);
3732 }
3733 NbKnotsResult = ideb;
3734 }
3735}
3736
3737//=======================================================================
3738//function : Resolution
3739//purpose :
3740// d
0d969553 3741// Let C(t) = SUM Ci Bi(t) a Bspline curve of degree d
7fd59977 3742// i = 1,n
0d969553 3743// with nodes tj for j = 1,n+d+1
7fd59977 3744//
3745//
3746// ' C1 - Ci-1 d-1
0d969553 3747// Then C (t) = SUM d * --------- Bi (t)
7fd59977 3748// i = 2,n ti+d - ti
3749//
3750// d-1
0d969553 3751// for the base of BSpline Bi (t) of degree d-1.
7fd59977 3752//
0d969553 3753// Consequently the upper bound of the norm of the derivative from C is :
7fd59977 3754//
3755//
3756// | Ci - Ci-1 |
3757// d * Max | --------- |
3758// i = 2,n | ti+d - ti |
3759//
3760// N(t)
0d969553 3761// In the rational case set C(t) = -----
7fd59977 3762// D(t)
3763//
3764//
3765// D(t) = SUM Di Bi(t)
3766// i=1,n
3767//
3768// N(t) = SUM Di * Ci Bi(t)
3769// i =1,n
3770//
3771// N'(t) - D'(t) C(t)
3772// C'(t) = -----------------------
3773// D(t)
3774//
3775//
3776// N'(t) - D'(t) C(t) =
3777//
3778// Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t)) d-1
3779// SUM d * ---------------------------------------- * Bi (t) =
3780// i=2,n ti+d - ti
3781//
3782//
3783// Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) d-1
3784// SUM SUM d * ----------------------------------- * Betaj(t) * Bi (t)
3785//i=2,n j=1,n ti+d - ti
3786//
3787//
3788//
3789// Dj Bj(t)
3790// Betaj(t) = --------
3791// D(t)
3792//
0d969553
Y
3793// Betaj(t) form a partition >= 0 of the entity with support
3794// tj, tj+d+1. Consequently if Rj = {j-d, ...., j+d+d+1}
3795// obtain an upper bound of the derivative of C by taking :
7fd59977 3796//
3797//
3798//
3799//
3800//
3801//
3802// Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)
3803// Max Max d * -----------------------------------
3804// j=1,n i dans Rj ti+d - ti
3805//
3806// --------------------------------------------------------
3807//
3808// Min Di
3809// i =1,n
3810//
3811//
3812//=======================================================================
3813
3814void BSplCLib::Resolution( Standard_Real& Poles,
3815 const Standard_Integer ArrayDimension,
3816 const Standard_Integer NumPoles,
3817 const TColStd_Array1OfReal& Weights,
3818 const TColStd_Array1OfReal& FlatKnots,
3819 const Standard_Integer Degree,
3820 const Standard_Real Tolerance3D,
3821 Standard_Real& UTolerance)
3822{
3823 Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim;
3824 Standard_Integer lower,upper,ii_minus,jj,ii_miDim;
3825 Standard_Integer Deg1 = Degree + 1;
3826 Standard_Integer Deg2 = (Degree << 1) + 1;
3827 Standard_Real value,factor,W,min_weights,inverse;
3828 Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3;
3829 Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3;
3830 Standard_Real wg_ii_index, wg_ii_minus;
3831 Standard_Real *PA,max_derivative;
3832 const Standard_Real * FK = &FlatKnots(FlatKnots.Lower());
3833 PA = &Poles;
3834 max_derivative = 0.0e0;
3835 num_poles = FlatKnots.Length() - Deg1;
3836 switch (ArrayDimension) {
3837 case 2 : {
3838 if (&Weights != NULL) {
3839 const Standard_Real * WG = &Weights(Weights.Lower());
3840 min_weights = WG[0];
3841
3842 for (ii = 1 ; ii < NumPoles ; ii++) {
3843 W = WG[ii];
3844 if (W < min_weights) min_weights = W;
3845 }
3846
3847 for (ii = 1 ; ii < num_poles ; ii++) {
3848 ii_index = ii % NumPoles;
3849 ii_inDim = ii_index << 1;
3850 ii_minus = (ii - 1) % NumPoles;
3851 ii_miDim = ii_minus << 1;
3852 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3853 pa_ii_inDim_1 = PA[ii_inDim];
3854 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3855 pa_ii_miDim_1 = PA[ii_miDim];
3856 wg_ii_index = WG[ii_index];
3857 wg_ii_minus = WG[ii_minus];
3858 inverse = FK[ii + Degree] - FK[ii];
3859 inverse = 1.0e0 / inverse;
3860 lower = ii - Deg1;
3861 if (lower < 0) lower = 0;
3862 upper = Deg2 + ii;
3863 if (upper > num_poles) upper = num_poles;
3864
3865 for (jj = lower ; jj < upper ; jj++) {
3866 jj_index = jj % NumPoles;
3867 jj_index = jj_index << 1;
3868 value = 0.0e0;
3869 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3870 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3871 if (factor < 0) factor = - factor;
3872 value += factor; jj_index++;
3873 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3874 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3875 if (factor < 0) factor = - factor;
3876 value += factor;
3877 value *= inverse;
3878 if (max_derivative < value) max_derivative = value;
3879 }
3880 }
3881 max_derivative /= min_weights;
3882 }
3883 else {
3884
3885 for (ii = 1 ; ii < num_poles ; ii++) {
3886 ii_index = ii % NumPoles;
3887 ii_index = ii_index << 1;
3888 ii_minus = (ii - 1) % NumPoles;
3889 ii_minus = ii_minus << 1;
3890 inverse = FK[ii + Degree] - FK[ii];
3891 inverse = 1.0e0 / inverse;
3892 value = 0.0e0;
3893 factor = PA[ii_index] - PA[ii_minus];
3894 if (factor < 0) factor = - factor;
3895 value += factor; ii_index++; ii_minus++;
3896 factor = PA[ii_index] - PA[ii_minus];
3897 if (factor < 0) factor = - factor;
3898 value += factor;
3899 value *= inverse;
3900 if (max_derivative < value) max_derivative = value;
3901 }
3902 }
3903 break;
3904 }
3905 case 3 : {
3906 if (&Weights != NULL) {
3907 const Standard_Real * WG = &Weights(Weights.Lower());
3908 min_weights = WG[0];
3909
3910 for (ii = 1 ; ii < NumPoles ; ii++) {
3911 W = WG[ii];
3912 if (W < min_weights) min_weights = W;
3913 }
3914
3915 for (ii = 1 ; ii < num_poles ; ii++) {
3916 ii_index = ii % NumPoles;
3917 ii_inDim = (ii_index << 1) + ii_index;
3918 ii_minus = (ii - 1) % NumPoles;
3919 ii_miDim = (ii_minus << 1) + ii_minus;
3920 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3921 pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
3922 pa_ii_inDim_2 = PA[ii_inDim];
3923 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3924 pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
3925 pa_ii_miDim_2 = PA[ii_miDim];
3926 wg_ii_index = WG[ii_index];
3927 wg_ii_minus = WG[ii_minus];
3928 inverse = FK[ii + Degree] - FK[ii];
3929 inverse = 1.0e0 / inverse;
3930 lower = ii - Deg1;
3931 if (lower < 0) lower = 0;
3932 upper = Deg2 + ii;
3933 if (upper > num_poles) upper = num_poles;
3934
3935 for (jj = lower ; jj < upper ; jj++) {
3936 jj_index = jj % NumPoles;
3937 jj_index = (jj_index << 1) + jj_index;
3938 value = 0.0e0;
3939 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3940 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3941 if (factor < 0) factor = - factor;
3942 value += factor; jj_index++;
3943 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3944 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3945 if (factor < 0) factor = - factor;
3946 value += factor; jj_index++;
3947 factor = (((PA[jj_index