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 <NCollection_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 //=======================================================================
73 //=======================================================================
75 void BSplCLib::Hunt (const Array1OfReal& XX,
76 const Standard_Real X,
77 Standard_Integer& Ilc)
79 // replaced by simple dichotomy (RLE)
81 const Standard_Real *px = &XX(Ilc);
88 Standard_Integer Ihi = XX.Upper();
95 while (Ihi - Ilc != 1) {
96 Im = (Ihi + Ilc) >> 1;
97 if (X > px[Im]) Ilc = Im;
102 //=======================================================================
103 //function : FirstUKnotIndex
105 //=======================================================================
107 Standard_Integer BSplCLib::FirstUKnotIndex (const Standard_Integer Degree,
108 const TColStd_Array1OfInteger& Mults)
110 Standard_Integer Index = Mults.Lower();
111 Standard_Integer SigmaMult = Mults(Index);
113 while (SigmaMult <= Degree) {
115 SigmaMult += Mults (Index);
120 //=======================================================================
121 //function : LastUKnotIndex
123 //=======================================================================
125 Standard_Integer BSplCLib::LastUKnotIndex (const Standard_Integer Degree,
126 const Array1OfInteger& Mults)
128 Standard_Integer Index = Mults.Upper();
129 Standard_Integer SigmaMult = Mults(Index);
131 while (SigmaMult <= Degree) {
133 SigmaMult += Mults.Value (Index);
138 //=======================================================================
139 //function : FlatIndex
141 //=======================================================================
143 Standard_Integer BSplCLib::FlatIndex
144 (const Standard_Integer Degree,
145 const Standard_Integer Index,
146 const TColStd_Array1OfInteger& Mults,
147 const Standard_Boolean Periodic)
149 Standard_Integer i, index = Index;
150 const Standard_Integer MLower = Mults.Lower();
151 const Standard_Integer *pmu = &Mults(MLower);
154 for (i = MLower + 1; i <= Index; i++)
159 index += pmu[MLower] - 1;
163 //=======================================================================
164 //function : LocateParameter
165 //purpose : Processing of nodes with multiplicities
166 //pmn 28-01-97 -> compute eventual of the period.
167 //=======================================================================
169 void BSplCLib::LocateParameter
170 (const Standard_Integer , //Degree,
171 const Array1OfReal& Knots,
172 const Array1OfInteger& , //Mults,
173 const Standard_Real U,
174 const Standard_Boolean IsPeriodic,
175 const Standard_Integer FromK1,
176 const Standard_Integer ToK2,
177 Standard_Integer& KnotIndex,
180 Standard_Real uf = 0, ul=1;
182 uf = Knots(Knots.Lower());
183 ul = Knots(Knots.Upper());
185 BSplCLib::LocateParameter(Knots,U,IsPeriodic,FromK1,ToK2,
186 KnotIndex,NewU, uf, ul);
189 //=======================================================================
190 //function : LocateParameter
191 //purpose : For plane nodes
192 // pmn 28-01-97 -> There is a need of the degre to calculate
193 // the eventual period
194 //=======================================================================
196 void BSplCLib::LocateParameter
197 (const Standard_Integer Degree,
198 const Array1OfReal& Knots,
199 const Standard_Real U,
200 const Standard_Boolean IsPeriodic,
201 const Standard_Integer FromK1,
202 const Standard_Integer ToK2,
203 Standard_Integer& KnotIndex,
207 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
209 Knots(Knots.Lower() + Degree),
210 Knots(Knots.Upper() - Degree));
212 BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2,
218 //=======================================================================
219 //function : LocateParameter
220 //purpose : Effective computation
221 // pmn 28-01-97 : Add limits of the period as input argument,
222 // as it is imposible to produce them at this level.
223 //=======================================================================
225 void BSplCLib::LocateParameter
226 (const TColStd_Array1OfReal& Knots,
227 const Standard_Real U,
228 const Standard_Boolean IsPeriodic,
229 const Standard_Integer FromK1,
230 const Standard_Integer ToK2,
231 Standard_Integer& KnotIndex,
233 const Standard_Real UFirst,
234 const Standard_Real ULast)
236 Standard_Integer First,Last;
245 Standard_Integer Last1 = Last - 1;
248 Standard_Real Period = ULast - UFirst;
250 while (NewU > ULast )
253 while (NewU < UFirst)
257 BSplCLib::Hunt (Knots, NewU, KnotIndex);
259 Standard_Real Eps = Epsilon(U);
261 if (Eps < 0) Eps = - Eps;
262 Standard_Integer KLower = Knots.Lower();
263 const Standard_Real *knots = &Knots(KLower);
265 if ( KnotIndex < Knots.Upper()) {
266 val = NewU - knots[KnotIndex + 1];
267 if (val < 0) val = - val;
268 // <= to be coherent with Segment where Eps corresponds to a bit of error.
269 if (val <= Eps) KnotIndex++;
271 if (KnotIndex < First) KnotIndex = First;
272 if (KnotIndex > Last1) KnotIndex = Last1;
274 if (KnotIndex != Last1) {
275 Standard_Real K1 = knots[KnotIndex];
276 Standard_Real K2 = knots[KnotIndex + 1];
278 if (val < 0) val = - val;
283 K2 = knots[KnotIndex + 1];
285 if (val < 0) val = - val;
290 //=======================================================================
291 //function : LocateParameter
292 //purpose : the index is recomputed only if out of range
293 //pmn 28-01-97 -> eventual computation of the period.
294 //=======================================================================
296 void BSplCLib::LocateParameter
297 (const Standard_Integer Degree,
298 const TColStd_Array1OfReal& Knots,
299 const TColStd_Array1OfInteger& Mults,
300 const Standard_Real U,
301 const Standard_Boolean Periodic,
302 Standard_Integer& KnotIndex,
305 Standard_Integer first,last;
308 first = Knots.Lower();
309 last = Knots.Upper();
312 first = FirstUKnotIndex(Degree,Mults);
313 last = LastUKnotIndex (Degree,Mults);
317 first = Knots.Lower() + Degree;
318 last = Knots.Upper() - Degree;
320 if ( KnotIndex < first || KnotIndex > last)
321 BSplCLib::LocateParameter(Knots, U, Periodic, first, last,
322 KnotIndex, NewU, Knots(first), Knots(last));
327 //=======================================================================
328 //function : MaxKnotMult
330 //=======================================================================
332 Standard_Integer BSplCLib::MaxKnotMult
333 (const Array1OfInteger& Mults,
334 const Standard_Integer FromK1,
335 const Standard_Integer ToK2)
337 Standard_Integer MLower = Mults.Lower();
338 const Standard_Integer *pmu = &Mults(MLower);
340 Standard_Integer MaxMult = pmu[FromK1];
342 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
343 if (MaxMult < pmu[i]) MaxMult = pmu[i];
348 //=======================================================================
349 //function : MinKnotMult
351 //=======================================================================
353 Standard_Integer BSplCLib::MinKnotMult
354 (const Array1OfInteger& Mults,
355 const Standard_Integer FromK1,
356 const Standard_Integer ToK2)
358 Standard_Integer MLower = Mults.Lower();
359 const Standard_Integer *pmu = &Mults(MLower);
361 Standard_Integer MinMult = pmu[FromK1];
363 for (Standard_Integer i = FromK1; i <= ToK2; i++) {
364 if (MinMult > pmu[i]) MinMult = pmu[i];
369 //=======================================================================
372 //=======================================================================
374 Standard_Integer BSplCLib::NbPoles(const Standard_Integer Degree,
375 const Standard_Boolean Periodic,
376 const TColStd_Array1OfInteger& Mults)
378 Standard_Integer i,sigma = 0;
379 Standard_Integer f = Mults.Lower();
380 Standard_Integer l = Mults.Upper();
381 const Standard_Integer * pmu = &Mults(f);
383 Standard_Integer Mf = pmu[f];
384 Standard_Integer Ml = pmu[l];
385 if (Mf <= 0) return 0;
386 if (Ml <= 0) return 0;
388 if (Mf > Degree) return 0;
389 if (Ml > Degree) return 0;
390 if (Mf != Ml ) return 0;
394 Standard_Integer Deg1 = Degree + 1;
395 if (Mf > Deg1) return 0;
396 if (Ml > Deg1) return 0;
397 sigma = Mf + Ml - Deg1;
400 for (i = f + 1; i < l; i++) {
401 if (pmu[i] <= 0 ) return 0;
402 if (pmu[i] > Degree) return 0;
408 //=======================================================================
409 //function : KnotSequenceLength
411 //=======================================================================
413 Standard_Integer BSplCLib::KnotSequenceLength
414 (const TColStd_Array1OfInteger& Mults,
415 const Standard_Integer Degree,
416 const Standard_Boolean Periodic)
418 Standard_Integer i,l = 0;
419 Standard_Integer MLower = Mults.Lower();
420 Standard_Integer MUpper = Mults.Upper();
421 const Standard_Integer * pmu = &Mults(MLower);
424 for (i = MLower; i <= MUpper; i++)
426 if (Periodic) l += 2 * (Degree + 1 - pmu[MLower]);
430 //=======================================================================
431 //function : KnotSequence
433 //=======================================================================
435 void BSplCLib::KnotSequence
436 (const TColStd_Array1OfReal& Knots,
437 const TColStd_Array1OfInteger& Mults,
438 TColStd_Array1OfReal& KnotSeq)
440 BSplCLib::KnotSequence(Knots,Mults,0,Standard_False,KnotSeq);
443 //=======================================================================
444 //function : KnotSequence
446 //=======================================================================
448 void BSplCLib::KnotSequence
449 (const TColStd_Array1OfReal& Knots,
450 const TColStd_Array1OfInteger& Mults,
451 const Standard_Integer Degree,
452 const Standard_Boolean Periodic,
453 TColStd_Array1OfReal& KnotSeq)
456 Standard_Integer Mult;
457 Standard_Integer MLower = Mults.Lower();
458 const Standard_Integer * pmu = &Mults(MLower);
460 Standard_Integer KLower = Knots.Lower();
461 Standard_Integer KUpper = Knots.Upper();
462 const Standard_Real * pkn = &Knots(KLower);
464 Standard_Integer M1 = Degree + 1 - pmu[MLower]; // for periodic
465 Standard_Integer i,j,index = Periodic ? M1 + 1 : 1;
467 for (i = KLower; i <= KUpper; i++) {
471 for (j = 1; j <= Mult; j++) {
477 Standard_Real period = pkn[KUpper] - pkn[KLower];
482 for (i = M1; i >= 1; i--) {
483 KnotSeq(i) = pkn[j] - period;
493 for (i = index; i <= KnotSeq.Upper(); i++) {
494 KnotSeq(i) = pkn[j] + period;
504 //=======================================================================
505 //function : KnotsLength
507 //=======================================================================
508 Standard_Integer BSplCLib::KnotsLength(const TColStd_Array1OfReal& SeqKnots,
509 // const Standard_Boolean Periodic)
510 const Standard_Boolean )
512 Standard_Integer sizeMult = 1;
513 Standard_Real val = SeqKnots(1);
514 for (Standard_Integer jj=2;
515 jj<=SeqKnots.Length();jj++)
517 // test on strict equality on nodes
518 if (SeqKnots(jj)!=val)
527 //=======================================================================
530 //=======================================================================
531 void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots,
532 TColStd_Array1OfReal &knots,
533 TColStd_Array1OfInteger &mult,
534 // const Standard_Boolean Periodic)
535 const Standard_Boolean )
537 Standard_Real val = SeqKnots(1);
538 Standard_Integer kk=1;
542 for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++)
544 // test on strict equality on nodes
545 if (SeqKnots(jj)!=val)
559 //=======================================================================
560 //function : KnotForm
562 //=======================================================================
564 BSplCLib_KnotDistribution BSplCLib::KnotForm
565 (const Array1OfReal& Knots,
566 const Standard_Integer FromK1,
567 const Standard_Integer ToK2)
569 Standard_Real DU0,DU1,Ui,Uj,Eps0,val;
570 BSplCLib_KnotDistribution KForm = BSplCLib_Uniform;
572 Standard_Integer KLower = Knots.Lower();
573 const Standard_Real * pkn = &Knots(KLower);
576 if (Ui < 0) Ui = - Ui;
577 Uj = pkn[FromK1 + 1];
578 if (Uj < 0) Uj = - Uj;
580 if (DU0 < 0) DU0 = - DU0;
581 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
582 Standard_Integer i = FromK1 + 1;
584 while (KForm != BSplCLib_NonUniform && i < ToK2) {
586 if (Ui < 0) Ui = - Ui;
589 if (Uj < 0) Uj = - Uj;
591 if (DU1 < 0) DU1 = - DU1;
593 if (val < 0) val = -val;
594 if (val > Eps0) KForm = BSplCLib_NonUniform;
596 Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0);
601 //=======================================================================
602 //function : MultForm
604 //=======================================================================
606 BSplCLib_MultDistribution BSplCLib::MultForm
607 (const Array1OfInteger& Mults,
608 const Standard_Integer FromK1,
609 const Standard_Integer ToK2)
611 Standard_Integer First,Last;
620 Standard_Integer MLower = Mults.Lower();
621 const Standard_Integer *pmu = &Mults(MLower);
623 Standard_Integer FirstMult = pmu[First];
624 BSplCLib_MultDistribution MForm = BSplCLib_Constant;
625 Standard_Integer i = First + 1;
626 Standard_Integer Mult = pmu[i];
628 // while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR????????
629 while (MForm != BSplCLib_NonConstant && i <= Last) {
630 if (i == First + 1) {
631 if (Mult != FirstMult) MForm = BSplCLib_QuasiConstant;
633 else if (i == Last) {
634 if (MForm == BSplCLib_QuasiConstant) {
635 if (FirstMult != pmu[i]) MForm = BSplCLib_NonConstant;
638 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
642 if (Mult != pmu[i]) MForm = BSplCLib_NonConstant;
650 //=======================================================================
651 //function : KnotAnalysis
653 //=======================================================================
655 void BSplCLib::KnotAnalysis (const Standard_Integer Degree,
656 const Standard_Boolean Periodic,
657 const TColStd_Array1OfReal& CKnots,
658 const TColStd_Array1OfInteger& CMults,
659 GeomAbs_BSplKnotDistribution& KnotForm,
660 Standard_Integer& MaxKnotMult)
662 KnotForm = GeomAbs_NonUniform;
664 BSplCLib_KnotDistribution KSet =
665 BSplCLib::KnotForm (CKnots, 1, CKnots.Length());
668 if (KSet == BSplCLib_Uniform) {
669 BSplCLib_MultDistribution MSet =
670 BSplCLib::MultForm (CMults, 1, CMults.Length());
672 case BSplCLib_NonConstant :
674 case BSplCLib_Constant :
675 if (CKnots.Length() == 2) {
676 KnotForm = GeomAbs_PiecewiseBezier;
679 if (CMults (1) == 1) KnotForm = GeomAbs_Uniform;
682 case BSplCLib_QuasiConstant :
683 if (CMults (1) == Degree + 1) {
684 Standard_Real M = CMults (2);
685 if (M == Degree ) KnotForm = GeomAbs_PiecewiseBezier;
686 else if (M == 1) KnotForm = GeomAbs_QuasiUniform;
692 Standard_Integer FirstKM =
693 Periodic ? CKnots.Lower() : BSplCLib::FirstUKnotIndex (Degree,CMults);
694 Standard_Integer LastKM =
695 Periodic ? CKnots.Upper() : BSplCLib::LastUKnotIndex (Degree,CMults);
697 if (LastKM - FirstKM != 1) {
698 Standard_Integer Multi;
699 for (Standard_Integer i = FirstKM + 1; i < LastKM; i++) {
701 MaxKnotMult = Max (MaxKnotMult, Multi);
706 //=======================================================================
707 //function : Reparametrize
709 //=======================================================================
711 void BSplCLib::Reparametrize
712 (const Standard_Real U1,
713 const Standard_Real U2,
716 Standard_Integer Lower = Knots.Lower();
717 Standard_Integer Upper = Knots.Upper();
718 Standard_Real UFirst = Min (U1, U2);
719 Standard_Real ULast = Max (U1, U2);
720 Standard_Real NewLength = ULast - UFirst;
721 BSplCLib_KnotDistribution KSet = BSplCLib::KnotForm (Knots, Lower, Upper);
722 if (KSet == BSplCLib_Uniform) {
723 Standard_Real DU = NewLength / (Upper - Lower);
724 Knots (Lower) = UFirst;
726 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
727 Knots (i) = Knots (i-1) + DU;
733 Standard_Real K1 = Knots (Lower);
734 Standard_Real Length = Knots (Upper) - Knots (Lower);
735 Knots (Lower) = UFirst;
737 for (Standard_Integer i = Lower + 1; i <= Upper; i++) {
739 Ratio = (K2 - K1) / Length;
740 Knots (i) = Knots (i-1) + (NewLength * Ratio);
743 Standard_Real Eps = Epsilon( Abs(Knots(i-1)) );
744 if (Knots(i) - Knots(i-1) <= Eps)
745 Knots(i) = NextAfter (Knots(i-1) + Eps, RealLast());
752 //=======================================================================
755 //=======================================================================
757 void BSplCLib::Reverse(TColStd_Array1OfReal& Knots)
759 Standard_Integer first = Knots.Lower();
760 Standard_Integer last = Knots.Upper();
761 Standard_Real kfirst = Knots(first);
762 Standard_Real klast = Knots(last);
763 Standard_Real tfirst = kfirst;
764 Standard_Real tlast = klast;
768 while (first <= last) {
769 tfirst += klast - Knots(last);
770 tlast -= Knots(first) - kfirst;
771 kfirst = Knots(first);
773 Knots(first) = tfirst;
780 //=======================================================================
783 //=======================================================================
785 void BSplCLib::Reverse(TColStd_Array1OfInteger& Mults)
787 Standard_Integer first = Mults.Lower();
788 Standard_Integer last = Mults.Upper();
789 Standard_Integer temp;
791 while (first < last) {
793 Mults(first) = Mults(last);
800 //=======================================================================
803 //=======================================================================
805 void BSplCLib::Reverse(TColStd_Array1OfReal& Weights,
806 const Standard_Integer L)
808 Standard_Integer i, l = L;
809 l = Weights.Lower()+(l-Weights.Lower())%(Weights.Upper()-Weights.Lower()+1);
811 TColStd_Array1OfReal temp(0,Weights.Length()-1);
813 for (i = Weights.Lower(); i <= l; i++)
814 temp(l-i) = Weights(i);
816 for (i = l+1; i <= Weights.Upper(); i++)
817 temp(l-Weights.Lower()+Weights.Upper()-i+1) = Weights(i);
819 for (i = Weights.Lower(); i <= Weights.Upper(); i++)
820 Weights(i) = temp(i-Weights.Lower());
823 //=======================================================================
824 //function : IsRational
826 //=======================================================================
828 Standard_Boolean BSplCLib::IsRational(const TColStd_Array1OfReal& Weights,
829 const Standard_Integer I1,
830 const Standard_Integer I2,
831 // const Standard_Real Epsi)
832 const Standard_Real )
834 Standard_Integer i, f = Weights.Lower(), l = Weights.Length();
835 Standard_Integer I3 = I2 - f;
836 const Standard_Real * WG = &Weights(f);
839 for (i = I1 - f; i < I3; i++) {
840 if (WG[f + (i % l)] != WG[f + ((i + 1) % l)]) return Standard_True;
842 return Standard_False ;
845 //=======================================================================
847 //purpose : evaluate point and derivatives
848 //=======================================================================
850 void BSplCLib::Eval(const Standard_Real U,
851 const Standard_Integer Degree,
852 Standard_Real& Knots,
853 const Standard_Integer Dimension,
854 Standard_Real& Poles)
856 Standard_Integer step,i,Dms,Dm1,Dpi,Sti;
857 Standard_Real X, Y, *poles, *knots = &Knots;
865 for (step = - 1; step < Dm1; step++) {
871 for (i = 0; i < Dms; i++) {
874 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
876 poles[0] *= X; poles[0] += Y * poles[1];
884 for (step = - 1; step < Dm1; step++) {
890 for (i = 0; i < Dms; i++) {
893 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
895 poles[0] *= X; poles[0] += Y * poles[2];
896 poles[1] *= X; poles[1] += Y * poles[3];
904 for (step = - 1; step < Dm1; step++) {
910 for (i = 0; i < Dms; i++) {
913 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
915 poles[0] *= X; poles[0] += Y * poles[3];
916 poles[1] *= X; poles[1] += Y * poles[4];
917 poles[2] *= X; poles[2] += Y * poles[5];
925 for (step = - 1; step < Dm1; step++) {
931 for (i = 0; i < Dms; i++) {
934 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
936 poles[0] *= X; poles[0] += Y * poles[4];
937 poles[1] *= X; poles[1] += Y * poles[5];
938 poles[2] *= X; poles[2] += Y * poles[6];
939 poles[3] *= X; poles[3] += Y * poles[7];
948 for (step = - 1; step < Dm1; step++) {
954 for (i = 0; i < Dms; i++) {
957 X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]);
960 for (k = 0; k < Dimension; k++) {
962 poles[k] += Y * poles[k + Dimension];
971 //=======================================================================
972 //function : BoorScheme
974 //=======================================================================
976 void BSplCLib::BoorScheme(const Standard_Real U,
977 const Standard_Integer Degree,
978 Standard_Real& Knots,
979 const Standard_Integer Dimension,
980 Standard_Real& Poles,
981 const Standard_Integer Depth,
982 const Standard_Integer Length)
985 // Compute the values
989 // for i = 0 to Depth,
990 // j = 0 to Length - i
992 // The Boor scheme is :
995 // P(i,j) = x * P(i-1,j) + (1-x) * P(i-1,j+1)
997 // where x = (knot(i+j+Degree) - U) / (knot(i+j+Degree) - knot(i+j))
1000 // The values are stored in the array Poles
1001 // They are alternatively written if the odd and even positions.
1003 // The successives contents of the array are
1004 // ***** means unitialised, l = Degree + Length
1006 // P(0,0) ****** P(0,1) ...... P(0,l-1) ******** P(0,l)
1007 // P(0,0) P(1,0) P(0,1) ...... P(0,l-1) P(1,l-1) P(0,l)
1008 // P(0,0) P(1,0) P(2,0) ...... P(2,l-1) P(1,l-1) P(0,l)
1011 Standard_Integer i,k,step;
1012 Standard_Real *knots = &Knots;
1013 Standard_Real *pole, *firstpole = &Poles - 2 * Dimension;
1014 // the steps of recursion
1016 for (step = 0; step < Depth; step++) {
1017 firstpole += Dimension;
1019 // compute the new row of poles
1021 for (i = step; i < Length; i++) {
1022 pole += 2 * Dimension;
1024 Standard_Real X = (knots[i+Degree-step] - U)
1025 / (knots[i+Degree-step] - knots[i]);
1026 Standard_Real Y = 1. - X;
1028 // P(i,j) = X * P(i-1,j) + (1-X) * P(i-1,j+1)
1030 for (k = 0; k < Dimension; k++)
1031 pole[k] = X * pole[k - Dimension] + Y * pole[k + Dimension];
1036 //=======================================================================
1037 //function : AntiBoorScheme
1039 //=======================================================================
1041 Standard_Boolean BSplCLib::AntiBoorScheme(const Standard_Real U,
1042 const Standard_Integer Degree,
1043 Standard_Real& Knots,
1044 const Standard_Integer Dimension,
1045 Standard_Real& Poles,
1046 const Standard_Integer Depth,
1047 const Standard_Integer Length,
1048 const Standard_Real Tolerance)
1050 // do the Boor scheme reverted.
1052 Standard_Integer i,k,step, half_length;
1053 Standard_Real *knots = &Knots;
1054 Standard_Real z,X,Y,*pole, *firstpole = &Poles + (Depth-1) * Dimension;
1056 // Test the special case length = 1
1057 // only verification of the central point
1060 X = (knots[Degree] - U) / (knots[Degree] - knots[0]);
1063 for (k = 0; k < Dimension; k++) {
1064 z = X * firstpole[k] + Y * firstpole[k+2*Dimension];
1065 if (Abs(z - firstpole[k+Dimension]) > Tolerance)
1066 return Standard_False;
1068 return Standard_True;
1072 // the steps of recursion
1074 for (step = Depth-1; step >= 0; step--) {
1075 firstpole -= Dimension;
1078 // first step from left to right
1080 for (i = step; i < Length-1; i++) {
1081 pole += 2 * Dimension;
1083 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1086 for (k = 0; k < Dimension; k++)
1087 pole[k+Dimension] = (pole[k] - X*pole[k-Dimension]) / Y;
1091 // second step from right to left
1092 pole += 4* Dimension;
1093 half_length = (Length - 1 + step) / 2 ;
1095 // only do half of the way from right to left
1096 // otherwise it start degenerating because of
1100 for (i = Length-1; i > half_length ; i--) {
1101 pole -= 2 * Dimension;
1104 X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]);
1107 for (k = 0; k < Dimension; k++) {
1108 z = (pole[k] - Y * pole[k+Dimension]) / X;
1109 if (Abs(z-pole[k-Dimension]) > Tolerance)
1110 return Standard_False;
1111 pole[k-Dimension] += z;
1112 pole[k-Dimension] /= 2.;
1116 return Standard_True;
1119 //=======================================================================
1120 //function : Derivative
1122 //=======================================================================
1124 void BSplCLib::Derivative(const Standard_Integer Degree,
1125 Standard_Real& Knots,
1126 const Standard_Integer Dimension,
1127 const Standard_Integer Length,
1128 const Standard_Integer Order,
1129 Standard_Real& Poles)
1131 Standard_Integer i,k,step,span = Degree;
1132 Standard_Real *knot = &Knots;
1134 for (step = 1; step <= Order; step++) {
1135 Standard_Real* pole = &Poles;
1137 for (i = step; i < Length; i++) {
1138 Standard_Real coef = - span / (knot[i+span] - knot[i]);
1140 for (k = 0; k < Dimension; k++) {
1141 pole[k] -= pole[k+Dimension];
1150 //=======================================================================
1153 //=======================================================================
1155 void BSplCLib::Bohm(const Standard_Real U,
1156 const Standard_Integer Degree,
1157 const Standard_Integer N,
1158 Standard_Real& Knots,
1159 const Standard_Integer Dimension,
1160 Standard_Real& Poles)
1162 // First phase independant of U, compute the poles of the derivatives
1163 Standard_Integer i,j,iDim,min,Dmi,DDmi,jDmi,Degm1;
1164 Standard_Real *knot = &Knots, *pole, coef, *tbis, *psav, *psDD, *psDDmDim;
1166 if (N < Degree) min = N;
1169 DDmi = (Degree << 1) + 1;
1170 switch (Dimension) {
1172 psDD = psav + Degree;
1173 psDDmDim = psDD - 1;
1175 for (i = 0; i < Degree; i++) {
1181 for (j = Degm1; j >= i; j--) {
1183 *pole -= *tbis; *pole /= (knot[jDmi] - knot[j]);
1188 // Second phase, dependant of U
1191 for (i = 0; i < Degree; i++) {
1197 for (j = i; j >= 0; j--) {
1198 *pole += coef * (*tbis);
1203 // multiply by the degrees
1208 for (i = 1; i <= min; i++) {
1209 *pole *= coef; pole++;
1216 psDD = psav + (Degree << 1);
1217 psDDmDim = psDD - 2;
1219 for (i = 0; i < Degree; i++) {
1225 for (j = Degm1; j >= i; j--) {
1227 coef = 1. / (knot[jDmi] - knot[j]);
1228 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1229 *pole -= *tbis; *pole *= coef;
1234 // Second phase, dependant of U
1237 for (i = 0; i < Degree; i++) {
1243 for (j = i; j >= 0; j--) {
1244 *pole += coef * (*tbis); pole++; tbis++;
1245 *pole += coef * (*tbis);
1250 // multiply by the degrees
1255 for (i = 1; i <= min; i++) {
1256 *pole *= coef; pole++;
1257 *pole *= coef; pole++;
1264 psDD = psav + (Degree << 1) + Degree;
1265 psDDmDim = psDD - 3;
1267 for (i = 0; i < Degree; i++) {
1273 for (j = Degm1; j >= i; j--) {
1275 coef = 1. / (knot[jDmi] - knot[j]);
1276 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1277 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1278 *pole -= *tbis; *pole *= coef;
1283 // Second phase, dependant of U
1286 for (i = 0; i < Degree; i++) {
1292 for (j = i; j >= 0; j--) {
1293 *pole += coef * (*tbis); pole++; tbis++;
1294 *pole += coef * (*tbis); pole++; tbis++;
1295 *pole += coef * (*tbis);
1300 // multiply by the degrees
1305 for (i = 1; i <= min; i++) {
1306 *pole *= coef; pole++;
1307 *pole *= coef; pole++;
1308 *pole *= coef; pole++;
1315 psDD = psav + (Degree << 2);
1316 psDDmDim = psDD - 4;
1318 for (i = 0; i < Degree; i++) {
1324 for (j = Degm1; j >= i; j--) {
1326 coef = 1. / (knot[jDmi] - knot[j]);
1327 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1328 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1329 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1330 *pole -= *tbis; *pole *= coef;
1335 // Second phase, dependant of U
1338 for (i = 0; i < Degree; i++) {
1344 for (j = i; j >= 0; j--) {
1345 *pole += coef * (*tbis); pole++; tbis++;
1346 *pole += coef * (*tbis); pole++; tbis++;
1347 *pole += coef * (*tbis); pole++; tbis++;
1348 *pole += coef * (*tbis);
1353 // multiply by the degrees
1358 for (i = 1; i <= min; i++) {
1359 *pole *= coef; pole++;
1360 *pole *= coef; pole++;
1361 *pole *= coef; pole++;
1362 *pole *= coef; pole++;
1370 Standard_Integer Dim2 = Dimension << 1;
1371 psDD = psav + Degree * Dimension;
1372 psDDmDim = psDD - Dimension;
1374 for (i = 0; i < Degree; i++) {
1380 for (j = Degm1; j >= i; j--) {
1382 coef = 1. / (knot[jDmi] - knot[j]);
1384 for (k = 0; k < Dimension; k++) {
1385 *pole -= *tbis; *pole *= coef; pole++; tbis++;
1391 // Second phase, dependant of U
1394 for (i = 0; i < Degree; i++) {
1397 tbis = pole + Dimension;
1400 for (j = i; j >= 0; j--) {
1402 for (k = 0; k < Dimension; k++) {
1403 *pole += coef * (*tbis); pole++; tbis++;
1409 // multiply by the degrees
1412 pole = psav + Dimension;
1414 for (i = 1; i <= min; i++) {
1416 for (k = 0; k < Dimension; k++) {
1417 *pole *= coef; pole++;
1426 //=======================================================================
1427 //function : BuildKnots
1429 //=======================================================================
1431 void BSplCLib::BuildKnots(const Standard_Integer Degree,
1432 const Standard_Integer Index,
1433 const Standard_Boolean Periodic,
1434 const TColStd_Array1OfReal& Knots,
1435 const TColStd_Array1OfInteger& Mults,
1438 Standard_Integer KLower = Knots.Lower();
1439 const Standard_Real * pkn = &Knots(KLower);
1441 Standard_Real *knot = &LK;
1442 if (&Mults == NULL) {
1445 Standard_Integer j = Index ;
1446 knot[0] = pkn[j]; j++;
1451 Standard_Integer j = Index - 1;
1452 knot[0] = pkn[j]; j++;
1453 knot[1] = pkn[j]; j++;
1454 knot[2] = pkn[j]; j++;
1459 Standard_Integer j = Index - 2;
1460 knot[0] = pkn[j]; j++;
1461 knot[1] = pkn[j]; j++;
1462 knot[2] = pkn[j]; j++;
1463 knot[3] = pkn[j]; j++;
1464 knot[4] = pkn[j]; j++;
1469 Standard_Integer j = Index - 3;
1470 knot[0] = pkn[j]; j++;
1471 knot[1] = pkn[j]; j++;
1472 knot[2] = pkn[j]; j++;
1473 knot[3] = pkn[j]; j++;
1474 knot[4] = pkn[j]; j++;
1475 knot[5] = pkn[j]; j++;
1476 knot[6] = pkn[j]; j++;
1481 Standard_Integer j = Index - 4;
1482 knot[0] = pkn[j]; j++;
1483 knot[1] = pkn[j]; j++;
1484 knot[2] = pkn[j]; j++;
1485 knot[3] = pkn[j]; j++;
1486 knot[4] = pkn[j]; j++;
1487 knot[5] = pkn[j]; j++;
1488 knot[6] = pkn[j]; j++;
1489 knot[7] = pkn[j]; j++;
1490 knot[8] = pkn[j]; j++;
1495 Standard_Integer j = Index - 5;
1496 knot[ 0] = pkn[j]; j++;
1497 knot[ 1] = pkn[j]; j++;
1498 knot[ 2] = pkn[j]; j++;
1499 knot[ 3] = pkn[j]; j++;
1500 knot[ 4] = pkn[j]; j++;
1501 knot[ 5] = pkn[j]; j++;
1502 knot[ 6] = pkn[j]; j++;
1503 knot[ 7] = pkn[j]; j++;
1504 knot[ 8] = pkn[j]; j++;
1505 knot[ 9] = pkn[j]; j++;
1506 knot[10] = pkn[j]; j++;
1511 Standard_Integer i,j;
1512 Standard_Integer Deg2 = Degree << 1;
1515 for (i = 0; i < Deg2; i++) {
1524 Standard_Integer Deg1 = Degree - 1;
1525 Standard_Integer KUpper = Knots.Upper();
1526 Standard_Integer MLower = Mults.Lower();
1527 Standard_Integer MUpper = Mults.Upper();
1528 const Standard_Integer * pmu = &Mults(MLower);
1530 Standard_Real dknot = 0;
1531 Standard_Integer ilow = Index , mlow = 0;
1532 Standard_Integer iupp = Index + 1, mupp = 0;
1533 Standard_Real loffset = 0., uoffset = 0.;
1534 Standard_Boolean getlow = Standard_True, getupp = Standard_True;
1536 dknot = pkn[KUpper] - pkn[KLower];
1537 if (iupp > MUpper) {
1542 // Find the knots around Index
1544 for (i = 0; i < Degree; i++) {
1547 if (mlow > pmu[ilow]) {
1550 getlow = (ilow >= MLower);
1551 if (Periodic && !getlow) {
1554 getlow = Standard_True;
1558 knot[Deg1 - i] = pkn[ilow] - loffset;
1562 if (mupp > pmu[iupp]) {
1565 getupp = (iupp <= MUpper);
1566 if (Periodic && !getupp) {
1569 getupp = Standard_True;
1573 knot[Degree + i] = pkn[iupp] + uoffset;
1579 //=======================================================================
1580 //function : PoleIndex
1582 //=======================================================================
1584 Standard_Integer BSplCLib::PoleIndex(const Standard_Integer Degree,
1585 const Standard_Integer Index,
1586 const Standard_Boolean Periodic,
1587 const TColStd_Array1OfInteger& Mults)
1589 Standard_Integer i, pindex = 0;
1591 for (i = Mults.Lower(); i <= Index; i++)
1594 pindex -= Mults(Mults.Lower());
1596 pindex -= Degree + 1;
1601 //=======================================================================
1602 //function : BuildBoor
1603 //purpose : builds the local array for boor
1604 //=======================================================================
1606 void BSplCLib::BuildBoor(const Standard_Integer Index,
1607 const Standard_Integer Length,
1608 const Standard_Integer Dimension,
1609 const TColStd_Array1OfReal& Poles,
1612 Standard_Real *poles = &LP;
1613 Standard_Integer i,k, ip = Poles.Lower() + Index * Dimension;
1615 for (i = 0; i < Length+1; i++) {
1617 for (k = 0; k < Dimension; k++) {
1618 poles[k] = Poles(ip);
1620 if (ip > Poles.Upper()) ip = Poles.Lower();
1622 poles += 2 * Dimension;
1626 //=======================================================================
1627 //function : BoorIndex
1629 //=======================================================================
1631 Standard_Integer BSplCLib::BoorIndex(const Standard_Integer Index,
1632 const Standard_Integer Length,
1633 const Standard_Integer Depth)
1635 if (Index <= Depth) return Index;
1636 if (Index <= Length) return 2 * Index - Depth;
1637 return Length + Index - Depth;
1640 //=======================================================================
1641 //function : GetPole
1643 //=======================================================================
1645 void BSplCLib::GetPole(const Standard_Integer Index,
1646 const Standard_Integer Length,
1647 const Standard_Integer Depth,
1648 const Standard_Integer Dimension,
1650 Standard_Integer& Position,
1651 TColStd_Array1OfReal& Pole)
1654 Standard_Real* pole = &LP + BoorIndex(Index,Length,Depth) * Dimension;
1656 for (k = 0; k < Dimension; k++) {
1657 Pole(Position) = pole[k];
1660 if (Position > Pole.Upper()) Position = Pole.Lower();
1663 //=======================================================================
1664 //function : PrepareInsertKnots
1666 //=======================================================================
1668 Standard_Boolean BSplCLib::PrepareInsertKnots
1669 (const Standard_Integer Degree,
1670 const Standard_Boolean Periodic,
1671 const TColStd_Array1OfReal& Knots,
1672 const TColStd_Array1OfInteger& Mults,
1673 const TColStd_Array1OfReal& AddKnots,
1674 const TColStd_Array1OfInteger& AddMults,
1675 Standard_Integer& NbPoles,
1676 Standard_Integer& NbKnots,
1677 const Standard_Real Tolerance,
1678 const Standard_Boolean Add)
1680 Standard_Boolean addflat = &AddMults == NULL;
1682 Standard_Integer first,last;
1684 first = Knots.Lower();
1685 last = Knots.Upper();
1688 first = FirstUKnotIndex(Degree,Mults);
1689 last = LastUKnotIndex(Degree,Mults);
1691 Standard_Real adeltaK1 = Knots(first)-AddKnots(AddKnots.Lower());
1692 Standard_Real adeltaK2 = AddKnots(AddKnots.Upper())-Knots(last);
1693 if (adeltaK1 > Tolerance) return Standard_False;
1694 if (adeltaK2 > Tolerance) return Standard_False;
1696 Standard_Integer sigma = 0, mult, amult;
1698 Standard_Integer k = Knots.Lower() - 1;
1699 Standard_Integer ak = AddKnots.Lower();
1701 if(Periodic && AddKnots.Length() > 1)
1703 //gka for case when segments was produced on full period only one knot
1704 //was added in the end of curve
1705 if(fabs(adeltaK1) <= Precision::PConfusion() &&
1706 fabs(adeltaK2) <= Precision::PConfusion())
1710 Standard_Real au,oldau = AddKnots(ak),Eps;
1712 while (ak <= AddKnots.Upper()) {
1714 if (au < oldau) return Standard_False;
1717 Eps = Max(Tolerance,Epsilon(au));
1719 while ((k < Knots.Upper()) && (Knots(k+1) - au <= Eps)) {
1725 if (addflat) amult = 1;
1726 else amult = Max(0,AddMults(ak));
1728 while ((ak < AddKnots.Upper()) &&
1729 (Abs(au - AddKnots(ak+1)) <= Eps)) {
1732 if (addflat) amult++;
1733 else amult += Max(0,AddMults(ak));
1738 if (Abs(au - Knots(k)) <= Eps) {
1739 // identic to existing knot
1742 if (mult + amult > Degree)
1743 amult = Max(0,Degree - mult);
1747 else if (amult > mult) {
1748 if (amult > Degree) amult = Degree;
1749 sigma += amult - mult;
1752 // on periodic curves if this is the last knot
1753 // the multiplicity is added twice to account for the first knot
1754 if (k == Knots.Upper() && Periodic) {
1758 sigma += amult - mult;
1763 // not identic to existing knot
1765 if (amult > Degree) amult = Degree;
1774 // count the last knots
1775 while (k < Knots.Upper()) {
1782 //for periodic B-Spline the requirement is that multiplicites of the first
1783 //and last knots must be equal (see Geom_BSplineCurve constructor for
1785 //respectively AddMults() must meet this requirement if AddKnots() contains
1786 //knot(s) coincident with first or last
1787 NbPoles = sigma - Mults(Knots.Upper());
1790 NbPoles = sigma - Degree - 1;
1793 return Standard_True;
1796 //=======================================================================
1798 //purpose : copy reals from an array to an other
1800 // NbValues are copied from OldPoles(OldFirst)
1801 // to NewPoles(NewFirst)
1803 // Periodicity is handled.
1804 // OldFirst and NewFirst are updated
1805 // to the position after the last copied pole.
1807 //=======================================================================
1809 static void Copy(const Standard_Integer NbPoles,
1810 Standard_Integer& OldFirst,
1811 const TColStd_Array1OfReal& OldPoles,
1812 Standard_Integer& NewFirst,
1813 TColStd_Array1OfReal& NewPoles)
1815 // reset the index in the range for periodicity
1817 OldFirst = OldPoles.Lower() +
1818 (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1);
1820 NewFirst = NewPoles.Lower() +
1821 (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1);
1826 for (i = 1; i <= NbPoles; i++) {
1827 NewPoles(NewFirst) = OldPoles(OldFirst);
1829 if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower();
1831 if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower();
1835 //=======================================================================
1836 //function : InsertKnots
1837 //purpose : insert an array of knots and multiplicities
1838 //=======================================================================
1840 void BSplCLib::InsertKnots
1841 (const Standard_Integer Degree,
1842 const Standard_Boolean Periodic,
1843 const Standard_Integer Dimension,
1844 const TColStd_Array1OfReal& Poles,
1845 const TColStd_Array1OfReal& Knots,
1846 const TColStd_Array1OfInteger& Mults,
1847 const TColStd_Array1OfReal& AddKnots,
1848 const TColStd_Array1OfInteger& AddMults,
1849 TColStd_Array1OfReal& NewPoles,
1850 TColStd_Array1OfReal& NewKnots,
1851 TColStd_Array1OfInteger& NewMults,
1852 const Standard_Real Tolerance,
1853 const Standard_Boolean Add)
1855 Standard_Boolean addflat = &AddMults == NULL;
1857 Standard_Integer i,k,mult,firstmult;
1858 Standard_Integer index,kn,curnk,curk;
1859 Standard_Integer p,np, curp, curnp, length, depth;
1861 Standard_Integer need;
1864 // -------------------
1865 // create local arrays
1866 // -------------------
1868 Standard_Real *knots = new Standard_Real[2*Degree];
1869 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
1871 //----------------------------
1872 // loop on the knots to insert
1873 //----------------------------
1875 curk = Knots.Lower()-1; // current position in Knots
1876 curnk = NewKnots.Lower()-1; // current position in NewKnots
1877 curp = Poles.Lower(); // current position in Poles
1878 curnp = NewPoles.Lower(); // current position in NewPoles
1880 // NewKnots, NewMults, NewPoles contains the current state of the curve
1882 // index is the first pole of the current curve for insertion schema
1884 if (Periodic) index = -Mults(Mults.Lower());
1885 else index = -Degree-1;
1887 // on Periodic curves the first knot and the last knot are inserted later
1888 // (they are the same knot)
1889 firstmult = 0; // multiplicity of the first-last knot for periodic
1892 // kn current knot to insert in AddKnots
1894 for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) {
1897 Eps = Max(Tolerance,Epsilon(u));
1899 //-----------------------------------
1900 // find the position in the old knots
1901 // and copy to the new knots
1902 //-----------------------------------
1904 while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) {
1906 NewKnots(curnk) = Knots(curk);
1907 index += NewMults(curnk) = Mults(curk);
1910 //-----------------------------------
1911 // Slice the knots and the mults
1912 // to the current size of the new curve
1913 //-----------------------------------
1915 i = curnk + Knots.Upper() - curk;
1916 TColStd_Array1OfReal nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i);
1917 TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i);
1919 //-----------------------------------
1920 // copy enough knots
1921 // to compute the insertion schema
1922 //-----------------------------------
1928 while (mult < Degree && k < Knots.Upper()) {
1930 nknots(i) = Knots(k);
1931 mult += nmults(i) = Mults(k);
1934 // copy knots at the end for periodic curve
1940 while (mult < Degree && i > curnk) {
1941 nknots(i) = Knots(k);
1942 mult += nmults(i) = Mults(k);
1946 nmults(nmults.Upper()) = nmults(nmults.Lower());
1951 //------------------------------------
1952 // do the boor scheme on the new curve
1953 // to insert the new knot
1954 //------------------------------------
1956 Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps);
1958 if (sameknot) length = Max(0,Degree - NewMults(curnk));
1959 else length = Degree;
1961 if (addflat) depth = 1;
1962 else depth = Min(Degree,AddMults(kn));
1966 if ((NewMults(curnk) + depth) > Degree)
1967 depth = Degree - NewMults(curnk);
1970 depth = Max(0,depth-NewMults(curnk));
1974 // on periodic curve the first and last knot are delayed to the end
1975 if (curk == Knots.Lower() || (curk == Knots.Upper())) {
1981 if (depth <= 0) continue;
1983 BuildKnots(Degree,curnk,Periodic,nknots,nmults,*knots);
1987 need = NewPoles.Lower()+(index+length+1)*Dimension - curnp;
1988 need = Min(need,Poles.Upper() - curp + 1);
1992 Copy(need,p,Poles,np,NewPoles);
1996 // slice the poles to the current number of poles in case of periodic
1997 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1999 BuildBoor(index,length,Dimension,npoles,*poles);
2000 BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length);
2002 //-------------------
2003 // copy the new poles
2004 //-------------------
2006 curnp += depth * Dimension; // number of poles is increased by depth
2007 TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
2008 np = NewKnots.Lower()+(index+1)*Dimension;
2010 for (i = 1; i <= length + depth; i++)
2011 GetPole(i,length,depth,Dimension,*poles,np,ThePoles);
2013 //-------------------
2015 //-------------------
2019 NewMults(curnk) += depth;
2023 NewKnots(curnk) = u;
2024 NewMults(curnk) = depth;
2028 //------------------------------
2029 // copy the last poles and knots
2030 //------------------------------
2032 Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles);
2034 while (curk < Knots.Upper()) {
2036 NewKnots(curnk) = Knots(curk);
2037 NewMults(curnk) = Mults(curk);
2040 //------------------------------
2041 // process the first-last knot
2042 // on periodic curves
2043 //------------------------------
2045 if (firstmult > 0) {
2046 curnk = NewKnots.Lower();
2047 if (NewMults(curnk) + firstmult > Degree) {
2048 firstmult = Degree - NewMults(curnk);
2050 if (firstmult > 0) {
2052 length = Degree - NewMults(curnk);
2055 BuildKnots(Degree,curnk,Periodic,NewKnots,NewMults,*knots);
2056 TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),
2058 NewPoles.Upper()-depth*Dimension);
2059 BuildBoor(0,length,Dimension,npoles,*poles);
2060 BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length);
2062 //---------------------------
2063 // copy the new poles
2064 // but rotate them with depth
2065 //---------------------------
2067 np = NewPoles.Lower();
2069 for (i = depth; i < length + depth; i++)
2070 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2072 np = NewPoles.Upper() - depth*Dimension + 1;
2074 for (i = 0; i < depth; i++)
2075 GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2077 NewMults(NewMults.Lower()) += depth;
2078 NewMults(NewMults.Upper()) += depth;
2081 // free local arrays
2086 //=======================================================================
2087 //function : RemoveKnot
2089 //=======================================================================
2091 Standard_Boolean BSplCLib::RemoveKnot
2092 (const Standard_Integer Index,
2093 const Standard_Integer Mult,
2094 const Standard_Integer Degree,
2095 const Standard_Boolean Periodic,
2096 const Standard_Integer Dimension,
2097 const TColStd_Array1OfReal& Poles,
2098 const TColStd_Array1OfReal& Knots,
2099 const TColStd_Array1OfInteger& Mults,
2100 TColStd_Array1OfReal& NewPoles,
2101 TColStd_Array1OfReal& NewKnots,
2102 TColStd_Array1OfInteger& NewMults,
2103 const Standard_Real Tolerance)
2105 Standard_Integer index,i,j,k,p,np;
2107 Standard_Integer TheIndex = Index;
2110 Standard_Integer first,last;
2112 first = Knots.Lower();
2113 last = Knots.Upper();
2116 first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1;
2117 last = BSplCLib::LastUKnotIndex(Degree,Mults) - 1;
2119 if (Index < first) return Standard_False;
2120 if (Index > last) return Standard_False;
2122 if ( Periodic && (Index == first)) TheIndex = last;
2124 Standard_Integer depth = Mults(TheIndex) - Mult;
2125 Standard_Integer length = Degree - Mult;
2127 // -------------------
2128 // create local arrays
2129 // -------------------
2131 Standard_Real *knots = new Standard_Real[4*Degree];
2132 Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
2135 // ------------------------------------
2136 // build the knots for anti Boor Scheme
2137 // ------------------------------------
2139 // the new sequence of knots
2140 // is obtained from the knots at Index-1 and Index
2142 BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,Mults,*knots);
2143 index = PoleIndex(Degree,TheIndex-1,Periodic,Mults);
2144 BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,Mults,knots[2*Degree]);
2148 for (i = 0; i < Degree - Mult; i++)
2149 knots[i] = knots[i+Mult];
2151 for (i = Degree-Mult; i < 2*Degree; i++)
2152 knots[i] = knots[2*Degree+i];
2155 // ------------------------------------
2156 // build the poles for anti Boor Scheme
2157 // ------------------------------------
2159 p = Poles.Lower()+index * Dimension;
2161 for (i = 0; i <= length + depth; i++) {
2162 j = Dimension * BoorIndex(i,length,depth);
2164 for (k = 0; k < Dimension; k++) {
2165 poles[j+k] = Poles(p+k);
2168 if (p > Poles.Upper()) p = Poles.Lower();
2176 Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots,
2178 depth,length,Tolerance);
2189 np = NewPoles.Lower();
2191 // unmodified poles before
2192 Copy((index+1)*Dimension,p,Poles,np,NewPoles);
2197 for (i = 1; i <= length; i++)
2198 BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles);
2199 p += (length + depth) * Dimension ;
2201 // unmodified poles after
2202 if (p != Poles.Lower()) {
2203 i = Poles.Upper() - p + 1;
2204 Copy(i,p,Poles,np,NewPoles);
2212 NewMults(TheIndex) = Mult;
2214 if (TheIndex == first) NewMults(last) = Mult;
2215 if (TheIndex == last) NewMults(first) = Mult;
2219 if (!Periodic || (TheIndex != first && TheIndex != last)) {
2221 for (i = Knots.Lower(); i < TheIndex; i++) {
2222 NewKnots(i) = Knots(i);
2223 NewMults(i) = Mults(i);
2226 for (i = TheIndex+1; i <= Knots.Upper(); i++) {
2227 NewKnots(i-1) = Knots(i);
2228 NewMults(i-1) = Mults(i);
2232 // The interesting case of a Periodic curve
2233 // where the first and last knot is removed.
2235 for (i = first; i < last-1; i++) {
2236 NewKnots(i) = Knots(i+1);
2237 NewMults(i) = Mults(i+1);
2239 NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first);
2240 NewMults(last-1) = NewMults(first);
2246 // free local arrays
2253 //=======================================================================
2254 //function : IncreaseDegreeCountKnots
2256 //=======================================================================
2258 Standard_Integer BSplCLib::IncreaseDegreeCountKnots
2259 (const Standard_Integer Degree,
2260 const Standard_Integer NewDegree,
2261 const Standard_Boolean Periodic,
2262 const TColStd_Array1OfInteger& Mults)
2264 if (Periodic) return Mults.Length();
2265 Standard_Integer f = FirstUKnotIndex(Degree,Mults);
2266 Standard_Integer l = LastUKnotIndex(Degree,Mults);
2267 Standard_Integer m,i,removed = 0, step = NewDegree - Degree;
2270 m = Degree + (f - i + 1) * step + 1;
2272 while (m > NewDegree+1) {
2274 m -= Mults(i) + step;
2277 if (m < NewDegree+1) removed--;
2280 m = Degree + (i - l + 1) * step + 1;
2282 while (m > NewDegree+1) {
2284 m -= Mults(i) + step;
2287 if (m < NewDegree+1) removed--;
2289 return Mults.Length() - removed;
2292 //=======================================================================
2293 //function : IncreaseDegree
2295 //=======================================================================
2297 void BSplCLib::IncreaseDegree
2298 (const Standard_Integer Degree,
2299 const Standard_Integer NewDegree,
2300 const Standard_Boolean Periodic,
2301 const Standard_Integer Dimension,
2302 const TColStd_Array1OfReal& Poles,
2303 const TColStd_Array1OfReal& Knots,
2304 const TColStd_Array1OfInteger& Mults,
2305 TColStd_Array1OfReal& NewPoles,
2306 TColStd_Array1OfReal& NewKnots,
2307 TColStd_Array1OfInteger& NewMults)
2309 // Degree elevation of a BSpline Curve
2311 // This algorithms loops on degree incrementation from Degree to NewDegree.
2312 // The variable curDeg is the current degree to increment.
2314 // Before degree incrementations a "working curve" is created.
2315 // It has the same knot, poles and multiplicities.
2317 // If the curve is periodic knots are added on the working curve before
2318 // and after the existing knots to make it a non-periodic curves.
2319 // The poles are also copied.
2321 // The first and last multiplicity of the working curve are set to Degree+1,
2322 // null poles are added if necessary.
2324 // Then the degree is incremented on the working curve.
2325 // The knots are unchanged but all multiplicities will be incremented.
2327 // Each degree incrementation is achieved by averaging curDeg+1 curves.
2329 // See : Degree elevation of B-spline curves
2330 // Hartmut PRAUTZSCH
2334 //-------------------------
2335 // create the working curve
2336 //-------------------------
2338 Standard_Integer i,k,f,l,m,pf,pl,firstknot;
2340 pf = 0; // number of null poles added at beginning
2341 pl = 0; // number of null poles added at end
2343 Standard_Integer nbwknots = Knots.Length();
2344 f = FirstUKnotIndex(Degree,Mults);
2345 l = LastUKnotIndex (Degree,Mults);
2348 // Periodic curves are transformed in non-periodic curves
2350 nbwknots += f - Mults.Lower();
2354 for (i = Mults.Lower(); i <= f; i++)
2357 nbwknots += Mults.Upper() - l;
2361 for (i = l; i <= Mults.Upper(); i++)
2365 // copy the knots and multiplicities
2366 TColStd_Array1OfReal wknots(1,nbwknots);
2367 TColStd_Array1OfInteger wmults(1,nbwknots);
2373 // copy the knots for a periodic curve
2374 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2377 for (k = l; k < Knots.Upper(); k++) {
2379 wknots(i) = Knots(k) - period;
2380 wmults(i) = Mults(k);
2383 for (k = Knots.Lower(); k <= Knots.Upper(); k++) {
2385 wknots(i) = Knots(k);
2386 wmults(i) = Mults(k);
2389 for (k = Knots.Lower()+1; k <= f; k++) {
2391 wknots(i) = Knots(k) + period;
2392 wmults(i) = Mults(k);
2396 // set the first and last mults to Degree+1
2397 // and add null poles
2399 pf += Degree + 1 - wmults(1);
2400 wmults(1) = Degree + 1;
2401 pl += Degree + 1 - wmults(nbwknots);
2402 wmults(nbwknots) = Degree + 1;
2404 //---------------------------
2405 // poles of the working curve
2406 //---------------------------
2408 Standard_Integer nbwpoles = 0;
2410 for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i);
2411 nbwpoles -= Degree + 1;
2413 // we provide space for degree elevation
2414 TColStd_Array1OfReal
2415 wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension);
2417 for (i = 1; i <= pf * Dimension; i++)
2422 for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) {
2423 wpoles(i) = Poles(k);
2425 if (k > Poles.Upper()) k = Poles.Lower();
2428 for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++)
2432 //-----------------------------------------------------------
2433 // Declare the temporary arrays used in degree incrementation
2434 //-----------------------------------------------------------
2436 Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree);
2437 // Arrays for storing the temporary curves
2438 TColStd_Array1OfReal tempc1(1,nbwp * Dimension);
2439 TColStd_Array1OfReal tempc2(1,nbwp * Dimension);
2441 // Array for storing the knots to insert
2442 TColStd_Array1OfReal iknots(1,nbwknots);
2444 // Arrays for receiving the knots after insertion
2445 TColStd_Array1OfReal nknots(1,nbwknots);
2449 //------------------------------
2450 // Loop on degree incrementation
2451 //------------------------------
2453 Standard_Integer step,curDeg;
2454 Standard_Integer nbp = nbwpoles;
2457 for (curDeg = Degree; curDeg < NewDegree; curDeg++) {
2459 nbp = nbwp; // current number of poles
2460 nbwp = nbp + nbwknots - 1; // new number of poles
2462 // For the averaging
2463 TColStd_Array1OfReal nwpoles(1,nbwp * Dimension);
2464 nwpoles.Init(0.0e0) ;
2467 for (step = 0; step <= curDeg; step++) {
2469 // Compute the bspline of rank step.
2471 // if not the first time, decrement the multiplicities back
2473 for (i = 1; i <= nbwknots; i++)
2477 // Poles are the current poles
2478 // but the poles congruent to step are duplicated.
2480 Standard_Integer offset = 0;
2482 for (i = 0; i < nbp; i++) {
2485 for (k = 0; k < Dimension; k++) {
2486 tempc1((offset-1)*Dimension+k+1) =
2487 wpoles(NewPoles.Lower()+i*Dimension+k);
2489 if (i % (curDeg+1) == step) {
2492 for (k = 0; k < Dimension; k++) {
2493 tempc1((offset-1)*Dimension+k+1) =
2494 wpoles(NewPoles.Lower()+i*Dimension+k);
2499 // Knots multiplicities are increased
2500 // For each knot where the sum of multiplicities is congruent to step
2502 Standard_Integer stepmult = step+1;
2503 Standard_Integer nbknots = 0;
2504 Standard_Integer smult = 0;
2506 for (k = 1; k <= nbwknots; k++) {
2508 if (smult >= stepmult) {
2509 // this knot is increased
2510 stepmult += curDeg+1;
2514 // this knot is inserted
2516 iknots(nbknots) = wknots(k);
2520 // the curve is obtained by inserting the knots
2521 // to raise the multiplicities
2523 // we build "slices" of the arrays to set the correct size
2525 TColStd_Array1OfReal aknots(iknots(1),1,nbknots);
2526 TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension);
2527 TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp * Dimension);
2528 // InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2529 // aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.));
2531 InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2532 aknots,NoMults(),ncurve,nknots,wmults,0.0);
2534 // add to the average
2536 for (i = 1; i <= nbwp * Dimension; i++)
2537 nwpoles(i) += ncurve(i);
2540 // add to the average
2542 for (i = 1; i <= nbwp * Dimension; i++)
2543 nwpoles(i) += tempc1(i);
2547 // The result is the average
2549 for (i = 1; i <= nbwp * Dimension; i++) {
2550 wpoles(i) = nwpoles(i) / (curDeg+1);
2558 // index in new knots of the first knot of the curve
2560 firstknot = Mults.Upper() - l + 1;
2564 // the new curve starts at index firstknot
2565 // so we must remove knots until the sum of multiplicities
2566 // from the first to the start is NewDegree+1
2568 // m is the current sum of multiplicities
2571 for (k = 1; k <= firstknot; k++)
2574 // compute the new first knot (k), pf will be the index of the first pole
2578 while (m > NewDegree+1) {
2583 if (m < NewDegree+1) {
2585 wmults(k) += m - NewDegree - 1;
2586 pf += m - NewDegree - 1;
2589 // on a periodic curve the knots start with firstknot
2595 for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) {
2596 NewKnots(i) = wknots(k);
2597 NewMults(i) = wmults(k);
2604 for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) {
2606 NewPoles(i) = wpoles(pf);
2610 //=======================================================================
2611 //function : PrepareUnperiodize
2613 //=======================================================================
2615 void BSplCLib::PrepareUnperiodize
2616 (const Standard_Integer Degree,
2617 const TColStd_Array1OfInteger& Mults,
2618 Standard_Integer& NbKnots,
2619 Standard_Integer& NbPoles)
2622 // initialize NbKnots and NbPoles
2623 NbKnots = Mults.Length();
2624 NbPoles = - Degree - 1;
2626 for (i = Mults.Lower(); i <= Mults.Upper(); i++)
2627 NbPoles += Mults(i);
2629 Standard_Integer sigma, k;
2630 // Add knots at the beginning of the curve to raise Multiplicities
2632 sigma = Mults(Mults.Lower());
2633 k = Mults.Upper() - 1;
2635 while ( sigma < Degree + 1) {
2637 NbPoles += Mults(k);
2641 // We must add exactly until Degree + 1 ->
2642 // Supress the excedent.
2643 if ( sigma > Degree + 1)
2644 NbPoles -= sigma - Degree - 1;
2646 // Add knots at the end of the curve to raise Multiplicities
2648 sigma = Mults(Mults.Upper());
2649 k = Mults.Lower() + 1;
2651 while ( sigma < Degree + 1) {
2653 NbPoles += Mults(k);
2657 // We must add exactly until Degree + 1 ->
2658 // Supress the excedent.
2659 if ( sigma > Degree + 1)
2660 NbPoles -= sigma - Degree - 1;
2663 //=======================================================================
2664 //function : Unperiodize
2666 //=======================================================================
2668 void BSplCLib::Unperiodize
2669 (const Standard_Integer Degree,
2670 const Standard_Integer , // Dimension,
2671 const TColStd_Array1OfInteger& Mults,
2672 const TColStd_Array1OfReal& Knots,
2673 const TColStd_Array1OfReal& Poles,
2674 TColStd_Array1OfInteger& NewMults,
2675 TColStd_Array1OfReal& NewKnots,
2676 TColStd_Array1OfReal& NewPoles)
2678 Standard_Integer sigma, k, index = 0;
2679 // evaluation of index : number of knots to insert before knot(1) to
2680 // raise sum of multiplicities to <Degree + 1>
2681 sigma = Mults(Mults.Lower());
2682 k = Mults.Upper() - 1;
2684 while ( sigma < Degree + 1) {
2690 Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2692 // set the 'interior' knots;
2694 for ( k = 1; k <= Knots.Length(); k++) {
2695 NewKnots ( k + index ) = Knots( k);
2696 NewMults ( k + index ) = Mults( k);
2699 // set the 'starting' knots;
2701 for ( k = 1; k <= index; k++) {
2702 NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period;
2703 NewMults( k) = NewMults( k + Knots.Length() - 1);
2705 NewMults( 1) -= sigma - Degree -1;
2707 // set the 'ending' knots;
2708 sigma = NewMults( index + Knots.Length() );
2710 for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) {
2711 NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period;
2712 NewMults( k) = NewMults( k - Knots.Length() + 1);
2713 sigma += NewMults( k - Knots.Length() + 1);
2715 NewMults(NewMults.Length()) -= sigma - Degree - 1;
2717 for ( k = 1 ; k <= NewPoles.Length(); k++) {
2718 NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1);
2722 //=======================================================================
2723 //function : PrepareTrimming
2725 //=======================================================================
2727 void BSplCLib::PrepareTrimming(const Standard_Integer Degree,
2728 const Standard_Boolean Periodic,
2729 const TColStd_Array1OfReal& Knots,
2730 const TColStd_Array1OfInteger& Mults,
2731 const Standard_Real U1,
2732 const Standard_Real U2,
2733 Standard_Integer& NbKnots,
2734 Standard_Integer& NbPoles)
2737 Standard_Real NewU1, NewU2;
2738 Standard_Integer index1 = 0, index2 = 0;
2740 // Eval index1, index2 : position of U1 and U2 in the Array Knots
2741 // such as Knots(index1-1) <= U1 < Knots(index1)
2742 // Knots(index2-1) <= U2 < Knots(index2)
2743 LocateParameter( Degree, Knots, Mults, U1, Periodic,
2744 Knots.Lower(), Knots.Upper(), index1, NewU1);
2745 LocateParameter( Degree, Knots, Mults, U2, Periodic,
2746 Knots.Lower(), Knots.Upper(), index2, NewU2);
2748 if ( Abs(Knots(index2) - U2) <= Epsilon( U1))
2752 NbKnots = index2 - index1 + 3;
2755 NbPoles = Degree + 1;
2757 for ( i = index1; i <= index2; i++)
2758 NbPoles += Mults(i);
2761 //=======================================================================
2762 //function : Trimming
2764 //=======================================================================
2766 void BSplCLib::Trimming(const Standard_Integer Degree,
2767 const Standard_Boolean Periodic,
2768 const Standard_Integer Dimension,
2769 const TColStd_Array1OfReal& Knots,
2770 const TColStd_Array1OfInteger& Mults,
2771 const TColStd_Array1OfReal& Poles,
2772 const Standard_Real U1,
2773 const Standard_Real U2,
2774 TColStd_Array1OfReal& NewKnots,
2775 TColStd_Array1OfInteger& NewMults,
2776 TColStd_Array1OfReal& NewPoles)
2778 Standard_Integer i, nbpoles, nbknots;
2779 Standard_Real kk[2];
2780 Standard_Integer mm[2];
2781 TColStd_Array1OfReal K( kk[0], 1, 2 );
2782 TColStd_Array1OfInteger M( mm[0], 1, 2 );
2784 K(1) = U1; K(2) = U2;
2785 mm[0] = mm[1] = Degree;
2786 if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, M,
2787 nbpoles, nbknots, Epsilon( U1), 0))
2788 Standard_OutOfRange::Raise();
2790 TColStd_Array1OfReal TempPoles(1, nbpoles*Dimension);
2791 TColStd_Array1OfReal TempKnots(1, nbknots);
2792 TColStd_Array1OfInteger TempMults(1, nbknots);
2795 // do not allow the multiplicities to Add : they must be less than Degree
2797 InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults,
2798 K, M, TempPoles, TempKnots, TempMults, Epsilon(U1),
2801 // find in TempPoles the index of the pole corresponding to U1
2802 Standard_Integer Kindex = 0, Pindex;
2803 Standard_Real NewU1;
2804 LocateParameter( Degree, TempKnots, TempMults, U1, Periodic,
2805 TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1);
2806 Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults);
2807 Pindex *= Dimension;
2809 for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i);
2811 for ( i = 1; i <= NewKnots.Length(); i++) {
2812 NewKnots(i) = TempKnots( Kindex+i-1);
2813 NewMults(i) = TempMults( Kindex+i-1);
2815 NewMults(1) = Min(Degree, NewMults(1)) + 1 ;
2816 NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ;
2819 //=======================================================================
2820 //function : Solves a LU factored Matrix
2822 //=======================================================================
2825 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2826 const Standard_Integer UpperBandWidth,
2827 const Standard_Integer LowerBandWidth,
2828 const Standard_Integer ArrayDimension,
2829 Standard_Real& Array)
2831 Standard_Integer ii,
2838 Standard_Real *PolesArray = &Array ;
2839 Standard_Real Inverse ;
2842 if (Matrix.LowerCol() != 1 ||
2843 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2848 for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow() ; ii++) {
2849 MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ?
2850 ii - LowerBandWidth : Matrix.LowerRow()) ;
2852 for ( jj = MinIndex ; jj < ii ; jj++) {
2854 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2855 PolesArray[(ii-1) * ArrayDimension + kk] +=
2856 PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2861 for (ii = Matrix.UpperRow() ; ii >= Matrix.LowerRow() ; ii--) {
2862 MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ?
2863 ii + UpperBandWidth : Matrix.UpperRow()) ;
2865 for (jj = MaxIndex ; jj > ii ; jj--) {
2867 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2868 PolesArray[(ii-1) * ArrayDimension + kk] -=
2869 PolesArray[(jj - 1) * ArrayDimension + kk] *
2870 Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2874 //fixing a bug PRO18577 to avoid divizion by zero
2876 Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ;
2877 Standard_Real Toler = 1.0e-16;
2878 if ( Abs(divizor) > Toler )
2879 Inverse = 1.0e0 / divizor ;
2882 // cout << " BSplCLib::SolveBandedSystem() : zero determinant " << endl;
2887 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2888 PolesArray[(ii-1) * ArrayDimension + kk] *= Inverse ;
2892 return (ReturnCode) ;
2895 //=======================================================================
2896 //function : Solves a LU factored Matrix
2898 //=======================================================================
2901 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
2902 const Standard_Integer UpperBandWidth,
2903 const Standard_Integer LowerBandWidth,
2904 const Standard_Boolean HomogeneousFlag,
2905 const Standard_Integer ArrayDimension,
2906 Standard_Real& Poles,
2907 Standard_Real& Weights)
2909 Standard_Integer ii,
2914 Standard_Real Inverse,
2915 *PolesArray = &Poles,
2916 *WeightsArray = &Weights ;
2918 if (Matrix.LowerCol() != 1 ||
2919 Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2923 if (HomogeneousFlag == Standard_False) {
2925 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) {
2927 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2928 PolesArray[ii * ArrayDimension + kk] *=
2934 BSplCLib::SolveBandedSystem(Matrix,
2939 if (ErrorCode != 0) {
2944 BSplCLib::SolveBandedSystem(Matrix,
2949 if (ErrorCode != 0) {
2953 if (HomogeneousFlag == Standard_False) {
2955 for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) {
2956 Inverse = 1.0e0 / WeightsArray[ii] ;
2958 for (kk = 0 ; kk < ArrayDimension ; kk++) {
2959 PolesArray[ii * ArrayDimension + kk] *= Inverse ;
2963 FINISH : return (ReturnCode) ;
2966 //=======================================================================
2967 //function : BuildSchoenbergPoints
2969 //=======================================================================
2971 void BSplCLib::BuildSchoenbergPoints(const Standard_Integer Degree,
2972 const TColStd_Array1OfReal& FlatKnots,
2973 TColStd_Array1OfReal& Parameters)
2975 Standard_Integer ii,
2977 Standard_Real Inverse ;
2978 Inverse = 1.0e0 / (Standard_Real)Degree ;
2980 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
2981 Parameters(ii) = 0.0e0 ;
2983 for (jj = 1 ; jj <= Degree ; jj++) {
2984 Parameters(ii) += FlatKnots(jj + ii) ;
2986 Parameters(ii) *= Inverse ;
2990 //=======================================================================
2991 //function : Interpolate
2993 //=======================================================================
2995 void BSplCLib::Interpolate(const Standard_Integer Degree,
2996 const TColStd_Array1OfReal& FlatKnots,
2997 const TColStd_Array1OfReal& Parameters,
2998 const TColStd_Array1OfInteger& ContactOrderArray,
2999 const Standard_Integer ArrayDimension,
3000 Standard_Real& Poles,
3001 Standard_Integer& InversionProblem)
3003 Standard_Integer ErrorCode,
3006 // Standard_Real *PolesArray = &Poles ;
3007 math_Matrix InterpolationMatrix(1, Parameters.Length(),
3008 1, 2 * Degree + 1) ;
3010 BSplCLib::BuildBSpMatrix(Parameters,
3014 InterpolationMatrix,
3017 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3020 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3024 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3027 BSplCLib::SolveBandedSystem(InterpolationMatrix,
3033 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
3036 //=======================================================================
3037 //function : Interpolate
3039 //=======================================================================
3041 void BSplCLib::Interpolate(const Standard_Integer Degree,
3042 const TColStd_Array1OfReal& FlatKnots,
3043 const TColStd_Array1OfReal& Parameters,
3044 const TColStd_Array1OfInteger& ContactOrderArray,
3045 const Standard_Integer ArrayDimension,
3046 Standard_Real& Poles,
3047 Standard_Real& Weights,
3048 Standard_Integer& InversionProblem)
3050 Standard_Integer ErrorCode,
3054 math_Matrix InterpolationMatrix(1, Parameters.Length(),
3055 1, 2 * Degree + 1) ;
3057 BSplCLib::BuildBSpMatrix(Parameters,
3061 InterpolationMatrix,
3064 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3067 BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3071 Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3074 BSplCLib::SolveBandedSystem(InterpolationMatrix,
3082 Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ;
3085 //=======================================================================
3086 //function : Evaluates a Bspline function : uses the ExtrapMode
3087 //purpose : the function is extrapolated using the Taylor expansion
3088 // of degree ExtrapMode[0] to the left and the Taylor
3089 // expansion of degree ExtrapMode[1] to the right
3090 // this evaluates the numerator by multiplying by the weights
3091 // and evaluating it but does not call RationalDerivatives after
3092 //=======================================================================
3095 (const Standard_Real Parameter,
3096 const Standard_Boolean PeriodicFlag,
3097 const Standard_Integer DerivativeRequest,
3098 Standard_Integer& ExtrapMode,
3099 const Standard_Integer Degree,
3100 const TColStd_Array1OfReal& FlatKnots,
3101 const Standard_Integer ArrayDimension,
3102 Standard_Real& Poles,
3103 Standard_Real& Weights,
3104 Standard_Real& PolesResults,
3105 Standard_Real& WeightsResults)
3107 Standard_Integer ii,
3116 ExtrapolatingFlag[2],
3120 FirstNonZeroBsplineIndex,
3121 LocalRequest = DerivativeRequest ;
3122 Standard_Real *PResultArray,
3130 PolesArray = &Poles ;
3131 WeightsArray = &Weights ;
3132 ExtrapModeArray = &ExtrapMode ;
3133 PResultArray = &PolesResults ;
3134 WResultArray = &WeightsResults ;
3135 LocalParameter = Parameter ;
3136 ExtrapolatingFlag[0] =
3137 ExtrapolatingFlag[1] = 0 ;
3139 // check if we are extrapolating to a degree which is smaller than
3140 // the degree of the Bspline
3143 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3145 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3146 LocalParameter -= Period ;
3149 while (LocalParameter < FlatKnots(2)) {
3150 LocalParameter += Period ;
3153 if (Parameter < FlatKnots(2) &&
3154 LocalRequest < ExtrapModeArray[0] &&
3155 ExtrapModeArray[0] < Degree) {
3156 LocalRequest = ExtrapModeArray[0] ;
3157 LocalParameter = FlatKnots(2) ;
3158 ExtrapolatingFlag[0] = 1 ;
3160 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3161 LocalRequest < ExtrapModeArray[1] &&
3162 ExtrapModeArray[1] < Degree) {
3163 LocalRequest = ExtrapModeArray[1] ;
3164 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3165 ExtrapolatingFlag[1] = 1 ;
3167 Delta = Parameter - LocalParameter ;
3168 if (LocalRequest >= Order) {
3169 LocalRequest = Degree ;
3172 Modulus = FlatKnots.Length() - Degree -1 ;
3175 Modulus = FlatKnots.Length() - Degree ;
3178 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3180 BSplCLib::EvalBsplineBasis(1,
3185 FirstNonZeroBsplineIndex,
3187 if (ErrorCode != 0) {
3191 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3195 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3196 Index1 = FirstNonZeroBsplineIndex ;
3198 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3199 PResultArray[Index + kk] = 0.0e0 ;
3201 WResultArray[Index] = 0.0e0 ;
3203 for (jj = 1 ; jj <= Order ; jj++) {
3205 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3206 PResultArray[Index + kk] +=
3207 PolesArray[(Index1-1) * ArrayDimension + kk]
3208 * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3210 WResultArray[Index2] += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3212 Index1 = Index1 % Modulus ;
3215 Index += ArrayDimension ;
3221 // store Taylor expansion in LocalRealArray
3223 NewRequest = DerivativeRequest ;
3224 if (NewRequest > Degree) {
3225 NewRequest = Degree ;
3227 NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3231 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3232 Index1 = FirstNonZeroBsplineIndex ;
3234 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3235 LocalRealArray[Index + kk] = 0.0e0 ;
3238 for (jj = 1 ; jj <= Order ; jj++) {
3240 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3241 LocalRealArray[Index + kk] +=
3242 PolesArray[(Index1-1)*ArrayDimension + kk] *
3243 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3245 Index1 = Index1 % Modulus ;
3249 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3250 LocalRealArray[Index + kk] *= Inverse ;
3252 Index += ArrayDimension ;
3253 Inverse /= (Standard_Real) ii ;
3255 PLib::EvalPolynomial(Delta,
3264 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3265 Index1 = FirstNonZeroBsplineIndex ;
3266 LocalRealArray[Index] = 0.0e0 ;
3268 for (jj = 1 ; jj <= Order ; jj++) {
3269 LocalRealArray[Index] +=
3270 WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3271 Index1 = Index1 % Modulus ;
3274 LocalRealArray[Index + kk] *= Inverse ;
3276 Inverse /= (Standard_Real) ii ;
3278 PLib::EvalPolynomial(Delta,
3288 //=======================================================================
3289 //function : Evaluates a Bspline function : uses the ExtrapMode
3290 //purpose : the function is extrapolated using the Taylor expansion
3291 // of degree ExtrapMode[0] to the left and the Taylor
3292 // expansion of degree ExtrapMode[1] to the right
3293 // WARNING : the array Results is supposed to have at least
3294 // (DerivativeRequest + 1) * ArrayDimension slots and the
3296 //=======================================================================
3299 (const Standard_Real Parameter,
3300 const Standard_Boolean PeriodicFlag,
3301 const Standard_Integer DerivativeRequest,
3302 Standard_Integer& ExtrapMode,
3303 const Standard_Integer Degree,
3304 const TColStd_Array1OfReal& FlatKnots,
3305 const Standard_Integer ArrayDimension,
3306 Standard_Real& Poles,
3307 Standard_Real& Results)
3309 Standard_Integer ii,
3317 ExtrapolatingFlag[2],
3321 FirstNonZeroBsplineIndex,
3322 LocalRequest = DerivativeRequest ;
3324 Standard_Real *ResultArray,
3331 PolesArray = &Poles ;
3332 ExtrapModeArray = &ExtrapMode ;
3333 ResultArray = &Results ;
3334 LocalParameter = Parameter ;
3335 ExtrapolatingFlag[0] =
3336 ExtrapolatingFlag[1] = 0 ;
3338 // check if we are extrapolating to a degree which is smaller than
3339 // the degree of the Bspline
3342 Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3344 while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3345 LocalParameter -= Period ;
3348 while (LocalParameter < FlatKnots(2)) {
3349 LocalParameter += Period ;
3352 if (Parameter < FlatKnots(2) &&
3353 LocalRequest < ExtrapModeArray[0] &&
3354 ExtrapModeArray[0] < Degree) {
3355 LocalRequest = ExtrapModeArray[0] ;
3356 LocalParameter = FlatKnots(2) ;
3357 ExtrapolatingFlag[0] = 1 ;
3359 if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3360 LocalRequest < ExtrapModeArray[1] &&
3361 ExtrapModeArray[1] < Degree) {
3362 LocalRequest = ExtrapModeArray[1] ;
3363 LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3364 ExtrapolatingFlag[1] = 1 ;
3366 Delta = Parameter - LocalParameter ;
3367 if (LocalRequest >= Order) {
3368 LocalRequest = Degree ;
3372 Modulus = FlatKnots.Length() - Degree -1 ;
3375 Modulus = FlatKnots.Length() - Degree ;
3378 BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3381 BSplCLib::EvalBsplineBasis(1,
3386 FirstNonZeroBsplineIndex,
3388 if (ErrorCode != 0) {
3392 if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3395 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3396 Index1 = FirstNonZeroBsplineIndex ;
3398 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3399 ResultArray[Index + kk] = 0.0e0 ;
3402 for (jj = 1 ; jj <= Order ; jj++) {
3404 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3405 ResultArray[Index + kk] +=
3406 PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3408 Index1 = Index1 % Modulus ;
3411 Index += ArrayDimension ;
3416 // store Taylor expansion in LocalRealArray
3418 NewRequest = DerivativeRequest ;
3419 if (NewRequest > Degree) {
3420 NewRequest = Degree ;
3422 NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3427 for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3428 Index1 = FirstNonZeroBsplineIndex ;
3430 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3431 LocalRealArray[Index + kk] = 0.0e0 ;
3434 for (jj = 1 ; jj <= Order ; jj++) {
3436 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3437 LocalRealArray[Index + kk] +=
3438 PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3440 Index1 = Index1 % Modulus ;
3444 for (kk = 0 ; kk < ArrayDimension ; kk++) {
3445 LocalRealArray[Index + kk] *= Inverse ;
3447 Index += ArrayDimension ;
3448 Inverse /= (Standard_Real) ii ;
3450 PLib::EvalPolynomial(Delta,
3460 //=======================================================================
3461 //function : TangExtendToConstraint
3462 //purpose : Extends a Bspline function using the tangency map
3466 //=======================================================================
3468 void BSplCLib::TangExtendToConstraint
3469 (const TColStd_Array1OfReal& FlatKnots,
3470 const Standard_Real C1Coefficient,
3471 const Standard_Integer NumPoles,
3472 Standard_Real& Poles,
3473 const Standard_Integer CDimension,
3474 const Standard_Integer CDegree,
3475 const TColStd_Array1OfReal& ConstraintPoint,
3476 const Standard_Integer Continuity,
3477 const Standard_Boolean After,
3478 Standard_Integer& NbPolesResult,
3479 Standard_Integer& NbKnotsResult,
3480 Standard_Real& KnotsResult,
3481 Standard_Real& PolesResult)
3484 if (CDegree<Continuity+1) {
3485 cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3488 Standard_Real * Padr = &Poles ;
3489 Standard_Real * KRadr = &KnotsResult ;
3490 Standard_Real * PRadr = &PolesResult ;
3492 ////////////////////////////////////////////////////////////////////////
3494 // 1. calculation of extension nD
3496 ////////////////////////////////////////////////////////////////////////
3499 Standard_Integer Csize = Continuity + 2;
3500 math_Matrix MatCoefs(1,Csize, 1,Csize);
3502 PLib::HermiteCoefficients(0, 1, // Limits
3503 Continuity, 0, // Orders of constraints
3507 PLib::HermiteCoefficients(0, 1, // Limits
3508 0, Continuity, // Orders of constraints
3513 // position at the node of connection
3514 Standard_Real Tbord ;
3516 Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3519 Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3521 Standard_Boolean periodic_flag = Standard_False ;
3522 Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3523 extrap_mode[0] = extrap_mode[1] = CDegree;
3524 TColStd_Array1OfReal EvalBS(1, CDimension * (derivative_request+1)) ;
3525 Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3526 BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3527 CDegree,FlatKnots,CDimension,Poles,*Eadr);
3529 // norm of the tangent at the node of connection
3530 math_Vector Tgte(1,CDimension);
3532 for (ipos=1;ipos<=CDimension;ipos++) {
3533 Tgte(ipos) = EvalBS(ipos+CDimension);
3535 Standard_Real L1=Tgte.Norm();
3538 // matrix of constraints
3539 math_Matrix Contraintes(1,Csize,1,CDimension);
3542 for (ipos=1;ipos<=CDimension;ipos++) {
3543 Contraintes(1,ipos) = EvalBS(ipos);
3544 Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3545 if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3546 if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3547 Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3552 for (ipos=1;ipos<=CDimension;ipos++) {
3553 Contraintes(1,ipos) = ConstraintPoint(ipos);
3554 Contraintes(2,ipos) = EvalBS(ipos);
3555 if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3556 if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3557 if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3561 // calculate the coefficients of extension
3562 Standard_Integer ii, jj, kk;
3563 TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3564 ExtraCoeffs.Init(0.);
3566 for (ii=1; ii<=Csize; ii++) {
3568 for (jj=1; jj<=Csize; jj++) {
3570 for (kk=1; kk<=CDimension; kk++) {
3571 ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3576 // calculate the poles of extension
3577 TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3578 Standard_Real * EPadr = &ExtrapPoles(1) ;
3579 PLib::CoefficientsPoles(CDimension,
3580 ExtraCoeffs, PLib::NoWeights(),
3581 ExtrapPoles, PLib::NoWeights());
3583 // calculate the nodes of extension with multiplicities
3584 TColStd_Array1OfReal ExtrapNoeuds(1,2);
3585 ExtrapNoeuds(1) = 0.;
3586 ExtrapNoeuds(2) = 1.;
3587 TColStd_Array1OfInteger ExtrapMults(1,2);
3588 ExtrapMults(1) = Csize;
3589 ExtrapMults(2) = Csize;
3591 // flat nodes of extension
3592 TColStd_Array1OfReal FK2(1, Csize*2);
3593 BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3595 // norm of the tangent at the connection point
3597 BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3598 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3601 BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3602 Csize-1,FK2,CDimension,*EPadr,*Eadr);
3605 for (ipos=1;ipos<=CDimension;ipos++) {
3606 Tgte(ipos) = EvalBS(ipos+CDimension);
3608 Standard_Real L2 = Tgte.Norm();
3610 // harmonisation of degrees
3611 TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3612 TColStd_Array1OfReal NewK2(1, 2);
3613 TColStd_Array1OfInteger NewM2(1, 2);
3614 if (Csize-1<CDegree) {
3615 BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3616 ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3620 NewP2 = ExtrapPoles;
3621 NewK2 = ExtrapNoeuds;
3622 NewM2 = ExtrapMults;
3625 // flat nodes of extension after harmonization of degrees
3626 TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3627 BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3630 ////////////////////////////////////////////////////////////////////////
3632 // 2. concatenation C0
3634 ////////////////////////////////////////////////////////////////////////
3636 // ratio of reparametrization
3637 Standard_Real Ratio=1, Delta;
3638 if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3641 if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3644 // do not touch the first BSpline
3645 Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3648 // do not touch the second BSpline
3649 Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower());
3652 // result of the concatenation
3653 Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1;
3654 Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1);
3655 TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension);
3656 TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2);
3659 Standard_Integer indNP, indP, indEP;
3662 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3664 for (jj=1; jj<=CDimension; jj++) {
3665 indNP = (ii-1)*CDimension+jj;
3666 indP = (ii-1)*CDimension+jj-1;
3667 indEP = (ii-NbP1)*CDimension+jj;
3668 if (ii<NbP1) NewPoles(indNP) = Padr[indP];
3669 else NewPoles(indNP) = NewP2(indEP);
3675 for (ii=1; ii<=NbP1+NbP2-1; ii++) {
3677 for (jj=1; jj<=CDimension; jj++) {
3678 indNP = (ii-1)*CDimension+jj;
3679 indEP = (ii-1)*CDimension+jj;
3680 indP = (ii-NbP2)*CDimension+jj-1;
3681 if (ii<NbP2) NewPoles(indNP) = NewP2(indEP);
3682 else NewPoles(indNP) = Padr[indP];
3689 // start with the nodes of the initial surface
3691 for (ii=1; ii<NbK1; ii++) {
3692 NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1);
3694 // continue with the reparameterized nodes of the extension
3696 for (ii=1; ii<=NbK2-CDegree-1; ii++) {
3697 NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta;
3701 // start with the reparameterized nodes of the extension
3703 for (ii=1; ii<NbK2-CDegree; ii++) {
3704 NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta;
3706 // continue with the nodes of the initial surface
3708 for (ii=2; ii<=NbK1; ii++) {
3709 NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1);
3714 ////////////////////////////////////////////////////////////////////////
3716 // 3. reduction of multiplicite at the node of connection
3718 ////////////////////////////////////////////////////////////////////////
3720 // number of separate nodes
3721 Standard_Integer KLength = 1;
3723 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3724 if (NewFlats(ii) != NewFlats(ii-1)) KLength++;
3727 // flat nodes --> nodes + multiplicities
3728 TColStd_Array1OfReal NewKnots (1, KLength);
3729 TColStd_Array1OfInteger NewMults (1, KLength);
3732 NewKnots(jj) = NewFlats(1);
3734 for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3735 if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++;
3738 NewKnots(jj) = NewFlats(ii);
3742 // reduction of multiplicity at the second or the last but one node
3743 Standard_Integer Index = 2, M = CDegree;
3744 if (After) Index = KLength-1;
3745 TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension);
3746 TColStd_Array1OfReal ResultKnots (1, KLength);
3747 TColStd_Array1OfInteger ResultMults (1, KLength);
3748 Standard_Real Tol = 1.e-6;
3749 Standard_Boolean Ok = Standard_True;
3751 while ( (M>CDegree-Continuity) && Ok) {
3752 Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension,
3753 NewPoles, NewKnots, NewMults,
3754 ResultPoles, ResultKnots, ResultMults, Tol);
3759 // number of poles of the concatenation
3760 NbPolesResult = NbP1 + NbP2 - 1;
3761 // the poles of the concatenation
3762 Standard_Integer PLength = NbPolesResult*CDimension;
3764 for (jj=1; jj<=PLength; jj++) {
3765 PRadr[jj-1] = NewPoles(jj);
3768 // flat nodes of the concatenation
3769 Standard_Integer ideb = 0;
3771 for (jj=0; jj<NewKnots.Length(); jj++) {
3772 for (ii=0; ii<NewMults(jj+1); ii++) {
3773 KRadr[ideb+ii] = NewKnots(jj+1);
3775 ideb += NewMults(jj+1);
3777 NbKnotsResult = ideb;
3781 // number of poles of the result
3782 NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M;
3783 // the poles of the result
3784 Standard_Integer PLength = NbPolesResult*CDimension;
3786 for (jj=0; jj<PLength; jj++) {
3787 PRadr[jj] = ResultPoles(jj+1);
3790 // flat nodes of the result
3791 Standard_Integer ideb = 0;
3793 for (jj=0; jj<ResultKnots.Length(); jj++) {
3794 for (ii=0; ii<ResultMults(jj+1); ii++) {
3795 KRadr[ideb+ii] = ResultKnots(jj+1);
3797 ideb += ResultMults(jj+1);
3799 NbKnotsResult = ideb;
3803 //=======================================================================
3804 //function : Resolution
3807 // Let C(t) = SUM Ci Bi(t) a Bspline curve of degree d
3809 // with nodes tj for j = 1,n+d+1
3813 // Then C (t) = SUM d * --------- Bi (t)
3814 // i = 2,n ti+d - ti
3817 // for the base of BSpline Bi (t) of degree d-1.
3819 // Consequently the upper bound of the norm of the derivative from C is :
3823 // d * Max | --------- |
3824 // i = 2,n | ti+d - ti |
3827 // In the rational case set C(t) = -----
3831 // D(t) = SUM Di Bi(t)
3834 // N(t) = SUM Di * Ci Bi(t)
3837 // N'(t) - D'(t) C(t)
3838 // C'(t) = -----------------------
3842 // N'(t) - D'(t) C(t) =
3844 // Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t)) d-1
3845 // SUM d * ---------------------------------------- * Bi (t) =
3849 // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) d-1
3850 // SUM SUM d * ----------------------------------- * Betaj(t) * Bi (t)
3851 //i=2,n j=1,n ti+d - ti
3856 // Betaj(t) = --------
3859 // Betaj(t) form a partition >= 0 of the entity with support
3860 // tj, tj+d+1. Consequently if Rj = {j-d, ...., j+d+d+1}
3861 // obtain an upper bound of the derivative of C by taking :
3868 // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)
3869 // Max Max d * -----------------------------------
3870 // j=1,n i dans Rj ti+d - ti
3872 // --------------------------------------------------------
3878 //=======================================================================
3880 void BSplCLib::Resolution( Standard_Real& Poles,
3881 const Standard_Integer ArrayDimension,
3882 const Standard_Integer NumPoles,
3883 const TColStd_Array1OfReal& Weights,
3884 const TColStd_Array1OfReal& FlatKnots,
3885 const Standard_Integer Degree,
3886 const Standard_Real Tolerance3D,
3887 Standard_Real& UTolerance)
3889 Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim;
3890 Standard_Integer lower,upper,ii_minus,jj,ii_miDim;
3891 Standard_Integer Deg1 = Degree + 1;
3892 Standard_Integer Deg2 = (Degree << 1) + 1;
3893 Standard_Real value,factor,W,min_weights,inverse;
3894 Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3;
3895 Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3;
3896 Standard_Real wg_ii_index, wg_ii_minus;
3897 Standard_Real *PA,max_derivative;
3898 const Standard_Real * FK = &FlatKnots(FlatKnots.Lower());
3900 max_derivative = 0.0e0;
3901 num_poles = FlatKnots.Length() - Deg1;
3902 switch (ArrayDimension) {
3904 if (&Weights != NULL) {
3905 const Standard_Real * WG = &Weights(Weights.Lower());
3906 min_weights = WG[0];
3908 for (ii = 1 ; ii < NumPoles ; ii++) {
3910 if (W < min_weights) min_weights = W;
3913 for (ii = 1 ; ii < num_poles ; ii++) {
3914 ii_index = ii % NumPoles;
3915 ii_inDim = ii_index << 1;
3916 ii_minus = (ii - 1) % NumPoles;
3917 ii_miDim = ii_minus << 1;
3918 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3919 pa_ii_inDim_1 = PA[ii_inDim];
3920 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3921 pa_ii_miDim_1 = PA[ii_miDim];
3922 wg_ii_index = WG[ii_index];
3923 wg_ii_minus = WG[ii_minus];
3924 inverse = FK[ii + Degree] - FK[ii];
3925 inverse = 1.0e0 / inverse;
3927 if (lower < 0) lower = 0;
3929 if (upper > num_poles) upper = num_poles;
3931 for (jj = lower ; jj < upper ; jj++) {
3932 jj_index = jj % NumPoles;
3933 jj_index = jj_index << 1;
3935 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3936 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3937 if (factor < 0) factor = - factor;
3938 value += factor; jj_index++;
3939 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3940 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3941 if (factor < 0) factor = - factor;
3944 if (max_derivative < value) max_derivative = value;
3947 max_derivative /= min_weights;
3951 for (ii = 1 ; ii < num_poles ; ii++) {
3952 ii_index = ii % NumPoles;
3953 ii_index = ii_index << 1;
3954 ii_minus = (ii - 1) % NumPoles;
3955 ii_minus = ii_minus << 1;
3956 inverse = FK[ii + Degree] - FK[ii];
3957 inverse = 1.0e0 / inverse;
3959 factor = PA[ii_index] - PA[ii_minus];
3960 if (factor < 0) factor = - factor;
3961 value += factor; ii_index++; ii_minus++;
3962 factor = PA[ii_index] - PA[ii_minus];
3963 if (factor < 0) factor = - factor;
3966 if (max_derivative < value) max_derivative = value;
3972 if (&Weights != NULL) {
3973 const Standard_Real * WG = &Weights(Weights.Lower());
3974 min_weights = WG[0];
3976 for (ii = 1 ; ii < NumPoles ; ii++) {
3978 if (W < min_weights) min_weights = W;
3981 for (ii = 1 ; ii < num_poles ; ii++) {
3982 ii_index = ii % NumPoles;
3983 ii_inDim = (ii_index << 1) + ii_index;
3984 ii_minus = (ii - 1) % NumPoles;
3985 ii_miDim = (ii_minus << 1) + ii_minus;
3986 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3987 pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
3988 pa_ii_inDim_2 = PA[ii_inDim];
3989 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3990 pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
3991 pa_ii_miDim_2 = PA[ii_miDim];
3992 wg_ii_index = WG[ii_index];
3993 wg_ii_minus = WG[ii_minus];
3994 inverse = FK[ii + Degree] - FK[ii];
3995 inverse = 1.0e0 / inverse;
3997 if (lower < 0) lower = 0;
3999 if (upper > num_poles) upper = num_poles;
4001 for (jj = lower ; jj < upper ; jj++) {
4002 jj_index = jj % NumPoles;
4003 jj_index = (jj_index << 1) + jj_index;
4005 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4006 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4007 if (factor < 0) factor = - factor;
4008 value += factor; jj_index++;
4009 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4010 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4011 if (factor < 0) factor = - factor;
4012 value += factor; jj_index++;
4013 factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4014 ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4015 if (factor < 0) factor = - factor;
4018 if (max_derivative < value) max_derivative = value;
4021 max_derivative /= min_weights;
4025 for (ii = 1 ; ii < num_poles ; ii++) {
4026 ii_index = ii % NumPoles;
4027 ii_index = (ii_index << 1) + ii_index;
4028 ii_minus = (ii - 1) % NumPoles;
4029 ii_minus = (ii_minus << 1) + ii_minus;
4030 inverse = FK[ii + Degree] - FK[ii];
4031 inverse = 1.0e0 / inverse;
4033 factor = PA[ii_index] - PA[ii_minus];
4034 if (factor < 0) factor = - factor;
4035 value += factor; ii_index++; ii_minus++;
4036 factor = PA[ii_index] - PA[ii_minus];
4037 if (factor < 0) factor = - factor;
4038 value += factor; ii_index++; ii_minus++;
4039 factor = PA[ii_index] - PA[ii_minus];
4040 if (factor < 0) factor = - factor;
4043 if (max_derivative < value) max_derivative = value;
4049 if (&Weights != NULL) {
4050 const Standard_Real * WG = &Weights(Weights.Lower());
4051 min_weights = WG[0];
4053 for (ii = 1 ; ii < NumPoles ; ii++) {
4055 if (W < min_weights) min_weights = W;
4058 for (ii = 1 ; ii < num_poles ; ii++) {
4059 ii_index = ii % NumPoles;
4060 ii_inDim = ii_index << 2;
4061 ii_minus = (ii - 1) % NumPoles;
4062 ii_miDim = ii_minus << 2;
4063 pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
4064 pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
4065 pa_ii_inDim_2 = PA[ii_inDim]; ii_inDim++;
4066 pa_ii_inDim_3 = PA[ii_inDim];
4067 pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
4068 pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
4069 pa_ii_miDim_2 = PA[ii_miDim]; ii_miDim++;
4070 pa_ii_miDim_3 = PA[ii_miDim];
4071 wg_ii_index = WG[ii_index];
4072 wg_ii_minus = WG[ii_minus];
4073 inverse = FK[ii + Degree] - FK[ii];
4074 inverse = 1.0e0 / inverse;
4076 if (lower < 0) lower = 0;
4078 if (upper > num_poles) upper = num_poles;
4080 for (jj = lower ; jj < upper ; jj++) {
4081 jj_index = jj % NumPoles;
4082 jj_index = jj_index << 2;
4084 factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4085 ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4086 if (factor < 0) factor = - factor;
4087 value += factor; jj_index++;
4088 factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4089 ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4090 if (factor < 0) factor = - factor;
4091 value += factor; jj_index++;
4092 factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4093 ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4094 if (factor < 0) factor = - factor;
4095 value += factor; jj_index++;
4096 factor = (((PA[jj_index] - pa_ii_inDim_3) * wg_ii_index) -
4097 ((PA[jj_index] - pa_ii_miDim_3) * wg_ii_minus));
4098 if (factor < 0) factor = - factor;
4101 if (max_derivative < value) max_derivative = value;
4104 max_derivative /= min_weights;
4108 for (ii = 1 ; ii < num_poles ; ii++) {
4109 ii_index = ii % NumPoles;
4110 ii_index = ii_index << 2;
4111 ii_minus = (ii - 1) % NumPoles;
4112 ii_minus = ii_minus << 2;
4113 inverse = FK[ii + Degree] - FK[ii];
4114 inverse = 1.0e0 / inverse;
4116 factor = PA[ii_index] - PA[ii_minus];
4117 if (factor < 0) factor = - factor;
4118 value += factor; ii_index++; ii_minus++;
4119 factor = PA[ii_index] - PA[ii_minus];
4120 if (factor < 0) factor = - factor;
4121 value += factor; ii_index++; ii_minus++;
4122 factor = PA[ii_index] - PA[ii_minus];
4123 if (factor < 0) factor = - factor;
4124 value += factor; ii_index++; ii_minus++;
4125 factor = PA[ii_index] - PA[ii_minus];
4126 if (factor < 0) factor = - factor;
4129 if (max_derivative < value) max_derivative = value;
4135 Standard_Integer kk;
4136 if (&Weights != NULL) {
4137 const Standard_Real * WG = &Weights(Weights.Lower());
4138 min_weights = WG[0];
4140 for (ii = 1 ; ii < NumPoles ; ii++) {
4142 if (W < min_weights) min_weights = W;
4145 for (ii = 1 ; ii < num_poles ; ii++) {
4146 ii_index = ii % NumPoles;
4147 ii_inDim = ii_index * ArrayDimension;
4148 ii_minus = (ii - 1) % NumPoles;
4149 ii_miDim = ii_minus * ArrayDimension;
4150 wg_ii_index = WG[ii_index];
4151 wg_ii_minus = WG[ii_minus];
4152 inverse = FK[ii + Degree] - FK[ii];
4153 inverse = 1.0e0 / inverse;
4155 if (lower < 0) lower = 0;
4157 if (upper > num_poles) upper = num_poles;
4159 for (jj = lower ; jj < upper ; jj++) {
4160 jj_index = jj % NumPoles;
4161 jj_index *= ArrayDimension;
4164 for (kk = 0 ; kk < ArrayDimension ; kk++) {
4165 factor = (((PA[jj_index + kk] - PA[ii_inDim + kk]) * wg_ii_index) -
4166 ((PA[jj_index + kk] - PA[ii_miDim + kk]) * wg_ii_minus));
4167 if (factor < 0) factor = - factor;
4171 if (max_derivative < value) max_derivative = value;
4174 max_derivative /= min_weights;
4178 for (ii = 1 ; ii < num_poles ; ii++) {
4179 ii_index = ii % NumPoles;
4180 ii_index *= ArrayDimension;
4181 ii_minus = (ii - 1) % NumPoles;
4182 ii_minus *= ArrayDimension;
4183 inverse = FK[ii + Degree] - FK[ii];
4184 inverse = 1.0e0 / inverse;
4187 for (kk = 0 ; kk < ArrayDimension ; kk++) {
4188 factor = PA[ii_index + kk] - PA[ii_minus + kk];
4189 if (factor < 0) factor = - factor;
4193 if (max_derivative < value) max_derivative = value;
4198 max_derivative *= Degree;
4199 if (max_derivative > RealSmall())
4200 UTolerance = Tolerance3D / max_derivative;
4202 UTolerance = Tolerance3D / RealSmall();
4205 //=======================================================================
4206 // function: FlatBezierKnots
4208 //=======================================================================
4210 // array of flat knots for bezier curve of maximum 25 degree
4211 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,
4212 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 };
4213 const Standard_Real& BSplCLib::FlatBezierKnots (const Standard_Integer Degree)
4215 Standard_OutOfRange_Raise_if (Degree < 1 || Degree > MaxDegree() || MaxDegree() != 25,
4216 "Bezier curve degree greater than maximal supported");
4218 return knots[25-Degree];