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