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