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