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
   9 // under the terms of the GNU Lesser General Public 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   ReturnCode,
3114   Order = Degree + 1,
3115   FirstNonZeroBsplineIndex,
3116   LocalRequest = DerivativeRequest ;
3117   Standard_Real  *PResultArray,
3118   *WResultArray,
3119   *PolesArray,
3120   *WeightsArray,
3121   LocalParameter,
3122   Period,
3123   Inverse,
3124   Delta ;
3125   PolesArray = &Poles     ;
3126   WeightsArray = &Weights ;
3127   ExtrapModeArray = &ExtrapMode ;
3128   PResultArray = &PolesResults ;
3129   WResultArray = &WeightsResults ;
3130   LocalParameter = Parameter ;
3131   ExtrapolatingFlag[0] =
3132     ExtrapolatingFlag[1] = 0 ;
3133   //
3134   // check if we are extrapolating to a degree which is smaller than
3135   // the degree of the Bspline
3136   //
3137   if (PeriodicFlag) {
3138     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3139
3140     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3141       LocalParameter -= Period ;
3142     }
3143
3144     while (LocalParameter < FlatKnots(2)) {
3145       LocalParameter +=  Period ;
3146     }
3147   }
3148   if (Parameter < FlatKnots(2) &&
3149       LocalRequest < ExtrapModeArray[0] &&
3150       ExtrapModeArray[0] < Degree) {
3151     LocalRequest = ExtrapModeArray[0] ;
3152     LocalParameter = FlatKnots(2) ;
3153     ExtrapolatingFlag[0] = 1 ;
3154   }
3155   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3156       LocalRequest < ExtrapModeArray[1]  &&
3157       ExtrapModeArray[1] < Degree) {
3158     LocalRequest = ExtrapModeArray[1] ;
3159     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3160     ExtrapolatingFlag[1] = 1 ;
3161   }
3162   Delta = Parameter - LocalParameter ;
3163   if (LocalRequest >= Order) {
3164     LocalRequest = Degree ;
3165   }
3166   if (PeriodicFlag) {
3167     Modulus = FlatKnots.Length() - Degree -1 ;
3168   }
3169   else {
3170     Modulus = FlatKnots.Length() - Degree ;
3171   }
3172
3173   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3174   ErrorCode =
3175     BSplCLib::EvalBsplineBasis(1,
3176                                LocalRequest,
3177                                Order,
3178                                FlatKnots,
3179                                LocalParameter,
3180                                FirstNonZeroBsplineIndex,
3181                                BsplineBasis) ;
3182   if (ErrorCode != 0) {
3183     ReturnCode = 1 ;
3184     goto FINISH ;
3185   }
3186   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3187     Index = 0 ;
3188     Index2 = 0 ;
3189
3190     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3191       Index1 = FirstNonZeroBsplineIndex ;
3192
3193       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3194         PResultArray[Index + kk] = 0.0e0 ;
3195       }
3196       WResultArray[Index] = 0.0e0 ;
3197
3198       for (jj = 1  ; jj <= Order ; jj++) {
3199
3200         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3201           PResultArray[Index + kk] +=
3202             PolesArray[(Index1-1) * ArrayDimension + kk]
3203               * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3204         }
3205         WResultArray[Index2]  += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3206
3207         Index1 = Index1 % Modulus ;
3208         Index1 += 1 ;
3209       }
3210       Index += ArrayDimension ;
3211       Index2 += 1 ;
3212     }
3213   }
3214   else {
3215     //
3216     //  store Taylor expansion in LocalRealArray
3217     //
3218     NewRequest = DerivativeRequest ;
3219     if (NewRequest > Degree) {
3220       NewRequest = Degree ;
3221     }
3222     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3223     Index = 0 ;
3224     Inverse = 1.0e0 ;
3225
3226     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3227       Index1 = FirstNonZeroBsplineIndex ;
3228
3229       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3230         LocalRealArray[Index + kk] = 0.0e0 ;
3231       }
3232
3233       for (jj = 1  ; jj <= Order ; jj++) {
3234
3235         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3236           LocalRealArray[Index + kk] +=
3237             PolesArray[(Index1-1)*ArrayDimension + kk] *
3238               WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3239         }
3240         Index1 = Index1 % Modulus ;
3241         Index1 += 1 ;
3242       }
3243
3244       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3245         LocalRealArray[Index + kk] *= Inverse ;
3246       }
3247       Index += ArrayDimension ;
3248       Inverse /= (Standard_Real) ii ;
3249     }
3250     PLib::EvalPolynomial(Delta,
3251                          NewRequest,
3252                          Degree,
3253                          ArrayDimension,
3254                          LocalRealArray[0],
3255                          PolesResults) ;
3256     Index = 0 ;
3257     Inverse = 1.0e0 ;
3258
3259     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3260       Index1 = FirstNonZeroBsplineIndex ;
3261       LocalRealArray[Index] = 0.0e0 ;
3262
3263       for (jj = 1  ; jj <= Order ; jj++) {
3264         LocalRealArray[Index] +=
3265           WeightsArray[Index1-1] * BsplineBasis(ii,jj) ;
3266         Index1 = Index1 % Modulus ;
3267         Index1 += 1 ;
3268       }
3269       LocalRealArray[Index + kk] *= Inverse ;
3270       Index += 1 ;
3271       Inverse /= (Standard_Real) ii ;
3272     }
3273     PLib::EvalPolynomial(Delta,
3274                          NewRequest,
3275                          Degree,
3276                          1,
3277                          LocalRealArray[0],
3278                          WeightsResults) ;
3279   }
3280   FINISH : ;
3281 }
3282
3283 //=======================================================================
3284 //function : Evaluates a Bspline function : uses the ExtrapMode
3285 //purpose  : the function is extrapolated using the Taylor expansion
3286 //           of degree ExtrapMode[0] to the left and the Taylor
3287 //           expansion of degree ExtrapMode[1] to the right
3288 // WARNING : the array Results is supposed to have at least
3289 // (DerivativeRequest + 1) * ArrayDimension slots and the
3290 //
3291 //=======================================================================
3292
3293 void  BSplCLib::Eval
3294 (const Standard_Real                   Parameter,
3295  const Standard_Boolean                PeriodicFlag,
3296  const Standard_Integer                DerivativeRequest,
3297  Standard_Integer&                     ExtrapMode,
3298  const Standard_Integer                Degree,
3299  const  TColStd_Array1OfReal&          FlatKnots,
3300  const Standard_Integer                ArrayDimension,
3301  Standard_Real&                        Poles,
3302  Standard_Real&                        Results)
3303 {
3304   Standard_Integer ii,
3305   jj,
3306   kk,
3307   Index,
3308   Index1,
3309   *ExtrapModeArray,
3310   Modulus,
3311   NewRequest,
3312   ExtrapolatingFlag[2],
3313   ErrorCode,
3314   ReturnCode,
3315   Order = Degree + 1,
3316   FirstNonZeroBsplineIndex,
3317   LocalRequest = DerivativeRequest ;
3318
3319   Standard_Real  *ResultArray,
3320   *PolesArray,
3321   LocalParameter,
3322   Period,
3323   Inverse,
3324   Delta ;
3325
3326   PolesArray = &Poles ;
3327   ExtrapModeArray = &ExtrapMode ;
3328   ResultArray = &Results ;
3329   LocalParameter = Parameter ;
3330   ExtrapolatingFlag[0] =
3331     ExtrapolatingFlag[1] = 0 ;
3332   //
3333   // check if we are extrapolating to a degree which is smaller than
3334   // the degree of the Bspline
3335   //
3336   if (PeriodicFlag) {
3337     Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ;
3338
3339     while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) {
3340       LocalParameter -= Period ;
3341     }
3342
3343     while (LocalParameter < FlatKnots(2)) {
3344       LocalParameter +=  Period ;
3345     }
3346   }
3347   if (Parameter < FlatKnots(2) &&
3348       LocalRequest < ExtrapModeArray[0] &&
3349       ExtrapModeArray[0] < Degree) {
3350     LocalRequest = ExtrapModeArray[0] ;
3351     LocalParameter = FlatKnots(2) ;
3352     ExtrapolatingFlag[0] = 1 ;
3353   }
3354   if (Parameter > FlatKnots(FlatKnots.Upper()-1) &&
3355       LocalRequest < ExtrapModeArray[1]  &&
3356       ExtrapModeArray[1] < Degree) {
3357     LocalRequest = ExtrapModeArray[1] ;
3358     LocalParameter = FlatKnots(FlatKnots.Upper()-1) ;
3359     ExtrapolatingFlag[1] = 1 ;
3360   }
3361   Delta = Parameter - LocalParameter ;
3362   if (LocalRequest >= Order) {
3363     LocalRequest = Degree ;
3364   }
3365
3366   if (PeriodicFlag) {
3367     Modulus = FlatKnots.Length() - Degree -1 ;
3368   }
3369   else {
3370     Modulus = FlatKnots.Length() - Degree ;
3371   }
3372
3373   BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order);
3374
3375   ErrorCode =
3376     BSplCLib::EvalBsplineBasis(1,
3377                                LocalRequest,
3378                                Order,
3379                                FlatKnots,
3380                                LocalParameter,
3381                                FirstNonZeroBsplineIndex,
3382                                BsplineBasis);
3383   if (ErrorCode != 0) {
3384     ReturnCode = 1 ;
3385     goto FINISH ;
3386   }
3387   if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) {
3388     Index = 0 ;
3389
3390     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3391       Index1 = FirstNonZeroBsplineIndex ;
3392
3393       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3394         ResultArray[Index + kk] = 0.0e0 ;
3395       }
3396
3397       for (jj = 1  ; jj <= Order ; jj++) {
3398
3399         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3400           ResultArray[Index + kk] +=
3401             PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3402         }
3403         Index1 = Index1 % Modulus ;
3404         Index1 += 1 ;
3405       }
3406       Index += ArrayDimension ;
3407     }
3408   }
3409   else {
3410     //
3411     //  store Taylor expansion in LocalRealArray
3412     //
3413     NewRequest = DerivativeRequest ;
3414     if (NewRequest > Degree) {
3415       NewRequest = Degree ;
3416     }
3417     NCollection_LocalArray<Standard_Real> LocalRealArray((LocalRequest + 1)*ArrayDimension);
3418
3419     Index = 0 ;
3420     Inverse = 1.0e0 ;
3421
3422     for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) {
3423       Index1 = FirstNonZeroBsplineIndex ;
3424
3425       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3426         LocalRealArray[Index + kk] = 0.0e0 ;
3427       }
3428
3429       for (jj = 1  ; jj <= Order ; jj++) {
3430
3431         for (kk = 0 ; kk < ArrayDimension ; kk++) {
3432           LocalRealArray[Index + kk] +=
3433             PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ;
3434         }
3435         Index1 = Index1 % Modulus ;
3436         Index1 += 1 ;
3437       }
3438
3439       for (kk = 0 ; kk < ArrayDimension ; kk++) {
3440         LocalRealArray[Index + kk] *= Inverse ;
3441       }
3442       Index += ArrayDimension ;
3443       Inverse /= (Standard_Real) ii ;
3444     }
3445     PLib::EvalPolynomial(Delta,
3446                          NewRequest,
3447                          Degree,
3448                          ArrayDimension,
3449                          LocalRealArray[0],
3450                          Results) ;
3451   }
3452   FINISH : ;
3453 }
3454
3455 //=======================================================================
3456 //function : TangExtendToConstraint
3457 //purpose  : Extends a Bspline function using the tangency map
3458 // WARNING :
3459 //
3460 //
3461 //=======================================================================
3462
3463 void  BSplCLib::TangExtendToConstraint
3464 (const  TColStd_Array1OfReal&          FlatKnots,
3465  const Standard_Real                   C1Coefficient,
3466  const Standard_Integer                NumPoles,
3467  Standard_Real&                        Poles,
3468  const Standard_Integer                CDimension,
3469  const Standard_Integer                CDegree,
3470  const  TColStd_Array1OfReal&          ConstraintPoint,
3471  const Standard_Integer                Continuity,
3472  const Standard_Boolean                After,
3473  Standard_Integer&                     NbPolesResult,
3474  Standard_Integer&                     NbKnotsResult,
3475  Standard_Real&                        KnotsResult,
3476  Standard_Real&                        PolesResult)
3477 {
3478 #if DEB
3479   if (CDegree<Continuity+1) {
3480     cout<<"The BSpline degree must be greater than the order of continuity"<<endl;
3481   }
3482 #endif
3483   Standard_Real * Padr = &Poles ;
3484   Standard_Real * KRadr = &KnotsResult ;
3485   Standard_Real * PRadr = &PolesResult ;
3486
3487 ////////////////////////////////////////////////////////////////////////
3488 //
3489 //    1. calculation of extension nD
3490 //
3491 ////////////////////////////////////////////////////////////////////////
3492
3493 //  Hermite matrix
3494   Standard_Integer Csize = Continuity + 2;
3495   math_Matrix  MatCoefs(1,Csize, 1,Csize);
3496   if (After) {
3497     PLib::HermiteCoefficients(0, 1,           // Limits
3498                               Continuity, 0,  // Orders of constraints
3499                               MatCoefs);
3500   }
3501   else {
3502     PLib::HermiteCoefficients(0, 1,           // Limits
3503                               0, Continuity,  // Orders of constraints
3504                               MatCoefs);
3505   }
3506
3507
3508 //  position at the node of connection
3509   Standard_Real Tbord ;
3510   if (After) {
3511     Tbord = FlatKnots(FlatKnots.Upper()-CDegree);
3512   }
3513   else {
3514     Tbord = FlatKnots(FlatKnots.Lower()+CDegree);
3515   }
3516   Standard_Boolean periodic_flag = Standard_False ;
3517   Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1);
3518   extrap_mode[0] = extrap_mode[1] = CDegree;
3519   TColStd_Array1OfReal  EvalBS(1, CDimension * (derivative_request+1)) ;
3520   Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
3521   BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
3522                   CDegree,FlatKnots,CDimension,Poles,*Eadr);
3523
3524 //  norm of the tangent at the node of connection
3525   math_Vector Tgte(1,CDimension);
3526
3527   for (ipos=1;ipos<=CDimension;ipos++) {
3528     Tgte(ipos) = EvalBS(ipos+CDimension);
3529   }
3530   Standard_Real L1=Tgte.Norm();
3531
3532
3533 //  matrix of constraints
3534   math_Matrix Contraintes(1,Csize,1,CDimension);
3535   if (After) {
3536
3537     for (ipos=1;ipos<=CDimension;ipos++) {
3538       Contraintes(1,ipos) = EvalBS(ipos);
3539       Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3540       if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3541       if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3542       Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos);
3543     }
3544   }
3545   else {
3546
3547     for (ipos=1;ipos<=CDimension;ipos++) {
3548       Contraintes(1,ipos) = ConstraintPoint(ipos);
3549       Contraintes(2,ipos) = EvalBS(ipos);
3550       if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension);
3551       if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2);
3552       if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3);
3553     }
3554   }
3555
3556 //  calculate the coefficients of extension
3557   Standard_Integer ii, jj, kk;
3558   TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension);
3559   ExtraCoeffs.Init(0.);
3560
3561   for (ii=1; ii<=Csize; ii++) {
3562
3563     for (jj=1; jj<=Csize; jj++) {
3564
3565       for (kk=1; kk<=CDimension; kk++) {
3566         ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk);
3567       }
3568     }
3569   }
3570
3571 //  calculate the poles of extension
3572   TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension);
3573   Standard_Real * EPadr = &ExtrapPoles(1) ;
3574   PLib::CoefficientsPoles(CDimension,
3575                           ExtraCoeffs,  PLib::NoWeights(),
3576                           ExtrapPoles,  PLib::NoWeights());
3577
3578 //  calculate the nodes of extension with multiplicities
3579   TColStd_Array1OfReal ExtrapNoeuds(1,2);
3580   ExtrapNoeuds(1) = 0.;
3581   ExtrapNoeuds(2) = 1.;
3582   TColStd_Array1OfInteger ExtrapMults(1,2);
3583   ExtrapMults(1) = Csize;
3584   ExtrapMults(2) = Csize;
3585
3586 // flat nodes of extension
3587   TColStd_Array1OfReal FK2(1, Csize*2);
3588   BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2);
3589
3590 //  norm of the tangent at the connection point
3591   if (After) {
3592     BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0],
3593                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3594   }
3595   else {
3596     BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0],
3597                   Csize-1,FK2,CDimension,*EPadr,*Eadr);
3598   }
3599
3600   for (ipos=1;ipos<=CDimension;ipos++) {
3601     Tgte(ipos) = EvalBS(ipos+CDimension);
3602   }
3603   Standard_Real L2 = Tgte.Norm();
3604
3605 //  harmonisation of degrees
3606   TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension);
3607   TColStd_Array1OfReal NewK2(1, 2);
3608   TColStd_Array1OfInteger NewM2(1, 2);
3609   if (Csize-1<CDegree) {
3610     BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension,
3611                              ExtrapPoles,ExtrapNoeuds,ExtrapMults,
3612                              NewP2,NewK2,NewM2);
3613   }
3614   else {
3615     NewP2 = ExtrapPoles;
3616     NewK2 = ExtrapNoeuds;
3617     NewM2 = ExtrapMults;
3618   }
3619
3620 //  flat nodes of extension after harmonization of degrees
3621   TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2);
3622   BSplCLib::KnotSequence(NewK2,NewM2,NewFK2);
3623
3624
3625 ////////////////////////////////////////////////////////////////////////
3626 //
3627 //    2.  concatenation C0
3628 //
3629 ////////////////////////////////////////////////////////////////////////
3630
3631 //  ratio of reparametrization
3632   Standard_Real Ratio=1, Delta;
3633   if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) {
3634     Ratio = L2 / L1;
3635   }
3636   if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1;
3637
3638   if (After) {
3639 //    do not touch the first BSpline
3640     Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper());
3641   }
3642   else {
3643 //    do not touch the second BSpline
3644     Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower());
3645   }
3646
3647 //  result of the concatenation
3648   Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1;
3649   Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1);
3650   TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension);
3651   TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2);
3652
3653 //  poles
3654   Standard_Integer indNP, indP, indEP;
3655   if (After) {
3656
3657     for (ii=1;  ii<=NbP1+NbP2-1; ii++) {
3658
3659       for (jj=1;  jj<=CDimension; jj++) {
3660         indNP = (ii-1)*CDimension+jj;
3661         indP = (ii-1)*CDimension+jj-1;
3662         indEP = (ii-NbP1)*CDimension+jj;
3663         if (ii<NbP1) NewPoles(indNP) =  Padr[indP];
3664         else NewPoles(indNP) = NewP2(indEP);
3665       }
3666     }
3667   }
3668   else {
3669
3670     for (ii=1;  ii<=NbP1+NbP2-1; ii++) {
3671
3672       for (jj=1;  jj<=CDimension; jj++) {
3673         indNP = (ii-1)*CDimension+jj;
3674         indEP = (ii-1)*CDimension+jj;
3675         indP = (ii-NbP2)*CDimension+jj-1;
3676         if (ii<NbP2) NewPoles(indNP) =  NewP2(indEP);
3677         else NewPoles(indNP) = Padr[indP];
3678       }
3679     }
3680   }
3681
3682 //  flat nodes
3683   if (After) {
3684 //    start with the nodes of the initial surface
3685
3686     for (ii=1; ii<NbK1; ii++) {
3687       NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1);
3688     }
3689 //    continue with the reparameterized nodes of the extension
3690
3691     for (ii=1; ii<=NbK2-CDegree-1; ii++) {
3692       NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta;
3693     }
3694   }
3695   else {
3696 //    start with the reparameterized nodes of the extension
3697
3698     for (ii=1; ii<NbK2-CDegree; ii++) {
3699       NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta;
3700     }
3701 //    continue with the nodes of the initial surface
3702
3703     for (ii=2; ii<=NbK1; ii++) {
3704       NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1);
3705     }
3706   }
3707
3708
3709 ////////////////////////////////////////////////////////////////////////
3710 //
3711 //    3.  reduction of multiplicite at the node of connection
3712 //
3713 ////////////////////////////////////////////////////////////////////////
3714
3715 //  number of separate nodes
3716   Standard_Integer KLength = 1;
3717
3718   for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3719     if (NewFlats(ii) != NewFlats(ii-1)) KLength++;
3720   }
3721
3722 //  flat nodes --> nodes + multiplicities
3723   TColStd_Array1OfReal NewKnots (1, KLength);
3724   TColStd_Array1OfInteger NewMults (1, KLength);
3725   NewMults.Init(1);
3726   jj = 1;
3727   NewKnots(jj) = NewFlats(1);
3728
3729   for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) {
3730     if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++;
3731     else {
3732       jj++;
3733       NewKnots(jj) = NewFlats(ii);
3734     }
3735   }
3736
3737 //  reduction of multiplicity at the second or the last but one node
3738   Standard_Integer Index = 2, M = CDegree;
3739   if (After) Index = KLength-1;
3740   TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension);
3741   TColStd_Array1OfReal ResultKnots (1, KLength);
3742   TColStd_Array1OfInteger ResultMults (1, KLength);
3743   Standard_Real Tol = 1.e-6;
3744   Standard_Boolean Ok = Standard_True;
3745
3746   while ( (M>CDegree-Continuity) && Ok) {
3747     Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension,
3748                     NewPoles, NewKnots, NewMults,
3749                     ResultPoles, ResultKnots, ResultMults, Tol);
3750     if (Ok) M--;
3751   }
3752
3753   if (M == CDegree) {
3754 //    number of poles of the concatenation
3755     NbPolesResult = NbP1 + NbP2 - 1;
3756 //    the poles of the concatenation
3757     Standard_Integer PLength = NbPolesResult*CDimension;
3758
3759     for (jj=1; jj<=PLength; jj++) {
3760       PRadr[jj-1] = NewPoles(jj);
3761     }
3762
3763 //    flat nodes of the concatenation
3764     Standard_Integer ideb = 0;
3765
3766     for (jj=0; jj<NewKnots.Length(); jj++) {
3767       for (ii=0; ii<NewMults(jj+1); ii++) {
3768         KRadr[ideb+ii] = NewKnots(jj+1);
3769       }
3770       ideb += NewMults(jj+1);
3771     }
3772     NbKnotsResult = ideb;
3773   }
3774
3775   else {
3776 //    number of poles of the result
3777     NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M;
3778 //    the poles of the result
3779     Standard_Integer PLength = NbPolesResult*CDimension;
3780
3781     for (jj=0; jj<PLength; jj++) {
3782       PRadr[jj] = ResultPoles(jj+1);
3783     }
3784
3785 //    flat nodes of the result
3786     Standard_Integer ideb = 0;
3787
3788     for (jj=0; jj<ResultKnots.Length(); jj++) {
3789       for (ii=0; ii<ResultMults(jj+1); ii++) {
3790         KRadr[ideb+ii] = ResultKnots(jj+1);
3791       }
3792       ideb += ResultMults(jj+1);
3793     }
3794     NbKnotsResult = ideb;
3795   }
3796 }
3797
3798 //=======================================================================
3799 //function : Resolution
3800 //purpose  :
3801 //                           d
3802 //  Let C(t) = SUM      Ci Bi(t)  a Bspline curve of degree d
3803 //            i = 1,n
3804 //  with nodes tj for j = 1,n+d+1
3805 //
3806 //
3807 //         '                    C1 - Ci-1   d-1
3808 //  Then C (t) = SUM     d *  ---------  Bi (t)
3809 //                i = 2,n      ti+d - ti
3810 //
3811 //                          d-1
3812 //  for the base of BSpline  Bi  (t) of degree d-1.
3813 //
3814 //  Consequently the upper bound of the norm of the derivative from C is :
3815 //
3816 //
3817 //                        |  Ci - Ci-1  |
3818 //          d *   Max     |  ---------  |
3819 //              i = 2,n |  ti+d - ti  |
3820 //
3821 //                                      N(t)
3822 //  In the rational case set    C(t) = -----
3823 //                                      D(t)
3824 //
3825 //
3826 //  D(t) =  SUM    Di Bi(t)
3827 //        i=1,n
3828 //
3829 //  N(t) =  SUM   Di * Ci Bi(t)
3830 //          i =1,n
3831 //
3832 //          N'(t)  -    D'(t) C(t)
3833 //   C'(t) = -----------------------
3834 //                   D(t)
3835 //
3836 //
3837 //   N'(t) - D'(t) C(t) =
3838 //
3839 //                     Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t))    d-1
3840 //      SUM   d *   ---------------------------------------- * Bi  (t)  =
3841 //        i=2,n                   ti+d   - ti
3842 //
3843 //
3844 //                   Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)                d-1
3845 // SUM   SUM     d * -----------------------------------  * Betaj(t) * Bi  (t)
3846 //i=2,n j=1,n               ti+d  - ti
3847 //
3848 //
3849 //
3850 //                 Dj Bj(t)
3851 //    Betaj(t) =   --------
3852 //                 D(t)
3853 //
3854 //  Betaj(t) form a partition >= 0 of the entity with support
3855 //  tj, tj+d+1. Consequently if Rj = {j-d, ....,  j+d+d+1}
3856 //  obtain an upper bound of the derivative of C by taking :
3857 //
3858 //
3859 //
3860 //
3861 //
3862 //
3863 //                         Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj)
3864 //   Max   Max       d  *  -----------------------------------
3865 // j=1,n  i dans Rj                   ti+d  - ti
3866 //
3867 //  --------------------------------------------------------
3868 //
3869 //               Min    Di
3870 //              i =1,n
3871 //
3872 //
3873 //=======================================================================
3874
3875 void BSplCLib::Resolution(      Standard_Real&        Poles,
3876                           const Standard_Integer      ArrayDimension,
3877                           const Standard_Integer      NumPoles,
3878                           const TColStd_Array1OfReal& Weights,
3879                           const TColStd_Array1OfReal& FlatKnots,
3880                           const Standard_Integer      Degree,
3881                           const Standard_Real         Tolerance3D,
3882                           Standard_Real&              UTolerance)
3883 {
3884   Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim;
3885   Standard_Integer lower,upper,ii_minus,jj,ii_miDim;
3886   Standard_Integer Deg1 = Degree + 1;
3887   Standard_Integer Deg2 = (Degree << 1) + 1;
3888   Standard_Real value,factor,W,min_weights,inverse;
3889   Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3;
3890   Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3;
3891   Standard_Real wg_ii_index, wg_ii_minus;
3892   Standard_Real *PA,max_derivative;
3893   const Standard_Real * FK = &FlatKnots(FlatKnots.Lower());
3894   PA = &Poles;
3895   max_derivative = 0.0e0;
3896   num_poles = FlatKnots.Length() - Deg1;
3897   switch (ArrayDimension) {
3898   case 2 : {
3899     if (&Weights != NULL) {
3900       const Standard_Real * WG = &Weights(Weights.Lower());
3901       min_weights = WG[0];
3902
3903       for (ii = 1 ; ii < NumPoles ; ii++) {
3904         W = WG[ii];
3905         if (W < min_weights) min_weights = W;
3906       }
3907
3908       for (ii = 1 ; ii < num_poles ; ii++) {
3909         ii_index = ii % NumPoles;
3910         ii_inDim = ii_index << 1;
3911         ii_minus = (ii - 1) % NumPoles;
3912         ii_miDim = ii_minus << 1;
3913         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3914         pa_ii_inDim_1 = PA[ii_inDim];
3915         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3916         pa_ii_miDim_1 = PA[ii_miDim];
3917         wg_ii_index   = WG[ii_index];
3918         wg_ii_minus   = WG[ii_minus];
3919         inverse = FK[ii + Degree] - FK[ii];
3920         inverse = 1.0e0 / inverse;
3921         lower = ii - Deg1;
3922         if (lower < 0) lower = 0;
3923         upper = Deg2 + ii;
3924         if (upper > num_poles) upper = num_poles;
3925
3926         for (jj = lower ; jj < upper ; jj++) {
3927           jj_index = jj % NumPoles;
3928           jj_index = jj_index << 1;
3929           value = 0.0e0;
3930           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
3931                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
3932           if (factor < 0) factor = - factor;
3933           value += factor; jj_index++;
3934           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
3935                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
3936           if (factor < 0) factor = - factor;
3937           value += factor;
3938           value *= inverse;
3939           if (max_derivative < value) max_derivative = value;
3940         }
3941       }
3942       max_derivative /= min_weights;
3943     }
3944     else {
3945
3946       for (ii = 1 ; ii < num_poles ; ii++) {
3947         ii_index = ii % NumPoles;
3948         ii_index = ii_index << 1;
3949         ii_minus = (ii - 1) % NumPoles;
3950         ii_minus = ii_minus << 1;
3951         inverse = FK[ii + Degree] - FK[ii];
3952         inverse = 1.0e0 / inverse;
3953         value = 0.0e0;
3954         factor = PA[ii_index] - PA[ii_minus];
3955         if (factor < 0) factor = - factor;
3956         value += factor; ii_index++; ii_minus++;
3957         factor = PA[ii_index] - PA[ii_minus];
3958         if (factor < 0) factor = - factor;
3959         value += factor;
3960         value *= inverse;
3961         if (max_derivative < value) max_derivative = value;
3962       }
3963     }
3964     break;
3965   }
3966   case 3 : {
3967     if (&Weights != NULL) {
3968       const Standard_Real * WG = &Weights(Weights.Lower());
3969       min_weights = WG[0];
3970
3971       for (ii = 1 ; ii < NumPoles ; ii++) {
3972         W = WG[ii];
3973         if (W < min_weights) min_weights = W;
3974       }
3975
3976       for (ii = 1 ; ii < num_poles ; ii++) {
3977         ii_index = ii % NumPoles;
3978         ii_inDim = (ii_index << 1) + ii_index;
3979         ii_minus = (ii - 1) % NumPoles;
3980         ii_miDim = (ii_minus << 1) + ii_minus;
3981         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
3982         pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
3983         pa_ii_inDim_2 = PA[ii_inDim];
3984         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
3985         pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
3986         pa_ii_miDim_2 = PA[ii_miDim];
3987         wg_ii_index   = WG[ii_index];
3988         wg_ii_minus   = WG[ii_minus];
3989         inverse = FK[ii + Degree] - FK[ii];
3990         inverse = 1.0e0 / inverse;
3991         lower = ii - Deg1;
3992         if (lower < 0) lower = 0;
3993         upper = Deg2 + ii;
3994         if (upper > num_poles) upper = num_poles;
3995
3996         for (jj = lower ; jj < upper ; jj++) {
3997           jj_index = jj % NumPoles;
3998           jj_index = (jj_index << 1) + jj_index;
3999           value = 0.0e0;
4000           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4001                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4002           if (factor < 0) factor = - factor;
4003           value += factor; jj_index++;
4004           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4005                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4006           if (factor < 0) factor = - factor;
4007           value += factor; jj_index++;
4008           factor  = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4009                      ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4010           if (factor < 0) factor = - factor;
4011           value += factor;
4012           value *= inverse;
4013           if (max_derivative < value) max_derivative = value;
4014         }
4015       }
4016       max_derivative /= min_weights;
4017     }
4018     else {
4019
4020       for (ii = 1 ; ii < num_poles ; ii++) {
4021         ii_index = ii % NumPoles;
4022         ii_index = (ii_index << 1) + ii_index;
4023         ii_minus = (ii - 1) % NumPoles;
4024         ii_minus = (ii_minus << 1) + ii_minus;
4025         inverse = FK[ii + Degree] - FK[ii];
4026         inverse = 1.0e0 / inverse;
4027         value = 0.0e0;
4028         factor = PA[ii_index] - PA[ii_minus];
4029         if (factor < 0) factor = - factor;
4030         value += factor; ii_index++; ii_minus++;
4031         factor = PA[ii_index] - PA[ii_minus];
4032         if (factor < 0) factor = - factor;
4033         value += factor; ii_index++; ii_minus++;
4034         factor = PA[ii_index] - PA[ii_minus];
4035         if (factor < 0) factor = - factor;
4036         value += factor;
4037         value *= inverse;
4038         if (max_derivative < value) max_derivative = value;
4039       }
4040     }
4041     break;
4042   }
4043   case 4 : {
4044     if (&Weights != NULL) {
4045       const Standard_Real * WG = &Weights(Weights.Lower());
4046       min_weights = WG[0];
4047
4048       for (ii = 1 ; ii < NumPoles ; ii++) {
4049         W = WG[ii];
4050         if (W < min_weights) min_weights = W;
4051       }
4052
4053       for (ii = 1 ; ii < num_poles ; ii++) {
4054         ii_index = ii % NumPoles;
4055         ii_inDim = ii_index << 2;
4056         ii_minus = (ii - 1) % NumPoles;
4057         ii_miDim = ii_minus << 2;
4058         pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++;
4059         pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++;
4060         pa_ii_inDim_2 = PA[ii_inDim]; ii_inDim++;
4061         pa_ii_inDim_3 = PA[ii_inDim];
4062         pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++;
4063         pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++;
4064         pa_ii_miDim_2 = PA[ii_miDim]; ii_miDim++;
4065         pa_ii_miDim_3 = PA[ii_miDim];
4066         wg_ii_index   = WG[ii_index];
4067         wg_ii_minus   = WG[ii_minus];
4068         inverse = FK[ii + Degree] - FK[ii];
4069         inverse = 1.0e0 / inverse;
4070         lower = ii - Deg1;
4071         if (lower < 0) lower = 0;
4072         upper = Deg2 + ii;
4073         if (upper > num_poles) upper = num_poles;
4074
4075         for (jj = lower ; jj < upper ; jj++) {
4076           jj_index = jj % NumPoles;
4077           jj_index = jj_index << 2;
4078           value = 0.0e0;
4079           factor  = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) -
4080                      ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus));
4081           if (factor < 0) factor = - factor;
4082           value += factor; jj_index++;
4083           factor  = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) -
4084                      ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus));
4085           if (factor < 0) factor = - factor;
4086           value += factor; jj_index++;
4087           factor  = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) -
4088                      ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus));
4089           if (factor < 0) factor = - factor;
4090           value += factor; jj_index++;
4091           factor  = (((PA[jj_index] - pa_ii_inDim_3) * wg_ii_index) -
4092                      ((PA[jj_index] - pa_ii_miDim_3) * wg_ii_minus));
4093           if (factor < 0) factor = - factor;
4094           value += factor;
4095           value *= inverse;
4096           if (max_derivative < value) max_derivative = value;
4097         }
4098       }
4099       max_derivative /= min_weights;
4100     }
4101     else {
4102
4103       for (ii = 1 ; ii < num_poles ; ii++) {
4104         ii_index = ii % NumPoles;
4105         ii_index = ii_index << 2;
4106         ii_minus = (ii - 1) % NumPoles;
4107         ii_minus = ii_minus << 2;
4108         inverse = FK[ii + Degree] - FK[ii];
4109         inverse = 1.0e0 / inverse;
4110         value = 0.0e0;
4111         factor = PA[ii_index] - PA[ii_minus];
4112         if (factor < 0) factor = - factor;
4113         value += factor; ii_index++; ii_minus++;
4114         factor = PA[ii_index] - PA[ii_minus];
4115         if (factor < 0) factor = - factor;
4116         value += factor; ii_index++; ii_minus++;
4117         factor = PA[ii_index] - PA[ii_minus];
4118         if (factor < 0) factor = - factor;
4119         value += factor; ii_index++; ii_minus++;
4120         factor = PA[ii_index] - PA[ii_minus];
4121         if (factor < 0) factor = - factor;
4122         value += factor;
4123         value *= inverse;
4124         if (max_derivative < value) max_derivative = value;
4125       }
4126     }
4127     break;
4128   }
4129     default : {
4130       Standard_Integer kk;
4131       if (&Weights != NULL) {
4132         const Standard_Real * WG = &Weights(Weights.Lower());
4133         min_weights = WG[0];
4134
4135         for (ii = 1 ; ii < NumPoles ; ii++) {
4136           W = WG[ii];
4137           if (W < min_weights) min_weights = W;
4138         }
4139
4140         for (ii = 1 ; ii < num_poles ; ii++) {
4141           ii_index  = ii % NumPoles;
4142           ii_inDim  = ii_index * ArrayDimension;
4143           ii_minus  = (ii - 1) % NumPoles;
4144           ii_miDim  = ii_minus * ArrayDimension;
4145           wg_ii_index   = WG[ii_index];
4146           wg_ii_minus   = WG[ii_minus];
4147           inverse = FK[ii + Degree] - FK[ii];
4148           inverse = 1.0e0 / inverse;
4149           lower = ii - Deg1;
4150           if (lower < 0) lower = 0;
4151           upper = Deg2 + ii;
4152           if (upper > num_poles) upper = num_poles;
4153
4154           for (jj = lower ; jj < upper ; jj++) {
4155             jj_index = jj % NumPoles;
4156             jj_index *= ArrayDimension;
4157             value = 0.0e0;
4158
4159             for (kk = 0 ; kk < ArrayDimension ; kk++) {
4160               factor  = (((PA[jj_index + kk] - PA[ii_inDim + kk]) * wg_ii_index) -
4161                          ((PA[jj_index + kk] - PA[ii_miDim + kk]) * wg_ii_minus));
4162               if (factor < 0) factor = - factor;
4163               value += factor;
4164             }
4165             value *= inverse;
4166             if (max_derivative < value) max_derivative = value;
4167           }
4168         }
4169         max_derivative /= min_weights;
4170       }
4171       else {
4172
4173         for (ii = 1 ; ii < num_poles ; ii++) {
4174           ii_index  = ii % NumPoles;
4175           ii_index *= ArrayDimension;
4176           ii_minus  = (ii - 1) % NumPoles;
4177           ii_minus *= ArrayDimension;
4178           inverse = FK[ii + Degree] - FK[ii];
4179           inverse = 1.0e0 / inverse;
4180           value = 0.0e0;
4181
4182           for (kk = 0 ; kk < ArrayDimension ; kk++) {
4183             factor = PA[ii_index + kk] - PA[ii_minus + kk];
4184             if (factor < 0) factor = - factor;
4185             value += factor;
4186           }
4187           value *= inverse;
4188           if (max_derivative < value) max_derivative = value;
4189         }
4190       }
4191     }
4192   }
4193   max_derivative *= Degree;
4194   if (max_derivative > RealSmall())
4195     UTolerance = Tolerance3D / max_derivative;
4196   else
4197     UTolerance = Tolerance3D / RealSmall();
4198 }
4199
4200 //=======================================================================
4201 // function: FlatBezierKnots
4202 // purpose :
4203 //=======================================================================
4204
4205 // array of flat knots for bezier curve of maximum 25 degree
4206 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,
4207                                          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 };
4208 const Standard_Real& BSplCLib::FlatBezierKnots (const Standard_Integer Degree)
4209 {
4210   Standard_OutOfRange_Raise_if (Degree < 1 || Degree > MaxDegree() || MaxDegree() != 25,
4211     "Bezier curve degree greater than maximal supported");
4212
4213   return knots[25-Degree];
4214 }