src/Geom/Geom_BSplineCurve.cxx

   1 // Created on: 1993-03-09
   2 // Created by: JCV
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //Avril 1991 : constructeurs + methodes de lecture.
  18 //Mai 1991   : revue des specifs + debut de realisation des classes tool =>
  19 //             implementation des methodes Set et calcul du point courant.
  20 //Juillet 1991 : voir egalement File Geom_BSplineCurve_1.cxx
  21 //Juin    1992 : mise a plat des valeurs nodales - amelioration des
  22 //               performances sur calcul du point courant
  23
  24 //RLE Aug 1993  Remove Swaps, Remove typedefs, Update BSplCLib
  25 //              debug periodic, IncreaseDegree
  26 // 21-Mar-95 : xab implemented cache
  27 // 14-Mar-96 : xab implemented MovePointAndTangent
  28 // 13-Oct-96 : pmn Bug dans SetPeriodic (PRO6088) et Segment (PRO6250)
  29
  30 #define No_Standard_OutOfRange
  31
  32 #include <Geom_BSplineCurve.ixx>
  33 #include <gp.hxx>
  34 #include <ElCLib.hxx>
  35 #include <BSplCLib.hxx>
  36 #include <BSplCLib_KnotDistribution.hxx>
  37 #include <BSplCLib_MultDistribution.hxx>
  38 #include <Standard_NotImplemented.hxx>
  39 #include <Standard_ConstructionError.hxx>
  40 #include <Standard_OutOfRange.hxx>
  41 #include <Standard_Real.hxx>
  42
  43 //=======================================================================
  44 //function : CheckCurveData
  45 //purpose  : Internal use only
  46 //=======================================================================
  47
  48 static void CheckCurveData
  49 (const TColgp_Array1OfPnt&         CPoles,
  50  const TColStd_Array1OfReal&       CKnots,
  51  const TColStd_Array1OfInteger&    CMults,
  52  const Standard_Integer            Degree,
  53  const Standard_Boolean            Periodic)
  54 {
  55   if (Degree < 1 || Degree > Geom_BSplineCurve::MaxDegree()) {
  56     Standard_ConstructionError::Raise();
  57   }
  58
  59   if (CPoles.Length() < 2)                Standard_ConstructionError::Raise();
  60   if (CKnots.Length() != CMults.Length()) Standard_ConstructionError::Raise();
  61
  62   for (Standard_Integer I = CKnots.Lower(); I < CKnots.Upper(); I++) {
  63     if (CKnots (I+1) - CKnots (I) <= Epsilon (Abs(CKnots (I)))) {
  64       Standard_ConstructionError::Raise();
  65     }
  66   }
  67
  68   if (CPoles.Length() != BSplCLib::NbPoles(Degree,Periodic,CMults))
  69     Standard_ConstructionError::Raise();
  70 }
  71
  72 //=======================================================================
  73 //function : Rational
  74 //purpose  : check rationality of an array of weights
  75 //=======================================================================
  76
  77 static Standard_Boolean Rational(const TColStd_Array1OfReal& W)
  78 {
  79   Standard_Integer i, n = W.Length();
  80   Standard_Boolean rat = Standard_False;
  81   for (i = 1; i < n; i++) {
  82     rat =  Abs(W(i) - W(i+1)) > gp::Resolution();
  83     if (rat) break;
  84   }
  85   return rat;
  86 }
  87
  88 //=======================================================================
  89 //function : Copy
  90 //purpose  :
  91 //=======================================================================
  92
  93 Handle(Geom_Geometry) Geom_BSplineCurve::Copy() const
  94 {
  95   Handle(Geom_BSplineCurve) C;
  96   if (IsRational())
  97     C = new Geom_BSplineCurve(poles->Array1(),
  98                               weights->Array1(),
  99                               knots->Array1(),
 100                               mults->Array1(),
 101                               deg,periodic);
 102   else
 103     C = new Geom_BSplineCurve(poles->Array1(),
 104                               knots->Array1(),
 105                               mults->Array1(),
 106                               deg,periodic);
 107   return C;
 108 }
 109
 110 //=======================================================================
 111 //function : Geom_BSplineCurve
 112 //purpose  :
 113 //=======================================================================
 114
 115 Geom_BSplineCurve::Geom_BSplineCurve
 116 (const TColgp_Array1OfPnt&       Poles,
 117  const TColStd_Array1OfReal&     Knots,
 118  const TColStd_Array1OfInteger&  Mults,
 119  const Standard_Integer          Degree,
 120  const Standard_Boolean          Periodic) :
 121  rational(Standard_False),
 122  periodic(Periodic),
 123  deg(Degree),
 124  maxderivinvok(Standard_False)
 125 {
 126   // check
 127
 128   CheckCurveData (Poles,
 129                   Knots,
 130                   Mults,
 131                   Degree,
 132                   Periodic);
 133
 134
 135   // copy arrays
 136
 137   poles =  new TColgp_HArray1OfPnt(1,Poles.Length());
 138   poles->ChangeArray1() = Poles;
 139
 140
 141   knots = new TColStd_HArray1OfReal(1,Knots.Length());
 142   knots->ChangeArray1() = Knots;
 143
 144   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
 145   mults->ChangeArray1() = Mults;
 146
 147   UpdateKnots();
 148   cachepoles = new TColgp_HArray1OfPnt(1,Degree + 1);
 149   parametercache = 0.0e0 ;
 150   spanlenghtcache = 0.0e0 ;
 151   spanindexcache = 0 ;
 152
 153 }
 154
 155 //=======================================================================
 156 //function : Geom_BSplineCurve
 157 //purpose  :
 158 //=======================================================================
 159
 160 Geom_BSplineCurve::Geom_BSplineCurve
 161 (const TColgp_Array1OfPnt&      Poles,
 162  const TColStd_Array1OfReal&    Weights,
 163  const TColStd_Array1OfReal&    Knots,
 164  const TColStd_Array1OfInteger& Mults,
 165  const Standard_Integer         Degree,
 166  const Standard_Boolean         Periodic,
 167  const Standard_Boolean         CheckRational)  :
 168  rational(Standard_True),
 169  periodic(Periodic),
 170  deg(Degree),
 171  maxderivinvok(Standard_False)
 172
 173 {
 174
 175   // check
 176
 177   CheckCurveData (Poles,
 178                   Knots,
 179                   Mults,
 180                   Degree,
 181                   Periodic);
 182
 183   if (Weights.Length() != Poles.Length())
 184     Standard_ConstructionError::Raise("Geom_BSplineCurve");
 185
 186   Standard_Integer i;
 187   for (i = Weights.Lower(); i <= Weights.Upper(); i++) {
 188     if (Weights(i) <= gp::Resolution())
 189       Standard_ConstructionError::Raise("Geom_BSplineCurve");
 190   }
 191
 192   // check really rational
 193   if (CheckRational)
 194     rational = Rational(Weights);
 195
 196   // copy arrays
 197
 198   poles =  new TColgp_HArray1OfPnt(1,Poles.Length());
 199   poles->ChangeArray1() = Poles;
 200   cachepoles = new TColgp_HArray1OfPnt(1,Degree + 1);
 201   if (rational) {
 202     weights =  new TColStd_HArray1OfReal(1,Weights.Length());
 203     weights->ChangeArray1() = Weights;
 204     cacheweights  = new TColStd_HArray1OfReal(1,Degree + 1);
 205   }
 206
 207   knots = new TColStd_HArray1OfReal(1,Knots.Length());
 208   knots->ChangeArray1() = Knots;
 209
 210   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
 211   mults->ChangeArray1() = Mults;
 212
 213   UpdateKnots();
 214   parametercache = 0.0e0 ;
 215   spanlenghtcache = 0.0e0 ;
 216   spanindexcache = 0 ;
 217 }
 218
 219 //=======================================================================
 220 //function : MaxDegree
 221 //purpose  :
 222 //=======================================================================
 223
 224 Standard_Integer Geom_BSplineCurve::MaxDegree ()
 225 {
 226   return BSplCLib::MaxDegree();
 227 }
 228
 229 //=======================================================================
 230 //function : IncreaseDegree
 231 //purpose  :
 232 //=======================================================================
 233
 234 void Geom_BSplineCurve::IncreaseDegree  (const Standard_Integer Degree)
 235 {
 236   if (Degree == deg) return;
 237
 238   if (Degree < deg || Degree > Geom_BSplineCurve::MaxDegree()) {
 239     Standard_ConstructionError::Raise();
 240   }
 241   Standard_Integer FromK1 = FirstUKnotIndex ();
 242   Standard_Integer ToK2   = LastUKnotIndex  ();
 243
 244   Standard_Integer Step   = Degree - deg;
 245
 246   Handle(TColgp_HArray1OfPnt) npoles = new
 247     TColgp_HArray1OfPnt(1,poles->Length() + Step * (ToK2-FromK1));
 248
 249   Standard_Integer nbknots = BSplCLib::IncreaseDegreeCountKnots
 250     (deg,Degree,periodic,mults->Array1());
 251
 252   Handle(TColStd_HArray1OfReal) nknots =
 253     new TColStd_HArray1OfReal(1,nbknots);
 254
 255   Handle(TColStd_HArray1OfInteger) nmults =
 256     new TColStd_HArray1OfInteger(1,nbknots);
 257
 258   Handle(TColStd_HArray1OfReal) nweights;
 259
 260   if (IsRational()) {
 261
 262     nweights = new TColStd_HArray1OfReal(1,npoles->Upper());
 263
 264     BSplCLib::IncreaseDegree
 265       (deg,Degree, periodic,
 266        poles->Array1(),weights->Array1(),
 267        knots->Array1(),mults->Array1(),
 268        npoles->ChangeArray1(),nweights->ChangeArray1(),
 269        nknots->ChangeArray1(),nmults->ChangeArray1());
 270   }
 271   else {
 272     BSplCLib::IncreaseDegree
 273       (deg,Degree, periodic,
 274        poles->Array1(),BSplCLib::NoWeights(),
 275        knots->Array1(),mults->Array1(),
 276        npoles->ChangeArray1(),
 277        *((TColStd_Array1OfReal*) NULL),
 278        nknots->ChangeArray1(),nmults->ChangeArray1());
 279   }
 280
 281   deg     = Degree;
 282   poles   = npoles;
 283   weights = nweights;
 284   knots   = nknots;
 285   mults   = nmults;
 286   UpdateKnots();
 287
 288 }
 289
 290 //=======================================================================
 291 //function : IncreaseMultiplicity
 292 //purpose  :
 293 //=======================================================================
 294
 295 void Geom_BSplineCurve::IncreaseMultiplicity  (const Standard_Integer Index,
 296                                                const Standard_Integer M)
 297 {
 298   TColStd_Array1OfReal k(1,1);
 299   k(1) = knots->Value(Index);
 300   TColStd_Array1OfInteger m(1,1);
 301   m(1) = M - mults->Value(Index);
 302   InsertKnots(k,m,Epsilon(1.),Standard_True);
 303 }
 304
 305 //=======================================================================
 306 //function : IncreaseMultiplicity
 307 //purpose  :
 308 //=======================================================================
 309
 310 void Geom_BSplineCurve::IncreaseMultiplicity  (const Standard_Integer I1,
 311                                                const Standard_Integer I2,
 312                                                const Standard_Integer M)
 313 {
 314   Handle(TColStd_HArray1OfReal)  tk = knots;
 315   TColStd_Array1OfReal k((knots->Array1())(I1),I1,I2);
 316   TColStd_Array1OfInteger m(I1,I2);
 317   Standard_Integer i;
 318   for (i = I1; i <= I2; i++)
 319     m(i) = M - mults->Value(i);
 320   InsertKnots(k,m,Epsilon(1.),Standard_True);
 321 }
 322
 323 //=======================================================================
 324 //function : IncrementMultiplicity
 325 //purpose  :
 326 //=======================================================================
 327
 328 void Geom_BSplineCurve::IncrementMultiplicity
 329 (const Standard_Integer I1,
 330  const Standard_Integer I2,
 331  const Standard_Integer Step)
 332 {
 333   Handle(TColStd_HArray1OfReal) tk = knots;
 334   TColStd_Array1OfReal    k((knots->Array1())(I1),I1,I2);
 335   TColStd_Array1OfInteger m(I1,I2) ;
 336   m.Init(Step);
 337   InsertKnots(k,m,Epsilon(1.),Standard_True);
 338 }
 339
 340 //=======================================================================
 341 //function : InsertKnot
 342 //purpose  :
 343 //=======================================================================
 344
 345 void Geom_BSplineCurve::InsertKnot
 346 (const Standard_Real U,
 347  const Standard_Integer M,
 348  const Standard_Real ParametricTolerance,
 349  const Standard_Boolean Add)
 350 {
 351   TColStd_Array1OfReal k(1,1);
 352   k(1) = U;
 353   TColStd_Array1OfInteger m(1,1);
 354   m(1) = M;
 355   InsertKnots(k,m,ParametricTolerance,Add);
 356 }
 357
 358 //=======================================================================
 359 //function : InsertKnots
 360 //purpose  :
 361 //=======================================================================
 362
 363 void  Geom_BSplineCurve::InsertKnots(const TColStd_Array1OfReal& Knots,
 364                                      const TColStd_Array1OfInteger& Mults,
 365                                      const Standard_Real Epsilon,
 366                                      const Standard_Boolean Add)
 367 {
 368   // Check and compute new sizes
 369   Standard_Integer nbpoles,nbknots;
 370
 371   if (!BSplCLib::PrepareInsertKnots(deg,periodic,
 372                                     knots->Array1(),mults->Array1(),
 373                                     Knots,Mults,nbpoles,nbknots,Epsilon,Add))
 374     Standard_ConstructionError::Raise("Geom_BSplineCurve::InsertKnots");
 375
 376   if (nbpoles == poles->Length()) return;
 377
 378   Handle(TColgp_HArray1OfPnt)      npoles = new TColgp_HArray1OfPnt(1,nbpoles);
 379   Handle(TColStd_HArray1OfReal)    nknots = knots;
 380   Handle(TColStd_HArray1OfInteger) nmults = mults;
 381
 382   if (nbknots != knots->Length()) {
 383     nknots = new TColStd_HArray1OfReal(1,nbknots);
 384     nmults = new TColStd_HArray1OfInteger(1,nbknots);
 385   }
 386
 387   if (rational) {
 388     Handle(TColStd_HArray1OfReal) nweights =
 389       new TColStd_HArray1OfReal(1,nbpoles);
 390     BSplCLib::InsertKnots(deg,periodic,
 391                           poles->Array1(), weights->Array1(),
 392                           knots->Array1(), mults->Array1(),
 393                           Knots, Mults,
 394                           npoles->ChangeArray1(), nweights->ChangeArray1(),
 395                           nknots->ChangeArray1(), nmults->ChangeArray1(),
 396                           Epsilon, Add);
 397     weights = nweights;
 398   }
 399   else {
 400     BSplCLib::InsertKnots(deg,periodic,
 401                           poles->Array1(), BSplCLib::NoWeights(),
 402                           knots->Array1(), mults->Array1(),
 403                           Knots, Mults,
 404                           npoles->ChangeArray1(),
 405                           *((TColStd_Array1OfReal*) NULL),
 406                           nknots->ChangeArray1(), nmults->ChangeArray1(),
 407                           Epsilon, Add);
 408   }
 409
 410   poles = npoles;
 411   knots = nknots;
 412   mults = nmults;
 413   UpdateKnots();
 414
 415 }
 416
 417 //=======================================================================
 418 //function : RemoveKnot
 419 //purpose  :
 420 //=======================================================================
 421
 422 Standard_Boolean  Geom_BSplineCurve::RemoveKnot(const Standard_Integer Index,
 423                                                 const Standard_Integer M,
 424                                                 const Standard_Real Tolerance)
 425 {
 426   if (M < 0) return Standard_True;
 427
 428   Standard_Integer I1  = FirstUKnotIndex ();
 429   Standard_Integer I2  = LastUKnotIndex  ();
 430
 431   if ( !periodic && (Index <= I1 || Index >= I2) ) {
 432     Standard_OutOfRange::Raise();
 433   }
 434   else if ( periodic  && (Index < I1 || Index > I2)) {
 435     Standard_OutOfRange::Raise();
 436   }
 437
 438   const TColgp_Array1OfPnt   & oldpoles   = poles->Array1();
 439
 440   Standard_Integer step = mults->Value(Index) - M;
 441   if (step <= 0) return Standard_True;
 442
 443   Handle(TColgp_HArray1OfPnt) npoles =
 444     new TColgp_HArray1OfPnt(1,oldpoles.Length()-step);
 445
 446   Handle(TColStd_HArray1OfReal)    nknots  = knots;
 447   Handle(TColStd_HArray1OfInteger) nmults  = mults;
 448
 449   if (M == 0) {
 450     nknots = new TColStd_HArray1OfReal(1,knots->Length()-1);
 451     nmults = new TColStd_HArray1OfInteger(1,knots->Length()-1);
 452   }
 453
 454   if (IsRational()) {
 455     Handle(TColStd_HArray1OfReal) nweights =
 456       new TColStd_HArray1OfReal(1,npoles->Length());
 457     if (!BSplCLib::RemoveKnot
 458         (Index, M, deg, periodic,
 459          poles->Array1(),weights->Array1(),
 460          knots->Array1(),mults->Array1(),
 461          npoles->ChangeArray1(), nweights->ChangeArray1(),
 462          nknots->ChangeArray1(),nmults->ChangeArray1(),
 463          Tolerance))
 464       return Standard_False;
 465     weights = nweights;
 466   }
 467   else {
 468     if (!BSplCLib::RemoveKnot
 469         (Index, M, deg, periodic,
 470          poles->Array1(), BSplCLib::NoWeights(),
 471          knots->Array1(),mults->Array1(),
 472          npoles->ChangeArray1(),
 473          *((TColStd_Array1OfReal*) NULL),
 474          nknots->ChangeArray1(),nmults->ChangeArray1(),
 475          Tolerance))
 476       return Standard_False;
 477   }
 478
 479   poles = npoles;
 480   knots = nknots;
 481   mults = nmults;
 482
 483   UpdateKnots();
 484   maxderivinvok = 0;
 485   return Standard_True;
 486 }
 487
 488 //=======================================================================
 489 //function : Reverse
 490 //purpose  :
 491 //=======================================================================
 492
 493 void Geom_BSplineCurve::Reverse ()
 494 {
 495   BSplCLib::Reverse(knots->ChangeArray1());
 496   BSplCLib::Reverse(mults->ChangeArray1());
 497   Standard_Integer last;
 498   if (periodic)
 499     last = flatknots->Upper() - deg - 1;
 500   else
 501     last = poles->Upper();
 502   BSplCLib::Reverse(poles->ChangeArray1(),last);
 503   if (rational)
 504     BSplCLib::Reverse(weights->ChangeArray1(),last);
 505   UpdateKnots();
 506 }
 507
 508 //=======================================================================
 509 //function : ReversedParameter
 510 //purpose  :
 511 //=======================================================================
 512
 513 Standard_Real Geom_BSplineCurve::ReversedParameter
 514 (const Standard_Real U) const
 515 {
 516   return (FirstParameter() + LastParameter() - U);
 517 }
 518
 519 //=======================================================================
 520 //function : Segment
 521 //purpose  :
 522 //=======================================================================
 523
 524 void Geom_BSplineCurve::Segment(const Standard_Real U1,
 525                                 const Standard_Real U2)
 526 {
 527   Standard_DomainError_Raise_if ( U2 < U1,
 528                                  "Geom_BSplineCurve::Segment");
 529
 530   Standard_Real NewU1, NewU2;
 531   Standard_Real U,DU=0,aDDU=0;
 532   Standard_Integer index;
 533   Standard_Boolean wasPeriodic = periodic;
 534
 535   TColStd_Array1OfReal    Knots(1,2);
 536   TColStd_Array1OfInteger Mults(1,2);
 537
 538   // define param ditance to keep (eap, Apr 18 2002, occ311)
 539   if (periodic) {
 540     Standard_Real Period = LastParameter() - FirstParameter();
 541     DU = U2 - U1;
 542     while (DU > Period)
 543       DU -= Period;
 544     if (DU <= Epsilon(Period))
 545       DU = Period;
 546     aDDU = DU;
 547   }
 548
 549   index = 0;
 550   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 551                             U1,periodic,knots->Lower(),knots->Upper(),
 552                             index,NewU1);
 553   index = 0;
 554   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 555                             U2,periodic,knots->Lower(),knots->Upper(),
 556                             index,NewU2);
 557
 558   //-- DBB
 559   Standard_Real aNu2 = NewU2;
 560   //-- DBB
 561
 562
 563   Knots( 1) = Min( NewU1, NewU2);
 564   Knots( 2) = Max( NewU1, NewU2);
 565   Mults( 1) = Mults( 2) = deg;
 566
 567   Standard_Real AbsUMax = Max(Abs(NewU1),Abs(NewU2));
 568
 569 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 Begin
 570   AbsUMax = Max(AbsUMax, Max(Abs(FirstParameter()),Abs(LastParameter())));
 571 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 End
 572
 573   Standard_Real Eps = 100. * Epsilon(AbsUMax);
 574
 575   InsertKnots( Knots, Mults, Eps);
 576
 577   if (periodic) { // set the origine at NewU1
 578     index = 0;
 579     BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 580                               U1,periodic,knots->Lower(),knots->Upper(),
 581                               index,U);
 582     // Test si l'insertion est Ok et decalage sinon.
 583     if ( Abs(knots->Value(index+1)-U) <= Eps) // <= pour etre homogene a InsertKnots
 584       index++;
 585     SetOrigin(index);
 586     SetNotPeriodic();
 587     NewU2 = NewU1 + DU;
 588   }
 589
 590   // compute index1 and index2 to set the new knots and mults
 591   Standard_Integer index1 = 0, index2 = 0;
 592   Standard_Integer FromU1 = knots->Lower();
 593   Standard_Integer ToU2   = knots->Upper();
 594   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 595                             NewU1,periodic,FromU1,ToU2,index1,U);
 596   if ( Abs(knots->Value(index1+1)-U) <= Eps)
 597     index1++;
 598
 599   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 600                             NewU2,periodic,FromU1,ToU2,index2,U);
 601   if ( Abs(knots->Value(index2+1)-U) <= Eps)
 602     index2++;
 603
 604   Standard_Integer nbknots = index2 - index1 + 1;
 605
 606   Handle(TColStd_HArray1OfReal)
 607     nknots = new TColStd_HArray1OfReal(1,nbknots);
 608   Handle(TColStd_HArray1OfInteger)
 609     nmults = new TColStd_HArray1OfInteger(1,nbknots);
 610
 611   // to restore changed U1
 612   if (DU > 0) // if was periodic
 613     DU = NewU1 - U1;
 614
 615   Standard_Integer i , k = 1;
 616   for ( i = index1; i<= index2; i++) {
 617     nknots->SetValue(k, knots->Value(i) - DU);
 618     nmults->SetValue(k, mults->Value(i));
 619     k++;
 620   }
 621   nmults->SetValue(      1, deg + 1);
 622   nmults->SetValue(nbknots, deg + 1);
 623
 624
 625   // compute index1 and index2 to set the new poles and weights
 626   Standard_Integer pindex1
 627     = BSplCLib::PoleIndex(deg,index1,periodic,mults->Array1());
 628   Standard_Integer pindex2
 629     = BSplCLib::PoleIndex(deg,index2,periodic,mults->Array1());
 630
 631   pindex1++;
 632   pindex2 = Min( pindex2+1, poles->Length());
 633
 634   Standard_Integer nbpoles  = pindex2 - pindex1 + 1;
 635
 636   Handle(TColStd_HArray1OfReal)
 637     nweights = new TColStd_HArray1OfReal(1,nbpoles);
 638   Handle(TColgp_HArray1OfPnt)
 639     npoles = new TColgp_HArray1OfPnt(1,nbpoles);
 640
 641   k = 1;
 642   if ( rational) {
 643     nweights = new TColStd_HArray1OfReal( 1, nbpoles);
 644     for ( i = pindex1; i <= pindex2; i++) {
 645       npoles->SetValue(k, poles->Value(i));
 646       nweights->SetValue(k, weights->Value(i));
 647       k++;
 648     }
 649   }
 650   else {
 651     for ( i = pindex1; i <= pindex2; i++) {
 652       npoles->SetValue(k, poles->Value(i));
 653       k++;
 654     }
 655   }
 656
 657   //-- DBB
 658   if( wasPeriodic ) {
 659     nknots->ChangeValue(nknots->Lower()) = U1;
 660     if( aNu2 < U2 ) {
 661       nknots->ChangeValue(nknots->Upper()) = U1 + aDDU;
 662     }
 663   }
 664   //-- DBB
 665
 666   knots = nknots;
 667   mults = nmults;
 668   poles = npoles;
 669   if ( rational)
 670     weights = nweights;
 671
 672   maxderivinvok = 0;
 673   UpdateKnots();
 674 }
 675
 676 //=======================================================================
 677 //function : SetKnot
 678 //purpose  :
 679 //=======================================================================
 680
 681 void Geom_BSplineCurve::SetKnot
 682 (const Standard_Integer Index,
 683  const Standard_Real K)
 684 {
 685   if (Index < 1 || Index > knots->Length())     Standard_OutOfRange::Raise();
 686   Standard_Real DK = Abs(Epsilon (K));
 687   if (Index == 1) {
 688     if (K >= knots->Value(2) - DK) Standard_ConstructionError::Raise();
 689   }
 690   else if (Index == knots->Length()) {
 691     if (K <= knots->Value (knots->Length()-1) + DK)  {
 692       Standard_ConstructionError::Raise();
 693     }
 694   }
 695   else {
 696     if (K <= knots->Value(Index-1) + DK ||
 697         K >= knots->Value(Index+1) - DK ) {
 698       Standard_ConstructionError::Raise();
 699     }
 700   }
 701   if (K != knots->Value (Index)) {
 702     knots->SetValue (Index, K);
 703     maxderivinvok = 0;
 704     UpdateKnots();
 705   }
 706 }
 707
 708 //=======================================================================
 709 //function : SetKnots
 710 //purpose  :
 711 //=======================================================================
 712
 713 void Geom_BSplineCurve::SetKnots
 714 (const TColStd_Array1OfReal& K)
 715 {
 716   CheckCurveData(poles->Array1(),K,mults->Array1(),deg,periodic);
 717   knots->ChangeArray1() = K;
 718   maxderivinvok = 0;
 719   UpdateKnots();
 720 }
 721
 722 //=======================================================================
 723 //function : SetKnot
 724 //purpose  :
 725 //=======================================================================
 726
 727 void Geom_BSplineCurve::SetKnot
 728 (const Standard_Integer Index,
 729  const Standard_Real K,
 730  const Standard_Integer M)
 731 {
 732   IncreaseMultiplicity (Index, M);
 733   SetKnot (Index, K);
 734 }
 735
 736 //=======================================================================
 737 //function : SetPeriodic
 738 //purpose  :
 739 //=======================================================================
 740
 741 void Geom_BSplineCurve::SetPeriodic ()
 742 {
 743   Standard_Integer first = FirstUKnotIndex();
 744   Standard_Integer last  = LastUKnotIndex();
 745
 746   Handle(TColStd_HArray1OfReal) tk = knots;
 747   TColStd_Array1OfReal    cknots((knots->Array1())(first),first,last);
 748   knots = new TColStd_HArray1OfReal(1,cknots.Length());
 749   knots->ChangeArray1() = cknots;
 750
 751   Handle(TColStd_HArray1OfInteger) tm = mults;
 752   TColStd_Array1OfInteger cmults((mults->Array1())(first),first,last);
 753   cmults(first) = cmults(last) = Min(deg, Max( cmults(first), cmults(last)));
 754   mults = new TColStd_HArray1OfInteger(1,cmults.Length());
 755   mults->ChangeArray1() = cmults;
 756
 757   // compute new number of poles;
 758   Standard_Integer nbp = BSplCLib::NbPoles(deg,Standard_True,cmults);
 759
 760   Handle(TColgp_HArray1OfPnt) tp = poles;
 761   TColgp_Array1OfPnt cpoles((poles->Array1())(1),1,nbp);
 762   poles = new TColgp_HArray1OfPnt(1,nbp);
 763   poles->ChangeArray1() = cpoles;
 764
 765   if (rational) {
 766     Handle(TColStd_HArray1OfReal) tw = weights;
 767     TColStd_Array1OfReal cweights((weights->Array1())(1),1,nbp);
 768     weights = new TColStd_HArray1OfReal(1,nbp);
 769     weights->ChangeArray1() = cweights;
 770   }
 771
 772   periodic = Standard_True;
 773
 774   maxderivinvok = 0;
 775   UpdateKnots();
 776 }
 777
 778 //=======================================================================
 779 //function : SetOrigin
 780 //purpose  :
 781 //=======================================================================
 782
 783 void Geom_BSplineCurve::SetOrigin(const Standard_Integer Index)
 784 {
 785   Standard_NoSuchObject_Raise_if( !periodic,
 786                                  "Geom_BSplineCurve::SetOrigin");
 787   Standard_Integer i,k;
 788   Standard_Integer first = FirstUKnotIndex();
 789   Standard_Integer last  = LastUKnotIndex();
 790
 791   Standard_DomainError_Raise_if( (Index < first) || (Index > last),
 792                                 "Geom_BSplineCurve::SetOrigine");
 793
 794   Standard_Integer nbknots = knots->Length();
 795   Standard_Integer nbpoles = poles->Length();
 796
 797   Handle(TColStd_HArray1OfReal) nknots =
 798     new TColStd_HArray1OfReal(1,nbknots);
 799   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
 800
 801   Handle(TColStd_HArray1OfInteger) nmults =
 802     new TColStd_HArray1OfInteger(1,nbknots);
 803   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
 804
 805   // set the knots and mults
 806   Standard_Real period = knots->Value(last) - knots->Value(first);
 807   k = 1;
 808   for ( i = Index; i <= last ; i++) {
 809     newknots(k) = knots->Value(i);
 810     newmults(k) = mults->Value(i);
 811     k++;
 812   }
 813   for ( i = first+1; i <= Index; i++) {
 814     newknots(k) = knots->Value(i) + period;
 815     newmults(k) = mults->Value(i);
 816     k++;
 817   }
 818
 819   Standard_Integer index = 1;
 820   for (i = first+1; i <= Index; i++)
 821     index += mults->Value(i);
 822
 823   // set the poles and weights
 824   Handle(TColgp_HArray1OfPnt) npoles =
 825     new TColgp_HArray1OfPnt(1,nbpoles);
 826   Handle(TColStd_HArray1OfReal) nweights =
 827     new TColStd_HArray1OfReal(1,nbpoles);
 828   TColgp_Array1OfPnt   & newpoles   = npoles->ChangeArray1();
 829   TColStd_Array1OfReal & newweights = nweights->ChangeArray1();
 830   first = poles->Lower();
 831   last  = poles->Upper();
 832   if ( rational) {
 833     k = 1;
 834     for ( i = index; i <= last; i++) {
 835       newpoles(k)   = poles->Value(i);
 836       newweights(k) = weights->Value(i);
 837       k++;
 838     }
 839     for ( i = first; i < index; i++) {
 840       newpoles(k)   = poles->Value(i);
 841       newweights(k) = weights->Value(i);
 842       k++;
 843     }
 844   }
 845   else {
 846     k = 1;
 847     for ( i = index; i <= last; i++) {
 848       newpoles(k) = poles->Value(i);
 849       k++;
 850     }
 851     for ( i = first; i < index; i++) {
 852       newpoles(k) = poles->Value(i);
 853       k++;
 854     }
 855   }
 856
 857   poles = npoles;
 858   knots = nknots;
 859   mults = nmults;
 860   if (rational)
 861     weights = nweights;
 862   maxderivinvok = 0;
 863   UpdateKnots();
 864 }
 865
 866 //=======================================================================
 867 //function : SetOrigin
 868 //purpose  :
 869 //=======================================================================
 870
 871 void Geom_BSplineCurve::SetOrigin(const Standard_Real U,
 872                                   const Standard_Real Tol)
 873 {
 874   Standard_NoSuchObject_Raise_if( !periodic,
 875                                  "Geom_BSplineCurve::SetOrigin");
 876   //U est il dans la period.
 877   Standard_Real uf = FirstParameter(), ul = LastParameter();
 878   Standard_Real u = U, period = ul - uf;
 879   while (Tol < (uf-u)) u += period;
 880   while (Tol > (ul-u)) u -= period;
 881
 882   if(Abs(U-u)>Tol) { //On reparametre la courbe
 883     Standard_Real delta = U-u;
 884     uf += delta;
 885     ul += delta;
 886     TColStd_Array1OfReal& kn = knots->ChangeArray1();
 887     Standard_Integer fk = kn.Lower(), lk = kn.Upper();
 888     for(Standard_Integer i = fk; i <= lk; i++){
 889       kn.ChangeValue(i) += delta;
 890     }
 891     UpdateKnots();
 892   }
 893   if(Abs(U-uf)<Tol) return;
 894
 895   TColStd_Array1OfReal& kn = knots->ChangeArray1();
 896   Standard_Integer fk = kn.Lower(), lk = kn.Upper(),ik=0;
 897   Standard_Real delta = RealLast();
 898   for(Standard_Integer i = fk; i<= lk; i++){
 899     Standard_Real dki = kn.Value(i)-U;
 900     if(Abs(dki)<Abs(delta)){
 901       ik = i;
 902       delta = dki;
 903     }
 904   }
 905   if(Abs(delta)>Tol){
 906     InsertKnot(U);
 907     if(delta < 0.) ik++;
 908   }
 909   SetOrigin(ik);
 910 }
 911
 912 //=======================================================================
 913 //function : SetNotPeriodic
 914 //purpose  :
 915 //=======================================================================
 916
 917 void Geom_BSplineCurve::SetNotPeriodic ()
 918 {
 919   if ( periodic) {
 920     Standard_Integer NbKnots, NbPoles;
 921     BSplCLib::PrepareUnperiodize( deg, mults->Array1(),NbKnots,NbPoles);
 922
 923     Handle(TColgp_HArray1OfPnt) npoles
 924       = new TColgp_HArray1OfPnt(1,NbPoles);
 925
 926     Handle(TColStd_HArray1OfReal) nknots
 927       = new TColStd_HArray1OfReal(1,NbKnots);
 928
 929     Handle(TColStd_HArray1OfInteger) nmults
 930       = new TColStd_HArray1OfInteger(1,NbKnots);
 931
 932     Handle(TColStd_HArray1OfReal) nweights;
 933
 934     if (IsRational()) {
 935
 936       nweights = new TColStd_HArray1OfReal(1,NbPoles);
 937
 938       BSplCLib::Unperiodize
 939         (deg,mults->Array1(),knots->Array1(),poles->Array1(),
 940          weights->Array1(),nmults->ChangeArray1(),
 941          nknots->ChangeArray1(),npoles->ChangeArray1(),
 942          nweights->ChangeArray1());
 943
 944     }
 945     else {
 946
 947       BSplCLib::Unperiodize
 948         (deg,mults->Array1(),knots->Array1(),poles->Array1(),
 949          BSplCLib::NoWeights(),nmults->ChangeArray1(),
 950          nknots->ChangeArray1(),npoles->ChangeArray1(),
 951          *((TColStd_Array1OfReal*) NULL));
 952
 953     }
 954     poles   = npoles;
 955     weights = nweights;
 956     mults   = nmults;
 957     knots   = nknots;
 958     periodic = Standard_False;
 959
 960     maxderivinvok = 0;
 961     UpdateKnots();
 962   }
 963 }
 964
 965 //=======================================================================
 966 //function : SetPole
 967 //purpose  :
 968 //=======================================================================
 969
 970 void Geom_BSplineCurve::SetPole
 971 (const Standard_Integer Index,
 972  const gp_Pnt& P)
 973 {
 974   if (Index < 1 || Index > poles->Length()) Standard_OutOfRange::Raise();
 975   poles->SetValue (Index, P);
 976   maxderivinvok = 0;
 977   InvalidateCache() ;
 978 }
 979
 980 //=======================================================================
 981 //function : SetPole
 982 //purpose  :
 983 //=======================================================================
 984
 985 void Geom_BSplineCurve::SetPole
 986 (const Standard_Integer Index,
 987  const gp_Pnt& P,
 988  const Standard_Real W)
 989 {
 990   SetPole(Index,P);
 991   SetWeight(Index,W);
 992 }
 993
 994 //=======================================================================
 995 //function : SetWeight
 996 //purpose  :
 997 //=======================================================================
 998
 999 void Geom_BSplineCurve::SetWeight
