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