345e3c6c1b11aa02a7d49cf4bd35242ac8a94d5c
[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_Real au,oldau = AddKnots(ak),Eps;
1706   
1707   while (ak <= AddKnots.Upper()) {
1708     au = AddKnots(ak);
1709     if (au < oldau) return Standard_False;
1710     oldau = au;
1711
1712     Eps = Max(Tolerance,Epsilon(au));
1713     
1714     while ((k < Knots.Upper()) && (Knots(k+1)  - au <= Eps)) {
1715       k++;
1716       NbKnots++;
1717       sigma += Mults(k);
1718     }
1719
1720     if (addflat) amult = 1;
1721     else         amult = Max(0,AddMults(ak));
1722     
1723     while ((ak < AddKnots.Upper()) &&
1724            (Abs(au - AddKnots(ak+1)) <= Eps)) {
1725       ak++;
1726       if (Add) {
1727         if (addflat) amult++;
1728         else         amult += Max(0,AddMults(ak));
1729       }
1730     }
1731     
1732     
1733     if (Abs(au - Knots(k)) <= Eps) {
1734       // identic to existing knot
1735       mult = Mults(k);
1736       if (Add) {
1737         if (mult + amult > Degree)
1738           amult = Max(0,Degree - mult);
1739         sigma += amult;
1740         
1741       }
1742       else if (amult > mult) {
1743         if (amult > Degree) amult = Degree;
1744         sigma += amult - mult;
1745       }
1746       /*
1747       // on periodic curves if this is the last knot
1748       // the multiplicity is added twice to account for the first knot
1749       if (k == Knots.Upper() && Periodic) {
1750         if (Add)
1751           sigma += amult;
1752         else
1753           sigma += amult - mult;
1754       }
1755       */
1756     }
1757     else {
1758       // not identic to existing knot
1759       if (amult > 0) {
1760         if (amult > Degree) amult = Degree;
1761         NbKnots++;
1762         sigma += amult;
1763       }
1764     }
1765     
1766     ak++;
1767   }
1768   
1769   // count the last knots
1770   while (k < Knots.Upper()) {
1771     k++;
1772     NbKnots++;
1773     sigma += Mults(k);
1774   }
1775
1776   if (Periodic) {
1777     //for periodic B-Spline the requirement is that multiplicites of the first
1778     //and last knots must be equal (see Geom_BSplineCurve constructor for
1779     //instance);
1780     //respectively AddMults() must meet this requirement if AddKnots() contains
1781     //knot(s) coincident with first or last
1782     NbPoles = sigma - Mults(Knots.Upper());
1783   }
1784   else {
1785     NbPoles = sigma - Degree - 1;
1786   }
1787  
1788   return Standard_True;
1789 }
1790
1791 //=======================================================================
1792 //function : Copy
1793 //purpose  : copy reals from an array to an other
1794 //        
1795 //   NbValues are copied from OldPoles(OldFirst)
1796 //                 to    NewPoles(NewFirst)
1797 //
1798 //   Periodicity is handled.
1799 //   OldFirst and NewFirst are updated 
1800 //   to the position after the last copied pole.
1801 //
1802 //=======================================================================
1803
1804 static void Copy(const Standard_Integer      NbPoles,
1805                  Standard_Integer&           OldFirst,
1806                  const TColStd_Array1OfReal& OldPoles,
1807                  Standard_Integer&           NewFirst,
1808                  TColStd_Array1OfReal&       NewPoles)
1809 {
1810   // reset the index in the range for periodicity
1811
1812   OldFirst = OldPoles.Lower() + 
1813     (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1);
1814
1815   NewFirst = NewPoles.Lower() + 
1816     (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1);
1817
1818   // copy
1819   Standard_Integer i;
1820
1821   for (i = 1; i <= NbPoles; i++) {
1822     NewPoles(NewFirst) = OldPoles(OldFirst);
1823     OldFirst++;
1824     if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower();
1825     NewFirst++;
1826     if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower();
1827   }
1828 }
1829                       
1830 //=======================================================================
1831 //function : InsertKnots
1832 //purpose  : insert an array of knots and multiplicities
1833 //=======================================================================
1834
1835 void BSplCLib::InsertKnots
1836 (const Standard_Integer         Degree, 
1837  const Standard_Boolean         Periodic,
1838  const Standard_Integer         Dimension, 
1839  const TColStd_Array1OfReal&    Poles,  
1840  const TColStd_Array1OfReal&    Knots,    
1841  const TColStd_Array1OfInteger& Mults, 
1842  const TColStd_Array1OfReal&    AddKnots,    
1843  const TColStd_Array1OfInteger& AddMults, 
1844  TColStd_Array1OfReal&          NewPoles,
1845  TColStd_Array1OfReal&          NewKnots,    
1846  TColStd_Array1OfInteger&       NewMults, 
1847  const Standard_Real            Tolerance,
1848  const Standard_Boolean         Add)
1849 {
1850   Standard_Boolean addflat  = &AddMults == NULL;
1851   
1852   Standard_Integer i,k,mult,firstmult;
1853   Standard_Integer index,kn,curnk,curk;
1854   Standard_Integer p,np, curp, curnp, length, depth;
1855   Standard_Real u;
1856   Standard_Integer need;
1857   Standard_Real Eps;
1858
1859   // -------------------
1860   // create local arrays
1861   // -------------------
1862
1863   Standard_Real *knots = new Standard_Real[2*Degree];
1864   Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
1865   
1866   //----------------------------
1867   // loop on the knots to insert
1868   //----------------------------
1869   
1870   curk   = Knots.Lower()-1;          // current position in Knots
1871   curnk  = NewKnots.Lower()-1;       // current position in NewKnots
1872   curp   = Poles.Lower();            // current position in Poles
1873   curnp  = NewPoles.Lower();         // current position in NewPoles
1874
1875   // NewKnots, NewMults, NewPoles contains the current state of the curve
1876
1877   // index is the first pole of the current curve for insertion schema
1878
1879   if (Periodic) index = -Mults(Mults.Lower());
1880   else          index = -Degree-1;
1881
1882   // on Periodic curves the first knot and the last knot are inserted later
1883   // (they are the same knot)
1884   firstmult = 0;  // multiplicity of the first-last knot for periodic
1885   
1886
1887   // kn current knot to insert in AddKnots
1888
1889   for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) {
1890     
1891     u = AddKnots(kn);
1892     Eps = Max(Tolerance,Epsilon(u));
1893     
1894     //-----------------------------------
1895     // find the position in the old knots
1896     // and copy to the new knots
1897     //-----------------------------------
1898     
1899     while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) {
1900       curk++; curnk++;
1901       NewKnots(curnk) = Knots(curk);
1902       index += NewMults(curnk) = Mults(curk);
1903     }
1904     
1905     //-----------------------------------
1906     // Slice the knots and the mults
1907     // to the current size of the new curve
1908     //-----------------------------------
1909
1910     i = curnk + Knots.Upper() - curk;
1911     TColStd_Array1OfReal    nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i);
1912     TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i);
1913
1914     //-----------------------------------
1915     // copy enough knots 
1916     // to compute the insertion schema
1917     //-----------------------------------
1918
1919     k = curk;
1920     i = curnk;
1921     mult = 0;
1922
1923     while (mult < Degree && k < Knots.Upper()) {
1924       k++; i++;
1925       nknots(i) = Knots(k);
1926       mult += nmults(i) = Mults(k);
1927     }
1928
1929     // copy knots at the end for periodic curve
1930     if (Periodic) {
1931       mult = 0;
1932       k = Knots.Upper();
1933       i = nknots.Upper();
1934
1935       while (mult < Degree && i > curnk) {
1936         nknots(i) = Knots(k);
1937         mult += nmults(i) = Mults(k);
1938         k--;
1939         i--;
1940       }
1941       nmults(nmults.Upper()) = nmults(nmults.Lower());
1942     }
1943
1944   
1945
1946     //------------------------------------
1947     // do the boor scheme on the new curve
1948     // to insert the new knot
1949     //------------------------------------
1950     
1951     Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps);
1952     
1953     if (sameknot) length = Max(0,Degree - NewMults(curnk));
1954     else          length = Degree;
1955     
1956     if (addflat) depth = 1;
1957     else         depth = Min(Degree,AddMults(kn));
1958
1959     if (sameknot) {
1960       if (Add) {
1961         if ((NewMults(curnk) + depth) > Degree)
1962           depth = Degree - NewMults(curnk);
1963       }
1964       else {
1965         depth = Max(0,depth-NewMults(curnk));
1966       }
1967
1968       if (Periodic) {
1969         // on periodic curve the first and last knot are delayed to the end
1970         if (curk == Knots.Lower() || (curk == Knots.Upper())) {
1971           firstmult += depth;
1972           depth = 0;
1973         }
1974       }
1975     }
1976     if (depth <= 0) continue;
1977     
1978     BuildKnots(Degree,curnk,Periodic,nknots,nmults,*knots);
1979
1980     // copy the poles
1981
1982     need   = NewPoles.Lower()+(index+length+1)*Dimension - curnp;
1983     need = Min(need,Poles.Upper() - curp + 1);
1984
1985     p = curp;
1986     np = curnp;
1987     Copy(need,p,Poles,np,NewPoles);
1988     curp  += need;
1989     curnp += need;
1990
1991     // slice the poles to the current number of poles in case of periodic
1992     TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
1993
1994     BuildBoor(index,length,Dimension,npoles,*poles);
1995     BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length);
1996     
1997     //-------------------
1998     // copy the new poles
1999     //-------------------
2000
2001     curnp += depth * Dimension; // number of poles is increased by depth
2002     TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1);
2003     np = NewKnots.Lower()+(index+1)*Dimension;
2004
2005     for (i = 1; i <= length + depth; i++)
2006       GetPole(i,length,depth,Dimension,*poles,np,ThePoles);
2007     
2008     //-------------------
2009     // insert the knot
2010     //-------------------
2011
2012     index += depth;
2013     if (sameknot) {
2014       NewMults(curnk) += depth;
2015     }
2016     else {
2017       curnk++;
2018       NewKnots(curnk) = u;
2019       NewMults(curnk) = depth;
2020     }
2021   }
2022   
2023   //------------------------------
2024   // copy the last poles and knots
2025   //------------------------------
2026   
2027   Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles);
2028   
2029   while (curk < Knots.Upper()) {
2030     curk++;  curnk++;
2031     NewKnots(curnk) = Knots(curk);
2032     NewMults(curnk) = Mults(curk);
2033   }
2034   
2035   //------------------------------
2036   // process the first-last knot 
2037   // on periodic curves
2038   //------------------------------
2039
2040   if (firstmult > 0) {
2041     curnk = NewKnots.Lower();
2042     if (NewMults(curnk) + firstmult > Degree) {
2043       firstmult = Degree - NewMults(curnk);
2044     }
2045     if (firstmult > 0) {
2046
2047       length = Degree - NewMults(curnk);
2048       depth  = firstmult;
2049
2050       BuildKnots(Degree,curnk,Periodic,NewKnots,NewMults,*knots);
2051       TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),
2052                                   NewPoles.Lower(),
2053                                   NewPoles.Upper()-depth*Dimension);
2054       BuildBoor(0,length,Dimension,npoles,*poles);
2055       BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length);
2056       
2057       //---------------------------
2058       // copy the new poles
2059       // but rotate them with depth
2060       //---------------------------
2061       
2062       np = NewPoles.Lower();
2063
2064       for (i = depth; i < length + depth; i++)
2065         GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2066
2067       np = NewPoles.Upper() - depth*Dimension + 1;
2068
2069       for (i = 0; i < depth; i++)
2070         GetPole(i,length,depth,Dimension,*poles,np,NewPoles);
2071       
2072       NewMults(NewMults.Lower()) += depth;
2073       NewMults(NewMults.Upper()) += depth;
2074     }
2075   }
2076   // free local arrays
2077   delete [] knots;
2078   delete [] poles;
2079 }
2080
2081 //=======================================================================
2082 //function : RemoveKnot
2083 //purpose  : 
2084 //=======================================================================
2085
2086 Standard_Boolean BSplCLib::RemoveKnot 
2087 (const Standard_Integer         Index,       
2088  const Standard_Integer         Mult,        
2089  const Standard_Integer         Degree,  
2090  const Standard_Boolean         Periodic,
2091  const Standard_Integer         Dimension,  
2092  const TColStd_Array1OfReal&    Poles,
2093  const TColStd_Array1OfReal&    Knots,  
2094  const TColStd_Array1OfInteger& Mults,
2095  TColStd_Array1OfReal&          NewPoles,
2096  TColStd_Array1OfReal&          NewKnots,  
2097  TColStd_Array1OfInteger&       NewMults,
2098  const Standard_Real            Tolerance)
2099 {
2100   Standard_Integer index,i,j,k,p,np;
2101
2102   Standard_Integer TheIndex = Index;
2103
2104   // protection
2105   Standard_Integer first,last;
2106   if (Periodic) {
2107     first = Knots.Lower();
2108     last  = Knots.Upper();
2109   }
2110   else {
2111     first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1;
2112     last  = BSplCLib::LastUKnotIndex(Degree,Mults) - 1;
2113   }
2114   if (Index < first) return Standard_False;
2115   if (Index > last)  return Standard_False;
2116
2117   if ( Periodic && (Index == first)) TheIndex = last;
2118
2119   Standard_Integer depth  = Mults(TheIndex) - Mult;
2120   Standard_Integer length = Degree - Mult;
2121
2122   // -------------------
2123   // create local arrays
2124   // -------------------
2125
2126   Standard_Real *knots = new Standard_Real[4*Degree];
2127   Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension];
2128   
2129
2130   // ------------------------------------
2131   // build the knots for anti Boor Scheme
2132   // ------------------------------------
2133
2134   // the new sequence of knots
2135   // is obtained from the knots at Index-1 and Index
2136   
2137   BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,Mults,*knots);
2138   index = PoleIndex(Degree,TheIndex-1,Periodic,Mults);
2139   BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,Mults,knots[2*Degree]);
2140
2141   index += Mult;
2142
2143   for (i = 0; i < Degree - Mult; i++)
2144     knots[i] = knots[i+Mult];
2145
2146   for (i = Degree-Mult; i < 2*Degree; i++)
2147     knots[i] = knots[2*Degree+i];
2148
2149
2150   // ------------------------------------
2151   // build the poles for anti Boor Scheme
2152   // ------------------------------------
2153
2154   p = Poles.Lower()+index * Dimension;
2155
2156   for (i = 0; i <= length + depth; i++) {
2157     j = Dimension * BoorIndex(i,length,depth);
2158
2159     for (k = 0; k < Dimension; k++) {
2160       poles[j+k] = Poles(p+k);
2161     }
2162     p += Dimension;
2163     if (p > Poles.Upper()) p = Poles.Lower();
2164   }
2165
2166
2167   // ----------------
2168   // Anti Boor Scheme
2169   // ----------------
2170
2171   Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots,
2172                                            Dimension,*poles,
2173                                            depth,length,Tolerance);
2174   
2175   // ----------------
2176   // copy the results
2177   // ----------------
2178
2179   if (result) {
2180
2181     // poles
2182
2183     p = Poles.Lower();
2184     np = NewPoles.Lower();
2185     
2186     // unmodified poles before
2187     Copy((index+1)*Dimension,p,Poles,np,NewPoles);
2188     
2189     
2190     // modified
2191
2192     for (i = 1; i <= length; i++)
2193       BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles);
2194     p += (length + depth) * Dimension ;
2195     
2196     // unmodified poles after
2197     if (p != Poles.Lower()) {
2198       i = Poles.Upper() - p + 1;
2199       Copy(i,p,Poles,np,NewPoles);
2200     }
2201
2202     // knots and mults
2203
2204     if (Mult > 0) {
2205       NewKnots = Knots;
2206       NewMults = Mults;
2207       NewMults(TheIndex) = Mult;
2208       if (Periodic) {
2209         if (TheIndex == first) NewMults(last)  = Mult;
2210         if (TheIndex == last)  NewMults(first) = Mult;
2211       }
2212     }
2213     else {
2214       if (!Periodic || (TheIndex != first && TheIndex != last)) {
2215
2216         for (i = Knots.Lower(); i < TheIndex; i++) {
2217           NewKnots(i) = Knots(i);
2218           NewMults(i) = Mults(i);
2219         }
2220
2221         for (i = TheIndex+1; i <= Knots.Upper(); i++) {
2222           NewKnots(i-1) = Knots(i);
2223           NewMults(i-1) = Mults(i);
2224         }
2225       }
2226       else {
2227         // The interesting case of a Periodic curve 
2228         // where the first and last knot is removed.
2229         
2230         for (i = first; i < last-1; i++) {
2231           NewKnots(i) = Knots(i+1);
2232           NewMults(i) = Mults(i+1);
2233         }
2234         NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first);
2235         NewMults(last-1) = NewMults(first);
2236       }
2237     }
2238   }
2239
2240
2241   // free local arrays
2242   delete [] knots;
2243   delete [] poles;
2244   
2245   return result;
2246 }
2247
2248 //=======================================================================
2249 //function : IncreaseDegreeCountKnots
2250 //purpose  : 
2251 //=======================================================================
2252
2253 Standard_Integer  BSplCLib::IncreaseDegreeCountKnots
2254 (const Standard_Integer Degree,
2255  const Standard_Integer NewDegree, 
2256  const Standard_Boolean Periodic, 
2257  const TColStd_Array1OfInteger& Mults)
2258 {
2259   if (Periodic) return Mults.Length();
2260   Standard_Integer f = FirstUKnotIndex(Degree,Mults);
2261   Standard_Integer l = LastUKnotIndex(Degree,Mults);
2262   Standard_Integer m,i,removed = 0, step = NewDegree - Degree;
2263   
2264   i = Mults.Lower();
2265   m = Degree + (f - i + 1) * step + 1;
2266
2267   while (m > NewDegree+1) {
2268     removed++;
2269     m -= Mults(i) + step;
2270     i++;
2271   }
2272   if (m < NewDegree+1) removed--;
2273
2274   i = Mults.Upper();
2275   m = Degree + (i - l + 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   return Mults.Length() - removed;
2285 }
2286
2287 //=======================================================================
2288 //function : IncreaseDegree
2289 //purpose  : 
2290 //=======================================================================
2291
2292 void BSplCLib::IncreaseDegree 
2293 (const Standard_Integer         Degree,
2294  const Standard_Integer         NewDegree,
2295  const Standard_Boolean         Periodic,
2296  const Standard_Integer         Dimension,
2297  const TColStd_Array1OfReal&    Poles,
2298  const TColStd_Array1OfReal&    Knots,
2299  const TColStd_Array1OfInteger& Mults,
2300  TColStd_Array1OfReal&          NewPoles,
2301  TColStd_Array1OfReal&          NewKnots,
2302  TColStd_Array1OfInteger&       NewMults)
2303
2304   // Degree elevation of a BSpline Curve
2305
2306   // This algorithms loops on degree incrementation from Degree to NewDegree.
2307   // The variable curDeg is the current degree to increment.
2308
2309   // Before degree incrementations a "working curve" is created.
2310   // It has the same knot, poles and multiplicities.
2311
2312   // If the curve is periodic knots are added on the working curve before
2313   // and after the existing knots to make it a non-periodic curves. 
2314   // The poles are also copied.
2315
2316   // The first and last multiplicity of the working curve are set to Degree+1,
2317   // null poles are  added if necessary.
2318
2319   // Then the degree is incremented on the working curve.
2320   // The knots are unchanged but all multiplicities will be incremented.
2321
2322   // Each degree incrementation is achieved by averaging curDeg+1 curves.
2323
2324   // See : Degree elevation of B-spline curves
2325   //       Hartmut PRAUTZSCH
2326   //       CAGD 1 (1984)
2327
2328
2329   //-------------------------
2330   // create the working curve
2331   //-------------------------
2332
2333   Standard_Integer i,k,f,l,m,pf,pl,firstknot;
2334
2335   pf = 0; // number of null poles added at beginning
2336   pl = 0; // number of null poles added at end
2337
2338   Standard_Integer nbwknots = Knots.Length();
2339   f = FirstUKnotIndex(Degree,Mults);
2340   l = LastUKnotIndex (Degree,Mults);
2341
2342   if (Periodic) {
2343     // Periodic curves are transformed in non-periodic curves
2344
2345     nbwknots += f - Mults.Lower();
2346
2347     pf = -Degree - 1;
2348
2349     for (i = Mults.Lower(); i <= f; i++)
2350       pf += Mults(i);
2351
2352     nbwknots += Mults.Upper() - l;
2353
2354     pl = -Degree - 1;
2355
2356     for (i = l; i <= Mults.Upper(); i++)
2357       pl += Mults(i);
2358   }
2359
2360   // copy the knots and multiplicities 
2361   TColStd_Array1OfReal    wknots(1,nbwknots);
2362   TColStd_Array1OfInteger wmults(1,nbwknots);
2363   if (!Periodic) {
2364     wknots  = Knots;
2365     wmults  = Mults;
2366   }
2367   else {
2368     // copy the knots for a periodic curve
2369     Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2370     i = 0;
2371
2372     for (k = l; k < Knots.Upper(); k++) {
2373       i++; 
2374       wknots(i) = Knots(k) - period;
2375       wmults(i) = Mults(k);
2376     }
2377
2378     for (k = Knots.Lower(); k <= Knots.Upper(); k++) {
2379       i++; 
2380       wknots(i) = Knots(k);
2381       wmults(i) = Mults(k);
2382     }
2383
2384     for (k = Knots.Lower()+1; k <= f; k++) {
2385       i++; 
2386       wknots(i) = Knots(k) + period;
2387       wmults(i) = Mults(k);
2388     }
2389   }
2390
2391   // set the first and last mults to Degree+1
2392   // and add null poles
2393
2394   pf += Degree + 1 - wmults(1);
2395   wmults(1) = Degree + 1;
2396   pl += Degree + 1 - wmults(nbwknots);
2397   wmults(nbwknots) = Degree + 1;
2398
2399   //---------------------------
2400   // poles of the working curve
2401   //---------------------------
2402
2403   Standard_Integer nbwpoles = 0;
2404
2405   for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i);
2406   nbwpoles -= Degree + 1;
2407
2408   // we provide space for degree elevation
2409   TColStd_Array1OfReal 
2410     wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension);
2411
2412   for (i = 1; i <= pf * Dimension; i++) 
2413     wpoles(i) = 0;
2414
2415   k = Poles.Lower();
2416
2417   for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) {
2418     wpoles(i) = Poles(k);
2419     k++;
2420     if (k > Poles.Upper()) k = Poles.Lower();
2421   }
2422
2423   for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++)
2424     wpoles(i) = 0;
2425   
2426   
2427   //-----------------------------------------------------------
2428   // Declare the temporary arrays used in degree incrementation
2429   //-----------------------------------------------------------
2430
2431   Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree);
2432   // Arrays for storing the temporary curves
2433   TColStd_Array1OfReal tempc1(1,nbwp * Dimension);
2434   TColStd_Array1OfReal tempc2(1,nbwp * Dimension);
2435
2436   // Array for storing the knots to insert
2437   TColStd_Array1OfReal iknots(1,nbwknots);
2438
2439   // Arrays for receiving the knots after insertion
2440   TColStd_Array1OfReal    nknots(1,nbwknots);
2441
2442
2443   
2444   //------------------------------
2445   // Loop on degree incrementation
2446   //------------------------------
2447
2448   Standard_Integer step,curDeg;
2449   Standard_Integer nbp = nbwpoles;
2450   nbwp = nbp;
2451
2452   for (curDeg = Degree; curDeg < NewDegree; curDeg++) {
2453     
2454     nbp  = nbwp;               // current number of poles
2455     nbwp = nbp + nbwknots - 1; // new number of poles
2456
2457     // For the averaging
2458     TColStd_Array1OfReal nwpoles(1,nbwp * Dimension);
2459     nwpoles.Init(0.0e0) ;
2460   
2461     
2462     for (step = 0; step <= curDeg; step++) {
2463     
2464       // Compute the bspline of rank step.
2465
2466       // if not the first time, decrement the multiplicities back
2467       if (step != 0) {
2468         for (i = 1; i <= nbwknots; i++)
2469           wmults(i)--;
2470       }
2471     
2472       // Poles are the current poles 
2473       // but the poles congruent to step are duplicated.
2474       
2475       Standard_Integer offset = 0;
2476
2477       for (i = 0; i < nbp; i++) {
2478         offset++;
2479
2480         for (k = 0; k < Dimension; k++) {
2481           tempc1((offset-1)*Dimension+k+1) = 
2482             wpoles(NewPoles.Lower()+i*Dimension+k);
2483         }
2484         if (i % (curDeg+1) == step) {
2485           offset++;
2486
2487           for (k = 0; k < Dimension; k++) {
2488             tempc1((offset-1)*Dimension+k+1) = 
2489               wpoles(NewPoles.Lower()+i*Dimension+k);
2490           }
2491         }
2492       }
2493         
2494       // Knots multiplicities are increased
2495       // For each knot where the sum of multiplicities is congruent to step
2496       
2497       Standard_Integer stepmult = step+1;
2498       Standard_Integer nbknots = 0;
2499       Standard_Integer smult = 0;
2500
2501       for (k = 1; k <= nbwknots; k++) {
2502         smult += wmults(k);
2503         if (smult  >= stepmult) {
2504           // this knot is increased
2505           stepmult += curDeg+1;
2506           wmults(k)++;
2507         }
2508         else {
2509           // this knot is inserted
2510           nbknots++;
2511           iknots(nbknots) = wknots(k);
2512         }
2513       }
2514       
2515       // the curve is obtained by inserting the knots
2516       // to raise the multiplicities
2517
2518       // we build "slices" of the arrays to set the correct size
2519       if (nbknots > 0) {
2520         TColStd_Array1OfReal aknots(iknots(1),1,nbknots);
2521         TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension);
2522         TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp   * Dimension);
2523 //      InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2524 //                  aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.));
2525
2526         InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults,
2527                     aknots,NoMults(),ncurve,nknots,wmults,0.0);
2528         
2529         // add to the average
2530
2531         for (i = 1; i <= nbwp * Dimension; i++)
2532           nwpoles(i) += ncurve(i);
2533       }
2534       else {
2535         // add to the average
2536
2537         for (i = 1; i <= nbwp * Dimension; i++)
2538           nwpoles(i) += tempc1(i);
2539       }
2540     }
2541     
2542     // The result is the average
2543
2544     for (i = 1; i <= nbwp * Dimension; i++) {
2545       wpoles(i) = nwpoles(i) / (curDeg+1);
2546     }
2547   }
2548   
2549   //-----------------
2550   // Copy the results
2551   //-----------------
2552
2553   // index in new knots of the first knot of the curve
2554   if (Periodic)
2555     firstknot = Mults.Upper() - l + 1;
2556   else 
2557     firstknot = f;
2558   
2559   // the new curve starts at index firstknot
2560   // so we must remove knots until the sum of multiplicities
2561   // from the first to the start is NewDegree+1
2562
2563   // m is the current sum of multiplicities
2564   m = 0;
2565
2566   for (k = 1; k <= firstknot; k++)
2567     m += wmults(k);
2568
2569   // compute the new first knot (k), pf will be the index of the first pole
2570   k = 1;
2571   pf = 0;
2572
2573   while (m > NewDegree+1) {
2574     k++;
2575     m  -= wmults(k);
2576     pf += wmults(k);
2577   }
2578   if (m < NewDegree+1) {
2579     k--;
2580     wmults(k) += m - NewDegree - 1;
2581     pf        += m - NewDegree - 1;
2582   }
2583
2584   // on a periodic curve the knots start with firstknot
2585   if (Periodic)
2586     k = firstknot;
2587
2588   // copy knots
2589
2590   for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) {
2591     NewKnots(i) = wknots(k);
2592     NewMults(i) = wmults(k);
2593     k++;
2594   }
2595
2596   // copy poles
2597   pf *= Dimension;
2598
2599   for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) {
2600     pf++;
2601     NewPoles(i) = wpoles(pf);
2602   }
2603 }
2604
2605 //=======================================================================
2606 //function : PrepareUnperiodize
2607 //purpose  : 
2608 //=======================================================================
2609
2610 void  BSplCLib::PrepareUnperiodize
2611 (const Standard_Integer         Degree, 
2612  const TColStd_Array1OfInteger& Mults, 
2613  Standard_Integer&        NbKnots, 
2614  Standard_Integer&        NbPoles)
2615 {
2616   Standard_Integer i;
2617   // initialize NbKnots and NbPoles
2618   NbKnots = Mults.Length();
2619   NbPoles = - Degree - 1;
2620
2621   for (i = Mults.Lower(); i <= Mults.Upper(); i++) 
2622     NbPoles += Mults(i);
2623
2624   Standard_Integer sigma, k;
2625   // Add knots at the beginning of the curve to raise Multiplicities 
2626   // to Degre + 1;
2627   sigma = Mults(Mults.Lower());
2628   k = Mults.Upper() - 1;
2629
2630   while ( sigma < Degree + 1) {
2631     sigma   += Mults(k);
2632     NbPoles += Mults(k);
2633     k--;
2634     NbKnots++;
2635   }
2636   // We must add exactly until Degree + 1 -> 
2637   //    Supress the excedent.
2638   if ( sigma > Degree + 1)
2639     NbPoles -= sigma - Degree - 1;
2640
2641   // Add knots at the end of the curve to raise Multiplicities 
2642   // to Degre + 1;
2643   sigma = Mults(Mults.Upper());
2644   k = Mults.Lower() + 1;
2645
2646   while ( sigma < Degree + 1) {
2647     sigma   += Mults(k);
2648     NbPoles += Mults(k);
2649     k++;
2650     NbKnots++;
2651   }
2652   // We must add exactly until Degree + 1 -> 
2653   //    Supress the excedent.
2654   if ( sigma > Degree + 1)
2655     NbPoles -= sigma - Degree - 1;
2656 }
2657
2658 //=======================================================================
2659 //function : Unperiodize
2660 //purpose  : 
2661 //=======================================================================
2662
2663 void  BSplCLib::Unperiodize
2664 (const Standard_Integer         Degree,
2665  const Standard_Integer         , // Dimension,
2666  const TColStd_Array1OfInteger& Mults,
2667  const TColStd_Array1OfReal&    Knots,
2668  const TColStd_Array1OfReal&    Poles,
2669  TColStd_Array1OfInteger& NewMults,
2670  TColStd_Array1OfReal&    NewKnots,
2671  TColStd_Array1OfReal&    NewPoles)
2672 {
2673   Standard_Integer sigma, k, index = 0;
2674   // evaluation of index : number of knots to insert before knot(1) to
2675   // raise sum of multiplicities to <Degree + 1>
2676   sigma = Mults(Mults.Lower());
2677   k = Mults.Upper() - 1;
2678
2679   while ( sigma < Degree + 1) {
2680     sigma   += Mults(k);
2681     k--;
2682     index++;
2683   }
2684
2685   Standard_Real    period = Knots(Knots.Upper()) - Knots(Knots.Lower());
2686
2687   // set the 'interior' knots;
2688
2689   for ( k = 1; k <= Knots.Length(); k++) {
2690     NewKnots ( k + index ) = Knots( k);
2691     NewMults ( k + index ) = Mults( k);
2692   }
2693   
2694   // set the 'starting' knots;
2695
2696   for ( k = 1; k <= index; k++) {
2697     NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period;
2698     NewMults( k) = NewMults( k + Knots.Length() - 1);
2699   }
2700   NewMults( 1) -= sigma - Degree -1;
2701   
2702   // set the 'ending' knots;
2703   sigma = NewMults( index + Knots.Length() );
2704
2705   for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) {
2706     NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period;
2707     NewMults( k) = NewMults( k - Knots.Length() + 1);
2708     sigma += NewMults( k - Knots.Length() + 1);
2709   }
2710   NewMults(NewMults.Length()) -= sigma - Degree - 1;
2711
2712   for ( k = 1 ; k <= NewPoles.Length(); k++) {
2713     NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1);
2714   }
2715 }
2716
2717 //=======================================================================
2718 //function : PrepareTrimming
2719 //purpose  : 
2720 //=======================================================================
2721
2722 void BSplCLib::PrepareTrimming(const Standard_Integer         Degree,
2723                                const Standard_Boolean         Periodic,
2724                                const TColStd_Array1OfReal&    Knots,
2725                                const TColStd_Array1OfInteger& Mults,
2726                                const Standard_Real            U1,
2727                                const Standard_Real            U2,
2728                                      Standard_Integer&        NbKnots,
2729                                      Standard_Integer&        NbPoles)
2730 {
2731   Standard_Integer i;
2732   Standard_Real NewU1, NewU2;
2733   Standard_Integer index1 = 0, index2 = 0;
2734
2735   // Eval index1, index2 : position of U1 and U2 in the Array Knots
2736   // such as Knots(index1-1) <= U1 < Knots(index1)
2737   //         Knots(index2-1) <= U2 < Knots(index2)
2738   LocateParameter( Degree, Knots, Mults, U1, Periodic,
2739                    Knots.Lower(), Knots.Upper(), index1, NewU1);
2740   LocateParameter( Degree, Knots, Mults, U2, Periodic,
2741                    Knots.Lower(), Knots.Upper(), index2, NewU2);
2742   index1++;
2743   if ( Abs(Knots(index2) - U2) <= Epsilon( U1))
2744     index2--;
2745
2746   // eval NbKnots:
2747   NbKnots = index2 - index1 + 3;
2748
2749   // eval NbPoles:
2750   NbPoles = Degree + 1;
2751
2752   for ( i = index1; i <= index2; i++) 
2753     NbPoles += Mults(i);
2754 }
2755
2756 //=======================================================================
2757 //function : Trimming
2758 //purpose  : 
2759 //=======================================================================
2760
2761 void BSplCLib::Trimming(const Standard_Integer         Degree,
2762                         const Standard_Boolean         Periodic,
2763                         const Standard_Integer         Dimension,
2764                         const TColStd_Array1OfReal&    Knots,
2765                         const TColStd_Array1OfInteger& Mults,
2766                         const TColStd_Array1OfReal&    Poles,
2767                         const Standard_Real            U1,
2768                         const Standard_Real            U2,
2769                               TColStd_Array1OfReal&    NewKnots,
2770                               TColStd_Array1OfInteger& NewMults,
2771                               TColStd_Array1OfReal&    NewPoles)
2772 {
2773   Standard_Integer i, nbpoles, nbknots;
2774   Standard_Real kk[2];
2775   Standard_Integer mm[2];
2776   TColStd_Array1OfReal    K( kk[0], 1, 2 );
2777   TColStd_Array1OfInteger M( mm[0], 1, 2 );
2778
2779   K(1) = U1;  K(2) = U2;
2780   mm[0] = mm[1] = Degree;
2781   if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, M, 
2782                           nbpoles, nbknots, Epsilon( U1), 0))
2783     Standard_OutOfRange::Raise();
2784
2785   TColStd_Array1OfReal    TempPoles(1, nbpoles*Dimension);
2786   TColStd_Array1OfReal    TempKnots(1, nbknots);
2787   TColStd_Array1OfInteger TempMults(1, nbknots);
2788
2789 //
2790 // do not allow the multiplicities to Add : they must be less than Degree
2791 //
2792   InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults,
2793               K, M, TempPoles, TempKnots, TempMults, Epsilon(U1),
2794               Standard_False);
2795
2796   // find in TempPoles the index of the pole corresponding to U1
2797   Standard_Integer Kindex = 0, Pindex;
2798   Standard_Real NewU1;
2799   LocateParameter( Degree, TempKnots, TempMults, U1, Periodic,
2800                    TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1);
2801   Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults);
2802   Pindex *= Dimension;
2803
2804   for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i);
2805
2806   for ( i = 1; i <= NewKnots.Length(); i++) {
2807     NewKnots(i) = TempKnots( Kindex+i-1);
2808     NewMults(i) = TempMults( Kindex+i-1);
2809   }
2810   NewMults(1) = Min(Degree, NewMults(1)) + 1 ;
2811   NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ;
2812 }
2813
2814 //=======================================================================
2815 //function : Solves a LU factored Matrix 
2816 //purpose  : 
2817 //=======================================================================
2818
2819 Standard_Integer 
2820 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
2821                             const Standard_Integer UpperBandWidth,
2822                             const Standard_Integer LowerBandWidth,
2823                             const Standard_Integer ArrayDimension,
2824                             Standard_Real&   Array) 
2825 {
2826   Standard_Integer ii,
2827   jj,
2828   kk,
2829   MinIndex,
2830   MaxIndex,
2831   ReturnCode = 0 ;
2832   
2833   Standard_Real   *PolesArray = &Array ;
2834   Standard_Real   Inverse ;
2835   
2836   
2837   if (Matrix.LowerCol() != 1 || 
2838       Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2839     ReturnCode = 1 ;
2840     goto FINISH ;
2841   }
2842   
2843   for (ii = Matrix.LowerRow() + 1; ii <=  Matrix.UpperRow() ; ii++) {
2844     MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ?
2845                 ii - LowerBandWidth : Matrix.LowerRow()) ;
2846     
2847     for ( jj = MinIndex  ; jj < ii  ; jj++) {
2848       
2849       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2850         PolesArray[(ii-1) * ArrayDimension + kk] += 
2851           PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2852       }
2853     }
2854   }
2855   
2856   for (ii = Matrix.UpperRow() ; ii >=  Matrix.LowerRow() ; ii--) {
2857     MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ? 
2858                 ii + UpperBandWidth : Matrix.UpperRow()) ;
2859     
2860     for (jj = MaxIndex  ; jj > ii ; jj--) {
2861       
2862       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2863         PolesArray[(ii-1)  * ArrayDimension + kk] -=
2864           PolesArray[(jj - 1) * ArrayDimension + kk] * 
2865             Matrix(ii, jj - ii + LowerBandWidth + 1) ;
2866       }
2867     }
2868     
2869     //fixing a bug PRO18577 to avoid divizion by zero
2870     
2871     Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ;
2872     Standard_Real Toler = 1.0e-16;
2873     if ( Abs(divizor) > Toler )
2874       Inverse = 1.0e0 / divizor ;
2875     else {
2876       Inverse = 1.0e0;
2877 //      cout << "  BSplCLib::SolveBandedSystem() : zero determinant " << endl;
2878       ReturnCode = 1;
2879       goto FINISH;
2880     }
2881         
2882     for (kk = 0 ; kk < ArrayDimension ; kk++) {
2883       PolesArray[(ii-1)  * ArrayDimension + kk] *=  Inverse ; 
2884     }
2885   }
2886   FINISH :
2887     return (ReturnCode) ;
2888 }
2889
2890 //=======================================================================
2891 //function : Solves a LU factored Matrix 
2892 //purpose  : 
2893 //=======================================================================
2894
2895 Standard_Integer 
2896 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
2897                             const Standard_Integer UpperBandWidth,
2898                             const Standard_Integer LowerBandWidth,
2899                             const Standard_Boolean HomogeneousFlag,
2900                             const Standard_Integer ArrayDimension,
2901                             Standard_Real&   Poles,
2902                             Standard_Real&   Weights) 
2903 {
2904   Standard_Integer ii,
2905   kk,
2906   ErrorCode = 0,
2907   ReturnCode = 0 ;
2908   
2909   Standard_Real   Inverse,
2910   *PolesArray   = &Poles,
2911   *WeightsArray = &Weights ;
2912   
2913   if (Matrix.LowerCol() != 1 || 
2914       Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) {
2915     ReturnCode = 1 ;
2916     goto FINISH ;
2917   }
2918   if (HomogeneousFlag == Standard_False) {
2919     
2920     for (ii = 0 ; ii <  Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) {
2921       
2922       for (kk = 0 ; kk < ArrayDimension ; kk++) {
2923         PolesArray[ii * ArrayDimension + kk] *=
2924           WeightsArray[ii] ;
2925       }
2926     }
2927   }
2928   ErrorCode = 
2929     BSplCLib::SolveBandedSystem(Matrix,
2930                                 UpperBandWidth,
2931                                 LowerBandWidth,
2932                                 ArrayDimension,
2933                                 Poles) ;
2934   if (ErrorCode != 0) {
2935     ReturnCode = 2 ;
2936     goto FINISH ;
2937   }
2938   ErrorCode = 
2939     BSplCLib::SolveBandedSystem(Matrix,
2940                                 UpperBandWidth,
2941                                 LowerBandWidth,
2942                                 1,
2943                                 Weights) ;
2944   if (ErrorCode != 0) {
2945     ReturnCode = 3 ;
2946     goto FINISH ;
2947   }
2948   if (HomogeneousFlag == Standard_False) {
2949
2950     for (ii = 0  ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) {
2951       Inverse = 1.0e0 / WeightsArray[ii] ;
2952       
2953       for (kk = 0  ; kk < ArrayDimension ; kk++) {
2954         PolesArray[ii * ArrayDimension + kk] *= Inverse ;
2955       }
2956     }
2957   }
2958   FINISH : return (ReturnCode) ;
2959 }
2960
2961 //=======================================================================
2962 //function : BuildSchoenbergPoints
2963 //purpose  : 
2964 //=======================================================================
2965
2966 void  BSplCLib::BuildSchoenbergPoints(const Standard_Integer         Degree,
2967                                       const TColStd_Array1OfReal&    FlatKnots,
2968                                       TColStd_Array1OfReal&          Parameters) 
2969 {
2970   Standard_Integer ii,
2971   jj ;
2972   Standard_Real Inverse ;
2973   Inverse = 1.0e0 / (Standard_Real)Degree ;
2974   
2975   for (ii = Parameters.Lower() ;   ii <= Parameters.Upper() ; ii++) {
2976     Parameters(ii) = 0.0e0 ;
2977     
2978     for (jj = 1 ; jj <= Degree ; jj++) {
2979       Parameters(ii) += FlatKnots(jj + ii) ;
2980     } 
2981     Parameters(ii) *= Inverse ; 
2982   }
2983 }
2984
2985 //=======================================================================
2986 //function : Interpolate
2987 //purpose  : 
2988 //=======================================================================
2989
2990 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
2991                             const TColStd_Array1OfReal&    FlatKnots,
2992                             const TColStd_Array1OfReal&    Parameters,
2993                             const TColStd_Array1OfInteger& ContactOrderArray,
2994                             const Standard_Integer         ArrayDimension,
2995                             Standard_Real&                 Poles,
2996                             Standard_Integer&              InversionProblem) 
2997 {
2998   Standard_Integer ErrorCode,
2999   UpperBandWidth,
3000   LowerBandWidth ;
3001 //  Standard_Real *PolesArray = &Poles ;
3002   math_Matrix InterpolationMatrix(1, Parameters.Length(),
3003                                   1, 2 * Degree + 1) ;
3004   ErrorCode =
3005   BSplCLib::BuildBSpMatrix(Parameters,
3006                            ContactOrderArray,
3007                            FlatKnots,
3008                            Degree,
3009                            InterpolationMatrix,
3010                            UpperBandWidth,
3011                            LowerBandWidth) ;
3012   Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3013
3014   ErrorCode =
3015   BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3016                            UpperBandWidth,
3017                            LowerBandWidth,
3018                            InversionProblem) ;
3019   Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3020
3021   ErrorCode  =
3022   BSplCLib::SolveBandedSystem(InterpolationMatrix,
3023                               UpperBandWidth,
3024                               LowerBandWidth,
3025                               ArrayDimension,
3026                               Poles) ;
3027
3028   Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate")  ;
3029
3030
3031 //=======================================================================
3032 //function : Interpolate
3033 //purpose  : 
3034 //=======================================================================
3035
3036 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
3037                             const TColStd_Array1OfReal&    FlatKnots,
3038                             const TColStd_Array1OfReal&    Parameters,
3039                             const TColStd_Array1OfInteger& ContactOrderArray,
3040                             const Standard_Integer         ArrayDimension,
3041                             Standard_Real&                 Poles,
3042                             Standard_Real&                 Weights,
3043                             Standard_Integer&              InversionProblem) 
3044 {
3045   Standard_Integer ErrorCode,
3046   UpperBandWidth,
3047   LowerBandWidth ;
3048
3049   math_Matrix InterpolationMatrix(1, Parameters.Length(),
3050                                   1, 2 * Degree + 1) ;
3051   ErrorCode =
3052   BSplCLib::BuildBSpMatrix(Parameters,
3053                            ContactOrderArray,
3054                            FlatKnots,
3055                            Degree,
3056                            InterpolationMatrix,
3057                            UpperBandWidth,
3058                            LowerBandWidth) ;
3059   Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3060
3061   ErrorCode =
3062   BSplCLib::FactorBandedMatrix(InterpolationMatrix,
3063                            UpperBandWidth,
3064                            LowerBandWidth,
3065                            InversionProblem) ;
3066   Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ;
3067
3068   ErrorCode  =
3069   BSplCLib::SolveBandedSystem(InterpolationMatrix,
3070                               UpperBandWidth,
3071                               LowerBandWidth,
3072                               Standard_False,
3073                               ArrayDimension,
3074                               Poles,
3075                               Weights) ;
3076
3077   Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate")  ;
3078 }
3079
3080 //=======================================================================
3081 //function : Evaluates a Bspline function : uses the ExtrapMode 
3082 //purpose  : the function is extrapolated using the Taylor expansion
3083 //           of degree ExtrapMode[0] to the left and the Taylor
3084 //           expansion of degree ExtrapMode[1] to the right 
3085 //  this evaluates the numerator by multiplying by the weights
3086 //  and evaluating it but does not call RationalDerivatives after 
3087 //=======================================================================
3088
3089 void  BSplCLib::Eval
3090 (const Standard_Real                   Parameter,
3091  const Standard_Boolean                PeriodicFlag,
3092  const Standard_Integer                DerivativeRequest,
3093  Standard_Integer&                     ExtrapMode,
3094  const Standard_Integer                Degree,
3095  const  TColStd_Array1OfReal&          FlatKnots, 
3096  const Standard_Integer                ArrayDimension,
3097  Standard_Real&                        Poles,
3098  Standard_Real&                        Weights,
3099  Standard_Real&                        PolesResults,
3100  Standard_Real&                        WeightsResults)
3101 {
3102   Standard_Integer ii,
3103   jj,
3104   kk=0,
3105   Index,
3106   Index1,
3107   Index2,
3108   *ExtrapModeArray,
3109   Modulus,
3110   NewRequest,
3111   ExtrapolatingFlag[2],
3112   ErrorCode,
3113   Order = Degree + 1,
3114   FirstNonZeroBsplineIndex,
3115   LocalRequest = DerivativeRequest ;
3116   Standard_Real  *PResultArray,
3117   *WResultArray,
3118   *PolesArray,
3119   *WeightsArray,
3120   LocalParameter,
3121   Period,
3122   Inverse,
3123   Delta ;
3124   PolesArray = &Poles     ;
3125   WeightsArray = &Weights ;
3126   ExtrapModeArray = &ExtrapMode ;
3127   PResultArray = &PolesResults ;
3128   WResultArray = &WeightsResults ;
3129   LocalParameter = Parameter ;
3130   ExtrapolatingFlag[0] = 
3131     ExtrapolatingFlag[1] = 0 ;
3132   //
3133   // check if we are extrapolating to a degree which is smaller than
3134   // the degree of the Bspline
3135   //
3136   if (PeriodicFlag) {
3137     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3138
3139     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3140       LocalParameter -= Period ;
3141     }
3142     
3143     while (LocalParameter < FlatKnots(2)) {
3144       LocalParameter +=  Period ;
3145     }
3146   }
3147   if (Parameter < FlatKnots(2) && 
3148       LocalRequest < ExtrapModeArray[0] &&
3149       ExtrapModeArray[0] < Degree) {
3150     LocalRequest = ExtrapModeArray[0] ;
3151     LocalParameter = FlatKnots(2) ;
3152     ExtrapolatingFlag[0] = 1 ;
3153   }
3154   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3155       LocalRequest < ExtrapModeArray[1]  &&
3156       ExtrapModeArray[1] < Degree) {
3157     LocalRequest = ExtrapModeArray[1] ;
3158     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3159     ExtrapolatingFlag[1] = 1 ;
3160   }
3161   Delta = Parameter - LocalParameter ;
3162   if (LocalRequest >= Order) {
3163     LocalRequest = Degree ;
3164   }
3165   if (PeriodicFlag) {
3166     Modulus = FlatKnots.Length() - Degree -1 ;
3167   }
3168   else {
3169     Modulus = FlatKnots.Length() - Degree ;
3170   }
3171
3172   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3173   ErrorCode =
3174     BSplCLib::EvalBsplineBasis(1,
3175                                LocalRequest,
3176                                Order,
3177                                FlatKnots,
3178                                LocalParameter,
3179                                FirstNonZeroBsplineIndex,
3180                                BsplineBasis) ;
3181   if (ErrorCode != 0) {
3182     goto FINISH ;
3183   }
3184   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3185     Index = 0 ;
3186     Index2 = 0 ;
3187
3188     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3189       Index1 = FirstNonZeroBsplineIndex ;
3190
3191       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3192         PResultArray[Index + kk] = 0.0e0 ;
3193       }
3194       WResultArray[Index] = 0.0e0 ;
3195
3196       for (jj = 1  ; jj <= Order ; jj++) {
3197         
3198         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3199           PResultArray[Index + kk] += 
3200             PolesArray[(Index1-1) * ArrayDimension + kk] 
3201               * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3202         }
3203         WResultArray[Index2]  += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3204         
3205         Index1 = Index1 % Modulus ;
3206         Index1 += 1 ;
3207       }
3208       Index += ArrayDimension ;
3209       Index2 += 1 ;
3210     }
3211   }
3212   else {
3213     // 
3214     //  store Taylor expansion in LocalRealArray
3215     //
3216     NewRequest = DerivativeRequest ;
3217     if (NewRequest > Degree) {
3218       NewRequest = Degree ;
3219     }
3220     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3221     Index = 0 ;
3222     Inverse = 1.0e0 ;
3223
3224     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3225       Index1 = FirstNonZeroBsplineIndex ;
3226       
3227       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3228         LocalRealArray[Index + kk] = 0.0e0 ;
3229       }
3230
3231       for (jj = 1  ; jj <= Order ; jj++) {
3232
3233         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3234           LocalRealArray[Index + kk] += 
3235             PolesArray[(Index1-1)*ArrayDimension + kk] * 
3236               WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3237         }
3238         Index1 = Index1 % Modulus ;
3239         Index1 += 1 ;
3240       }
3241
3242       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3243         LocalRealArray[Index + kk] *= Inverse ;
3244       }
3245       Index += ArrayDimension ;
3246       Inverse /= (Standard_Real) ii ;
3247     }
3248     PLib::EvalPolynomial(Delta,
3249                          NewRequest,
3250                          Degree,
3251                          ArrayDimension,
3252                          LocalRealArray[0],
3253                          PolesResults) ;
3254     Index = 0 ;
3255     Inverse = 1.0e0 ;
3256
3257     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3258       Index1 = FirstNonZeroBsplineIndex ;
3259       LocalRealArray[Index] = 0.0e0 ;
3260
3261       for (jj = 1  ; jj <= Order ; jj++) {
3262         LocalRealArray[Index] += 
3263           WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3264         Index1 = Index1 % Modulus ;
3265         Index1 += 1 ;
3266       }
3267       LocalRealArray[Index + kk] *= Inverse ;
3268       Index += 1 ;
3269       Inverse /= (Standard_Real) ii ;
3270     }
3271     PLib::EvalPolynomial(Delta,
3272                          NewRequest,
3273                          Degree,
3274                          1,
3275                          LocalRealArray[0],
3276                          WeightsResults) ;
3277   }
3278   FINISH : ;
3279 }
3280
3281 //=======================================================================
3282 //function : Evaluates a Bspline function : uses the ExtrapMode 
3283 //purpose  : the function is extrapolated using the Taylor expansion
3284 //           of degree ExtrapMode[0] to the left and the Taylor
3285 //           expansion of degree ExtrapMode[1] to the right 
3286 // WARNING : the array Results is supposed to have at least 
3287 // (DerivativeRequest + 1) * ArrayDimension slots and the 
3288 // 
3289 //=======================================================================
3290
3291 void  BSplCLib::Eval
3292 (const Standard_Real                   Parameter,
3293  const Standard_Boolean                PeriodicFlag,
3294  const Standard_Integer                DerivativeRequest,
3295  Standard_Integer&                     ExtrapMode,
3296  const Standard_Integer                Degree,
3297  const  TColStd_Array1OfReal&          FlatKnots, 
3298  const Standard_Integer                ArrayDimension,
3299  Standard_Real&                        Poles,
3300  Standard_Real&                        Results) 
3301 {
3302   Standard_Integer ii,
3303   jj,
3304   kk,
3305   Index,
3306   Index1,
3307   *ExtrapModeArray,
3308   Modulus,
3309   NewRequest,
3310   ExtrapolatingFlag[2],
3311   ErrorCode,
3312   Order = Degree + 1,
3313   FirstNonZeroBsplineIndex,
3314   LocalRequest = DerivativeRequest ;
3315
3316   Standard_Real  *ResultArray,
3317   *PolesArray,
3318   LocalParameter,
3319   Period,
3320   Inverse,
3321   Delta ;
3322          
3323   PolesArray = &Poles ;
3324   ExtrapModeArray = &ExtrapMode ;
3325   ResultArray = &Results ;  
3326   LocalParameter = Parameter ;
3327   ExtrapolatingFlag[0] = 
3328     ExtrapolatingFlag[1] = 0 ;
3329   //
3330   // check if we are extrapolating to a degree which is smaller than
3331   // the degree of the Bspline
3332   //
3333   if (PeriodicFlag) {
3334     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3335
3336     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3337       LocalParameter -= Period ;
3338     }
3339
3340     while (LocalParameter < FlatKnots(2)) {
3341       LocalParameter +=  Period ;
3342     }
3343   }
3344   if (Parameter < FlatKnots(2) && 
3345       LocalRequest < ExtrapModeArray[0] &&
3346       ExtrapModeArray[0] < Degree) {
3347     LocalRequest = ExtrapModeArray[0] ;
3348     LocalParameter = FlatKnots(2) ;
3349     ExtrapolatingFlag[0] = 1 ;
3350   }
3351   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3352       LocalRequest < ExtrapModeArray[1]  &&
3353       ExtrapModeArray[1] < Degree) {
3354     LocalRequest = ExtrapModeArray[1] ;
3355     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3356     ExtrapolatingFlag[1] = 1 ;
3357   }
3358   Delta = Parameter - LocalParameter ;
3359   if (LocalRequest >= Order) {
3360     LocalRequest = Degree ;
3361   }
3362   
3363   if (PeriodicFlag) {
3364     Modulus = FlatKnots.Length() - Degree -1 ;
3365   }
3366   else {
3367     Modulus = FlatKnots.Length() - Degree ;
3368   }
3369   
3370   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3371   
3372   ErrorCode =
3373     BSplCLib::EvalBsplineBasis(1,
3374                                LocalRequest,
3375                                Order,
3376                                FlatKnots,
3377                                LocalParameter,
3378                                FirstNonZeroBsplineIndex,
3379                                BsplineBasis);
3380   if (ErrorCode != 0) {
3381     goto FINISH ;
3382   }
3383   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3384     Index = 0 ;
3385     
3386     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3387       Index1 = FirstNonZeroBsplineIndex ;
3388       
3389       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3390         ResultArray[Index + kk] = 0.0e0 ;
3391       }
3392
3393       for (jj = 1  ; jj <= Order ; jj++) {
3394         
3395         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3396           ResultArray[Index + kk] += 
3397             PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3398         }
3399         Index1 = Index1 % Modulus ;
3400         Index1 += 1 ;
3401       }
3402       Index += ArrayDimension ;
3403     }
3404   }
3405   else {
3406     // 
3407     //  store Taylor expansion in LocalRealArray
3408     //
3409     NewRequest = DerivativeRequest ;
3410     if (NewRequest > Degree) {
3411       NewRequest = Degree ;
3412     }
3413     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3414
3415     Index = 0 ;
3416     Inverse = 1.0e0 ;
3417
3418     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3419       Index1 = FirstNonZeroBsplineIndex ;
3420       
3421       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3422         LocalRealArray[Index + kk] = 0.0e0 ;
3423       }
3424
3425       for (jj = 1  ; jj <= Order ; jj++) {
3426
3427         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3428           LocalRealArray[Index + kk] += 
3429             PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3430         }
3431         Index1 = Index1 % Modulus ;
3432         Index1 += 1 ;
3433       }
3434
3435       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3436         LocalRealArray[Index + kk] *= Inverse ;
3437       }
3438       Index += ArrayDimension ;
3439       Inverse /= (Standard_Real) ii ;
3440     }
3441     PLib::EvalPolynomial(Delta,
3442                          NewRequest,
3443                          Degree,
3444                          ArrayDimension,
3445                          LocalRealArray[0],
3446                          Results) ;
3447   }
3448   FINISH : ;
3449 }
3450
3451 //=======================================================================
3452 //function : TangExtendToConstraint 
3453 //purpose  : Extends a Bspline function using the tangency map
3454 // WARNING :  
3455 //  
3456 // 
3457 //=======================================================================
3458
3459 void  BSplCLib::TangExtendToConstraint
3460 (const  TColStd_Array1OfReal&          FlatKnots, 
3461  const Standard_Real                   C1Coefficient,
3462  const Standard_Integer                NumPoles,
3463  Standard_Real&                        Poles,
3464  const Standard_Integer                CDimension,
3465  const Standard_Integer                CDegree,
3466  const  TColStd_Array1OfReal&          ConstraintPoint, 
3467  const Standard_Integer                Continuity,
3468  const Standard_Boolean                After,
3469  Standard_Integer&                     NbPolesResult,
3470  Standard_Integer&                     NbKnotsResult,
3471  Standard_Real&                        KnotsResult, 
3472  Standard_Real&                        PolesResult) 
3473 {
3474 #if DEB
3475   if (CDegree<Continuity+1) {
3476     cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3477   }
3478 #endif
3479   Standard_Real * Padr = &Poles ;
3480   Standard_Real * KRadr = &KnotsResult ;
3481   Standard_Real * PRadr = &PolesResult ;
3482
3483 ////////////////////////////////////////////////////////////////////////
3484 //
3485 //    1. calculation of extension nD
3486 //
3487 ////////////////////////////////////////////////////////////////////////
3488
3489 //  Hermite matrix
3490   Standard_Integer Csize = Continuity + 2;
3491   math_Matrix  MatCoefs(1,Csize, 1,Csize);
3492   if (After) {
3493     PLib::HermiteCoefficients(0, 1,           // Limits 
3494                               Continuity, 0,  // Orders of constraints
3495                               MatCoefs);
3496   }
3497   else {
3498     PLib::HermiteCoefficients(0, 1,           // Limits 
3499                               0, Continuity,  // Orders of constraints
3500                               MatCoefs);    
3501   }
3502
3503
3504 //  position at the node of connection
3505   Standard_Real Tbord ;
3506   if (After) {
3507     Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3508   }
3509   else {
3510     Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3511   }
3512   Standard_Boolean periodic_flag = Standard_False ;
3513   Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3514   extrap_mode[0] = extrap_mode[1] = CDegree;
3515   TColStd_Array1OfReal  EvalBS(1, CDimension * (derivative_request+1)) ; 
3516   Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3517   BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3518                   CDegree,FlatKnots,CDimension,Poles,*Eadr);
3519
3520 //  norm of the tangent at the node of connection
3521   math_Vector Tgte(1,CDimension);
3522
3523   for (ipos=1;ipos<=CDimension;ipos++) {
3524     Tgte(ipos) = EvalBS(ipos+CDimension);
3525   }
3526   Standard_Real L1=Tgte.Norm();
3527
3528
3529 //  matrix of constraints
3530   math_Matrix Contraintes(1,Csize,1,CDimension);
3531   if (After) {
3532
3533     for (ipos=1;ipos<=CDimension;ipos++) {
3534       Contraintes(1,ipos) = EvalBS(ipos);
3535       Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3536       if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3537       if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3538       Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3539     }
3540   }
3541   else {
3542
3543     for (ipos=1;ipos<=CDimension;ipos++) {
3544       Contraintes(1,ipos) = ConstraintPoint(ipos);
3545       Contraintes(2,ipos) = EvalBS(ipos);
3546       if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3547       if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3548       if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3549     }
3550   }
3551
3552 //  calculate the coefficients of extension
3553   Standard_Integer ii, jj, kk;
3554   TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3555   ExtraCoeffs.Init(0.);
3556
3557   for (ii=1; ii<=Csize; ii++) {
3558
3559     for (jj=1; jj<=Csize; jj++) {
3560
3561       for (kk=1; kk<=CDimension; kk++) {
3562         ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3563       }
3564     }
3565   }
3566
3567 //  calculate the poles of extension
3568   TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3569   Standard_Real * EPadr = &ExtrapPoles(1) ;
3570   PLib::CoefficientsPoles(CDimension,
3571                           ExtraCoeffs,  PLib::NoWeights(),
3572                           ExtrapPoles,  PLib::NoWeights());
3573
3574 //  calculate the nodes of extension with multiplicities
3575   TColStd_Array1OfReal ExtrapNoeuds(1,2);
3576   ExtrapNoeuds(1) = 0.;
3577   ExtrapNoeuds(2) = 1.;
3578   TColStd_Array1OfInteger ExtrapMults(1,2);
3579   ExtrapMults(1) = Csize;
3580   ExtrapMults(2) = Csize;
3581
3582 // flat nodes of extension
3583   TColStd_Array1OfReal FK2(1, Csize*2);
3584   BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3585
3586 //  norm of the tangent at the connection point 
3587   if (After) {
3588     BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3589                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3590   }
3591   else {
3592     BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3593                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3594   }
3595
3596   for (ipos=1;ipos<=CDimension;ipos++) {
3597     Tgte(ipos) = EvalBS(ipos+CDimension);
3598   }
3599   Standard_Real L2 = Tgte.Norm();
3600
3601 //  harmonisation of degrees
3602   TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3603   TColStd_Array1OfReal NewK2(1, 2);
3604   TColStd_Array1OfInteger NewM2(1, 2);
3605   if (Csize-1<CDegree) {
3606     BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3607                              ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3608                              NewP2,NewK2,NewM2);
3609   }
3610   else {
3611     NewP2 = ExtrapPoles;
3612     NewK2 = ExtrapNoeuds;
3613     NewM2 = ExtrapMults;
3614   }
3615
3616 //  flat nodes of extension after harmonization of degrees
3617   TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3618   BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3619
3620
3621 ////////////////////////////////////////////////////////////////////////
3622 //
3623 //    2.  concatenation C0
3624 //
3625 ////////////////////////////////////////////////////////////////////////
3626
3627 //  ratio of reparametrization
3628   Standard_Real Ratio=1, Delta;
3629   if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3630     Ratio = L2 / L1;
3631   }
3632   if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3633
3634   if (After) {
3635 //    do not touch the first BSpline
3636     Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3637   }
3638   else {
3639 //    do not touch the second BSpline
3640     Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower());
3641   }
3642
3643 //  result of the concatenation
3644   Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1;
3645   Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1);
3646   TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension);
3647   TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2);
3648
3649 //  poles
3650   Standard_Integer indNP, indP, indEP;
3651   if (After) {
3652
3653     for (ii=1;  ii<=NbP1+NbP2-1; ii++) {
3654
3655       for (jj=1;  jj<=CDimension; jj++) {
3656         indNP = (ii-1)*CDimension+jj;
3657         indP = (ii-1)*CDimension+jj-1;
3658         indEP = (ii-NbP1)*CDimension+jj;
3659         if (ii<NbP1) NewPoles(indNP) =  Padr[indP];
3660         else NewPoles(indNP) = NewP2(indEP);
3661       }
3662     }
3663   }
3664   else {
3665
3666     for (ii=1;  ii<=NbP1+NbP2-1; ii++) {
3667
3668       for (jj=1;  jj<=CDimension; jj++) {
3669         indNP = (ii-1)*CDimension+jj;
3670         indEP = (ii-1)*CDimension+jj;
3671         indP = (ii-NbP2)*CDimension+jj-1;
3672         if (ii<NbP2) NewPoles(indNP) =  NewP2(indEP);
3673         else NewPoles(indNP) = Padr[indP];
3674       }
3675     }
3676   }
3677
3678 //  flat nodes 
3679   if (After) {
3680 //    start with the nodes of the initial surface
3681
3682     for (ii=1; ii<NbK1; ii++) {
3683       NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1);
3684     }
3685 //    continue with the reparameterized nodes of the extension
3686
3687     for (ii=1; ii<=NbK2-CDegree-1; ii++) {
3688       NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta;
3689     }
3690   }
3691   else {
3692 //    start with the reparameterized nodes of the extension
3693
3694     for (ii=1; ii<NbK2-CDegree; ii++) {
3695       NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta;
3696     }
3697 //    continue with the nodes of the initial surface
3698
3699     for (ii=2; ii<=NbK1; ii++) {
3700       NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1);
3701     }
3702   }
3703
3704
3705 ////////////////////////////////////////////////////////////////////////
3706 //
3707 //    3.  reduction of multiplicite at the node of connection
3708 //
3709 ////////////////////////////////////////////////////////////////////////
3710
3711 //  number of separate nodes
3712   Standard_Integer KLength = 1;
3713
3714   for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3715     if (NewFlats(ii) != NewFlats(ii-1)) KLength++;
3716   }
3717
3718 //  flat nodes --> nodes + multiplicities
3719   TColStd_Array1OfReal NewKnots (1, KLength);
3720   TColStd_Array1OfInteger NewMults (1, KLength);
3721   NewMults.Init(1);
3722   jj = 1;
3723   NewKnots(jj) = NewFlats(1);
3724
3725   for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3726     if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++;
3727     else {
3728       jj++;
3729       NewKnots(jj) = NewFlats(ii);
3730     }
3731   }
3732
3733 //  reduction of multiplicity at the second or the last but one node
3734   Standard_Integer Index = 2, M = CDegree;
3735   if (After) Index = KLength-1;
3736   TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension);
3737   TColStd_Array1OfReal ResultKnots (1, KLength);
3738   TColStd_Array1OfInteger ResultMults (1, KLength);
3739   Standard_Real Tol = 1.e-6;
3740   Standard_Boolean Ok = Standard_True;
3741
3742   while ( (M>CDegree-Continuity) && Ok) {
3743     Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension,
3744                     NewPoles, NewKnots, NewMults,
3745                     ResultPoles, ResultKnots, ResultMults, Tol);
3746     if (Ok) M--;
3747   }
3748
3749   if (M == CDegree) {
3750 //    number of poles of the concatenation
3751     NbPolesResult = NbP1 + NbP2 - 1;
3752 //    the poles of the concatenation
3753     Standard_Integer PLength = NbPolesResult*CDimension;
3754
3755     for (jj=1; jj<=PLength; jj++) {
3756       PRadr[jj-1] = NewPoles(jj);
3757     }
3758   
3759 //    flat nodes of the concatenation
3760     Standard_Integer ideb = 0;
3761
3762     for (jj=0; jj<NewKnots.Length(); jj++) {
3763       for (ii=0; ii<NewMults(jj+1); ii++) {
3764         KRadr[ideb+ii] = NewKnots(jj+1);
3765       }
3766       ideb += NewMults(jj+1);
3767     }
3768     NbKnotsResult = ideb;
3769   }
3770
3771   else {
3772 //    number of poles of the result
3773     NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M;
3774 //    the poles of the result
3775     Standard_Integer PLength = NbPolesResult*CDimension;
3776
3777     for (jj=0; jj<PLength; jj++) {
3778       PRadr[jj] = ResultPoles(jj+1);
3779     }
3780   
3781 //    flat nodes of the result
3782     Standard_Integer ideb = 0;
3783
3784     for (jj=0; jj<ResultKnots.Length(); jj++) {
3785       for (ii=0; ii<ResultMults(jj+1); ii++) {
3786         KRadr[ideb+ii] = ResultKnots(jj+1);
3787       }
3788       ideb += ResultMults(jj+1);
3789     }
3790     NbKnotsResult = ideb;
3791   }
3792 }
3793
3794 //=======================================================================
3795 //function : Resolution
3796 //purpose  : 
3797 //                           d
3798 //  Let C(t) = SUM      Ci Bi(t)  a Bspline curve of degree d  
3799 //            i = 1,n      
3800 //  with nodes tj for j = 1,n+d+1 
3801 //
3802 //
3803 //         '                    C1 - Ci-1   d-1
3804 //  Then C (t) = SUM     d *  ---------  Bi (t) 
3805 //                i = 2,n      ti+d - ti
3806 //
3807 //                          d-1
3808 //  for the base of BSpline  Bi  (t) of degree d-1.
3809 //
3810 //  Consequently the upper bound of the norm of the derivative from C is :
3811 //
3812 //
3813 //                        |  Ci - Ci-1  |
3814 //          d *   Max     |  ---------  |
3815 //              i = 2,n |  ti+d - ti  |
3816 //     
3817 //                                      N(t) 
3818 //  In the rational case set    C(t) = -----
3819 //                                      D(t) 
3820 //
3821 //  
3822 //  D(t) =  SUM    Di Bi(t) 
3823 //        i=1,n
3824 //
3825 //  N(t) =  SUM   Di * Ci Bi(t) 
3826 //          i =1,n
3827 //
3828 //          N'(t)  -    D'(t) C(t) 
3829 //   C'(t) = -----------------------
3830 //                   D(t)
3831 //
3832 //                                   
3833 //   N'(t) - D'(t) C(t) = 
3834 //      
3835 //                     Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t))    d-1
3836 //      SUM   d *   ---------------------------------------- * Bi  (t)  =
3837 //        i=2,n                   ti+d   - ti
3838 //
3839 //    
3840 //                   Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)                d-1
3841 // SUM   SUM     d * -----------------------------------  * Betaj(t) * Bi  (t) 
3842 //i=2,n j=1,n               ti+d  - ti  
3843 //  
3844 //
3845 //
3846 //                 Dj Bj(t) 
3847 //    Betaj(t) =   --------
3848 //                 D(t) 
3849 //
3850 //  Betaj(t) form a partition >= 0 of the entity with support
3851 //  tj, tj+d+1. Consequently if Rj = {j-d, ....,  j+d+d+1} 
3852 //  obtain an upper bound of the derivative of C by taking :
3853 //
3854 //
3855 //
3856 //
3857 //
3858 //    
3859 //                         Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) 
3860 //   Max   Max       d  *  -----------------------------------  
3861 // j=1,n  i dans Rj                   ti+d  - ti  
3862 //
3863 //  --------------------------------------------------------
3864 //
3865 //               Min    Di
3866 //              i =1,n
3867 //  
3868 //
3869 //=======================================================================
3870
3871 void BSplCLib::Resolution(      Standard_Real&        Poles,
3872                           const Standard_Integer      ArrayDimension,
3873                           const Standard_Integer      NumPoles,
3874                           const TColStd_Array1OfReal& Weights,
3875                           const TColStd_Array1OfReal& FlatKnots,
3876                           const Standard_Integer      Degree,
3877                           const Standard_Real         Tolerance3D,
3878                           Standard_Real&              UTolerance) 
3879 {
3880   Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim;
3881   Standard_Integer lower,upper,ii_minus,jj,ii_miDim;
3882   Standard_Integer Deg1 = Degree + 1;
3883   Standard_Integer Deg2 = (Degree << 1) + 1;
3884   Standard_Real value,factor,W,min_weights,inverse;
3885   Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3;
3886   Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3;
3887   Standard_Real wg_ii_index, wg_ii_minus;
3888   Standard_Real *PA,max_derivative;
3889   const Standard_Real * FK = &FlatKnots(FlatKnots.Lower());
3890   PA = &Poles;
3891   max_derivative = 0.0e0;
3892   num_poles = FlatKnots.Length() - Deg1;
3893   switch (ArrayDimension) {
3894   case 2 : {
3895     if (&Weights != NULL) {
3896       const Standard_Real * WG = &Weights(Weights.Lower());
3897       min_weights = WG[0];
3898       
3899       for (ii = 1 ; ii < NumPoles ; ii++) {
3900         W = WG[ii];
3901         if (W < min_weights) min_weights = W;
3902       }
3903       
3904       for (ii = 1 ; ii < num_poles ; ii++) {
3905         ii_index = ii % NumPoles;
3906         ii_inDim = ii_index << 1;
3907         ii_minus = (ii - 1) % NumPoles;
3908         ii_miDim = ii_minus << 1;
3909         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3910         pa_ii_inDim_1 = PA[ii_inDim];
3911         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3912         pa_ii_miDim_1 = PA[ii_miDim];
3913         wg_ii_index   = WG[ii_index];
3914         wg_ii_minus   = WG[ii_minus];
3915         inverse = FK[ii + Degree] - FK[ii];
3916         inverse = 1.0e0 / inverse;
3917         lower = ii - Deg1;
3918         if (lower < 0) lower = 0;
3919         upper = Deg2 + ii;
3920         if (upper > num_poles) upper = num_poles;
3921         
3922         for (jj = lower ; jj < upper ; jj++) {
3923           jj_index = jj % NumPoles;
3924           jj_index = jj_index << 1;
3925           value = 0.0e0;
3926           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3927                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3928           if (factor < 0) factor = - factor;
3929           value += factor; jj_index++;
3930           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3931                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3932           if (factor < 0) factor = - factor;
3933           value += factor;
3934           value *= inverse;
3935           if (max_derivative < value) max_derivative = value;
3936         }
3937       }
3938       max_derivative /= min_weights;
3939     }
3940     else {
3941       
3942       for (ii = 1 ; ii < num_poles ; ii++) {
3943         ii_index = ii % NumPoles;
3944         ii_index = ii_index << 1;
3945         ii_minus = (ii - 1) % NumPoles;
3946         ii_minus = ii_minus << 1;
3947         inverse = FK[ii + Degree] - FK[ii];
3948         inverse = 1.0e0 / inverse;
3949         value = 0.0e0;
3950         factor = PA[ii_index] - PA[ii_minus];
3951         if (factor < 0) factor = - factor;
3952         value += factor; ii_index++; ii_minus++;
3953         factor = PA[ii_index] - PA[ii_minus];
3954         if (factor < 0) factor = - factor;
3955         value += factor;
3956         value *= inverse;
3957         if (max_derivative < value) max_derivative = value;
3958       }
3959     }
3960     break;
3961   }
3962   case 3 : {
3963     if (&Weights != NULL) {
3964       const Standard_Real * WG = &Weights(Weights.Lower());
3965       min_weights = WG[0];
3966       
3967       for (ii = 1 ; ii < NumPoles ; ii++) {
3968         W = WG[ii];
3969         if (W < min_weights) min_weights = W;
3970       }
3971       
3972       for (ii = 1 ; ii < num_poles ; ii++) {
3973         ii_index = ii % NumPoles;
3974         ii_inDim = (ii_index << 1) + ii_index;
3975         ii_minus = (ii - 1) % NumPoles;
3976         ii_miDim = (ii_minus << 1) + ii_minus;
3977         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3978         pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
3979         pa_ii_inDim_2 = PA[ii_inDim];
3980         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3981         pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
3982         pa_ii_miDim_2 = PA[ii_miDim];
3983         wg_ii_index   = WG[ii_index];
3984         wg_ii_minus   = WG[ii_minus];
3985         inverse = FK[ii + Degree] - FK[ii];
3986         inverse = 1.0e0 / inverse;
3987         lower = ii - Deg1;
3988         if (lower < 0) lower = 0;
3989         upper = Deg2 + ii;
3990         if (upper > num_poles) upper = num_poles;
3991         
3992         for (jj = lower ; jj < upper ; jj++) {
3993           jj_index = jj % NumPoles;
3994           jj_index = (jj_index << 1) + jj_index;
3995           value = 0.0e0;
3996           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3997                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3998           if (factor < 0) factor = - factor;
3999           value += factor; jj_index++;
4000           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4001                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4002           if (factor < 0) factor = - factor;
4003           value += factor; jj_index++;
4004           factor  = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4005                      ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4006           if (factor < 0) factor = - factor;
4007           value += factor;
4008           value *= inverse;
4009           if (max_derivative < value) max_derivative = value;
4010         }
4011       }
4012       max_derivative /= min_weights;
4013     }
4014     else {
4015       
4016       for (ii = 1 ; ii < num_poles ; ii++) {
4017         ii_index = ii % NumPoles;
4018         ii_index = (ii_index << 1) + ii_index;
4019         ii_minus = (ii - 1) % NumPoles;
4020         ii_minus = (ii_minus << 1) + ii_minus;
4021         inverse = FK[ii + Degree] - FK[ii];
4022         inverse = 1.0e0 / inverse;
4023         value = 0.0e0;
4024         factor = PA[ii_index] - PA[ii_minus];
4025         if (factor < 0) factor = - factor;
4026         value += factor; ii_index++; ii_minus++;
4027         factor = PA[ii_index] - PA[ii_minus];
4028         if (factor < 0) factor = - factor;
4029         value += factor; ii_index++; ii_minus++;
4030         factor = PA[ii_index] - PA[ii_minus];
4031         if (factor < 0) factor = - factor;
4032         value += factor;
4033         value *= inverse;
4034         if (max_derivative < value) max_derivative = value;
4035       }
4036     }
4037     break;
4038   }
4039   case 4 : {
4040     if (&Weights != NULL) {
4041       const Standard_Real * WG = &Weights(Weights.Lower());
4042       min_weights = WG[0];
4043       
4044       for (ii = 1 ; ii < NumPoles ; ii++) {
4045         W = WG[ii];
4046         if (W < min_weights) min_weights = W;
4047       }
4048       
4049       for (ii = 1 ; ii < num_poles ; ii++) {
4050         ii_index = ii % NumPoles;
4051         ii_inDim = ii_index << 2;
4052         ii_minus = (ii - 1) % NumPoles;
4053         ii_miDim = ii_minus << 2;
4054         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
4055         pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
4056         pa_ii_inDim_2 = PA[ii_inDim]; ii_inDim++;
4057         pa_ii_inDim_3 = PA[ii_inDim];
4058         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
4059         pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
4060         pa_ii_miDim_2 = PA[ii_miDim]; ii_miDim++;
4061         pa_ii_miDim_3 = PA[ii_miDim];
4062         wg_ii_index   = WG[ii_index];
4063         wg_ii_minus   = WG[ii_minus];
4064         inverse = FK[ii + Degree] - FK[ii];
4065         inverse = 1.0e0 / inverse;
4066         lower = ii - Deg1;
4067         if (lower < 0) lower = 0;
4068         upper = Deg2 + ii;
4069         if (upper > num_poles) upper = num_poles;
4070         
4071         for (jj = lower ; jj < upper ; jj++) {
4072           jj_index = jj % NumPoles;
4073           jj_index = jj_index << 2;
4074           value = 0.0e0;
4075           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4076                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4077           if (factor < 0) factor = - factor;
4078           value += factor; jj_index++;
4079           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4080                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4081           if (factor < 0) factor = - factor;
4082           value += factor; jj_index++;
4083           factor  = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4084                      ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4085           if (factor < 0) factor = - factor;
4086           value += factor; jj_index++;
4087           factor  = (((PA[jj_index] - pa_ii_inDim_3) * wg_ii_index) -
4088                      ((PA[jj_index] - pa_ii_miDim_3) * wg_ii_minus));
4089           if (factor < 0) factor = - factor;
4090           value += factor;
4091           value *= inverse;
4092           if (max_derivative < value) max_derivative = value;
4093         }
4094       }
4095       max_derivative /= min_weights;
4096     }
4097     else {
4098       
4099       for (ii = 1 ; ii < num_poles ; ii++) {
4100         ii_index = ii % NumPoles;
4101         ii_index = ii_index << 2;
4102         ii_minus = (ii - 1) % NumPoles;
4103         ii_minus = ii_minus << 2;
4104         inverse = FK[ii + Degree] - FK[ii];
4105         inverse = 1.0e0 / inverse;
4106         value = 0.0e0;
4107         factor = PA[ii_index] - PA[ii_minus];
4108         if (factor < 0) factor = - factor;
4109         value += factor; ii_index++; ii_minus++;
4110         factor = PA[ii_index] - PA[ii_minus];
4111         if (factor < 0) factor = - factor;
4112         value += factor; ii_index++; ii_minus++;
4113         factor = PA[ii_index] - PA[ii_minus];
4114         if (factor < 0) factor = - factor;
4115         value += factor; ii_index++; ii_minus++;
4116         factor = PA[ii_index] - PA[ii_minus];
4117         if (factor < 0) factor = - factor;
4118         value += factor;
4119         value *= inverse;
4120         if (max_derivative < value) max_derivative = value;
4121       }
4122     }
4123     break;
4124   }
4125     default : {
4126       Standard_Integer kk;
4127       if (&Weights != NULL) {
4128         const Standard_Real * WG = &Weights(Weights.Lower());
4129         min_weights = WG[0];
4130         
4131         for (ii = 1 ; ii < NumPoles ; ii++) {
4132           W = WG[ii];
4133           if (W < min_weights) min_weights = W;
4134         }
4135         
4136         for (ii = 1 ; ii < num_poles ; ii++) {
4137           ii_index  = ii % NumPoles;
4138           ii_inDim  = ii_index * ArrayDimension;
4139           ii_minus  = (ii - 1) % NumPoles;
4140           ii_miDim  = ii_minus * ArrayDimension;
4141           wg_ii_index   = WG[ii_index];
4142           wg_ii_minus   = WG[ii_minus];
4143           inverse = FK[ii + Degree] - FK[ii];
4144           inverse = 1.0e0 / inverse;
4145           lower = ii - Deg1;
4146           if (lower < 0) lower = 0;
4147           upper = Deg2 + ii;
4148           if (upper > num_poles) upper = num_poles;
4149           
4150           for (jj = lower ; jj < upper ; jj++) {
4151             jj_index = jj % NumPoles;
4152             jj_index *= ArrayDimension;
4153             value = 0.0e0;
4154             
4155             for (kk = 0 ; kk < ArrayDimension ; kk++) {
4156               factor  = (((PA[jj_index + kk] - PA[ii_inDim + kk]) * wg_ii_index) -
4157                          ((PA[jj_index + kk] - PA[ii_miDim + kk]) * wg_ii_minus));
4158               if (factor < 0) factor = - factor;
4159               value += factor;
4160             }
4161             value *= inverse;
4162             if (max_derivative < value) max_derivative = value;
4163           }
4164         }
4165         max_derivative /= min_weights;
4166       }
4167       else {
4168         
4169         for (ii = 1 ; ii < num_poles ; ii++) {
4170           ii_index  = ii % NumPoles;
4171           ii_index *= ArrayDimension;
4172           ii_minus  = (ii - 1) % NumPoles;
4173           ii_minus *= ArrayDimension;
4174           inverse = FK[ii + Degree] - FK[ii];
4175           inverse = 1.0e0 / inverse;
4176           value = 0.0e0;
4177           
4178           for (kk = 0 ; kk < ArrayDimension ; kk++) {
4179             factor = PA[ii_index + kk] - PA[ii_minus + kk];
4180             if (factor < 0) factor = - factor;
4181             value += factor;
4182           }
4183           value *= inverse;
4184           if (max_derivative < value) max_derivative = value;
4185         }
4186       }
4187     }
4188   }
4189   max_derivative *= Degree;
4190   if (max_derivative > RealSmall())
4191     UTolerance = Tolerance3D / max_derivative; 
4192   else
4193     UTolerance = Tolerance3D / RealSmall();
4194 }
4195
4196 //=======================================================================
4197 // function: FlatBezierKnots
4198 // purpose :
4199 //=======================================================================
4200
4201 // array of flat knots for bezier curve of maximum 25 degree
4202 static const Standard_Real knots[52] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
4203                                          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
4204 const Standard_Real& BSplCLib::FlatBezierKnots (const Standard_Integer Degree)
4205 {
4206   Standard_OutOfRange_Raise_if (Degree < 1 || Degree > MaxDegree() || MaxDegree() != 25,
4207     "Bezier curve degree greater than maximal supported");
4208
4209   return knots[25-Degree];
4210 }