1000 (const Standard_Integer Index,
1001  const Standard_Real W)
1002 {
1003   if (Index < 1 || Index > poles->Length())   Standard_OutOfRange::Raise();
1004
1005   if (W <= gp::Resolution ())     Standard_ConstructionError::Raise();
1006
1007
1008   Standard_Boolean rat = IsRational() || (Abs(W - 1.) > gp::Resolution());
1009
1010   if ( rat) {
1011     if (rat && !IsRational()) {
1012       weights = new TColStd_HArray1OfReal(1,poles->Length());
1013       weights->Init(1.);
1014     }
1015
1016     weights->SetValue (Index, W);
1017
1018     if (IsRational()) {
1019       rat = Rational(weights->Array1());
1020       if (!rat) weights.Nullify();
1021     }
1022
1023     rational = !weights.IsNull();
1024   }
1025   maxderivinvok = 0;
1026   InvalidateCache() ;
1027 }
1028
1029 //=======================================================================
1030 //function : MovePoint
1031 //purpose  :
1032 //=======================================================================
1033
1034 void Geom_BSplineCurve::MovePoint(const Standard_Real U,
1035                                   const gp_Pnt& P,
1036                                   const Standard_Integer Index1,
1037                                   const Standard_Integer Index2,
1038                                   Standard_Integer& FirstModifiedPole,
1039                                   Standard_Integer& LastmodifiedPole)
1040 {
1041   if (Index1 < 1 || Index1 > poles->Length() ||
1042       Index2 < 1 || Index2 > poles->Length() || Index1 > Index2) {
1043     Standard_OutOfRange::Raise();
1044   }
1045   TColgp_Array1OfPnt npoles(1, poles->Length());
1046   gp_Pnt P0;
1047   D0(U, P0);
1048   gp_Vec Displ(P0, P);
1049   BSplCLib::MovePoint(U, Displ, Index1, Index2, deg, rational, poles->Array1(),
1050                       weights->Array1(), flatknots->Array1(),
1051                       FirstModifiedPole, LastmodifiedPole, npoles);
1052   if (FirstModifiedPole) {
1053     poles->ChangeArray1() = npoles;
1054     maxderivinvok = 0;
1055     InvalidateCache() ;
1056   }
1057 }
1058
1059 //=======================================================================
1060 //function : MovePointAndTangent
1061 //purpose  :
1062 //=======================================================================
1063
1064 void Geom_BSplineCurve::
1065 MovePointAndTangent(const Standard_Real U,
1066                     const gp_Pnt&       P,
1067                     const gp_Vec&       Tangent,
1068                     const Standard_Real    Tolerance,
1069                     const Standard_Integer StartingCondition,
1070                     const Standard_Integer EndingCondition,
1071                     Standard_Integer&      ErrorStatus)
1072 {
1073   Standard_Integer ii ;
1074   if (IsPeriodic()) {
1075     //
1076     // for the time being do not deal with periodic curves
1077     //
1078     SetNotPeriodic() ;
1079   }
1080   TColgp_Array1OfPnt new_poles(1, poles->Length());
1081   gp_Pnt P0;
1082
1083
1084   gp_Vec delta_derivative;
1085   D1(U, P0,
1086      delta_derivative) ;
1087   gp_Vec delta(P0, P);
1088   for (ii = 1 ; ii <= 3 ; ii++) {
1089     delta_derivative.SetCoord(ii,
1090                               Tangent.Coord(ii)- delta_derivative.Coord(ii)) ;
1091   }
1092   BSplCLib::MovePointAndTangent(U,
1093                                 delta,
1094                                 delta_derivative,
1095                                 Tolerance,
1096                                 deg,
1097                                 rational,
1098                                 StartingCondition,
1099                                 EndingCondition,
1100                                 poles->Array1(),
1101                                 weights->Array1(),
1102                                 flatknots->Array1(),
1103                                 new_poles,
1104                                 ErrorStatus) ;
1105   if (!ErrorStatus) {
1106     poles->ChangeArray1() = new_poles;
1107     maxderivinvok = 0;
1108     InvalidateCache() ;
1109   }
1110 }
1111
1112 //=======================================================================
1113 //function : UpdateKnots
1114 //purpose  :
1115 //=======================================================================
1116
1117 void Geom_BSplineCurve::UpdateKnots()
1118 {
1119   rational = !weights.IsNull();
1120
1121   Standard_Integer MaxKnotMult = 0;
1122   BSplCLib::KnotAnalysis (deg,
1123                 periodic,
1124                 knots->Array1(),
1125                 mults->Array1(),
1126                 knotSet, MaxKnotMult);
1127
1128   if (knotSet == GeomAbs_Uniform && !periodic)  {
1129     flatknots = knots;
1130   }
1131   else {
1132     flatknots = new TColStd_HArray1OfReal
1133       (1, BSplCLib::KnotSequenceLength(mults->Array1(),deg,periodic));
1134
1135     BSplCLib::KnotSequence (knots->Array1(),
1136                             mults->Array1(),
1137                             deg,periodic,
1138                             flatknots->ChangeArray1());
1139   }
1140
1141   if (MaxKnotMult == 0)  smooth = GeomAbs_CN;
1142   else {
1143     switch (deg - MaxKnotMult) {
1144     case 0 :   smooth = GeomAbs_C0;   break;
1145     case 1 :   smooth = GeomAbs_C1;   break;
1146     case 2 :   smooth = GeomAbs_C2;   break;
1147     case 3 :   smooth = GeomAbs_C3;   break;
1148       default :  smooth = GeomAbs_C3;   break;
1149     }
1150   }
1151   InvalidateCache() ;
1152 }
1153
1154 //=======================================================================
1155 //function : Invalidate the Cache
1156 //purpose  : as the name says
1157 //=======================================================================
1158
1159 void Geom_BSplineCurve::InvalidateCache()
1160 {
1161   validcache = 0 ;
1162 }
1163
1164 //=======================================================================
1165 //function : check if the Cache is valid
1166 //purpose  : as the name says
1167 //=======================================================================
1168
1169 Standard_Boolean Geom_BSplineCurve::IsCacheValid
1170 (const Standard_Real  U)  const
1171 {
1172   //Roman Lygin 26.12.08, performance improvements
1173   //1. avoided using NewParameter = (U - parametercache) / spanlenghtcache
1174   //to check against [0, 1), as division is CPU consuming
1175   //2. minimized use of if, as branching is also CPU consuming
1176   Standard_Real aDelta = U - parametercache;
1177
1178   return ( validcache &&
1179       (aDelta >= 0.0e0) &&
1180       ((aDelta < spanlenghtcache) || (spanindexcache == flatknots->Upper() - deg)) );
1181 }
1182
1183 //=======================================================================
1184 //function : Normalizes the parameters if the curve is periodic
1185 //purpose  : that is compute the cache so that it is valid
1186 //=======================================================================
1187
1188 void Geom_BSplineCurve::PeriodicNormalization(Standard_Real&  Parameter) const
1189 {
1190   Standard_Real Period ;
1191
1192   if (periodic) {
1193     Period = flatknots->Value(flatknots->Upper() - deg) - flatknots->Value (deg + 1) ;
1194     while (Parameter > flatknots->Value(flatknots->Upper()-deg)) {
1195       Parameter -= Period ;
1196     }
1197     while (Parameter < flatknots->Value((deg + 1))) {
1198       Parameter +=  Period ;
1199     }
1200   }
1201 }
1202
1203 //=======================================================================
1204 //function : Validate the Cache
1205 //purpose  : that is compute the cache so that it is valid
1206 //=======================================================================
1207
1208 void Geom_BSplineCurve::ValidateCache(const Standard_Real  Parameter)
1209 {
1210   Standard_Real NewParameter ;
1211   Standard_Integer LocalIndex = 0 ;
1212   //
1213   // check if the degree did not change
1214   //
1215   if (cachepoles->Upper() < deg + 1) {
1216     cachepoles = new TColgp_HArray1OfPnt(1,deg + 1);
1217     if (rational) {
1218       cacheweights  = new TColStd_HArray1OfReal(1,deg + 1);
1219     }
1220   }
1221   BSplCLib::LocateParameter(deg,
1222                             (flatknots->Array1()),
1223                             (BSplCLib::NoMults()),
1224                             Parameter,
1225                             periodic,
1226                             LocalIndex,
1227                             NewParameter);
1228   spanindexcache = LocalIndex ;
1229   if (Parameter == flatknots->Value(LocalIndex + 1)) {
1230
1231     LocalIndex += 1 ;
1232     parametercache = flatknots->Value(LocalIndex) ;
1233     if (LocalIndex == flatknots->Upper() - deg) {
1234       //
1235       // for the last span if the parameter is outside of
1236       // the domain of the curve than use the last knot
1237       // and normalize with the last span Still set the
1238       // spanindexcache to flatknots->Upper() - deg so that
1239       // the IsCacheValid will know for sure we are extending
1240       // the Bspline
1241       //
1242
1243       spanlenghtcache = flatknots->Value(LocalIndex - 1) - parametercache ;
1244     }
1245     else {
1246       spanlenghtcache = flatknots->Value(LocalIndex + 1) - parametercache ;
1247     }
1248   }
1249   else {
1250     parametercache = flatknots->Value(LocalIndex) ;
1251     spanlenghtcache = flatknots->Value(LocalIndex + 1) - parametercache ;
1252   }
1253
1254   if  (rational) {
1255     BSplCLib::BuildCache(parametercache,
1256                          spanlenghtcache,
1257                          periodic,
1258                          deg,
1259                          (flatknots->Array1()),
1260                          poles->Array1(),
1261                          weights->Array1(),
1262                          cachepoles->ChangeArray1(),
1263                          cacheweights->ChangeArray1()) ;
1264   }
1265   else {
1266     BSplCLib::BuildCache(parametercache,
1267                          spanlenghtcache,
1268                          periodic,
1269                          deg,
1270                          (flatknots->Array1()),
1271                          poles->Array1(),
1272                          *((TColStd_Array1OfReal*) NULL),
1273                          cachepoles->ChangeArray1(),
1274                          *((TColStd_Array1OfReal*) NULL)) ;
1275   }
1276   validcache = 1 ;
1277 }
1278