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