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