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