src/BSplCLib/BSplCLib.cxx

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