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