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