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