1 // Created on: 1991-08-09
3 // Copyright (c) 1991-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
22 // Modified RLE 9 Sep 1993
23 // pmn : modified 28-01-97 : fixed a mistake in LocateParameter (PRO6973)
24 // pmn : modified 4-11-96 : fixed a mistake in BuildKnots (PRO6124)
25 // pmn : modified 28-Jun-96 : fixed a mistake in AntiBoorScheme
26 // xab : modified 15-Jun-95 : fixed a mistake in IsRational
27 // xab : modified 15-Mar-95 : removed Epsilon comparison in IsRational
28 // added RationalDerivatives.
29 // xab : 30-Mar-95 : fixed coupling with lti in RationalDerivatives
30 // xab : 15-Mar-96 : fixed a typo in Eval with extrapolation
31 // jct : 15-Apr-97 : added TangExtendToConstraint
32 // jct : 24-Apr-97 : correction on computation of Tbord and NewFlatKnots
33 // in TangExtendToConstraint; Continuity can be equal to 0
35 #include <BSplCLib.ixx>
37 #include <PLib_LocalArray.hxx>
38 #include <Precision.hxx>
39 #include <Standard_NotImplemented.hxx>
43 typedef TColgp_Array1OfPnt Array1OfPnt;
44 typedef TColStd_Array1OfReal Array1OfReal;
45 typedef TColStd_Array1OfInteger Array1OfInteger;
47 //=======================================================================
48 //class : BSplCLib_LocalMatrix
49 //purpose: Auxiliary class optimizing creation of matrix buffer for
50 // evaluation of bspline (using stack allocation for main matrix)
51 //=======================================================================
53 class BSplCLib_LocalMatrix : public math_Matrix
56 BSplCLib_LocalMatrix (Standard_Integer DerivativeRequest, Standard_Integer Order)
57 : math_Matrix (myBuffer, 1, DerivativeRequest + 1, 1, Order)
59 Standard_OutOfRange_Raise_if (DerivativeRequest > BSplCLib::MaxDegree() ||
60 Order > BSplCLib::MaxDegree()+1 || BSplCLib::MaxDegree() > 25,
61 "BSplCLib: bspline degree is greater than maximum supported");
65 // local buffer, to be sufficient for addressing by index [Degree+1][Degree+1]
66 // (see math_Matrix implementation)
67 Standard_Real myBuffer[27*27];
70 typedef PLib_LocalArray BSplCLib_LocalArray;
72 //=======================================================================
75 //=======================================================================
77 void BSplCLib::Hunt (const Array1OfReal& XX,
78 const Standard_Real X,
79 Standard_Integer& Ilc)
81 // replaced by simple dichotomy (RLE)
83 const Standard_Real *px = &XX(Ilc);
90 Standard_Integer Ihi = XX.Upper();
97 while (Ihi - Ilc != 1) {
98 Im = (Ihi + Ilc) >> 1;
99 if (X > px[Im]) Ilc = Im;
104 //=======================================================================
105 //function : FirstUKnotIndex
107 //=======================================================================
109 Standard_Integer BSplCLib::FirstUKnotIndex (const Standard_Integer Degree,
110 const TColStd_Array1OfInteger& Mults)
112 Standard_Integer Index = Mults.Lower();
113 Standard_Integer SigmaMult = Mults(Index);
115 while (SigmaMult <= Degree) {
117 SigmaMult += Mults (Index);
122 //=======================================================================
123 //function : LastUKnotIndex
125 //=======================================================================
127 Standard_Integer BSplCLib::LastUKnotIndex (const Standard_Integer Degree,
128 const Array1OfInteger& Mults)
130 Standard_Integer Index = Mults.Upper();
131 Standard_Integer SigmaMult = Mults(Index);
133 while (SigmaMult <= Degree) {
135 SigmaMult += Mults.Value (Index);
140 //=======================================================================
141 //function : FlatIndex
143 //=======================================================================
145 Standard_Integer BSplCLib::FlatIndex
146 (const Standard_Integer Degree,
147 const Standard_Integer Index,
148 const TColStd_Array1OfInteger& Mults,
149 const Standard_Boolean Periodic)
151 Standard_Integer i, index = Index;
152 const Standard_Integer MLower = Mults.Lower();
153 const Standard_Integer *pmu = &Mults(MLower);
156 for (i = MLower + 1; i <= Index; i++)
161 index += pmu[MLower] - 1;
165 //=======================================================================
166 //function : LocateParameter
167 //purpose : Processing of nodes with multiplicities
168 //pmn 28-01-97 -> compute eventual of the period.
169 //=======================================================================
171 void BSplCLib::LocateParameter
172 (const Standard_Integer , //Degree,
173 const Array1OfReal& Knots,
174 const Array1OfInteger& , //Mults,
175 const Standard_Real U,
176 const Standard_Boolean IsPeriodic,
177 const Standard_Integer FromK1,
178 const Standard_Integer ToK2,
179 Standard_Integer& KnotIndex,
182 Standard_Real uf = 0, ul=1;
184 uf = Knots(Knots.Lower());
185 ul = Knots(Knots.Upper());
187 BSplCLib::LocateParameter(Knots,U,IsPeriodic,FromK1,ToK2,
188 KnotIndex,NewU, uf, ul);
191 //=======================================================================
192 //function : LocateParameter
193 //purpose : For plane nodes
194 // pmn 28-01-97 -> There is a need of the degre to calculate
195 // the eventual period
196 //=======================================================================
198 void BSplCLib::LocateParameter
199 (const Standard_Integer Degree,
200 const Array1OfReal& Knots,
201 const Standard_Real U,
202 const Standard_Boolean IsPeriodic,
203 const Standard_Integer FromK1,
204 const Standard_Integer ToK2,
205 Standard_Integer& KnotIndex,
209 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
211 Knots(Knots.Lower() + Degree),
212 Knots(Knots.Upper() - Degree));
214 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
220 //=======================================================================
221 //function : LocateParameter
222 //purpose : Effective computation
223 // pmn 28-01-97 : Add limits of the period as input argument,
224 // as it is imposible to produce them at this level.
225 //=======================================================================
227 void BSplCLib::LocateParameter
228 (const TColStd_Array1OfReal& Knots,
229 const Standard_Real U,
230 const Standard_Boolean IsPeriodic,
231 const Standard_Integer FromK1,
232 const Standard_Integer ToK2,
233 Standard_Integer& KnotIndex,
235 const Standard_Real UFirst,
236 const Standard_Real ULast)
238 Standard_Integer First,Last;
247 Standard_Integer Last1 = Last - 1;
250 Standard_Real Period = ULast - UFirst;
252 while (NewU > ULast )
255 while (NewU < UFirst)
259 BSplCLib::Hunt (Knots, NewU, KnotIndex);
261 Standard_Real Eps = Epsilon(U);
263 if (Eps < 0) Eps = - Eps;
264 Standard_Integer KLower = Knots.Lower();
265 const Standard_Real *knots = &Knots(KLower);
267 if ( KnotIndex < Knots.Upper()) {
268 val = NewU - knots[KnotIndex + 1];
269 if (val < 0) val = - val;
270 // <= to be coherent with Segment where Eps corresponds to a bit of error.
271 if (val <= Eps) KnotIndex++;
273 if (KnotIndex < First) KnotIndex = First;
274 if (KnotIndex > Last1) KnotIndex = Last1;
276 if (KnotIndex != Last1) {
277 Standard_Real K1 = knots[KnotIndex];
278 Standard_Real K2 = knots[KnotIndex + 1];
280 if (val < 0) val = - val;
285 K2 = knots[KnotIndex + 1];
287 if (val < 0) val = - val;
292 //=======================================================================
293 //function : LocateParameter
294 //purpose : the index is recomputed only if out of range
295 //pmn 28-01-97 -> eventual computation of the period.
296 //=======================================================================
298 void BSplCLib::LocateParameter
299 (const Standard_Integer Degree,
300 const TColStd_Array1OfReal& Knots,
301 const TColStd_Array1OfInteger& Mults,
302 const Standard_Real U,
303 const Standard_Boolean Periodic,
304 Standard_Integer& KnotIndex,
307 Standard_Integer first,last;
310 first = Knots.Lower();
311 last = Knots.Upper();
314 first = FirstUKnotIndex(Degree,Mults);
315 last = LastUKnotIndex (Degree,Mults);
319 first = Knots.Lower() + Degree;
320 last = Knots.Upper() - Degree;
322 if ( KnotIndex < first || KnotIndex > last)
323 BSplCLib::LocateParameter(Knots, U, Periodic, first, last,
324 KnotIndex, NewU, Knots(first), Knots(last));
329 //=======================================================================
330 //function : MaxKnotMult
332 //=======================================================================
334 Standard_Integer BSplCLib::MaxKnotMult
335 (const Array1OfInteger& Mults,
336 const Standard_Integer FromK1,
337 const Standard_Integer ToK2)
339 Standard_Integer MLower = Mults.Lower();
340 const Standard_Integer *pmu = &Mults(MLower);
342 Standard_Integer MaxMult = pmu[FromK1];
344 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
345 if (MaxMult < pmu[i]) MaxMult = pmu[i];
350 //=======================================================================
351 //function : MinKnotMult
353 //=======================================================================
355 Standard_Integer BSplCLib::MinKnotMult
356 (const Array1OfInteger& Mults,
357 const Standard_Integer FromK1,
358 const Standard_Integer ToK2)
360 Standard_Integer MLower = Mults.Lower();
361 const Standard_Integer *pmu = &Mults(MLower);
363 Standard_Integer MinMult = pmu[FromK1];
365 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
366 if (MinMult > pmu[i]) MinMult = pmu[i];
371 //=======================================================================
374 //=======================================================================
376 Standard_Integer BSplCLib::NbPoles(const Standard_Integer Degree,
377 const Standard_Boolean Periodic,
378 const TColStd_Array1OfInteger& Mults)
380 Standard_Integer i,sigma = 0;
381 Standard_Integer f = Mults.Lower();
382 Standard_Integer l = Mults.Upper();
383 const Standard_Integer * pmu = &Mults(f);
385 Standard_Integer Mf = pmu[f];
386 Standard_Integer Ml = pmu[l];
387 if (Mf <= 0) return 0;
388 if (Ml <= 0) return 0;
390 if (Mf > Degree) return 0;
391 if (Ml > Degree) return 0;
392 if (Mf != Ml ) return 0;
396 Standard_Integer Deg1 = Degree + 1;
397 if (Mf > Deg1) return 0;
398 if (Ml > Deg1) return 0;
399 sigma = Mf + Ml - Deg1;
402 for (i = f + 1; i < l; i++) {
403 if (pmu[i] <= 0 ) return 0;
404 if (pmu[i] > Degree) return 0;
410 //=======================================================================
411 //function : KnotSequenceLength
413 //=======================================================================
415 Standard_Integer BSplCLib::KnotSequenceLength
416 (const TColStd_Array1OfInteger& Mults,
417 const Standard_Integer Degree,
418 const Standard_Boolean Periodic)
420 Standard_Integer i,l = 0;
421 Standard_Integer MLower = Mults.Lower();
422 Standard_Integer MUpper = Mults.Upper();
423 const Standard_Integer * pmu = &Mults(MLower);
426 for (i = MLower; i <= MUpper; i++)
428 if (Periodic) l += 2 * (Degree + 1 - pmu[MLower]);
432 //=======================================================================
433 //function : KnotSequence
435 //=======================================================================
437 void BSplCLib::KnotSequence
438 (const TColStd_Array1OfReal& Knots,
439 const TColStd_Array1OfInteger& Mults,
440 TColStd_Array1OfReal& KnotSeq)
442 BSplCLib::KnotSequence(Knots,Mults,0,Standard_False,KnotSeq);
445 //=======================================================================
446 //function : KnotSequence
448 //=======================================================================
450 void BSplCLib::KnotSequence
451 (const TColStd_Array1OfReal& Knots,
452 const TColStd_Array1OfInteger& Mults,
453 const Standard_Integer Degree,
454 const Standard_Boolean Periodic,
455 TColStd_Array1OfReal& KnotSeq)
458 Standard_Integer Mult;
459 Standard_Integer MLower = Mults.Lower();
460 const Standard_Integer * pmu = &Mults(MLower);
462 Standard_Integer KLower = Knots.Lower();
463 Standard_Integer KUpper = Knots.Upper();
464 const Standard_Real * pkn = &Knots(KLower);
466 Standard_Integer M1 = Degree + 1 - pmu[MLower]; // for periodic
467 Standard_Integer i,j,index = Periodic ? M1 + 1 : 1;
469 for (i = KLower; i <= KUpper; i++) {
473 for (j = 1; j <= Mult; j++) {
479 Standard_Real period = pkn[KUpper] - pkn[KLower];
484 for (i = M1; i >= 1; i--) {
485 KnotSeq(i) = pkn[j] - period;
495 for (i = index; i <= KnotSeq.Upper(); i++) {
496 KnotSeq(i) = pkn[j] + period;
506 //=======================================================================
507 //function : KnotsLength
509 //=======================================================================
510 Standard_Integer BSplCLib::KnotsLength(const TColStd_Array1OfReal& SeqKnots,
511 // const Standard_Boolean Periodic)
512 const Standard_Boolean )
514 Standard_Integer sizeMult = 1;
515 Standard_Real val = SeqKnots(1);
516 for (Standard_Integer jj=2;
517 jj<=SeqKnots.Length();jj++)
519 // test on strict equality on nodes
520 if (SeqKnots(jj)!=val)
529 //=======================================================================
532 //=======================================================================
533 void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots,
534 TColStd_Array1OfReal &knots,
535 TColStd_Array1OfInteger &mult,
536 // const Standard_Boolean Periodic)
537 const Standard_Boolean )
539 Standard_Real val = SeqKnots(1);
540 Standard_Integer kk=1;
544 for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++)
546 // test on strict equality on nodes
547 if (SeqKnots(jj)!=val)
561 //=======================================================================
562 //function : KnotForm
564 //=======================================================================
566 BSplCLib_KnotDistribution BSplCLib::KnotForm
567 (const Array1OfReal& Knots,
568 const Standard_Integer FromK1,
569 const Standard_Integer ToK2)
571 Standard_Real DU0,DU1,Ui,Uj,Eps0,val;
572 BSplCLib_KnotDistribution KForm = BSplCLib_Uniform;
574 Standard_Integer KLower = Knots.Lower();
575 const Standard_Real * pkn = &Knots(KLower);
578 if (Ui < 0) Ui = - Ui;
579 Uj = pkn[FromK1 + 1];
580 if (Uj < 0) Uj = - Uj;
582 if (DU0 < 0) DU0 = - DU0;
583 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
584 Standard_Integer i = FromK1 + 1;
586 while (KForm != BSplCLib_NonUniform && i < ToK2) {
588 if (Ui < 0) Ui = - Ui;
591 if (Uj < 0) Uj = - Uj;
593 if (DU1 < 0) DU1 = - DU1;
595 if (val < 0) val = -val;
596 if (val > Eps0) KForm = BSplCLib_NonUniform;
598 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
603 //=======================================================================
604 //function : MultForm
606 //=======================================================================
608 BSplCLib_MultDistribution BSplCLib::MultForm
609 (const Array1OfInteger& Mults,
610 const Standard_Integer FromK1,
611 const Standard_Integer ToK2)
613 Standard_Integer First,Last;
622 Standard_Integer MLower = Mults.Lower();
623 const Standard_Integer *pmu = &Mults(MLower);
625 Standard_Integer FirstMult = pmu[First];
626 BSplCLib_MultDistribution MForm = BSplCLib_Constant;
627 Standard_Integer i = First + 1;
628 Standard_Integer Mult = pmu[i];
630 // while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR????????
631 while (MForm != BSplCLib_NonConstant && i <= Last) {
632 if (i == First + 1) {
633 if (Mult != FirstMult) MForm = BSplCLib_QuasiConstant;
635 else if (i == Last) {
636 if (MForm == BSplCLib_QuasiConstant) {
637 if (FirstMult != pmu[i]) MForm = BSplCLib_NonConstant;
640 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
644 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
652 //=======================================================================
653 //function : Reparametrize
655 //=======================================================================
657 void BSplCLib::Reparametrize
658 (const Standard_Real U1,
659 const Standard_Real U2,
662 Standard_Integer Lower = Knots.Lower();
663 Standard_Integer Upper = Knots.Upper();
664 Standard_Real UFirst = Min (U1, U2);
665 Standard_Real ULast = Max (U1, U2);
666 Standard_Real NewLength = ULast - UFirst;
667 BSplCLib_KnotDistribution KSet = BSplCLib::KnotForm (Knots, Lower, Upper);
668 if (KSet == BSplCLib_Uniform) {
669 Standard_Real DU = NewLength / (Upper - Lower);
670 Knots (Lower) = UFirst;
672 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
673 Knots (i) = Knots (i-1) + DU;
679 Standard_Real K1 = Knots (Lower);
680 Standard_Real Length = Knots (Upper) - Knots (Lower);
681 Knots (Lower) = UFirst;
683 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
685 Ratio = (K2 - K1) / Length;
686 Knots (i) = Knots (i-1) + (NewLength * Ratio);
689 Standard_Real Eps = Epsilon( Abs(Knots(i-1)) );
690 if (Knots(i) - Knots(i-1) <= Eps)
691 Knots(i) = NextAfter (Knots(i-1) + Eps, RealLast());
698 //=======================================================================
701 //=======================================================================
703 void BSplCLib::Reverse(TColStd_Array1OfReal& Knots)
705 Standard_Integer first = Knots.Lower();
706 Standard_Integer last = Knots.Upper();
707 Standard_Real kfirst = Knots(first);
708 Standard_Real klast = Knots(last);
709 Standard_Real tfirst = kfirst;
710 Standard_Real tlast = klast;
714 while (first <= last) {
715 tfirst += klast - Knots(last);
716 tlast -= Knots(first) - kfirst;
717 kfirst = Knots(first);
719 Knots(first) = tfirst;
726 //=======================================================================
729 //=======================================================================
731 void BSplCLib::Reverse(TColStd_Array1OfInteger& Mults)
733 Standard_Integer first = Mults.Lower();
734 Standard_Integer last = Mults.Upper();
735 Standard_Integer temp;
737 while (first < last) {
739 Mults(first) = Mults(last);
746 //=======================================================================
749 //=======================================================================
751 void BSplCLib::Reverse(TColStd_Array1OfReal& Weights,
752 const Standard_Integer L)
754 Standard_Integer i, l = L;
755 l = Weights.Lower()+(l-Weights.Lower())%(Weights.Upper()-Weights.Lower()+1);
757 TColStd_Array1OfReal temp(0,Weights.Length()-1);
759 for (i = Weights.Lower(); i <= l; i++)
760 temp(l-i) = Weights(i);
762 for (i = l+1; i <= Weights.Upper(); i++)
763 temp(l-Weights.Lower()+Weights.Upper()-i+1) = Weights(i);
765 for (i = Weights.Lower(); i <= Weights.Upper(); i++)
766 Weights(i) = temp(i-Weights.Lower());
769 //=======================================================================
770 //function : IsRational
772 //=======================================================================
774 Standard_Boolean BSplCLib::IsRational(const TColStd_Array1OfReal& Weights,
775 const Standard_Integer I1,
776 const Standard_Integer I2,
777 // const Standard_Real Epsi)
778 const Standard_Real )
780 Standard_Integer i, f = Weights.Lower(), l = Weights.Length();
781 Standard_Integer I3 = I2 - f;
782 const Standard_Real * WG = &Weights(f);
785 for (i = I1 - f; i < I3; i++) {
786 if (WG[f + (i % l)] != WG[f + ((i + 1) % l)]) return Standard_True;
788 return Standard_False ;
791 //=======================================================================
793 //purpose : evaluate point and derivatives
794 //=======================================================================
796 void BSplCLib::Eval(const Standard_Real U,
797 const Standard_Integer Degree,
798 Standard_Real& Knots,
799 const Standard_Integer Dimension,
800 Standard_Real& Poles)
802 Standard_Integer step,i,Dms,Dm1,Dpi,Sti;
803 Standard_Real X, Y, *poles, *knots = &Knots;
811 for (step = - 1; step < Dm1; step++) {
817 for (i = 0; i < Dms; i++) {
820 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
822 poles[0] *= X; poles[0] += Y * poles[1];
830 for (step = - 1; step < Dm1; step++) {
836 for (i = 0; i < Dms; i++) {
839 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
841 poles[0] *= X; poles[0] += Y * poles[2];
842 poles[1] *= X; poles[1] += Y * poles[3];
850 for (step = - 1; step < Dm1; step++) {
856 for (i = 0; i < Dms; i++) {
859 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
861 poles[0] *= X; poles[0] += Y * poles[3];
862 poles[1] *= X; poles[1] += Y * poles[4];
863 poles[2] *= X; poles[2] += Y * poles[5];
871 for (step = - 1; step < Dm1; step++) {
877 for (i = 0; i < Dms; i++) {
880 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
882 poles[0] *= X; poles[0] += Y * poles[4];
883 poles[1] *= X; poles[1] += Y * poles[5];
884 poles[2] *= X; poles[2] += Y * poles[6];
885 poles[3] *= X; poles[3] += Y * poles[7];
894 for (step = - 1; step < Dm1; step++) {
900 for (i = 0; i < Dms; i++) {
903 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
906 for (k = 0; k < Dimension; k++) {
908 poles[k] += Y * poles[k + Dimension];
917 //=======================================================================
918 //function : BoorScheme
920 //=======================================================================
922 void BSplCLib::BoorScheme(const Standard_Real U,
923 const Standard_Integer Degree,
924 Standard_Real& Knots,
925 const Standard_Integer Dimension,
926 Standard_Real& Poles,
927 const Standard_Integer Depth,
928 const Standard_Integer Length)
931 // Compute the values
935 // for i = 0 to Depth,
936 // j = 0 to Length - i
938 // The Boor scheme is :
941 // P(i,j) = x * P(i-1,j) + (1-x) * P(i-1,j+1)
943 // where x = (knot(i+j+Degree) - U) / (knot(i+j+Degree) - knot(i+j))
946 // The values are stored in the array Poles
947 // They are alternatively written if the odd and even positions.
949 // The successives contents of the array are
950 // ***** means unitialised, l = Degree + Length
952 // P(0,0) ****** P(0,1) ...... P(0,l-1) ******** P(0,l)
953 // P(0,0) P(1,0) P(0,1) ...... P(0,l-1) P(1,l-1) P(0,l)
954 // P(0,0) P(1,0) P(2,0) ...... P(2,l-1) P(1,l-1) P(0,l)
957 Standard_Integer i,k,step;
958 Standard_Real *knots = &Knots;
959 Standard_Real *pole, *firstpole = &Poles - 2 * Dimension;
960 // the steps of recursion
962 for (step = 0; step < Depth; step++) {
963 firstpole += Dimension;
965 // compute the new row of poles
967 for (i = step; i < Length; i++) {
968 pole += 2 * Dimension;
970 Standard_Real X = (knots[i+Degree-step] - U)
971 / (knots[i+Degree-step] - knots[i]);
972 Standard_Real Y = 1. - X;
974 // P(i,j) = X * P(i-1,j) + (1-X) * P(i-1,j+1)
976 for (k = 0; k < Dimension; k++)
977 pole[k] = X * pole[k - Dimension] + Y * pole[k + Dimension];
982 //=======================================================================
983 //function : AntiBoorScheme
985 //=======================================================================
987 Standard_Boolean BSplCLib::AntiBoorScheme(const Standard_Real U,
988 const Standard_Integer Degree,
989 Standard_Real& Knots,
990 const Standard_Integer Dimension,
991 Standard_Real& Poles,
992 const Standard_Integer Depth,
993 const Standard_Integer Length,
994 const Standard_Real Tolerance)
996 // do the Boor scheme reverted.
998 Standard_Integer i,k,step, half_length;
999 Standard_Real *knots = &Knots;
1000 Standard_Real z,X,Y,*pole, *firstpole = &Poles + (Depth-1) * Dimension;
1002 // Test the special case length = 1
1003 // only verification of the central point
1006 X = (knots[Degree] - U) / (knots[Degree] - knots[0]);
1009 for (k = 0; k < Dimension; k++) {
1010 z = X * firstpole[k] + Y * firstpole[k+2*Dimension];
1011 if (Abs(z - firstpole[k+Dimension]) > Tolerance)
1012 return Standard_False;
1014 return Standard_True;
1018 // the steps of recursion
1020 for (step = Depth-1; step >= 0; step--) {
1021 firstpole -= Dimension;
1024 // first step from left to right
1026 for (i = step; i < Length-1; i++) {
1027 pole += 2 * Dimension;
1029 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1032 for (k = 0; k < Dimension; k++)
1033 pole[k+Dimension] = (pole[k] - X*pole[k-Dimension]) / Y;
1037 // second step from right to left
1038 pole += 4* Dimension;
1039 half_length = (Length - 1 + step) / 2 ;
1041 // only do half of the way from right to left
1042 // otherwise it start degenerating because of
1046 for (i = Length-1; i > half_length ; i--) {
1047 pole -= 2 * Dimension;
1050 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1053 for (k = 0; k < Dimension; k++) {
1054 z = (pole[k] - Y * pole[k+Dimension]) / X;
1055 if (Abs(z-pole[k-Dimension]) > Tolerance)
1056 return Standard_False;
1057 pole[k-Dimension] += z;
1058 pole[k-Dimension] /= 2.;
1062 return Standard_True;
1065 //=======================================================================
1066 //function : Derivative
1068 //=======================================================================
1070 void BSplCLib::Derivative(const Standard_Integer Degree,
1071 Standard_Real& Knots,
1072 const Standard_Integer Dimension,
1073 const Standard_Integer Length,
1074 const Standard_Integer Order,
1075 Standard_Real& Poles)
1077 Standard_Integer i,k,step,span = Degree;
1078 Standard_Real *knot = &Knots;
1080 for (step = 1; step <= Order; step++) {
1081 Standard_Real* pole = &Poles;
1083 for (i = step; i < Length; i++) {
1084 Standard_Real coef = - span / (knot[i+span] - knot[i]);
1086 for (k = 0; k < Dimension; k++) {
1087 pole[k] -= pole[k+Dimension];
1096 //=======================================================================
1099 //=======================================================================
1101 void BSplCLib::Bohm(const Standard_Real U,
1102 const Standard_Integer Degree,
1103 const Standard_Integer N,
1104 Standard_Real& Knots,
1105 const Standard_Integer Dimension,
1106 Standard_Real& Poles)
1108 // First phase independant of U, compute the poles of the derivatives
1109 Standard_Integer i,j,iDim,min,Dmi,DDmi,jDmi,Degm1;
1110 Standard_Real *knot = &Knots, *pole, coef, *tbis, *psav, *psDD, *psDDmDim;
1112 if (N < Degree) min = N;
1115 DDmi = (Degree << 1) + 1;
1116 switch (Dimension) {
1118 psDD = psav + Degree;
1119 psDDmDim = psDD - 1;
1121 for (i = 0; i < Degree; i++) {
1127 for (j = Degm1; j >= i; j--) {
1129 *pole -= *tbis; *pole /= (knot[jDmi] - knot[j]);
1134 // Second phase, dependant of U
1137 for (i = 0; i < Degree; i++) {
1143 for (j = i; j >= 0; j--) {
1144 *pole += coef * (*tbis);
1149 // multiply by the degrees
1154 for (i = 1; i <= min; i++) {
1155 *pole *= coef; pole++;
1162 psDD = psav + (Degree << 1);
1163 psDDmDim = psDD - 2;
1165 for (i = 0; i < Degree; i++) {
1171 for (j = Degm1; j >= i; j--) {
1173 coef = 1. / (knot[jDmi] - knot[j]);
1174 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1175 *pole -= *tbis; *pole *= coef;
1180 // Second phase, dependant of U
1183 for (i = 0; i < Degree; i++) {
1189 for (j = i; j >= 0; j--) {
1190 *pole += coef * (*tbis); pole++; tbis++;
1191 *pole += coef * (*tbis);
1196 // multiply by the degrees
1201 for (i = 1; i <= min; i++) {
1202 *pole *= coef; pole++;
1203 *pole *= coef; pole++;
1210 psDD = psav + (Degree << 1) + Degree;
1211 psDDmDim = psDD - 3;
1213 for (i = 0; i < Degree; i++) {
1219 for (j = Degm1; j >= i; j--) {
1221 coef = 1. / (knot[jDmi] - knot[j]);
1222 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1223 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1224 *pole -= *tbis; *pole *= coef;
1229 // Second phase, dependant of U
1232 for (i = 0; i < Degree; i++) {
1238 for (j = i; j >= 0; j--) {
1239 *pole += coef * (*tbis); pole++; tbis++;
1240 *pole += coef * (*tbis); pole++; tbis++;
1241 *pole += coef * (*tbis);
1246 // multiply by the degrees
1251 for (i = 1; i <= min; i++) {
1252 *pole *= coef; pole++;
1253 *pole *= coef; pole++;
1254 *pole *= coef; pole++;
1261 psDD = psav + (Degree << 2);
1262 psDDmDim = psDD - 4;
1264 for (i = 0; i < Degree; i++) {
1270 for (j = Degm1; j >= i; j--) {
1272 coef = 1. / (knot[jDmi] - knot[j]);
1273 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1274 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1275 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1276 *pole -= *tbis; *pole *= coef;
1281 // Second phase, dependant of U
1284 for (i = 0; i < Degree; i++) {
1290 for (j = i; j >= 0; j--) {
1291 *pole += coef * (*tbis); pole++; tbis++;
1292 *pole += coef * (*tbis); pole++; tbis++;
1293 *pole += coef * (*tbis); pole++; tbis++;
1294 *pole += coef * (*tbis);
1299 // multiply by the degrees
1304 for (i = 1; i <= min; i++) {
1305 *pole *= coef; pole++;
1306 *pole *= coef; pole++;
1307 *pole *= coef; pole++;
1308 *pole *= coef; pole++;
1316 Standard_Integer Dim2 = Dimension << 1;
1317 psDD = psav + Degree * Dimension;
1318 psDDmDim = psDD - Dimension;
1320 for (i = 0; i < Degree; i++) {
1326 for (j = Degm1; j >= i; j--) {
1328 coef = 1. / (knot[jDmi] - knot[j]);
1330 for (k = 0; k < Dimension; k++) {
1331 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1337 // Second phase, dependant of U
1340 for (i = 0; i < Degree; i++) {
1343 tbis = pole + Dimension;
1346 for (j = i; j >= 0; j--) {
1348 for (k = 0; k < Dimension; k++) {
1349 *pole += coef * (*tbis); pole++; tbis++;
1355 // multiply by the degrees
1358 pole = psav + Dimension;
1360 for (i = 1; i <= min; i++) {
1362 for (k = 0; k < Dimension; k++) {
1363 *pole *= coef; pole++;
1372 //=======================================================================
1373 //function : BuildKnots
1375 //=======================================================================
1377 void BSplCLib::BuildKnots(const Standard_Integer Degree,
1378 const Standard_Integer Index,
1379 const Standard_Boolean Periodic,
1380 const TColStd_Array1OfReal& Knots,
1381 const TColStd_Array1OfInteger& Mults,
1384 Standard_Integer KLower = Knots.Lower();
1385 const Standard_Real * pkn = &Knots(KLower);
1387 Standard_Real *knot = &LK;
1388 if (&Mults == NULL) {
1391 Standard_Integer j = Index ;
1392 knot[0] = pkn[j]; j++;
1397 Standard_Integer j = Index - 1;
1398 knot[0] = pkn[j]; j++;
1399 knot[1] = pkn[j]; j++;
1400 knot[2] = pkn[j]; j++;
1405 Standard_Integer j = Index - 2;
1406 knot[0] = pkn[j]; j++;
1407 knot[1] = pkn[j]; j++;
1408 knot[2] = pkn[j]; j++;
1409 knot[3] = pkn[j]; j++;
1410 knot[4] = pkn[j]; j++;
1415 Standard_Integer j = Index - 3;
1416 knot[0] = pkn[j]; j++;
1417 knot[1] = pkn[j]; j++;
1418 knot[2] = pkn[j]; j++;
1419 knot[3] = pkn[j]; j++;
1420 knot[4] = pkn[j]; j++;
1421 knot[5] = pkn[j]; j++;
1422 knot[6] = pkn[j]; j++;
1427 Standard_Integer j = Index - 4;
1428 knot[0] = pkn[j]; j++;
1429 knot[1] = pkn[j]; j++;
1430 knot[2] = pkn[j]; j++;
1431 knot[3] = pkn[j]; j++;
1432 knot[4] = pkn[j]; j++;
1433 knot[5] = pkn[j]; j++;
1434 knot[6] = pkn[j]; j++;
1435 knot[7] = pkn[j]; j++;
1436 knot[8] = pkn[j]; j++;
1441 Standard_Integer j = Index - 5;
1442 knot[ 0] = pkn[j]; j++;
1443 knot[ 1] = pkn[j]; j++;
1444 knot[ 2] = pkn[j]; j++;
1445 knot[ 3] = pkn[j]; j++;
1446 knot[ 4] = pkn[j]; j++;
1447 knot[ 5] = pkn[j]; j++;
1448 knot[ 6] = pkn[j]; j++;
1449 knot[ 7] = pkn[j]; j++;
1450 knot[ 8] = pkn[j]; j++;
1451 knot[ 9] = pkn[j]; j++;
1452 knot[10] = pkn[j]; j++;
1457 Standard_Integer i,j;
1458 Standard_Integer Deg2 = Degree << 1;
1461 for (i = 0; i < Deg2; i++) {
1470 Standard_Integer Deg1 = Degree - 1;
1471 Standard_Integer KUpper = Knots.Upper();
1472 Standard_Integer MLower = Mults.Lower();
1473 Standard_Integer MUpper = Mults.Upper();
1474 const Standard_Integer * pmu = &Mults(MLower);
1476 Standard_Real dknot = 0;
1477 Standard_Integer ilow = Index , mlow = 0;
1478 Standard_Integer iupp = Index + 1, mupp = 0;
1479 Standard_Real loffset = 0., uoffset = 0.;
1480 Standard_Boolean getlow = Standard_True, getupp = Standard_True;
1482 dknot = pkn[KUpper] - pkn[KLower];
1483 if (iupp > MUpper) {
1488 // Find the knots around Index
1490 for (i = 0; i < Degree; i++) {
1493 if (mlow > pmu[ilow]) {
1496 getlow = (ilow >= MLower);
1497 if (Periodic && !getlow) {
1500 getlow = Standard_True;
1504 knot[Deg1 - i] = pkn[ilow] - loffset;
1508 if (mupp > pmu[iupp]) {
1511 getupp = (iupp <= MUpper);
1512 if (Periodic && !getupp) {
1515 getupp = Standard_True;
1519 knot[Degree + i] = pkn[iupp] + uoffset;
1525 //=======================================================================
1526 //function : PoleIndex
1528 //=======================================================================
1530 Standard_Integer BSplCLib::PoleIndex(const Standard_Integer Degree,
1531 const Standard_Integer Index,
1532 const Standard_Boolean Periodic,
1533 const TColStd_Array1OfInteger& Mults)
1535 Standard_Integer i, pindex = 0;
1537 for (i = Mults.Lower(); i <= Index; i++)
1540 pindex -= Mults(Mults.Lower());
1542 pindex -= Degree + 1;
1547 //=======================================================================
1548 //function : BuildBoor
1549 //purpose : builds the local array for boor
1550 //=======================================================================
1552 void BSplCLib::BuildBoor(const Standard_Integer Index,
1553 const Standard_Integer Length,
1554 const Standard_Integer Dimension,
1555 const TColStd_Array1OfReal& Poles,
1558 Standard_Real *poles = &LP;
1559 Standard_Integer i,k, ip = Poles.Lower() + Index * Dimension;
1561 for (i = 0; i < Length+1; i++) {
1563 for (k = 0; k < Dimension; k++) {
1564 poles[k] = Poles(ip);
1566 if (ip > Poles.Upper()) ip = Poles.Lower();
1568 poles += 2 * Dimension;
1572 //=======================================================================
1573 //function : BoorIndex
1575 //=======================================================================
1577 Standard_Integer BSplCLib::BoorIndex(const Standard_Integer Index,
1578 const Standard_Integer Length,
1579 const Standard_Integer Depth)
1581 if (Index <= Depth) return Index;
1582 if (Index <= Length) return 2 * Index - Depth;
1583 return Length + Index - Depth;
1586 //=======================================================================
1587 //function : GetPole
1589 //=======================================================================
1591 void BSplCLib::GetPole(const Standard_Integer Index,
1592 const Standard_Integer Length,
1593 const Standard_Integer Depth,
1594 const Standard_Integer Dimension,
1596 Standard_Integer& Position,
1597 TColStd_Array1OfReal& Pole)
1600 Standard_Real* pole = &LP + BoorIndex(Index,Length,Depth) * Dimension;
1602 for (k = 0; k < Dimension; k++) {
1603 Pole(Position) = pole[k];
1606 if (Position > Pole.Upper()) Position = Pole.Lower();
1609 //=======================================================================
1610 //function : PrepareInsertKnots
1612 //=======================================================================
1614 Standard_Boolean BSplCLib::PrepareInsertKnots
1615 (const Standard_Integer Degree,
1616 const Standard_Boolean Periodic,
1617 const TColStd_Array1OfReal& Knots,
1618 const TColStd_Array1OfInteger& Mults,
1619 const TColStd_Array1OfReal& AddKnots,
1620 const TColStd_Array1OfInteger& AddMults,
1621 Standard_Integer& NbPoles,
1622 Standard_Integer& NbKnots,
1623 const Standard_Real Tolerance,
1624 const Standard_Boolean Add)
1626 Standard_Boolean addflat = &AddMults == NULL;
1628 Standard_Integer first,last;
1630 first = Knots.Lower();
1631 last = Knots.Upper();
1634 first = FirstUKnotIndex(Degree,Mults);
1635 last = LastUKnotIndex(Degree,Mults);
1637 Standard_Real adeltaK1 = Knots(first)-AddKnots(AddKnots.Lower());
1638 Standard_Real adeltaK2 = AddKnots(AddKnots.Upper())-Knots(last);
1639 if (adeltaK1 > Tolerance) return Standard_False;
1640 if (adeltaK2 > Tolerance) return Standard_False;
1642 Standard_Integer sigma = 0, mult, amult;
1644 Standard_Integer k = Knots.Lower() - 1;
1645 Standard_Integer ak = AddKnots.Lower();
1647 if(Periodic && AddKnots.Length() > 1)
1649 //gka for case when segments was produced on full period only one knot
1650 //was added in the end of curve
1651 if(fabs(adeltaK1) <= Precision::PConfusion() &&
1652 fabs(adeltaK2) <= Precision::PConfusion())
1656 Standard_Real au,oldau = AddKnots(ak),Eps;
1658 while (ak <= AddKnots.Upper()) {
1660 if (au < oldau) return Standard_False;
1663 Eps = Max(Tolerance,Epsilon(au));
1665 while ((k < Knots.Upper()) && (Knots(k+1) - au <= Eps)) {
1671 if (addflat) amult = 1;
1672 else amult = Max(0,AddMults(ak));
1674 while ((ak < AddKnots.Upper()) &&
1675 (Abs(au - AddKnots(ak+1)) <= Eps)) {
1678 if (addflat) amult++;
1679 else amult += Max(0,AddMults(ak));
1684 if (Abs(au - Knots(k)) <= Eps) {
1685 // identic to existing knot
1688 if (mult + amult > Degree)
1689 amult = Max(0,Degree - mult);
1693 else if (amult > mult) {
1694 if (amult > Degree) amult = Degree;
1695 sigma += amult - mult;
1698 // on periodic curves if this is the last knot
1699 // the multiplicity is added twice to account for the first knot
1700 if (k == Knots.Upper() && Periodic) {
1704 sigma += amult - mult;
1709 // not identic to existing knot
1711 if (amult > Degree) amult = Degree;
1720 // count the last knots
1721 while (k < Knots.Upper()) {
1728 //for periodic B-Spline the requirement is that multiplicites of the first
1729 //and last knots must be equal (see Geom_BSplineCurve constructor for
1731 //respectively AddMults() must meet this requirement if AddKnots() contains
1732 //knot(s) coincident with first or last
1733 NbPoles = sigma - Mults(Knots.Upper());
1736 NbPoles = sigma - Degree - 1;
1739 return Standard_True;
1742 //=======================================================================
1744 //purpose : copy reals from an array to an other
1746 // NbValues are copied from OldPoles(OldFirst)
1747 // to NewPoles(NewFirst)
1749 // Periodicity is handled.
1750 // OldFirst and NewFirst are updated
1751 // to the position after the last copied pole.
1753 //=======================================================================
1755 static void Copy(const Standard_Integer NbPoles,
1756 Standard_Integer& OldFirst,
1757 const TColStd_Array1OfReal& OldPoles,
1758 Standard_Integer& NewFirst,
1759 TColStd_Array1OfReal& NewPoles)
1761 // reset the index in the range for periodicity
1763 OldFirst = OldPoles.Lower() +
1764 (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1);
1766 NewFirst = NewPoles.Lower() +
1767 (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1);
1772 for (i = 1; i <= NbPoles; i++) {
1773 NewPoles(NewFirst) = OldPoles(OldFirst);
1775 if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower();
1777 if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower();
1781 //=======================================================================
1782 //function : InsertKnots
1783 //purpose : insert an array of knots and multiplicities
1784 //=======================================================================
1786 void BSplCLib::InsertKnots
1787 (const Standard_Integer Degree,
1788 const Standard_Boolean Periodic,
1789 const Standard_Integer Dimension,
1790 const TColStd_Array1OfReal& Poles,
1791 const TColStd_Array1OfReal& Knots,
1792 const TColStd_Array1OfInteger& Mults,
1793 const TColStd_Array1OfReal& AddKnots,
1794 const TColStd_Array1OfInteger& AddMults,
1795 TColStd_Array1OfReal& NewPoles,
1796 TColStd_Array1OfReal& NewKnots,
1797 TColStd_Array1OfInteger& NewMults,
1798 const Standard_Real Tolerance,
1799 const Standard_Boolean Add)
1801 Standard_Boolean addflat = &AddMults == NULL;
1803 Standard_Integer i,k,mult,firstmult;
1804 Standard_Integer index,kn,curnk,curk;
1805 Standard_Integer p,np, curp, curnp, length, depth;
1807 Standard_Integer need;
1810 // -------------------
1811 // create local arrays
1812 // -------------------
1814 Standard_Real *knots = new Standard_Real[2*Degree];
1815 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
1817 //----------------------------
1818 // loop on the knots to insert
1819 //----------------------------
1821 curk = Knots.Lower()-1; // current position in Knots
1822 curnk = NewKnots.Lower()-1; // current position in NewKnots
1823 curp = Poles.Lower(); // current position in Poles
1824 curnp = NewPoles.Lower(); // current position in NewPoles
1826 // NewKnots, NewMults, NewPoles contains the current state of the curve
1828 // index is the first pole of the current curve for insertion schema
1830 if (Periodic) index = -Mults(Mults.Lower());
1831 else index = -Degree-1;
1833 // on Periodic curves the first knot and the last knot are inserted later
1834 // (they are the same knot)
1835 firstmult = 0; // multiplicity of the first-last knot for periodic
1838 // kn current knot to insert in AddKnots
1840 for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) {
1843 Eps = Max(Tolerance,Epsilon(u));
1845 //-----------------------------------
1846 // find the position in the old knots
1847 // and copy to the new knots
1848 //-----------------------------------
1850 while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) {
1852 NewKnots(curnk) = Knots(curk);
1853 index += NewMults(curnk) = Mults(curk);
1856 //-----------------------------------
1857 // Slice the knots and the mults
1858 // to the current size of the new curve
1859 //-----------------------------------
1861 i = curnk + Knots.Upper() - curk;
1862 TColStd_Array1OfReal nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i);
1863 TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i);
1865 //-----------------------------------
1866 // copy enough knots
1867 // to compute the insertion schema
1868 //-----------------------------------
1874 while (mult < Degree && k < Knots.Upper()) {
1876 nknots(i) = Knots(k);
1877 mult += nmults(i) = Mults(k);
1880 // copy knots at the end for periodic curve
1886 while (mult < Degree && i > curnk) {
1887 nknots(i) = Knots(k);
1888 mult += nmults(i) = Mults(k);
1892 nmults(nmults.Upper()) = nmults(nmults.Lower());
1897 //------------------------------------
1898 // do the boor scheme on the new curve
1899 // to insert the new knot
1900 //------------------------------------
1902 Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps);
1904 if (sameknot) length = Max(0,Degree - NewMults(curnk));
1905 else length = Degree;
1907 if (addflat) depth = 1;
1908 else depth = Min(Degree,AddMults(kn));
1912 if ((NewMults(curnk) + depth) > Degree)
1913 depth = Degree - NewMults(curnk);
1916 depth = Max(0,depth-NewMults(curnk));
1920 // on periodic curve the first and last knot are delayed to the end
1921 if (curk == Knots.Lower() || (curk == Knots.Upper())) {
1927 if (depth <= 0) continue;
1929 BuildKnots(Degree,curnk,Periodic,nknots,nmults,*knots);
1933 need = NewPoles.Lower()+(index+length+1)*Dimension - curnp;
1934 need = Min(need,Poles.Upper() - curp + 1);
1938 Copy(need,p,Poles,np,NewPoles);
1942 // slice the poles to the current number of poles in case of periodic
1943 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1945 BuildBoor(index,length,Dimension,npoles,*poles);
1946 BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length);
1948 //-------------------
1949 // copy the new poles
1950 //-------------------
1952 curnp += depth * Dimension; // number of poles is increased by depth
1953 TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1954 np = NewKnots.Lower()+(index+1)*Dimension;
1956 for (i = 1; i <= length + depth; i++)
1957 GetPole(i,length,depth,Dimension,*poles,np,ThePoles);
1959 //-------------------
1961 //-------------------
1965 NewMults(curnk) += depth;
1969 NewKnots(curnk) = u;
1970 NewMults(curnk) = depth;
1974 //------------------------------
1975 // copy the last poles and knots
1976 //------------------------------
1978 Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles);
1980 while (curk < Knots.Upper()) {
1982 NewKnots(curnk) = Knots(curk);
1983 NewMults(curnk) = Mults(curk);
1986 //------------------------------
1987 // process the first-last knot
1988 // on periodic curves
1989 //------------------------------
1991 if (firstmult > 0) {
1992 curnk = NewKnots.Lower();
1993 if (NewMults(curnk) + firstmult > Degree) {
1994 firstmult = Degree - NewMults(curnk);
1996 if (firstmult > 0) {
1998 length = Degree - NewMults(curnk);
2001 BuildKnots(Degree,curnk,Periodic,NewKnots,NewMults,*knots);
2002 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),
2004 NewPoles.Upper()-depth*Dimension);
2005 BuildBoor(0,length,Dimension,npoles,*poles);
2006 BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length);
2008 //---------------------------
2009 // copy the new poles
2010 // but rotate them with depth
2011 //---------------------------
2013 np = NewPoles.Lower();
2015 for (i = depth; i < length + depth; i++)
2016 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2018 np = NewPoles.Upper() - depth*Dimension + 1;
2020 for (i = 0; i < depth; i++)
2021 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2023 NewMults(NewMults.Lower()) += depth;
2024 NewMults(NewMults.Upper()) += depth;
2027 // free local arrays
2032 //=======================================================================
2033 //function : RemoveKnot
2035 //=======================================================================
2037 Standard_Boolean BSplCLib::RemoveKnot
2038 (const Standard_Integer Index,
2039 const Standard_Integer Mult,
2040 const Standard_Integer Degree,
2041 const Standard_Boolean Periodic,
2042 const Standard_Integer Dimension,
2043 const TColStd_Array1OfReal& Poles,
2044 const TColStd_Array1OfReal& Knots,
2045 const TColStd_Array1OfInteger& Mults,
2046 TColStd_Array1OfReal& NewPoles,
2047 TColStd_Array1OfReal& NewKnots,
2048 TColStd_Array1OfInteger& NewMults,
2049 const Standard_Real Tolerance)
2051 Standard_Integer index,i,j,k,p,np;
2053 Standard_Integer TheIndex = Index;
2056 Standard_Integer first,last;
2058 first = Knots.Lower();
2059 last = Knots.Upper();
2062 first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1;
2063 last = BSplCLib::LastUKnotIndex(Degree,Mults) - 1;
2065 if (Index < first) return Standard_False;
2066 if (Index > last) return Standard_False;
2068 if ( Periodic && (Index == first)) TheIndex = last;
2070 Standard_Integer depth = Mults(TheIndex) - Mult;
2071 Standard_Integer length = Degree - Mult;
2073 // -------------------
2074 // create local arrays
2075 // -------------------
2077 Standard_Real *knots = new Standard_Real[4*Degree];
2078 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
2081 // ------------------------------------
2082 // build the knots for anti Boor Scheme
2083 // ------------------------------------
2085 // the new sequence of knots
2086 // is obtained from the knots at Index-1 and Index
2088 BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,Mults,*knots);
2089 index = PoleIndex(Degree,TheIndex-1,Periodic,Mults);
2090 BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,Mults,knots[2*Degree]);
2094 for (i = 0; i < Degree - Mult; i++)
2095 knots[i] = knots[i+Mult];
2097 for (i = Degree-Mult; i < 2*Degree; i++)
2098 knots[i] = knots[2*Degree+i];
2101 // ------------------------------------
2102 // build the poles for anti Boor Scheme
2103 // ------------------------------------
2105 p = Poles.Lower()+index * Dimension;
2107 for (i = 0; i <= length + depth; i++) {
2108 j = Dimension * BoorIndex(i,length,depth);
2110 for (k = 0; k < Dimension; k++) {
2111 poles[j+k] = Poles(p+k);
2114 if (p > Poles.Upper()) p = Poles.Lower();
2122 Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots,
2124 depth,length,Tolerance);
2135 np = NewPoles.Lower();
2137 // unmodified poles before
2138 Copy((index+1)*Dimension,p,Poles,np,NewPoles);
2143 for (i = 1; i <= length; i++)
2144 BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles);
2145 p += (length + depth) * Dimension ;
2147 // unmodified poles after
2148 if (p != Poles.Lower()) {
2149 i = Poles.Upper() - p + 1;
2150 Copy(i,p,Poles,np,NewPoles);
2158 NewMults(TheIndex) = Mult;
2160 if (TheIndex == first) NewMults(last) = Mult;
2161 if (TheIndex == last) NewMults(first) = Mult;
2165 if (!Periodic || (TheIndex != first && TheIndex != last)) {
2167 for (i = Knots.Lower(); i < TheIndex; i++) {
2168 NewKnots(i) = Knots(i);
2169 NewMults(i) = Mults(i);
2172 for (i = TheIndex+1; i <= Knots.Upper(); i++) {
2173 NewKnots(i-1) = Knots(i);
2174 NewMults(i-1) = Mults(i);
2178 // The interesting case of a Periodic curve
2179 // where the first and last knot is removed.
2181 for (i = first; i < last-1; i++) {
2182 NewKnots(i) = Knots(i+1);
2183 NewMults(i) = Mults(i+1);
2185 NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first);
2186 NewMults(last-1) = NewMults(first);
2192 // free local arrays
2199 //=======================================================================
2200 //function : IncreaseDegreeCountKnots
2202 //=======================================================================
2204 Standard_Integer BSplCLib::IncreaseDegreeCountKnots
2205 (const Standard_Integer Degree,
2206 const Standard_Integer NewDegree,
2207 const Standard_Boolean Periodic,
2208 const TColStd_Array1OfInteger& Mults)
2210 if (Periodic) return Mults.Length();
2211 Standard_Integer f = FirstUKnotIndex(Degree,Mults);
2212 Standard_Integer l = LastUKnotIndex(Degree,Mults);
2213 Standard_Integer m,i,removed = 0, step = NewDegree - Degree;
2216 m = Degree + (f - i + 1) * step + 1;
2218 while (m > NewDegree+1) {
2220 m -= Mults(i) + step;
2223 if (m < NewDegree+1) removed--;
2226 m = Degree + (i - l + 1) * step + 1;
2228 while (m > NewDegree+1) {
2230 m -= Mults(i) + step;
2233 if (m < NewDegree+1) removed--;
2235 return Mults.Length() - removed;
2238 //=======================================================================
2239 //function : IncreaseDegree
2241 //=======================================================================
2243 void BSplCLib::IncreaseDegree
2244 (const Standard_Integer Degree,
2245 const Standard_Integer NewDegree,
2246 const Standard_Boolean Periodic,
2247 const Standard_Integer Dimension,
2248 const TColStd_Array1OfReal& Poles,
2249 const TColStd_Array1OfReal& Knots,
2250 const TColStd_Array1OfInteger& Mults,
2251 TColStd_Array1OfReal& NewPoles,
2252 TColStd_Array1OfReal& NewKnots,
2253 TColStd_Array1OfInteger& NewMults)
2255 // Degree elevation of a BSpline Curve
2257 // This algorithms loops on degree incrementation from Degree to NewDegree.
2258 // The variable curDeg is the current degree to increment.
2260 // Before degree incrementations a "working curve" is created.
2261 // It has the same knot, poles and multiplicities.
2263 // If the curve is periodic knots are added on the working curve before
2264 // and after the existing knots to make it a non-periodic curves.
2265 // The poles are also copied.
2267 // The first and last multiplicity of the working curve are set to Degree+1,
2268 // null poles are added if necessary.
2270 // Then the degree is incremented on the working curve.
2271 // The knots are unchanged but all multiplicities will be incremented.
2273 // Each degree incrementation is achieved by averaging curDeg+1 curves.
2275 // See : Degree elevation of B-spline curves
2276 // Hartmut PRAUTZSCH
2280 //-------------------------
2281 // create the working curve
2282 //-------------------------
2284 Standard_Integer i,k,f,l,m,pf,pl,firstknot;
2286 pf = 0; // number of null poles added at beginning
2287 pl = 0; // number of null poles added at end
2289 Standard_Integer nbwknots = Knots.Length();
2290 f = FirstUKnotIndex(Degree,Mults);
2291 l = LastUKnotIndex (Degree,Mults);
2294 // Periodic curves are transformed in non-periodic curves
2296 nbwknots += f - Mults.Lower();
2300 for (i = Mults.Lower(); i <= f; i++)
2303 nbwknots += Mults.Upper() - l;
2307 for (i = l; i <= Mults.Upper(); i++)
2311 // copy the knots and multiplicities
2312 TColStd_Array1OfReal wknots(1,nbwknots);
2313 TColStd_Array1OfInteger wmults(1,nbwknots);
2319 // copy the knots for a periodic curve
2320 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2323 for (k = l; k < Knots.Upper(); k++) {
2325 wknots(i) = Knots(k) - period;
2326 wmults(i) = Mults(k);
2329 for (k = Knots.Lower(); k <= Knots.Upper(); k++) {
2331 wknots(i) = Knots(k);
2332 wmults(i) = Mults(k);
2335 for (k = Knots.Lower()+1; k <= f; k++) {
2337 wknots(i) = Knots(k) + period;
2338 wmults(i) = Mults(k);
2342 // set the first and last mults to Degree+1
2343 // and add null poles
2345 pf += Degree + 1 - wmults(1);
2346 wmults(1) = Degree + 1;
2347 pl += Degree + 1 - wmults(nbwknots);
2348 wmults(nbwknots) = Degree + 1;
2350 //---------------------------
2351 // poles of the working curve
2352 //---------------------------
2354 Standard_Integer nbwpoles = 0;
2356 for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i);
2357 nbwpoles -= Degree + 1;
2359 // we provide space for degree elevation
2360 TColStd_Array1OfReal
2361 wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension);
2363 for (i = 1; i <= pf * Dimension; i++)
2368 for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) {
2369 wpoles(i) = Poles(k);
2371 if (k > Poles.Upper()) k = Poles.Lower();
2374 for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++)
2378 //-----------------------------------------------------------
2379 // Declare the temporary arrays used in degree incrementation
2380 //-----------------------------------------------------------
2382 Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree);
2383 // Arrays for storing the temporary curves
2384 TColStd_Array1OfReal tempc1(1,nbwp * Dimension);
2385 TColStd_Array1OfReal tempc2(1,nbwp * Dimension);
2387 // Array for storing the knots to insert
2388 TColStd_Array1OfReal iknots(1,nbwknots);
2390 // Arrays for receiving the knots after insertion
2391 TColStd_Array1OfReal nknots(1,nbwknots);
2395 //------------------------------
2396 // Loop on degree incrementation
2397 //------------------------------
2399 Standard_Integer step,curDeg;
2400 Standard_Integer nbp = nbwpoles;
2403 for (curDeg = Degree; curDeg < NewDegree; curDeg++) {
2405 nbp = nbwp; // current number of poles
2406 nbwp = nbp + nbwknots - 1; // new number of poles
2408 // For the averaging
2409 TColStd_Array1OfReal nwpoles(1,nbwp * Dimension);
2410 nwpoles.Init(0.0e0) ;
2413 for (step = 0; step <= curDeg; step++) {
2415 // Compute the bspline of rank step.
2417 // if not the first time, decrement the multiplicities back
2419 for (i = 1; i <= nbwknots; i++)
2423 // Poles are the current poles
2424 // but the poles congruent to step are duplicated.
2426 Standard_Integer offset = 0;
2428 for (i = 0; i < nbp; i++) {
2431 for (k = 0; k < Dimension; k++) {
2432 tempc1((offset-1)*Dimension+k+1) =
2433 wpoles(NewPoles.Lower()+i*Dimension+k);
2435 if (i % (curDeg+1) == step) {
2438 for (k = 0; k < Dimension; k++) {
2439 tempc1((offset-1)*Dimension+k+1) =
2440 wpoles(NewPoles.Lower()+i*Dimension+k);
2445 // Knots multiplicities are increased
2446 // For each knot where the sum of multiplicities is congruent to step
2448 Standard_Integer stepmult = step+1;
2449 Standard_Integer nbknots = 0;
2450 Standard_Integer smult = 0;
2452 for (k = 1; k <= nbwknots; k++) {
2454 if (smult >= stepmult) {
2455 // this knot is increased
2456 stepmult += curDeg+1;
2460 // this knot is inserted
2462 iknots(nbknots) = wknots(k);
2466 // the curve is obtained by inserting the knots
2467 // to raise the multiplicities
2469 // we build "slices" of the arrays to set the correct size
2471 TColStd_Array1OfReal aknots(iknots(1),1,nbknots);
2472 TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension);
2473 TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp * Dimension);
2474 // InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2475 // aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.));
2477 InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2478 aknots,NoMults(),ncurve,nknots,wmults,0.0);
2480 // add to the average
2482 for (i = 1; i <= nbwp * Dimension; i++)
2483 nwpoles(i) += ncurve(i);
2486 // add to the average
2488 for (i = 1; i <= nbwp * Dimension; i++)
2489 nwpoles(i) += tempc1(i);
2493 // The result is the average
2495 for (i = 1; i <= nbwp * Dimension; i++) {
2496 wpoles(i) = nwpoles(i) / (curDeg+1);
2504 // index in new knots of the first knot of the curve
2506 firstknot = Mults.Upper() - l + 1;
2510 // the new curve starts at index firstknot
2511 // so we must remove knots until the sum of multiplicities
2512 // from the first to the start is NewDegree+1
2514 // m is the current sum of multiplicities
2517 for (k = 1; k <= firstknot; k++)
2520 // compute the new first knot (k), pf will be the index of the first pole
2524 while (m > NewDegree+1) {
2529 if (m < NewDegree+1) {
2531 wmults(k) += m - NewDegree - 1;
2532 pf += m - NewDegree - 1;
2535 // on a periodic curve the knots start with firstknot
2541 for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) {
2542 NewKnots(i) = wknots(k);
2543 NewMults(i) = wmults(k);
2550 for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) {
2552 NewPoles(i) = wpoles(pf);
2556 //=======================================================================
2557 //function : PrepareUnperiodize
2559 //=======================================================================
2561 void BSplCLib::PrepareUnperiodize
2562 (const Standard_Integer Degree,
2563 const TColStd_Array1OfInteger& Mults,
2564 Standard_Integer& NbKnots,
2565 Standard_Integer& NbPoles)
2568 // initialize NbKnots and NbPoles
2569 NbKnots = Mults.Length();
2570 NbPoles = - Degree - 1;
2572 for (i = Mults.Lower(); i <= Mults.Upper(); i++)
2573 NbPoles += Mults(i);
2575 Standard_Integer sigma, k;
2576 // Add knots at the beginning of the curve to raise Multiplicities
2578 sigma = Mults(Mults.Lower());
2579 k = Mults.Upper() - 1;
2581 while ( sigma < Degree + 1) {
2583 NbPoles += Mults(k);
2587 // We must add exactly until Degree + 1 ->
2588 // Supress the excedent.
2589 if ( sigma > Degree + 1)
2590 NbPoles -= sigma - Degree - 1;
2592 // Add knots at the end of the curve to raise Multiplicities
2594 sigma = Mults(Mults.Upper());
2595 k = Mults.Lower() + 1;
2597 while ( sigma < Degree + 1) {
2599 NbPoles += Mults(k);
2603 // We must add exactly until Degree + 1 ->
2604 // Supress the excedent.
2605 if ( sigma > Degree + 1)
2606 NbPoles -= sigma - Degree - 1;
2609 //=======================================================================
2610 //function : Unperiodize
2612 //=======================================================================
2614 void BSplCLib::Unperiodize
2615 (const Standard_Integer Degree,
2616 const Standard_Integer , // Dimension,
2617 const TColStd_Array1OfInteger& Mults,
2618 const TColStd_Array1OfReal& Knots,
2619 const TColStd_Array1OfReal& Poles,
2620 TColStd_Array1OfInteger& NewMults,
2621 TColStd_Array1OfReal& NewKnots,
2622 TColStd_Array1OfReal& NewPoles)
2624 Standard_Integer sigma, k, index = 0;
2625 // evaluation of index : number of knots to insert before knot(1) to
2626 // raise sum of multiplicities to <Degree + 1>
2627 sigma = Mults(Mults.Lower());
2628 k = Mults.Upper() - 1;
2630 while ( sigma < Degree + 1) {
2636 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2638 // set the 'interior' knots;
2640 for ( k = 1; k <= Knots.Length(); k++) {
2641 NewKnots ( k + index ) = Knots( k);
2642 NewMults ( k + index ) = Mults( k);
2645 // set the 'starting' knots;
2647 for ( k = 1; k <= index; k++) {
2648 NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period;
2649 NewMults( k) = NewMults( k + Knots.Length() - 1);
2651 NewMults( 1) -= sigma - Degree -1;
2653 // set the 'ending' knots;
2654 sigma = NewMults( index + Knots.Length() );
2656 for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) {
2657 NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period;
2658 NewMults( k) = NewMults( k - Knots.Length() + 1);
2659 sigma += NewMults( k - Knots.Length() + 1);
2661 NewMults(NewMults.Length()) -= sigma - Degree - 1;
2663 for ( k = 1 ; k <= NewPoles.Length(); k++) {
2664 NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1);
2668 //=======================================================================
2669 //function : PrepareTrimming
2671 //=======================================================================
2673 void BSplCLib::PrepareTrimming(const Standard_Integer Degree,
2674 const Standard_Boolean Periodic,
2675 const TColStd_Array1OfReal& Knots,
2676 const TColStd_Array1OfInteger& Mults,
2677 const Standard_Real U1,
2678 const Standard_Real U2,
2679 Standard_Integer& NbKnots,
2680 Standard_Integer& NbPoles)
2683 Standard_Real NewU1, NewU2;
2684 Standard_Integer index1 = 0, index2 = 0;
2686 // Eval index1, index2 : position of U1 and U2 in the Array Knots
2687 // such as Knots(index1-1) <= U1 < Knots(index1)
2688 // Knots(index2-1) <= U2 < Knots(index2)
2689 LocateParameter( Degree, Knots, Mults, U1, Periodic,
2690 Knots.Lower(), Knots.Upper(), index1, NewU1);
2691 LocateParameter( Degree, Knots, Mults, U2, Periodic,
2692 Knots.Lower(), Knots.Upper(), index2, NewU2);
2694 if ( Abs(Knots(index2) - U2) <= Epsilon( U1))
2698 NbKnots = index2 - index1 + 3;
2701 NbPoles = Degree + 1;
2703 for ( i = index1; i <= index2; i++)
2704 NbPoles += Mults(i);
2707 //=======================================================================
2708 //function : Trimming
2710 //=======================================================================
2712 void BSplCLib::Trimming(const Standard_Integer Degree,
2713 const Standard_Boolean Periodic,
2714 const Standard_Integer Dimension,
2715 const TColStd_Array1OfReal& Knots,
2716 const TColStd_Array1OfInteger& Mults,
2717 const TColStd_Array1OfReal& Poles,
2718 const Standard_Real U1,
2719 const Standard_Real U2,
2720 TColStd_Array1OfReal& NewKnots,
2721 TColStd_Array1OfInteger& NewMults,
2722 TColStd_Array1OfReal& NewPoles)
2724 Standard_Integer i, nbpoles, nbknots;
2725 Standard_Real kk[2];
2726 Standard_Integer mm[2];
2727 TColStd_Array1OfReal K( kk[0], 1, 2 );
2728 TColStd_Array1OfInteger M( mm[0], 1, 2 );
2730 K(1) = U1; K(2) = U2;
2731 mm[0] = mm[1] = Degree;
2732 if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, M,
2733 nbpoles, nbknots, Epsilon( U1), 0))
2734 Standard_OutOfRange::Raise();
2736 TColStd_Array1OfReal TempPoles(1, nbpoles*Dimension);
2737 TColStd_Array1OfReal TempKnots(1, nbknots);
2738 TColStd_Array1OfInteger TempMults(1, nbknots);
2741 // do not allow the multiplicities to Add : they must be less than Degree
2743 InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults,
2744 K, M, TempPoles, TempKnots, TempMults, Epsilon(U1),
2747 // find in TempPoles the index of the pole corresponding to U1
2748 Standard_Integer Kindex = 0, Pindex;
2749 Standard_Real NewU1;
2750 LocateParameter( Degree, TempKnots, TempMults, U1, Periodic,
2751 TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1);
2752 Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults);
2753 Pindex *= Dimension;
2755 for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i);
2757 for ( i = 1; i <= NewKnots.Length(); i++) {
2758 NewKnots(i) = TempKnots( Kindex+i-1);
2759 NewMults(i) = TempMults( Kindex+i-1);
2761 NewMults(1) = Min(Degree, NewMults(1)) + 1 ;
2762 NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ;
2765 //=======================================================================
2766 //function : Solves a LU factored Matrix
2768 //=======================================================================
2771 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2772 const Standard_Integer UpperBandWidth,
2773 const Standard_Integer LowerBandWidth,
2774 const Standard_Integer ArrayDimension,
2775 Standard_Real& Array)
2777 Standard_Integer ii,
2784 Standard_Real *PolesArray = &Array ;
2785 Standard_Real Inverse ;
2788 if (Matrix.LowerCol() != 1 ||
2789 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2794 for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow() ; ii++) {
2795 MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ?
2796 ii - LowerBandWidth : Matrix.LowerRow()) ;
2798 for ( jj = MinIndex ; jj < ii ; jj++) {
2800 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2801 PolesArray[(ii-1) * ArrayDimension + kk] +=
2802 PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2807 for (ii = Matrix.UpperRow() ; ii >= Matrix.LowerRow() ; ii--) {
2808 MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ?
2809 ii + UpperBandWidth : Matrix.UpperRow()) ;
2811 for (jj = MaxIndex ; jj > ii ; jj--) {
2813 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2814 PolesArray[(ii-1) * ArrayDimension + kk] -=
2815 PolesArray[(jj - 1) * ArrayDimension + kk] *
2816 Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2820 //fixing a bug PRO18577 to avoid divizion by zero
2822 Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ;
2823 Standard_Real Toler = 1.0e-16;
2824 if ( Abs(divizor) > Toler )
2825 Inverse = 1.0e0 / divizor ;
2828 // cout << " BSplCLib::SolveBandedSystem() : zero determinant " << endl;
2833 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2834 PolesArray[(ii-1) * ArrayDimension + kk] *= Inverse ;
2838 return (ReturnCode) ;
2841 //=======================================================================
2842 //function : Solves a LU factored Matrix
2844 //=======================================================================
2847 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2848 const Standard_Integer UpperBandWidth,
2849 const Standard_Integer LowerBandWidth,
2850 const Standard_Boolean HomogeneousFlag,
2851 const Standard_Integer ArrayDimension,
2852 Standard_Real& Poles,
2853 Standard_Real& Weights)
2855 Standard_Integer ii,
2860 Standard_Real Inverse,
2861 *PolesArray = &Poles,
2862 *WeightsArray = &Weights ;
2864 if (Matrix.LowerCol() != 1 ||
2865 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2869 if (HomogeneousFlag == Standard_False) {
2871 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) {
2873 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2874 PolesArray[ii * ArrayDimension + kk] *=
2880 BSplCLib::SolveBandedSystem(Matrix,
2885 if (ErrorCode != 0) {
2890 BSplCLib::SolveBandedSystem(Matrix,
2895 if (ErrorCode != 0) {
2899 if (HomogeneousFlag == Standard_False) {
2901 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) {
2902 Inverse = 1.0e0 / WeightsArray[ii] ;
2904 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2905 PolesArray[ii * ArrayDimension + kk] *= Inverse ;
2909 FINISH : return (ReturnCode) ;
2912 //=======================================================================
2913 //function : BuildSchoenbergPoints
2915 //=======================================================================
2917 void BSplCLib::BuildSchoenbergPoints(const Standard_Integer Degree,
2918 const TColStd_Array1OfReal& FlatKnots,
2919 TColStd_Array1OfReal& Parameters)
2921 Standard_Integer ii,
2923 Standard_Real Inverse ;
2924 Inverse = 1.0e0 / (Standard_Real)Degree ;
2926 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
2927 Parameters(ii) = 0.0e0 ;
2929 for (jj = 1 ; jj <= Degree ; jj++) {
2930 Parameters(ii) += FlatKnots(jj + ii) ;
2932 Parameters(ii) *= Inverse ;
2936 //=======================================================================
2937 //function : Interpolate
2939 //=======================================================================
2941 void BSplCLib::Interpolate(const Standard_Integer Degree,
2942 const TColStd_Array1OfReal& FlatKnots,
2943 const TColStd_Array1OfReal& Parameters,
2944 const TColStd_Array1OfInteger& ContactOrderArray,
2945 const Standard_Integer ArrayDimension,
2946 Standard_Real& Poles,
2947 Standard_Integer& InversionProblem)
2949 Standard_Integer ErrorCode,
2952 // Standard_Real *PolesArray = &Poles ;
2953 math_Matrix InterpolationMatrix(1, Parameters.Length(),
2954 1, 2 * Degree + 1) ;
2956 BSplCLib::BuildBSpMatrix(Parameters,
2960 InterpolationMatrix,
2963 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
2966 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
2970 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
2973 BSplCLib::SolveBandedSystem(InterpolationMatrix,
2979 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
2982 //=======================================================================
2983 //function : Interpolate
2985 //=======================================================================
2987 void BSplCLib::Interpolate(const Standard_Integer Degree,
2988 const TColStd_Array1OfReal& FlatKnots,
2989 const TColStd_Array1OfReal& Parameters,
2990 const TColStd_Array1OfInteger& ContactOrderArray,
2991 const Standard_Integer ArrayDimension,
2992 Standard_Real& Poles,
2993 Standard_Real& Weights,
2994 Standard_Integer& InversionProblem)
2996 Standard_Integer ErrorCode,
3000 math_Matrix InterpolationMatrix(1, Parameters.Length(),
3001 1, 2 * Degree + 1) ;
3003 BSplCLib::BuildBSpMatrix(Parameters,
3007 InterpolationMatrix,
3010 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3013 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3017 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3020 BSplCLib::SolveBandedSystem(InterpolationMatrix,
3028 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
3031 //=======================================================================
3032 //function : Evaluates a Bspline function : uses the ExtrapMode
3033 //purpose : the function is extrapolated using the Taylor expansion
3034 // of degree ExtrapMode[0] to the left and the Taylor
3035 // expansion of degree ExtrapMode[1] to the right
3036 // this evaluates the numerator by multiplying by the weights
3037 // and evaluating it but does not call RationalDerivatives after
3038 //=======================================================================
3041 (const Standard_Real Parameter,
3042 const Standard_Boolean PeriodicFlag,
3043 const Standard_Integer DerivativeRequest,
3044 Standard_Integer& ExtrapMode,
3045 const Standard_Integer Degree,
3046 const TColStd_Array1OfReal& FlatKnots,
3047 const Standard_Integer ArrayDimension,
3048 Standard_Real& Poles,
3049 Standard_Real& Weights,
3050 Standard_Real& PolesResults,
3051 Standard_Real& WeightsResults)
3053 Standard_Integer ii,
3062 ExtrapolatingFlag[2],
3066 FirstNonZeroBsplineIndex,
3067 LocalRequest = DerivativeRequest ;
3068 Standard_Real *PResultArray,
3076 PolesArray = &Poles ;
3077 WeightsArray = &Weights ;
3078 ExtrapModeArray = &ExtrapMode ;
3079 PResultArray = &PolesResults ;
3080 WResultArray = &WeightsResults ;
3081 LocalParameter = Parameter ;
3082 ExtrapolatingFlag[0] =
3083 ExtrapolatingFlag[1] = 0 ;
3085 // check if we are extrapolating to a degree which is smaller than
3086 // the degree of the Bspline
3089 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3091 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3092 LocalParameter -= Period ;
3095 while (LocalParameter < FlatKnots(2)) {
3096 LocalParameter += Period ;
3099 if (Parameter < FlatKnots(2) &&
3100 LocalRequest < ExtrapModeArray[0] &&
3101 ExtrapModeArray[0] < Degree) {
3102 LocalRequest = ExtrapModeArray[0] ;
3103 LocalParameter = FlatKnots(2) ;
3104 ExtrapolatingFlag[0] = 1 ;
3106 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3107 LocalRequest < ExtrapModeArray[1] &&
3108 ExtrapModeArray[1] < Degree) {
3109 LocalRequest = ExtrapModeArray[1] ;
3110 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3111 ExtrapolatingFlag[1] = 1 ;
3113 Delta = Parameter - LocalParameter ;
3114 if (LocalRequest >= Order) {
3115 LocalRequest = Degree ;
3118 Modulus = FlatKnots.Length() - Degree -1 ;
3121 Modulus = FlatKnots.Length() - Degree ;
3124 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3126 BSplCLib::EvalBsplineBasis(1,
3131 FirstNonZeroBsplineIndex,
3133 if (ErrorCode != 0) {
3137 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3141 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3142 Index1 = FirstNonZeroBsplineIndex ;
3144 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3145 PResultArray[Index + kk] = 0.0e0 ;
3147 WResultArray[Index] = 0.0e0 ;
3149 for (jj = 1 ; jj <= Order ; jj++) {
3151 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3152 PResultArray[Index + kk] +=
3153 PolesArray[(Index1-1) * ArrayDimension + kk]
3154 * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3156 WResultArray[Index2] += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3158 Index1 = Index1 % Modulus ;
3161 Index += ArrayDimension ;
3167 // store Taylor expansion in LocalRealArray
3169 NewRequest = DerivativeRequest ;
3170 if (NewRequest > Degree) {
3171 NewRequest = Degree ;
3173 BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension);
3177 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3178 Index1 = FirstNonZeroBsplineIndex ;
3180 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3181 LocalRealArray[Index + kk] = 0.0e0 ;
3184 for (jj = 1 ; jj <= Order ; jj++) {
3186 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3187 LocalRealArray[Index + kk] +=
3188 PolesArray[(Index1-1)*ArrayDimension + kk] *
3189 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3191 Index1 = Index1 % Modulus ;
3195 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3196 LocalRealArray[Index + kk] *= Inverse ;
3198 Index += ArrayDimension ;
3199 Inverse /= (Standard_Real) ii ;
3201 PLib::EvalPolynomial(Delta,
3210 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3211 Index1 = FirstNonZeroBsplineIndex ;
3212 LocalRealArray[Index] = 0.0e0 ;
3214 for (jj = 1 ; jj <= Order ; jj++) {
3215 LocalRealArray[Index] +=
3216 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3217 Index1 = Index1 % Modulus ;
3220 LocalRealArray[Index + kk] *= Inverse ;
3222 Inverse /= (Standard_Real) ii ;
3224 PLib::EvalPolynomial(Delta,
3234 //=======================================================================
3235 //function : Evaluates a Bspline function : uses the ExtrapMode
3236 //purpose : the function is extrapolated using the Taylor expansion
3237 // of degree ExtrapMode[0] to the left and the Taylor
3238 // expansion of degree ExtrapMode[1] to the right
3239 // WARNING : the array Results is supposed to have at least
3240 // (DerivativeRequest + 1) * ArrayDimension slots and the
3242 //=======================================================================
3245 (const Standard_Real Parameter,
3246 const Standard_Boolean PeriodicFlag,
3247 const Standard_Integer DerivativeRequest,
3248 Standard_Integer& ExtrapMode,
3249 const Standard_Integer Degree,
3250 const TColStd_Array1OfReal& FlatKnots,
3251 const Standard_Integer ArrayDimension,
3252 Standard_Real& Poles,
3253 Standard_Real& Results)
3255 Standard_Integer ii,
3263 ExtrapolatingFlag[2],
3267 FirstNonZeroBsplineIndex,
3268 LocalRequest = DerivativeRequest ;
3270 Standard_Real *ResultArray,
3277 PolesArray = &Poles ;
3278 ExtrapModeArray = &ExtrapMode ;
3279 ResultArray = &Results ;
3280 LocalParameter = Parameter ;
3281 ExtrapolatingFlag[0] =
3282 ExtrapolatingFlag[1] = 0 ;
3284 // check if we are extrapolating to a degree which is smaller than
3285 // the degree of the Bspline
3288 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3290 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3291 LocalParameter -= Period ;
3294 while (LocalParameter < FlatKnots(2)) {
3295 LocalParameter += Period ;
3298 if (Parameter < FlatKnots(2) &&
3299 LocalRequest < ExtrapModeArray[0] &&
3300 ExtrapModeArray[0] < Degree) {
3301 LocalRequest = ExtrapModeArray[0] ;
3302 LocalParameter = FlatKnots(2) ;
3303 ExtrapolatingFlag[0] = 1 ;
3305 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3306 LocalRequest < ExtrapModeArray[1] &&
3307 ExtrapModeArray[1] < Degree) {
3308 LocalRequest = ExtrapModeArray[1] ;
3309 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3310 ExtrapolatingFlag[1] = 1 ;
3312 Delta = Parameter - LocalParameter ;
3313 if (LocalRequest >= Order) {
3314 LocalRequest = Degree ;
3318 Modulus = FlatKnots.Length() - Degree -1 ;
3321 Modulus = FlatKnots.Length() - Degree ;
3324 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3327 BSplCLib::EvalBsplineBasis(1,
3332 FirstNonZeroBsplineIndex,
3334 if (ErrorCode != 0) {
3338 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3341 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3342 Index1 = FirstNonZeroBsplineIndex ;
3344 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3345 ResultArray[Index + kk] = 0.0e0 ;
3348 for (jj = 1 ; jj <= Order ; jj++) {
3350 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3351 ResultArray[Index + kk] +=
3352 PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3354 Index1 = Index1 % Modulus ;
3357 Index += ArrayDimension ;
3362 // store Taylor expansion in LocalRealArray
3364 NewRequest = DerivativeRequest ;
3365 if (NewRequest > Degree) {
3366 NewRequest = Degree ;
3368 BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension);
3373 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3374 Index1 = FirstNonZeroBsplineIndex ;
3376 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3377 LocalRealArray[Index + kk] = 0.0e0 ;
3380 for (jj = 1 ; jj <= Order ; jj++) {
3382 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3383 LocalRealArray[Index + kk] +=
3384 PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3386 Index1 = Index1 % Modulus ;
3390 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3391 LocalRealArray[Index + kk] *= Inverse ;
3393 Index += ArrayDimension ;
3394 Inverse /= (Standard_Real) ii ;
3396 PLib::EvalPolynomial(Delta,
3406 //=======================================================================
3407 //function : TangExtendToConstraint
3408 //purpose : Extends a Bspline function using the tangency map
3412 //=======================================================================
3414 void BSplCLib::TangExtendToConstraint
3415 (const TColStd_Array1OfReal& FlatKnots,
3416 const Standard_Real C1Coefficient,
3417 const Standard_Integer NumPoles,
3418 Standard_Real& Poles,
3419 const Standard_Integer CDimension,
3420 const Standard_Integer CDegree,
3421 const TColStd_Array1OfReal& ConstraintPoint,
3422 const Standard_Integer Continuity,
3423 const Standard_Boolean After,
3424 Standard_Integer& NbPolesResult,
3425 Standard_Integer& NbKnotsResult,
3426 Standard_Real& KnotsResult,
3427 Standard_Real& PolesResult)
3430 if (CDegree<Continuity+1) {
3431 cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3434 Standard_Real * Padr = &Poles ;
3435 Standard_Real * KRadr = &KnotsResult ;
3436 Standard_Real * PRadr = &PolesResult ;
3438 ////////////////////////////////////////////////////////////////////////
3440 // 1. calculation of extension nD
3442 ////////////////////////////////////////////////////////////////////////
3445 Standard_Integer Csize = Continuity + 2;
3446 math_Matrix MatCoefs(1,Csize, 1,Csize);
3448 PLib::HermiteCoefficients(0, 1, // Limits
3449 Continuity, 0, // Orders of constraints
3453 PLib::HermiteCoefficients(0, 1, // Limits
3454 0, Continuity, // Orders of constraints
3459 // position at the node of connection
3460 Standard_Real Tbord ;
3462 Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3465 Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3467 Standard_Boolean periodic_flag = Standard_False ;
3468 Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3469 extrap_mode[0] = extrap_mode[1] = CDegree;
3470 TColStd_Array1OfReal EvalBS(1, CDimension * (derivative_request+1)) ;
3471 Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3472 BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3473 CDegree,FlatKnots,CDimension,Poles,*Eadr);
3475 // norm of the tangent at the node of connection
3476 math_Vector Tgte(1,CDimension);
3478 for (ipos=1;ipos<=CDimension;ipos++) {
3479 Tgte(ipos) = EvalBS(ipos+CDimension);
3481 Standard_Real L1=Tgte.Norm();
3484 // matrix of constraints
3485 math_Matrix Contraintes(1,Csize,1,CDimension);
3488 for (ipos=1;ipos<=CDimension;ipos++) {
3489 Contraintes(1,ipos) = EvalBS(ipos);
3490 Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3491 if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3492 if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3493 Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3498 for (ipos=1;ipos<=CDimension;ipos++) {
3499 Contraintes(1,ipos) = ConstraintPoint(ipos);
3500 Contraintes(2,ipos) = EvalBS(ipos);
3501 if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3502 if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3503 if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3507 // calculate the coefficients of extension
3508 Standard_Integer ii, jj, kk;
3509 TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3510 ExtraCoeffs.Init(0.);
3512 for (ii=1; ii<=Csize; ii++) {
3514 for (jj=1; jj<=Csize; jj++) {
3516 for (kk=1; kk<=CDimension; kk++) {
3517 ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3522 // calculate the poles of extension
3523 TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3524 Standard_Real * EPadr = &ExtrapPoles(1) ;
3525 PLib::CoefficientsPoles(CDimension,
3526 ExtraCoeffs, PLib::NoWeights(),
3527 ExtrapPoles, PLib::NoWeights());
3529 // calculate the nodes of extension with multiplicities
3530 TColStd_Array1OfReal ExtrapNoeuds(1,2);
3531 ExtrapNoeuds(1) = 0.;
3532 ExtrapNoeuds(2) = 1.;
3533 TColStd_Array1OfInteger ExtrapMults(1,2);
3534 ExtrapMults(1) = Csize;
3535 ExtrapMults(2) = Csize;
3537 // flat nodes of extension
3538 TColStd_Array1OfReal FK2(1, Csize*2);
3539 BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3541 // norm of the tangent at the connection point
3543 BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3544 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3547 BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3548 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3551 for (ipos=1;ipos<=CDimension;ipos++) {
3552 Tgte(ipos) = EvalBS(ipos+CDimension);
3554 Standard_Real L2 = Tgte.Norm();
3556 // harmonisation of degrees
3557 TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3558 TColStd_Array1OfReal NewK2(1, 2);
3559 TColStd_Array1OfInteger NewM2(1, 2);
3560 if (Csize-1<CDegree) {
3561 BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3562 ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3566 NewP2 = ExtrapPoles;
3567 NewK2 = ExtrapNoeuds;
3568 NewM2 = ExtrapMults;
3571 // flat nodes of extension after harmonization of degrees
3572 TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3573 BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3576 ////////////////////////////////////////////////////////////////////////
3578 // 2. concatenation C0
3580 ////////////////////////////////////////////////////////////////////////
3582 // ratio of reparametrization
3583 Standard_Real Ratio=1, Delta;
3584 if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3587 if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3590 // do not touch the first BSpline
3591 Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3594 // do not touch the second BSpline
3595 Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower());
3598 // result of the concatenation
3599 Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1;
3600 Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1);
3601 TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension);
3602 TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2);
3605 Standard_Integer indNP, indP, indEP;
3608 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3610 for (jj=1; jj<=CDimension; jj++) {
3611 indNP = (ii-1)*CDimension+jj;
3612 indP = (ii-1)*CDimension+jj-1;
3613 indEP = (ii-NbP1)*CDimension+jj;
3614 if (ii<NbP1) NewPoles(indNP) = Padr[indP];
3615 else NewPoles(indNP) = NewP2(indEP);
3621 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3623 for (jj=1; jj<=CDimension; jj++) {
3624 indNP = (ii-1)*CDimension+jj;
3625 indEP = (ii-1)*CDimension+jj;
3626 indP = (ii-NbP2)*CDimension+jj-1;
3627 if (ii<NbP2) NewPoles(indNP) = NewP2(indEP);
3628 else NewPoles(indNP) = Padr[indP];
3635 // start with the nodes of the initial surface
3637 for (ii=1; ii<NbK1; ii++) {
3638 NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1);
3640 // continue with the reparameterized nodes of the extension
3642 for (ii=1; ii<=NbK2-CDegree-1; ii++) {
3643 NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta;
3647 // start with the reparameterized nodes of the extension
3649 for (ii=1; ii<NbK2-CDegree; ii++) {
3650 NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta;
3652 // continue with the nodes of the initial surface
3654 for (ii=2; ii<=NbK1; ii++) {
3655 NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1);
3660 ////////////////////////////////////////////////////////////////////////
3662 // 3. reduction of multiplicite at the node of connection
3664 ////////////////////////////////////////////////////////////////////////
3666 // number of separate nodes
3667 Standard_Integer KLength = 1;
3669 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3670 if (NewFlats(ii) != NewFlats(ii-1)) KLength++;
3673 // flat nodes --> nodes + multiplicities
3674 TColStd_Array1OfReal NewKnots (1, KLength);
3675 TColStd_Array1OfInteger NewMults (1, KLength);
3678 NewKnots(jj) = NewFlats(1);
3680 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3681 if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++;
3684 NewKnots(jj) = NewFlats(ii);
3688 // reduction of multiplicity at the second or the last but one node
3689 Standard_Integer Index = 2, M = CDegree;
3690 if (After) Index = KLength-1;
3691 TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension);
3692 TColStd_Array1OfReal ResultKnots (1, KLength);
3693 TColStd_Array1OfInteger ResultMults (1, KLength);
3694 Standard_Real Tol = 1.e-6;
3695 Standard_Boolean Ok = Standard_True;
3697 while ( (M>CDegree-Continuity) && Ok) {
3698 Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension,
3699 NewPoles, NewKnots, NewMults,
3700 ResultPoles, ResultKnots, ResultMults, Tol);
3705 // number of poles of the concatenation
3706 NbPolesResult = NbP1 + NbP2 - 1;
3707 // the poles of the concatenation
3708 Standard_Integer PLength = NbPolesResult*CDimension;
3710 for (jj=1; jj<=PLength; jj++) {
3711 PRadr[jj-1] = NewPoles(jj);
3714 // flat nodes of the concatenation
3715 Standard_Integer ideb = 0;
3717 for (jj=0; jj<NewKnots.Length(); jj++) {
3718 for (ii=0; ii<NewMults(jj+1); ii++) {
3719 KRadr[ideb+ii] = NewKnots(jj+1);
3721 ideb += NewMults(jj+1);
3723 NbKnotsResult = ideb;
3727 // number of poles of the result
3728 NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M;
3729 // the poles of the result
3730 Standard_Integer PLength = NbPolesResult*CDimension;
3732 for (jj=0; jj<PLength; jj++) {
3733 PRadr[jj] = ResultPoles(jj+1);
3736 // flat nodes of the result
3737 Standard_Integer ideb = 0;
3739 for (jj=0; jj<ResultKnots.Length(); jj++) {
3740 for (ii=0; ii<ResultMults(jj+1); ii++) {
3741 KRadr[ideb+ii] = ResultKnots(jj+1);
3743 ideb += ResultMults(jj+1);
3745 NbKnotsResult = ideb;
3749 //=======================================================================
3750 //function : Resolution
3753 // Let C(t) = SUM Ci Bi(t) a Bspline curve of degree d
3755 // with nodes tj for j = 1,n+d+1
3759 // Then C (t) = SUM d * --------- Bi (t)
3760 // i = 2,n ti+d - ti
3763 // for the base of BSpline Bi (t) of degree d-1.
3765 // Consequently the upper bound of the norm of the derivative from C is :
3769 // d * Max | --------- |
3770 // i = 2,n | ti+d - ti |
3773 // In the rational case set C(t) = -----
3777 // D(t) = SUM Di Bi(t)
3780 // N(t) = SUM Di * Ci Bi(t)
3783 // N'(t) - D'(t) C(t)
3784 // C'(t) = -----------------------
3788 // N'(t) - D'(t) C(t) =
3790 // Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t)) d-1
3791 // SUM d * ---------------------------------------- * Bi (t) =
3795 // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) d-1
3796 // SUM SUM d * ----------------------------------- * Betaj(t) * Bi (t)
3797 //i=2,n j=1,n ti+d - ti
3802 // Betaj(t) = --------
3805 // Betaj(t) form a partition >= 0 of the entity with support
3806 // tj, tj+d+1. Consequently if Rj = {j-d, ...., j+d+d+1}
3807 // obtain an upper bound of the derivative of C by taking :
3814 // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)
3815 // Max Max d * -----------------------------------
3816 // j=1,n i dans Rj ti+d - ti
3818 // --------------------------------------------------------
3824 //=======================================================================
3826 void BSplCLib::Resolution( Standard_Real& Poles,
3827 const Standard_Integer ArrayDimension,
3828 const Standard_Integer NumPoles,
3829 const TColStd_Array1OfReal& Weights,
3830 const TColStd_Array1OfReal& FlatKnots,
3831 const Standard_Integer Degree,
3832 const Standard_Real Tolerance3D,
3833 Standard_Real& UTolerance)
3835 Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim;
3836 Standard_Integer lower,upper,ii_minus,jj,ii_miDim;
3837 Standard_Integer Deg1 = Degree + 1;
3838 Standard_Integer Deg2 = (Degree << 1) + 1;
3839 Standard_Real value,factor,W,min_weights,inverse;
3840 Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3;
3841 Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3;
3842 Standard_Real wg_ii_index, wg_ii_minus;
3843 Standard_Real *PA,max_derivative;
3844 const Standard_Real * FK = &FlatKnots(FlatKnots.Lower());
3846 max_derivative = 0.0e0;
3847 num_poles = FlatKnots.Length() - Deg1;
3848 switch (ArrayDimension) {
3850 if (&Weights != NULL) {
3851 const Standard_Real * WG = &Weights(Weights.Lower());
3852 min_weights = WG[0];
3854 for (ii = 1 ; ii < NumPoles ; ii++) {
3856 if (W < min_weights) min_weights = W;
3859 for (ii = 1 ; ii < num_poles ; ii++) {
3860 ii_index = ii % NumPoles;
3861 ii_inDim = ii_index << 1;
3862 ii_minus = (ii - 1) % NumPoles;
3863 ii_miDim = ii_minus << 1;
3864 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3865 pa_ii_inDim_1 = PA[ii_inDim];
3866 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3867 pa_ii_miDim_1 = PA[ii_miDim];
3868 wg_ii_index = WG[ii_index];
3869 wg_ii_minus = WG[ii_minus];
3870 inverse = FK[ii + Degree] - FK[ii];
3871 inverse = 1.0e0 / inverse;
3873 if (lower < 0) lower = 0;
3875 if (upper > num_poles) upper = num_poles;
3877 for (jj = lower ; jj < upper ; jj++) {
3878 jj_index = jj % NumPoles;
3879 jj_index = jj_index << 1;
3881 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3882 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3883 if (factor < 0) factor = - factor;
3884 value += factor; jj_index++;
3885 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3886 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3887 if (factor < 0) factor = - factor;
3890 if (max_derivative < value) max_derivative = value;
3893 max_derivative /= min_weights;
3897 for (ii = 1 ; ii < num_poles ; ii++) {
3898 ii_index = ii % NumPoles;
3899 ii_index = ii_index << 1;
3900 ii_minus = (ii - 1) % NumPoles;
3901 ii_minus = ii_minus << 1;
3902 inverse = FK[ii + Degree] - FK[ii];
3903 inverse = 1.0e0 / inverse;
3905 factor = PA[ii_index] - PA[ii_minus];
3906 if (factor < 0) factor = - factor;
3907 value += factor; ii_index++; ii_minus++;
3908 factor = PA[ii_index] - PA[ii_minus];
3909 if (factor < 0) factor = - factor;
3912 if (max_derivative < value) max_derivative = value;
3918 if (&Weights != NULL) {
3919 const Standard_Real * WG = &Weights(Weights.Lower());
3920 min_weights = WG[0];
3922 for (ii = 1 ; ii < NumPoles ; ii++) {
3924 if (W < min_weights) min_weights = W;
3927 for (ii = 1 ; ii < num_poles ; ii++) {
3928 ii_index = ii % NumPoles;
3929 ii_inDim = (ii_index << 1) + ii_index;
3930 ii_minus = (ii - 1) % NumPoles;
3931 ii_miDim = (ii_minus << 1) + ii_minus;
3932 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3933 pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
3934 pa_ii_inDim_2 = PA[ii_inDim];
3935 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3936 pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
3937 pa_ii_miDim_2 = PA[ii_miDim];
3938 wg_ii_index = WG[ii_index];
3939 wg_ii_minus = WG[ii_minus];
3940 inverse = FK[ii + Degree] - FK[ii];
3941 inverse = 1.0e0 / inverse;
3943 if (lower < 0) lower = 0;
3945 if (upper > num_poles) upper = num_poles;
3947 for (jj = lower ; jj < upper ; jj++) {
3948 jj_index = jj % NumPoles;
3949 jj_index = (jj_index << 1) + jj_index;
3951 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3952 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3953 if (factor < 0) factor = - factor;
3954 value += factor; jj_index++;
3955 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3956 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3957 if (factor < 0) factor = - factor;
3958 value += factor; jj_index++;
3959 factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
3960 ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
3961 if (factor < 0) factor = - factor;
3964 if (max_derivative < value) max_derivative = value;
3967 max_derivative /= min_weights;
3971 for (ii = 1 ; ii < num_poles ; ii++) {
3972 ii_index = ii % NumPoles;
3973 ii_index = (ii_index << 1) + ii_index;
3974 ii_minus = (ii - 1) % NumPoles;
3975 ii_minus = (ii_minus << 1) + ii_minus;
3976 inverse = FK[ii + Degree] - FK[ii];
3977 inverse = 1.0e0 / inverse;
3979 factor = PA[ii_index] - PA[ii_minus];
3980 if (factor < 0) factor = - factor;
3981 value += factor; ii_index++; ii_minus++;
3982 factor = PA[ii_index] - PA[ii_minus];
3983 if (factor < 0) factor = - factor;
3984 value += factor; ii_index++; ii_minus++;
3985 factor = PA[ii_index] - PA[ii_minus];
3986 if (factor < 0) factor = - factor;
3989 if (max_derivative < value) max_derivative = value;
3995 if (&Weights != NULL) {
3996 const Standard_Real * WG = &Weights(Weights.Lower());
3997 min_weights = WG[0];
3999 for (ii = 1 ; ii < NumPoles ; ii++) {
4001 if (W < min_weights) min_weights = W;
4004 for (ii = 1 ; ii < num_poles ; ii++) {
4005 ii_index = ii % NumPoles;
4006 ii_inDim = ii_index << 2;
4007 ii_minus = (ii - 1) % NumPoles;
4008 ii_miDim = ii_minus << 2;
4009 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
4010 pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
4011 pa_ii_inDim_2 = PA[ii_inDim]; ii_inDim++;
4012 pa_ii_inDim_3 = PA[ii_inDim];
4013 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
4014 pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
4015 pa_ii_miDim_2 = PA[ii_miDim]; ii_miDim++;
4016 pa_ii_miDim_3 = PA[ii_miDim];
4017 wg_ii_index = WG[ii_index];
4018 wg_ii_minus = WG[ii_minus];
4019 inverse = FK[ii + Degree] - FK[ii];
4020 inverse = 1.0e0 / inverse;
4022 if (lower < 0) lower = 0;
4024 if (upper > num_poles) upper = num_poles;
4026 for (jj = lower ; jj < upper ; jj++) {
4027 jj_index = jj % NumPoles;
4028 jj_index = jj_index << 2;
4030 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4031 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4032 if (factor < 0) factor = - factor;
4033 value += factor; jj_index++;
4034 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4035 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4036 if (factor < 0) factor = - factor;
4037 value += factor; jj_index++;
4038 factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4039 ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4040 if (factor < 0) factor = - factor;
4041 value += factor; jj_index++;
4042 factor = (((PA[jj_index] - pa_ii_inDim_3) * wg_ii_index) -
4043 ((PA[jj_index] - pa_ii_miDim_3) * wg_ii_minus));
4044 if (factor < 0) factor = - factor;
4047 if (max_derivative < value) max_derivative = value;
4050 max_derivative /= min_weights;
4054 for (ii = 1 ; ii < num_poles ; ii++) {
4055 ii_index = ii % NumPoles;
4056 ii_index = ii_index << 2;
4057 ii_minus = (ii - 1) % NumPoles;
4058 ii_minus = ii_minus << 2;
4059 inverse = FK[ii + Degree] - FK[ii];
4060 inverse = 1.0e0 / inverse;
4062 factor = PA[ii_index] - PA[ii_minus];
4063 if (factor < 0) factor = - factor;
4064 value += factor; ii_index++; ii_minus++;
4065 factor = PA[ii_index] - PA[ii_minus];
4066 if (factor < 0) factor = - factor;
4067 value += factor; ii_index++; ii_minus++;
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;
4075 if (max_derivative < value) max_derivative = value;
4081 Standard_Integer kk;
4082 if (&Weights != NULL) {
4083 const Standard_Real * WG = &Weights(Weights.Lower());
4084 min_weights = WG[0];
4086 for (ii = 1 ; ii < NumPoles ; ii++) {
4088 if (W < min_weights) min_weights = W;
4091 for (ii = 1 ; ii < num_poles ; ii++) {
4092 ii_index = ii % NumPoles;
4093 ii_inDim = ii_index * ArrayDimension;
4094 ii_minus = (ii - 1) % NumPoles;
4095 ii_miDim = ii_minus * ArrayDimension;
4096 wg_ii_index = WG[ii_index];
4097 wg_ii_minus = WG[ii_minus];
4098 inverse = FK[ii + Degree] - FK[ii];
4099 inverse = 1.0e0 / inverse;
4101 if (lower < 0) lower = 0;
4103 if (upper > num_poles) upper = num_poles;
4105 for (jj = lower ; jj < upper ; jj++) {
4106 jj_index = jj % NumPoles;
4107 jj_index *= ArrayDimension;
4110 for (kk = 0 ; kk < ArrayDimension ; kk++) {
4111 factor = (((PA[jj_index + kk] - PA[ii_inDim + kk]) * wg_ii_index) -
4112 ((PA[jj_index + kk] - PA[ii_miDim + kk]) * wg_ii_minus));
4113 if (factor < 0) factor = - factor;
4117 if (max_derivative < value) max_derivative = value;
4120 max_derivative /= min_weights;
4124 for (ii = 1 ; ii < num_poles ; ii++) {
4125 ii_index = ii % NumPoles;
4126 ii_index *= ArrayDimension;
4127 ii_minus = (ii - 1) % NumPoles;
4128 ii_minus *= ArrayDimension;
4129 inverse = FK[ii + Degree] - FK[ii];
4130 inverse = 1.0e0 / inverse;
4133 for (kk = 0 ; kk < ArrayDimension ; kk++) {
4134 factor = PA[ii_index + kk] - PA[ii_minus + kk];
4135 if (factor < 0) factor = - factor;
4139 if (max_derivative < value) max_derivative = value;
4144 max_derivative *= Degree;
4145 if (max_derivative > RealSmall())
4146 UTolerance = Tolerance3D / max_derivative;
4148 UTolerance = Tolerance3D / RealSmall();
4151 //=======================================================================
4152 // function: FlatBezierKnots
4154 //=======================================================================
4156 // array of flat knots for bezier curve of maximum 25 degree
4157 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,
4158 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 };
4159 const Standard_Real& BSplCLib::FlatBezierKnots (const Standard_Integer Degree)
4161 Standard_OutOfRange_Raise_if (Degree < 1 || Degree > MaxDegree() || MaxDegree() != 25,
4162 "Bezier curve degree greater than maximal supported");
4164 return knots[25-Degree];