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