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