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