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