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