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