src/Law/Law_BSpline.cxx

   1 // Created on: 1995-10-20
   2 // Created by: Laurent BOURESCHE
   3 // Copyright (c) 1995-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 // Cut and past sauvage depuis Geom!?!?
  18 // 14-Mar-96 : xab implemented MovePointAndTangent
  19 // 03-02-97 : pmn ->LocateU sur Periodic (PRO6963),
  20 //            bon appel a LocateParameter (PRO6973) et mise en conformite avec
  21 //            le cdl de LocateU, lorsque U est un noeud (PRO6988)
  22
  23 #include <BSplCLib.hxx>
  24 #include <BSplCLib_KnotDistribution.hxx>
  25 #include <BSplCLib_MultDistribution.hxx>
  26 #include <gp.hxx>
  27 #include <Law_BSpline.hxx>
  28 #include <Standard_ConstructionError.hxx>
  29 #include <Standard_DimensionError.hxx>
  30 #include <Standard_DomainError.hxx>
  31 #include <Standard_NoSuchObject.hxx>
  32 #include <Standard_NotImplemented.hxx>
  33 #include <Standard_OutOfRange.hxx>
  34 #include <Standard_RangeError.hxx>
  35 #include <Standard_Type.hxx>
  36
  37 #define  POLES    (poles->Array1())
  38 #define  KNOTS    (knots->Array1())
  39 #define  FKNOTS   (flatknots->Array1())
  40 #define  FMULTS   (BSplCLib::NoMults())
  41
  42 //=======================================================================
  43 //function : SetPoles
  44 //purpose  :
  45 //=======================================================================
  46
  47 static void SetPoles(const TColStd_Array1OfReal& Poles,
  48                      const TColStd_Array1OfReal& Weights,
  49                      TColStd_Array1OfReal&       FP)
  50 {
  51   Standard_Integer i,j = FP.Lower();
  52   for (i = Poles.Lower(); i <= Poles.Upper(); i++) {
  53     Standard_Real w = Weights(i);
  54     FP(j) = Poles(i) * w;
  55     j++;
  56     FP(j) = w;
  57     j++;
  58   }
  59 }
  60
  61
  62 //=======================================================================
  63 //function : GetPoles
  64 //purpose  :
  65 //=======================================================================
  66
  67 static void GetPoles(const TColStd_Array1OfReal& FP,
  68                      TColStd_Array1OfReal&       Poles,
  69                      TColStd_Array1OfReal&       Weights)
  70
  71 {
  72   Standard_Integer i,j = FP.Lower();
  73   for (i = Poles.Lower(); i <= Poles.Upper(); i++) {
  74     Standard_Real w = FP(j+1);
  75     Weights(i) = w;
  76     Poles(i) = FP(j) /w;
  77     j+=2;
  78   }
  79 }
  80
  81
  82 //=======================================================================
  83 //function : CheckCurveData
  84 //purpose  : Internal use only
  85 //=======================================================================
  86
  87 static void CheckCurveData
  88 (const TColStd_Array1OfReal&       CPoles,
  89  const TColStd_Array1OfReal&       CKnots,
  90  const TColStd_Array1OfInteger&    CMults,
  91  const Standard_Integer            Degree,
  92  const Standard_Boolean            Periodic)
  93 {
  94   if (Degree < 1 || Degree > Law_BSpline::MaxDegree()) {
  95     Standard_ConstructionError::Raise();
  96   }
  97
  98   if (CPoles.Length() < 2)                Standard_ConstructionError::Raise();
  99   if (CKnots.Length() != CMults.Length()) Standard_ConstructionError::Raise();
 100
 101   for (Standard_Integer I = CKnots.Lower(); I < CKnots.Upper(); I++) {
 102     if (CKnots (I+1) - CKnots (I) <= Epsilon (Abs(CKnots (I)))) {
 103       Standard_ConstructionError::Raise();
 104     }
 105   }
 106
 107   if (CPoles.Length() != BSplCLib::NbPoles(Degree,Periodic,CMults))
 108     Standard_ConstructionError::Raise();
 109 }
 110
 111
 112 //=======================================================================
 113 //function : KnotAnalysis
 114 //purpose  : Internal use only
 115 //=======================================================================
 116
 117 static void KnotAnalysis
 118 (const Standard_Integer           Degree,
 119  const Standard_Boolean           Periodic,
 120  const TColStd_Array1OfReal&      CKnots,
 121  const TColStd_Array1OfInteger&   CMults,
 122  GeomAbs_BSplKnotDistribution&    KnotForm,
 123  Standard_Integer&                MaxKnotMult)
 124 {
 125   KnotForm = GeomAbs_NonUniform;
 126
 127   BSplCLib_KnotDistribution KSet =
 128     BSplCLib::KnotForm (CKnots, 1, CKnots.Length());
 129
 130
 131   if (KSet == BSplCLib_Uniform) {
 132     BSplCLib_MultDistribution MSet =
 133       BSplCLib::MultForm (CMults, 1, CMults.Length());
 134     switch (MSet) {
 135     case BSplCLib_NonConstant   :
 136       break;
 137     case BSplCLib_Constant      :
 138       if (CKnots.Length() == 2) {
 139         KnotForm = GeomAbs_PiecewiseBezier;
 140       }
 141       else {
 142         if (CMults (1) == 1)  KnotForm = GeomAbs_Uniform;
 143       }
 144       break;
 145     case BSplCLib_QuasiConstant :
 146       if (CMults (1) == Degree + 1) {
 147         Standard_Real M = CMults (2);
 148         if (M == Degree )   KnotForm = GeomAbs_PiecewiseBezier;
 149         else if  (M == 1)   KnotForm = GeomAbs_QuasiUniform;
 150       }
 151       break;
 152     }
 153   }
 154
 155   Standard_Integer FirstKM =
 156     Periodic ? CKnots.Lower() :  BSplCLib::FirstUKnotIndex (Degree,CMults);
 157   Standard_Integer LastKM =
 158     Periodic ? CKnots.Upper() :  BSplCLib::LastUKnotIndex (Degree,CMults);
 159   MaxKnotMult = 0;
 160   if (LastKM - FirstKM != 1) {
 161     Standard_Integer Multi;
 162     for (Standard_Integer i = FirstKM + 1; i < LastKM; i++) {
 163       Multi = CMults (i);
 164       MaxKnotMult = Max (MaxKnotMult, Multi);
 165     }
 166   }
 167 }
 168
 169
 170 //=======================================================================
 171 //function : Rational
 172 //purpose  : check rationality of an array of weights
 173 //=======================================================================
 174
 175 static Standard_Boolean Rational(const TColStd_Array1OfReal& W)
 176 {
 177   Standard_Integer i, n = W.Length();
 178   Standard_Boolean rat = Standard_False;
 179   for (i = 1; i < n; i++) {
 180     rat =  Abs(W(i) - W(i+1)) > gp::Resolution();
 181     if (rat) break;
 182   }
 183   return rat;
 184 }
 185
 186 //=======================================================================
 187 //function : Copy
 188 //purpose  :
 189 //=======================================================================
 190
 191 Handle(Law_BSpline) Law_BSpline::Copy() const
 192 {
 193   Handle(Law_BSpline) C;
 194   if (IsRational())
 195     C = new Law_BSpline(poles->Array1(),
 196                         weights->Array1(),
 197                         knots->Array1(),
 198                         mults->Array1(),
 199                         deg,periodic);
 200   else
 201     C = new Law_BSpline(poles->Array1(),
 202                         knots->Array1(),
 203                         mults->Array1(),
 204                         deg,periodic);
 205   return C;
 206 }
 207
 208
 209
 210 //=======================================================================
 211 //function : Law_BSpline
 212 //purpose  :
 213 //=======================================================================
 214
 215 Law_BSpline::Law_BSpline
 216 (const TColStd_Array1OfReal&       Poles,
 217  const TColStd_Array1OfReal&     Knots,
 218  const TColStd_Array1OfInteger&  Mults,
 219  const Standard_Integer          Degree,
 220  const Standard_Boolean          Periodic) :
 221  rational(Standard_False),periodic(Periodic), deg(Degree)
 222 {
 223   // check
 224
 225   CheckCurveData (Poles,
 226                   Knots,
 227                   Mults,
 228                   Degree,
 229                   Periodic);
 230
 231
 232   // copy arrays
 233
 234   poles =  new TColStd_HArray1OfReal(1,Poles.Length());
 235   poles->ChangeArray1() = Poles;
 236
 237
 238   knots = new TColStd_HArray1OfReal(1,Knots.Length());
 239   knots->ChangeArray1() = Knots;
 240
 241   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
 242   mults->ChangeArray1() = Mults;
 243
 244   UpdateKnots();
 245 }
 246
 247
 248
 249 //=======================================================================
 250 //function : Law_BSpline
 251 //purpose  :
 252 //=======================================================================
 253
 254 Law_BSpline::Law_BSpline
 255 (const TColStd_Array1OfReal&      Poles,
 256  const TColStd_Array1OfReal&    Weights,
 257  const TColStd_Array1OfReal&    Knots,
 258  const TColStd_Array1OfInteger& Mults,
 259  const Standard_Integer         Degree,
 260  const Standard_Boolean         Periodic)  :
 261  rational(Standard_True),  periodic(Periodic), deg(Degree)
 262
 263 {
 264
 265   // check
 266
 267   CheckCurveData (Poles,
 268                   Knots,
 269                   Mults,
 270                   Degree,
 271                   Periodic);
 272
 273   if (Weights.Length() != Poles.Length())
 274     Standard_ConstructionError::Raise("Law_BSpline");
 275
 276   Standard_Integer i;
 277   for (i = Weights.Lower(); i <= Weights.Upper(); i++) {
 278     if (Weights(i) <= gp::Resolution())
 279       Standard_ConstructionError::Raise("Law_BSpline");
 280   }
 281
 282   // check really rational
 283   rational = Rational(Weights);
 284
 285   // copy arrays
 286
 287   poles =  new TColStd_HArray1OfReal(1,Poles.Length());
 288   poles->ChangeArray1() = Poles;
 289   if (rational) {
 290     weights =  new TColStd_HArray1OfReal(1,Weights.Length());
 291     weights->ChangeArray1() = Weights;
 292   }
 293
 294   knots = new TColStd_HArray1OfReal(1,Knots.Length());
 295   knots->ChangeArray1() = Knots;
 296
 297   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
 298   mults->ChangeArray1() = Mults;
 299
 300   UpdateKnots();
 301 }
 302
 303
 304 //=======================================================================
 305 //function : MaxDegree
 306 //purpose  :
 307 //=======================================================================
 308
 309 Standard_Integer Law_BSpline::MaxDegree ()
 310 {
 311   return BSplCLib::MaxDegree();
 312 }
 313
 314
 315 //=======================================================================
 316 //function : IncreaseDegree
 317 //purpose  :
 318 //=======================================================================
 319
 320 void Law_BSpline::IncreaseDegree  (const Standard_Integer Degree)
 321 {
 322   if (Degree == deg) return;
 323
 324   if (Degree < deg || Degree > Law_BSpline::MaxDegree()) {
 325     Standard_ConstructionError::Raise();
 326   }
 327
 328   Standard_Integer FromK1 = FirstUKnotIndex ();
 329   Standard_Integer ToK2   = LastUKnotIndex  ();
 330
 331   Standard_Integer Step   = Degree - deg;
 332
 333   Handle(TColStd_HArray1OfReal) npoles = new
 334     TColStd_HArray1OfReal(1,poles->Length() + Step * (ToK2-FromK1));
 335
 336   Standard_Integer nbknots = BSplCLib::IncreaseDegreeCountKnots
 337     (deg,Degree,periodic,mults->Array1());
 338
 339   Handle(TColStd_HArray1OfReal) nknots =
 340     new TColStd_HArray1OfReal(1,nbknots);
 341
 342   Handle(TColStd_HArray1OfInteger) nmults =
 343     new TColStd_HArray1OfInteger(1,nbknots);
 344
 345   Handle(TColStd_HArray1OfReal) nweights;
 346
 347   if (IsRational()) {
 348     nweights = new TColStd_HArray1OfReal(1,npoles->Upper());
 349     TColStd_Array1OfReal adimpol(1,2*poles->Upper());
 350     SetPoles(poles->Array1(),weights->Array1(),adimpol);
 351     TColStd_Array1OfReal adimnpol(1,2*npoles->Upper());
 352     BSplCLib::IncreaseDegree
 353       (deg,Degree, periodic,2,adimpol,
 354        knots->Array1(),mults->Array1(),adimnpol,
 355        nknots->ChangeArray1(),nmults->ChangeArray1());
 356     GetPoles(adimnpol,npoles->ChangeArray1(),nweights->ChangeArray1());
 357   }
 358   else {
 359     BSplCLib::IncreaseDegree
 360       (deg,Degree, periodic,1,poles->Array1(),
 361        knots->Array1(),mults->Array1(),npoles->ChangeArray1(),
 362        nknots->ChangeArray1(),nmults->ChangeArray1());
 363   }
 364
 365   deg     = Degree;
 366   poles   = npoles;
 367   weights = nweights;
 368   knots   = nknots;
 369   mults   = nmults;
 370   UpdateKnots();
 371
 372 }
 373
 374
 375 //=======================================================================
 376 //function : IncreaseMultiplicity
 377 //purpose  :
 378 //=======================================================================
 379
 380 void Law_BSpline::IncreaseMultiplicity  (const Standard_Integer Index,
 381                                          const Standard_Integer M)
 382 {
 383   TColStd_Array1OfReal k(1,1);
 384   k(1) = knots->Value(Index);
 385   TColStd_Array1OfInteger m(1,1);
 386   m(1) = M - mults->Value(Index);
 387   InsertKnots(k,m,Epsilon(1.));
 388 }
 389
 390
 391 //=======================================================================
 392 //function : IncreaseMultiplicity
 393 //purpose  :
 394 //=======================================================================
 395
 396 void Law_BSpline::IncreaseMultiplicity  (const Standard_Integer I1,
 397                                          const Standard_Integer I2,
 398                                          const Standard_Integer M)
 399 {
 400   Handle(TColStd_HArray1OfReal)  tk = knots;
 401   TColStd_Array1OfReal k((knots->Array1())(I1),I1,I2);
 402   TColStd_Array1OfInteger m(I1,I2);
 403   Standard_Integer i;
 404   for (i = I1; i <= I2; i++)
 405     m(i) = M - mults->Value(i);
 406   InsertKnots(k,m,Epsilon(1.));
 407 }
 408
 409 //=======================================================================
 410 //function : IncrementMultiplicity
 411 //purpose  :
 412 //=======================================================================
 413
 414 void Law_BSpline::IncrementMultiplicity
 415 (const Standard_Integer I1,
 416  const Standard_Integer I2,
 417  const Standard_Integer Step)
 418 {
 419   Handle(TColStd_HArray1OfReal) tk = knots;
 420   TColStd_Array1OfReal    k((knots->Array1())(I1),I1,I2);
 421   TColStd_Array1OfInteger m(I1,I2) ;
 422   m.Init(Step);
 423   InsertKnots(k,m,Epsilon(1.));
 424 }
 425
 426
 427 //=======================================================================
 428 //function : InsertKnot
 429 //purpose  :
 430 //=======================================================================
 431
 432 void Law_BSpline::InsertKnot
 433 (const Standard_Real U,
 434  const Standard_Integer M,
 435  const Standard_Real ParametricTolerance,
 436  const Standard_Boolean Add)
 437 {
 438   TColStd_Array1OfReal k(1,1);
 439   k(1) = U;
 440   TColStd_Array1OfInteger m(1,1);
 441   m(1) = M;
 442   InsertKnots(k,m,ParametricTolerance,Add);
 443 }
 444
 445 //=======================================================================
 446 //function : InsertKnots
 447 //purpose  :
 448 //=======================================================================
 449
 450 void  Law_BSpline::InsertKnots(const TColStd_Array1OfReal& Knots,
 451                                const TColStd_Array1OfInteger& Mults,
 452                                const Standard_Real Epsilon,
 453                                const Standard_Boolean Add)
 454 {
 455   // Check and compute new sizes
 456   Standard_Integer nbpoles,nbknots;
 457
 458   if (!BSplCLib::PrepareInsertKnots(deg,periodic,
 459                                     knots->Array1(),mults->Array1(),
 460                                     Knots,Mults,nbpoles,nbknots,Epsilon,Add))
 461     Standard_ConstructionError::Raise("Law_BSpline::InsertKnots");
 462
 463   if (nbpoles == poles->Length()) return;
 464
 465   Handle(TColStd_HArray1OfReal) npoles = new TColStd_HArray1OfReal(1,nbpoles);
 466   Handle(TColStd_HArray1OfReal) nknots = knots;
 467   Handle(TColStd_HArray1OfInteger) nmults = mults;
 468
 469   if (nbknots != knots->Length()) {
 470     nknots = new TColStd_HArray1OfReal(1,nbknots);
 471     nmults = new TColStd_HArray1OfInteger(1,nbknots);
 472   }
 473
 474   if (rational) {
 475     Handle(TColStd_HArray1OfReal) nweights =
 476       new TColStd_HArray1OfReal(1,nbpoles);
 477     TColStd_Array1OfReal adimpol(1,2*poles->Upper());
 478     SetPoles(poles->Array1(),weights->Array1(),adimpol);
 479     TColStd_Array1OfReal adimnpol(1,2*npoles->Upper());
 480     BSplCLib::InsertKnots(deg,periodic,2,adimpol,
 481                           knots->Array1(), mults->Array1(),
 482                           Knots, Mults,adimnpol,
 483                           nknots->ChangeArray1(), nmults->ChangeArray1(),
 484                           Epsilon, Add);
 485     GetPoles(adimnpol,npoles->ChangeArray1(),nweights->ChangeArray1());
 486     weights = nweights;
 487   }
 488   else {
 489     BSplCLib::InsertKnots(deg,periodic,1,poles->Array1(),
 490                           knots->Array1(), mults->Array1(),
 491                           Knots, Mults,
 492                           npoles->ChangeArray1(),
 493                           nknots->ChangeArray1(), nmults->ChangeArray1(),
 494                           Epsilon, Add);
 495   }
 496
 497   poles = npoles;
 498   knots = nknots;
 499   mults = nmults;
 500   UpdateKnots();
 501
 502 }
 503
 504 //=======================================================================
 505 //function : RemoveKnot
 506 //purpose  :
 507 //=======================================================================
 508
 509 Standard_Boolean  Law_BSpline::RemoveKnot(const Standard_Integer Index,
 510                                           const Standard_Integer M,
 511                                           const Standard_Real Tolerance)
 512 {
 513   if (M < 0) return Standard_True;
 514
 515   Standard_Integer I1  = FirstUKnotIndex ();
 516   Standard_Integer I2  = LastUKnotIndex  ();
 517
 518   if ( !periodic && (Index <= I1 || Index >= I2) ) {
 519     Standard_OutOfRange::Raise();
 520   }
 521   else if ( periodic  && (Index < I1 || Index > I2)) {
 522     Standard_OutOfRange::Raise();
 523   }
 524
 525   const TColStd_Array1OfReal   & oldpoles   = poles->Array1();
 526   Standard_Integer step = mults->Value(Index) - M;
 527   if (step <= 0) return Standard_True;
 528
 529   Handle(TColStd_HArray1OfReal) npoles =
 530     new TColStd_HArray1OfReal(1,oldpoles.Length()-step);
 531
 532   Handle(TColStd_HArray1OfReal)    nknots  = knots;
 533   Handle(TColStd_HArray1OfInteger) nmults  = mults;
 534
 535   if (M == 0) {
 536     nknots = new TColStd_HArray1OfReal(1,knots->Length()-1);
 537     nmults = new TColStd_HArray1OfInteger(1,knots->Length()-1);
 538   }
 539
 540   if (IsRational()) {
 541     Handle(TColStd_HArray1OfReal) nweights =
 542       new TColStd_HArray1OfReal(1,npoles->Length());
 543     TColStd_Array1OfReal adimpol(1,2*poles->Upper());
 544     SetPoles(poles->Array1(),weights->Array1(),adimpol);
 545     TColStd_Array1OfReal adimnpol(1,2*npoles->Upper());
 546     if (!BSplCLib::RemoveKnot(Index, M, deg, periodic,2,adimpol,
 547                               knots->Array1(),mults->Array1(),adimnpol,
 548                               nknots->ChangeArray1(),
 549                               nmults->ChangeArray1(),Tolerance))
 550       return Standard_False;
 551     GetPoles(adimnpol,npoles->ChangeArray1(),nweights->ChangeArray1());
 552     weights = nweights;
 553   }
 554   else {
 555     if (!BSplCLib::RemoveKnot(Index, M, deg, periodic,1,poles->Array1(),
 556                               knots->Array1(),mults->Array1(),
 557                               npoles->ChangeArray1(), nknots->ChangeArray1(),
 558                               nmults->ChangeArray1(),Tolerance))
 559       return Standard_False;
 560   }
 561
 562   poles = npoles;
 563   knots = nknots;
 564   mults = nmults;
 565
 566   UpdateKnots();
 567   return Standard_True;
 568 }
 569
 570 //----------------------------------------------------------------------
 571 //----------------------------------------------------------------------
 572
 573 # if 0
 574
 575 --- methodes otees du CDL -> spec trop vagues : on ne sait pas ou rajouter
 576 le noeud
 577
 578 //=======================================================================
 579 //function : InsertPoleAfter
 580 //purpose  :
 581 //=======================================================================
 582
 583 void Law_BSpline::InsertPoleAfter
 584 (const Standard_Integer Index,
 585  const Standard_Real& P)
 586 {
 587   InsertPoleAfter(Index,P,1.);
 588 }
 589
 590
 591
 592 //=======================================================================
 593 //function : InsertPoleAfter
 594 //purpose  :
 595 //=======================================================================
 596
 597 void Law_BSpline::InsertPoleAfter
 598 (const Standard_Integer Index,
 599  const Standard_Real& P,
 600  const Standard_Real Weight)
 601 {
 602   if (Index < 0 || Index > poles->Length())  Standard_OutOfRange::Raise();
 603
 604   if (Weight <= gp::Resolution())     Standard_ConstructionError::Raise();
 605
 606
 607   // find the spans which are modified with the inserting pole
 608   //   --> evaluate NewKnot & KnotIndex : Value of the new knot to insert.
 609   Standard_Integer KnotIndex, k, sigma;
 610   Standard_Real    NewKnot;
 611
 612   if (periodic) {
 613     sigma = 0;
 614     k     = 1;
 615     while ( sigma < Index) {
 616       sigma += mults->Value(k);
 617       k++;
 618     }
 619     KnotIndex = k - 1;
 620     NewKnot   = ( knots->Value(KnotIndex) + knots->Value(KnotIndex+1)) / 2.;
 621   }
 622   else {
 623     sigma = 0;
 624     k     = 1;
 625     while ( sigma < Index) {
 626       sigma += mults->Value(k);
 627       k++;
 628     }
 629     Standard_Integer first = k - 1;
 630     sigma -= Index;
 631     while ( sigma < (deg+1)) {
 632       sigma += mults->Value(k);
 633       k++;
 634     }
 635     Standard_Integer last = k - 1;
 636
 637     KnotIndex = first + (( last - first) / 2);
 638     NewKnot = ( knots->Value(KnotIndex) + knots->Value(KnotIndex+1)) / 2.;
 639   }
 640
 641   Standard_Integer nbknots = knots->Length();
 642   Handle(TColStd_HArray1OfReal) nknots =
 643     new TColStd_HArray1OfReal(1,nbknots+1);
 644   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
 645   Handle(TColStd_HArray1OfInteger) nmults =
 646     new TColStd_HArray1OfInteger(1,nbknots+1);
 647   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
 648
 649   // insert the knot
 650
 651   Standard_Integer i;
 652   for ( i = 1; i<= KnotIndex; i++) {
 653     newknots(i) = knots->Value(i);
 654     newmults(i) = mults->Value(i);
 655   }
 656   newknots(KnotIndex+1) = NewKnot;
 657   newmults(KnotIndex+1) = 1;
 658   for ( i = KnotIndex+1; i <= nbknots; i++) {
 659     newknots(i+1) = knots->Value(i);
 660     newmults(i+1) = mults->Value(i);
 661   }
 662
 663   Standard_Integer nbpoles = poles->Length();
 664   Handle(TColStd_HArray1OfReal) npoles =
 665     new TColStd_HArray1OfReal(1,nbpoles+1);
 666   TColStd_Array1OfReal& newpoles = npoles->ChangeArray1();
 667
 668   // insert the pole
 669
 670   for (i = 1; i <= Index; i++)
 671     newpoles(i) = poles->Value(i);
 672
 673   newpoles(Index+1) = P;
 674
 675   for (i = Index+1; i <= nbpoles; i++)
 676     newpoles(i+1) = poles->Value(i);
 677
 678   // insert the weight
 679
 680   Handle(TColStd_HArray1OfReal) nweights;
 681   Standard_Boolean rat = IsRational() || Abs(Weight-1.) > gp::Resolution();
 682
 683   if (rat) {
 684     nweights = new TColStd_HArray1OfReal(1,nbpoles+1);
 685     TColStd_Array1OfReal& newweights = nweights->ChangeArray1();
 686
 687     for (i = 1; i <= Index; i++)
 688       if (IsRational())
 689         newweights(i) = weights->Value(i);
 690       else
 691         newweights(i) = 1.;
 692
 693     newweights(Index+1) = Weight;
 694
 695     for (i = Index+1; i <= nbpoles; i++)
 696       if (IsRational())
 697         newweights(i+1) = weights->Value(i);
 698       else
 699         newweights(i+1) = 1.;
 700   }
 701
 702   poles   = npoles;
 703   weights = nweights;
 704   knots   = nknots;
 705   mults   = nmults;
 706   UpdateKnots();
 707 }
 708
 709
 710 //=======================================================================
 711 //function : InsertPoleBefore
 712 //purpose  :
 713 //=======================================================================
 714
 715 void Law_BSpline::InsertPoleBefore
 716 (const Standard_Integer Index,
 717  const Standard_Real& P )
 718 {
 719   InsertPoleAfter(Index-1,P,1.);
 720 }
 721
 722
 723
 724 //=======================================================================
 725 //function : InsertPoleBefore
 726 //purpose  :
 727 //=======================================================================
 728
 729 void Law_BSpline::InsertPoleBefore
 730 (const Standard_Integer Index,
 731  const Standard_Real& P,
 732  const Standard_Real Weight)
 733 {
 734   InsertPoleAfter(Index-1,P,Weight);
 735 }
 736
 737 //=======================================================================
 738 //function : RemovePole
 739 //purpose  :
 740 //=======================================================================
 741
 742 void Law_BSpline::RemovePole
 743 (const Standard_Integer Index)
 744 {
 745   if (Index < 1 || Index > poles->Length())  Standard_OutOfRange::Raise();
 746
 747   if (poles->Length() <= 2)           Standard_ConstructionError::Raise();
 748
 749   if (knotSet == GeomAbs_NonUniform || knotSet == GeomAbs_PiecewiseBezier)
 750     Standard_ConstructionError::Raise();
 751
 752   Standard_Integer i;
 753   Handle(TColStd_HArray1OfReal) nknots =
 754     new TColStd_HArray1OfReal(1,knots->Length()-1);
 755   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
 756
 757   Handle(TColStd_HArray1OfInteger) nmults =
 758     new TColStd_HArray1OfInteger(1,mults->Length()-1);
 759   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
 760
 761   for (i = 1; i < newknots.Length(); i++) {
 762     newknots (i) = knots->Value (i);
 763     newmults (i) = 1;
 764   }
 765   newmults(1) = mults->Value(1);
 766   newknots(newknots.Upper()) = knots->Value (knots->Upper());
 767   newmults(newmults.Upper()) = mults->Value (mults->Upper());
 768
 769
 770   Handle(TColStd_HArray1OfReal) npoles =
 771     new TColStd_HArray1OfReal(1, poles->Upper()-1);
 772   TColStd_Array1OfReal& newpoles = npoles->ChangeArray1();
 773
 774   for (i = 1; i < Index; i++)
 775     newpoles(i) = poles->Value(i);
 776   for (i = Index; i < newpoles.Length(); i++)
 777     newpoles(i) = poles->Value(i+1);
 778
 779   Handle(TColStd_HArray1OfReal) nweights;
 780   if (IsRational()) {
 781     nweights = new TColStd_HArray1OfReal(1,newpoles.Length());
 782     TColStd_Array1OfReal& newweights = nweights->ChangeArray1();
 783     for (i = 1; i < Index; i++)
 784       newweights(i) = weights->Value(i);
 785     for (i = Index; i < newweights.Length(); i++)
 786       newweights(i) = weights->Value(i+1);
 787   }
 788
 789   poles   = npoles;
 790   weights = nweights;
 791   knots   = nknots;
 792   mults   = nmults;
 793   UpdateKnots();
 794 }
 795
 796
 797
 798
 799 #endif
 800
 801
 802 //=======================================================================
 803 //function : Reverse
 804 //purpose  :
 805 //=======================================================================
 806
 807 void Law_BSpline::Reverse ()
 808 {
 809   BSplCLib::Reverse(knots->ChangeArray1());
 810   BSplCLib::Reverse(mults->ChangeArray1());
 811   Standard_Integer last;
 812   if (periodic)
 813     last = flatknots->Upper() - deg - 1;
 814   else
 815     last = poles->Upper();
 816   BSplCLib::Reverse(poles->ChangeArray1(),last);
 817   if (rational)
 818     BSplCLib::Reverse(weights->ChangeArray1(),last);
 819   UpdateKnots();
 820 }
 821
 822
 823 //=======================================================================
 824 //function : ReversedParameter
 825 //purpose  :
 826 //=======================================================================
 827
 828 Standard_Real Law_BSpline::ReversedParameter
 829 (const Standard_Real U) const
 830 {
 831   return (FirstParameter() + LastParameter() - U);
 832 }
 833
 834
 835 //=======================================================================
 836 //function : Segment
 837 //purpose  :
 838 //=======================================================================
 839
 840 void Law_BSpline::Segment(const Standard_Real U1,
 841                           const Standard_Real U2)
 842 {
 843   Standard_DomainError_Raise_if ( U2 < U1,
 844                                  "Law_BSpline::Segment");
 845   Standard_Real Eps = Epsilon(Max(Abs(U1),Abs(U2)));
 846   Standard_Real delta = U2 - U1;
 847
 848   Standard_Real NewU1, NewU2;
 849   Standard_Real U;
 850   Standard_Integer index;
 851
 852   TColStd_Array1OfReal    Knots(1,2);
 853   TColStd_Array1OfInteger Mults(1,2);
 854
 855   index = 0;
 856   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 857                             U1,periodic,knots->Lower(),knots->Upper(),
 858                             index,NewU1);
 859   index = 0;
 860   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 861                             U2,periodic,knots->Lower(),knots->Upper(),
 862                             index,NewU2);
 863   Knots( 1) = Min( NewU1, NewU2);
 864   Knots( 2) = Max( NewU1, NewU2);
 865   Mults( 1) = Mults( 2) = deg;
 866   InsertKnots( Knots, Mults, Eps);
 867
 868   if (periodic) { // set the origine at NewU1
 869     Standard_Integer index0 = 0;
 870     BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 871                               U1,periodic,knots->Lower(),knots->Upper(),
 872                               index0,U);
 873     if ( Abs(knots->Value(index0+1)-U) < Eps)
 874       index0++;
 875     SetOrigin(index0);
 876     SetNotPeriodic();
 877   }
 878
 879   // compute index1 and index2 to set the new knots and mults
 880   Standard_Integer index1 = 0, index2 = 0;
 881   Standard_Integer FromU1 = knots->Lower();
 882   Standard_Integer ToU2   = knots->Upper();
 883   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 884                             NewU1,periodic,FromU1,ToU2,index1,U);
 885   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
 886                             NewU1 + delta,periodic,FromU1,ToU2,index2,U);
 887   if ( Abs(knots->Value(index2+1)-U) < Eps)
 888     index2++;
 889
 890   Standard_Integer nbknots = index2 - index1 + 1;
 891
 892   Handle(TColStd_HArray1OfReal)
 893     nknots = new TColStd_HArray1OfReal(1,nbknots);
 894   Handle(TColStd_HArray1OfInteger)
 895     nmults = new TColStd_HArray1OfInteger(1,nbknots);
 896
 897   Standard_Integer i , k = 1;
 898   for ( i = index1; i<= index2; i++) {
 899     nknots->SetValue(k, knots->Value(i));
 900     nmults->SetValue(k, mults->Value(i));
 901     k++;
 902   }
 903   nmults->SetValue(      1, deg + 1);
 904   nmults->SetValue(nbknots, deg + 1);
 905
 906
 907   // compute index1 and index2 to set the new poles and weights
 908   Standard_Integer pindex1
 909     = BSplCLib::PoleIndex(deg,index1,periodic,mults->Array1());
 910   Standard_Integer pindex2
 911     = BSplCLib::PoleIndex(deg,index2,periodic,mults->Array1());
 912
 913   pindex1++;
 914   pindex2 = Min( pindex2+1, poles->Length());
 915
 916   Standard_Integer nbpoles  = pindex2 - pindex1 + 1;
 917
 918   Handle(TColStd_HArray1OfReal)
 919     nweights = new TColStd_HArray1OfReal(1,nbpoles);
 920   Handle(TColStd_HArray1OfReal)
 921     npoles = new TColStd_HArray1OfReal(1,nbpoles);
 922
 923   k = 1;
 924   if ( rational) {
 925     nweights = new TColStd_HArray1OfReal( 1, nbpoles);
 926     for ( i = pindex1; i <= pindex2; i++) {
 927       npoles->SetValue(k, poles->Value(i));
 928       nweights->SetValue(k, weights->Value(i));
 929       k++;
 930     }
 931   }
 932   else {
 933     for ( i = pindex1; i <= pindex2; i++) {
 934       npoles->SetValue(k, poles->Value(i));
 935       k++;
 936     }
 937   }
 938
 939   knots = nknots;
 940   mults = nmults;
 941   poles = npoles;
 942   if ( rational)
 943     weights = nweights;
 944
 945   UpdateKnots();
 946 }
 947
 948
 949 //=======================================================================
 950 //function : SetKnot
 951 //purpose  :
 952 //=======================================================================
 953
 954 void Law_BSpline::SetKnot
 955 (const Standard_Integer Index,
 956  const Standard_Real K)
 957 {
 958   if (Index < 1 || Index > knots->Length())     Standard_OutOfRange::Raise();
 959   Standard_Real DK = Abs(Epsilon (K));
 960   if (Index == 1) {
 961     if (K >= knots->Value(2) - DK) Standard_ConstructionError::Raise();
 962   }
 963   else if (Index == knots->Length()) {
 964     if (K <= knots->Value (knots->Length()-1) + DK)  {
 965       Standard_ConstructionError::Raise();
 966     }
 967   }
 968   else {
 969     if (K <= knots->Value(Index-1) + DK ||
 970         K >= knots->Value(Index+1) - DK ) {
 971       Standard_ConstructionError::Raise();
 972     }
 973   }
 974   if (K != knots->Value (Index)) {
 975     knots->SetValue (Index, K);
 976     UpdateKnots();
 977   }
 978 }
 979
 980
 981 //=======================================================================
 982 //function : SetKnots
 983 //purpose  :
 984 //=======================================================================
 985
 986 void Law_BSpline::SetKnots
 987 (const TColStd_Array1OfReal& K)
 988 {
 989   CheckCurveData(poles->Array1(),K,mults->Array1(),deg,periodic);
 990   knots->ChangeArray1() = K;
 991   UpdateKnots();
 992 }
 993
 994
 995 //=======================================================================
 996 //function : SetKnot
 997 //purpose  :
 998 //=======================================================================
 999
1000 void Law_BSpline::SetKnot
1001 (const Standard_Integer Index,
1002  const Standard_Real K,
1003  const Standard_Integer M)
1004 {
1005   IncreaseMultiplicity (Index, M);
1006   SetKnot (Index, K);
1007 }
1008
1009
1010 //=======================================================================
1011 //function : SetPeriodic
1012 //purpose  :
1013 //=======================================================================
1014
1015 void Law_BSpline::SetPeriodic ()
1016 {
1017   Standard_Integer first = FirstUKnotIndex();
1018   Standard_Integer last  = LastUKnotIndex();
1019
1020   Handle(TColStd_HArray1OfReal) tk = knots;
1021   TColStd_Array1OfReal    cknots((knots->Array1())(first),first,last);
1022   knots = new TColStd_HArray1OfReal(1,cknots.Length());
1023   knots->ChangeArray1() = cknots;
1024
1025   Handle(TColStd_HArray1OfInteger) tm = mults;
1026   TColStd_Array1OfInteger cmults((mults->Array1())(first),first,last);
1027   cmults(first) = cmults(last) = Max( cmults(first), cmults(last));
1028   mults = new TColStd_HArray1OfInteger(1,cmults.Length());
1029   mults->ChangeArray1() = cmults;
1030
1031   // compute new number of poles;
1032   Standard_Integer nbp = BSplCLib::NbPoles(deg,Standard_True,cmults);
1033
1034   Handle(TColStd_HArray1OfReal) tp = poles;
1035   TColStd_Array1OfReal cpoles((poles->Array1())(1),1,nbp);
1036   poles = new TColStd_HArray1OfReal(1,nbp);
1037   poles->ChangeArray1() = cpoles;
1038
1039   if (rational) {
1040     Handle(TColStd_HArray1OfReal) tw = weights;
1041     TColStd_Array1OfReal cweights((weights->Array1())(1),1,nbp);
1042     weights = new TColStd_HArray1OfReal(1,nbp);
1043     weights->ChangeArray1() = cweights;
1044   }
1045
1046   periodic = Standard_True;
1047
1048   UpdateKnots();
1049 }
1050
1051
1052 //=======================================================================
1053 //function : SetOrigin
1054 //purpose  :
1055 //=======================================================================
1056
1057 void Law_BSpline::SetOrigin(const Standard_Integer Index)
1058 {
1059   Standard_NoSuchObject_Raise_if( !periodic,
1060                                  "Law_BSpline::SetOrigin");
1061   Standard_Integer i,k;
1062   Standard_Integer first = FirstUKnotIndex();
1063   Standard_Integer last  = LastUKnotIndex();
1064
1065   Standard_DomainError_Raise_if( (Index < first) || (Index > last),
1066                                 "Law_BSpline::SetOrigine");
1067
1068   Standard_Integer nbknots = knots->Length();
1069   Standard_Integer nbpoles = poles->Length();
1070
1071   Handle(TColStd_HArray1OfReal) nknots =
1072     new TColStd_HArray1OfReal(1,nbknots);
1073   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
1074
1075   Handle(TColStd_HArray1OfInteger) nmults =
1076     new TColStd_HArray1OfInteger(1,nbknots);
1077   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
1078
1079   // set the knots and mults
1080   Standard_Real period = knots->Value(last) - knots->Value(first);
1081   k = 1;
1082   for ( i = Index; i <= last ; i++) {
1083     newknots(k) = knots->Value(i);
1084     newmults(k) = mults->Value(i);
1085     k++;
1086   }
1087   for ( i = first+1; i <= Index; i++) {
1088     newknots(k) = knots->Value(i) + period;
1089     newmults(k) = mults->Value(i);
1090     k++;
1091   }
1092
1093   Standard_Integer index = 1;
1094   for (i = first+1; i <= Index; i++)
1095     index += mults->Value(i);
1096
1097   // set the poles and weights
1098   Handle(TColStd_HArray1OfReal) npoles =
1099     new TColStd_HArray1OfReal(1,nbpoles);
1100   Handle(TColStd_HArray1OfReal) nweights =
1101     new TColStd_HArray1OfReal(1,nbpoles);
1102   TColStd_Array1OfReal   & newpoles   = npoles->ChangeArray1();
1103   TColStd_Array1OfReal & newweights = nweights->ChangeArray1();
1104   first = poles->Lower();
1105   last  = poles->Upper();
1106   if ( rational) {
1107     k = 1;
1108     for ( i = index; i <= last; i++) {
1109       newpoles(k)   = poles->Value(i);
1110       newweights(k) = weights->Value(i);
1111       k++;
1112     }
1113     for ( i = first; i < index; i++) {
1114       newpoles(k)   = poles->Value(i);
1115       newweights(k) = weights->Value(i);
1116       k++;
1117     }
1118   }
1119   else {
1120     k = 1;
1121     for ( i = index; i <= last; i++) {
1122       newpoles(k) = poles->Value(i);
1123       k++;
1124     }
1125     for ( i = first; i < index; i++) {
1126       newpoles(k) = poles->Value(i);
1127       k++;
1128     }
1129   }
1130
1131   poles = npoles;
1132   knots = nknots;
1133   mults = nmults;
1134   if (rational)
1135     weights = nweights;
1136   UpdateKnots();
1137 }
1138
1139
1140
1141 //=======================================================================
1142 //function : SetNotPeriodic
1143 //purpose  :
1144 //=======================================================================
1145
1146 void Law_BSpline::SetNotPeriodic ()
1147 {
1148   if ( periodic) {
1149     Standard_Integer NbKnots, NbPoles;
1150     BSplCLib::PrepareUnperiodize( deg, mults->Array1(),NbKnots,NbPoles);
1151
1152     Handle(TColStd_HArray1OfReal) npoles
1153       = new TColStd_HArray1OfReal(1,NbPoles);
1154
1155     Handle(TColStd_HArray1OfReal) nknots
1156       = new TColStd_HArray1OfReal(1,NbKnots);
1157
1158     Handle(TColStd_HArray1OfInteger) nmults
1159       = new TColStd_HArray1OfInteger(1,NbKnots);
1160
1161     Handle(TColStd_HArray1OfReal) nweights;
1162
1163     if (IsRational()) {
1164
1165       nweights = new TColStd_HArray1OfReal(1,NbPoles);
1166
1167       TColStd_Array1OfReal adimpol(1,2*poles->Upper());
1168       SetPoles(poles->Array1(),weights->Array1(),adimpol);
1169       TColStd_Array1OfReal adimnpol(1,2*npoles->Upper());
1170       BSplCLib::Unperiodize
1171         (deg,1,mults->Array1(),knots->Array1(),adimpol,
1172          nmults->ChangeArray1(),nknots->ChangeArray1(),
1173          adimnpol);
1174       GetPoles(adimnpol,npoles->ChangeArray1(),nweights->ChangeArray1());
1175     }
1176     else {
1177
1178       BSplCLib::Unperiodize(deg,1,mults->Array1(),knots->Array1(),
1179                             poles->Array1(),nmults->ChangeArray1(),
1180                             nknots->ChangeArray1(),npoles->ChangeArray1());
1181
1182     }
1183     poles   = npoles;
1184     weights = nweights;
1185     mults   = nmults;
1186     knots   = nknots;
1187     periodic = Standard_False;
1188
1189     UpdateKnots();
1190   }
1191 }
1192
1193
1194 //=======================================================================
1195 //function : SetPole
1196 //purpose  :
1197 //=======================================================================
1198
1199 void Law_BSpline::SetPole
1200 (const Standard_Integer Index,
1201  const Standard_Real P)
1202 {
1203   if (Index < 1 || Index > poles->Length()) Standard_OutOfRange::Raise();
1204   poles->SetValue (Index, P);
1205 }
1206
1207
1208 //=======================================================================
1209 //function : SetPole
1210 //purpose  :
1211 //=======================================================================
1212
1213 void Law_BSpline::SetPole
1214 (const Standard_Integer Index,
1215  const Standard_Real P,
1216  const Standard_Real W)
1217 {
1218   SetPole(Index,P);
1219   SetWeight(Index,W);
1220 }
1221
1222 //=======================================================================
1223 //function : SetWeight
1224 //purpose  :
1225 //=======================================================================
1226
1227 void Law_BSpline::SetWeight
1228 (const Standard_Integer Index,
1229  const Standard_Real W)
1230 {
1231   if (Index < 1 || Index > poles->Length())   Standard_OutOfRange::Raise();
1232
1233   if (W <= gp::Resolution ())     Standard_ConstructionError::Raise();
1234
1235
1236   Standard_Boolean rat = IsRational() || (Abs(W - 1.) > gp::Resolution());
1237
1238   if ( rat) {
1239     if (rat && !IsRational())
1240       weights = new TColStd_HArray1OfReal(1,poles->Length(),1.);
1241
1242     weights->SetValue (Index, W);
1243
1244     if (IsRational()) {
1245       rat = Rational(weights->Array1());
1246       if (!rat) weights.Nullify();
1247     }
1248
1249     rational = !weights.IsNull();
1250   }
1251 }
1252
1253
1254 //=======================================================================
1255 //function : UpdateKnots
1256 //purpose  :
1257 //=======================================================================
1258
1259
1260 void Law_BSpline::UpdateKnots()
1261 {
1262
1263   rational = !weights.IsNull();
1264
1265   Standard_Integer MaxKnotMult = 0;
1266   KnotAnalysis (deg,
1267                 periodic,
1268                 knots->Array1(),
1269                 mults->Array1(),
1270                 knotSet, MaxKnotMult);
1271
1272   if (knotSet == GeomAbs_Uniform && !periodic)  {
1273     flatknots = knots;
1274   }
1275   else {
1276     flatknots = new TColStd_HArray1OfReal
1277       (1, BSplCLib::KnotSequenceLength(mults->Array1(),deg,periodic));
1278
1279     BSplCLib::KnotSequence (knots->Array1(),
1280                             mults->Array1(),
1281                             deg,periodic,
1282                             flatknots->ChangeArray1());
1283   }
1284
1285   if (MaxKnotMult == 0)  smooth = GeomAbs_CN;
1286   else {
1287     switch (deg - MaxKnotMult) {
1288     case 0:   smooth = GeomAbs_C0;   break;
1289     case 1:   smooth = GeomAbs_C1;   break;
1290     case 2:   smooth = GeomAbs_C2;   break;
1291     case 3:   smooth = GeomAbs_C3;   break;
1292     default:  smooth = GeomAbs_C3;   break;
1293     }
1294   }
1295 }
1296
1297
1298 //=======================================================================
1299 //function : Normalizes the parameters if the curve is periodic
1300 //purpose  : that is compute the cache so that it is valid
1301 //=======================================================================
1302
1303 void Law_BSpline::PeriodicNormalization(Standard_Real&  Parameter) const
1304 {
1305   Standard_Real Period ;
1306
1307   if (periodic){
1308     Period = flatknots->Value(flatknots->Upper() - deg) -
1309       flatknots->Value (deg + 1) ;
1310     while (Parameter > flatknots->Value(flatknots->Upper()-deg)){
1311       Parameter -= Period ;
1312     }
1313     while (Parameter < flatknots->Value((deg + 1))){
1314       Parameter +=  Period ;
1315     }
1316   }
1317 }
1318
1319
1320 //=======================================================================
1321 //function : IsCN
1322 //purpose  :
1323 //=======================================================================
1324
1325 Standard_Boolean Law_BSpline::IsCN ( const Standard_Integer N) const
1326 {
1327   Standard_RangeError_Raise_if
1328     (N < 0, "Law_BSpline::IsCN");
1329
1330   switch (smooth) {
1331   case GeomAbs_CN : return Standard_True;
1332   case GeomAbs_C0 : return N <= 0;
1333   case GeomAbs_G1 : return N <= 0;
1334   case GeomAbs_C1 : return N <= 1;
1335   case GeomAbs_G2 : return N <= 1;
1336   case GeomAbs_C2 : return N <= 2;
1337   case GeomAbs_C3 :
1338     return N <= 3 ? Standard_True :
1339            N <= deg - BSplCLib::MaxKnotMult (mults->Array1(), mults->Lower() + 1, mults->Upper() - 1);
1340   default:
1341     return Standard_False;
1342   }
1343 }
1344
1345
1346
1347 //=======================================================================
1348 //function : IsClosed
1349 //purpose  :
1350 //=======================================================================
1351
1352 Standard_Boolean Law_BSpline::IsClosed () const
1353 { return (Abs(StartPoint()-EndPoint())) <= gp::Resolution (); }
1354
1355
1356
1357 //=======================================================================
1358 //function : IsPeriodic
1359 //purpose  :
1360 //=======================================================================
1361
1362 Standard_Boolean Law_BSpline::IsPeriodic () const
1363 { return periodic; }
1364
1365 //=======================================================================
1366 //function : Continuity
1367 //purpose  :
1368 //=======================================================================
1369
1370 GeomAbs_Shape Law_BSpline::Continuity () const
1371 { return smooth; }
1372
1373 //=======================================================================
1374 //function : Degree
1375 //purpose  :
1376 //=======================================================================
1377
1378 Standard_Integer Law_BSpline::Degree () const
1379 { return deg; }
1380
1381
1382 //=======================================================================
1383 //function : Value
1384 //purpose  :
1385 //=======================================================================
1386
1387 Standard_Real Law_BSpline::Value(const Standard_Real U)const
1388 {
1389   Standard_Real  P;
1390   D0(U,P);
1391   return P;
1392 }
1393
1394
1395 //=======================================================================
1396 //function : D0
1397 //purpose  :
1398 //=======================================================================
1399
1400 void Law_BSpline::D0 (const Standard_Real U,
1401                       Standard_Real& P)  const
1402 {
1403   Standard_Real  NewU = U ;
1404   PeriodicNormalization(NewU) ;
1405   if (rational) {
1406     BSplCLib::D0(NewU,0,deg,periodic,POLES, weights->Array1(),FKNOTS,FMULTS,P);
1407   }
1408   else {
1409     BSplCLib::D0(NewU,0,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,P);
1410   }
1411 }
1412
1413
1414
1415
1416 //=======================================================================
1417 //function : D1
1418 //purpose  :
1419 //=======================================================================
1420
1421 void Law_BSpline::D1 (const Standard_Real U,
1422                       Standard_Real& P,
1423                       Standard_Real& V1) const
1424 {
1425   Standard_Real  NewU = U ;
1426   PeriodicNormalization(NewU) ;
1427   if (rational) {
1428     BSplCLib::D1(NewU,0,deg,periodic,POLES, weights->Array1(),FKNOTS,FMULTS,
1429                  P,V1) ;
1430   }
1431   else {
1432     BSplCLib::D1(NewU,0,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,
1433                  P,V1) ;
1434   }
1435 }
1436
1437
1438
1439 //=======================================================================
1440 //function : D2
1441 //purpose  :
1442 //=======================================================================
1443
1444 void Law_BSpline::D2(const Standard_Real U ,
1445                      Standard_Real& P ,
1446                      Standard_Real& V1,
1447                      Standard_Real& V2 ) const
1448 {
1449   Standard_Real  NewU = U ;
1450   PeriodicNormalization(NewU) ;
1451   if (rational) {
1452     BSplCLib::D2(NewU,0,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,
1453                  P, V1, V2) ;
1454   }
1455   else {
1456     BSplCLib::D2(NewU,0,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,
1457                  P, V1, V2) ;
1458   }
1459 }
1460
1461
1462
1463 //=======================================================================
1464 //function : D3
1465 //purpose  :
1466 //=======================================================================
1467
1468 void Law_BSpline::D3(const Standard_Real U ,
1469                      Standard_Real& P ,
1470                      Standard_Real& V1,
1471                      Standard_Real& V2,
1472                      Standard_Real& V3 ) const
1473 {
1474   Standard_Real  NewU = U ;
1475   PeriodicNormalization(NewU) ;
1476   if (rational) {
1477     BSplCLib::D3(NewU,0,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,
1478                  P, V1, V2, V3) ;
1479   }
1480   else {
1481     BSplCLib::D3(NewU,0,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,
1482                  P, V1, V2, V3) ;
1483   }
1484 }
1485
1486
1487
1488 //=======================================================================
1489 //function : DN
1490 //purpose  :
1491 //=======================================================================
1492
1493 Standard_Real Law_BSpline::DN(const Standard_Real    U,
1494                               const Standard_Integer N ) const
1495 {
1496   Standard_Real V;
1497   if (rational) {
1498     BSplCLib::DN(U,N,0,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,V);
1499   }
1500   else {
1501     BSplCLib::DN(U,N,0,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,V);
1502   }
1503   return V;
1504 }
1505
1506
1507
1508 //=======================================================================
1509 //function : EndPoint
1510 //purpose  :
1511 //=======================================================================
1512
1513 Standard_Real Law_BSpline::EndPoint () const
1514 {
1515   if (mults->Value (knots->Upper ()) == deg + 1)
1516     return poles->Value (poles->Upper());
1517   else
1518     return Value(LastParameter());
1519 }
1520
1521
1522 //=======================================================================
1523 //function : FirstUKnotIndex
1524 //purpose  :
1525 //=======================================================================
1526
1527 Standard_Integer Law_BSpline::FirstUKnotIndex () const
1528 {
1529   if (periodic) return 1;
1530   else return BSplCLib::FirstUKnotIndex (deg, mults->Array1());
1531 }
1532
1533
1534 //=======================================================================
1535 //function : FirstParameter
1536 //purpose  :
1537 //=======================================================================
1538
1539 Standard_Real Law_BSpline::FirstParameter () const
1540 {
1541   return flatknots->Value (deg+1);
1542 }
1543
1544
1545 //=======================================================================
1546 //function : Knot
1547 //purpose  :
1548 //=======================================================================
1549
1550 Standard_Real Law_BSpline::Knot (const Standard_Integer Index) const
1551 {
1552   Standard_OutOfRange_Raise_if
1553     (Index < 1 || Index > knots->Length(), "Law_BSpline::Knot");
1554   return knots->Value (Index);
1555 }
1556
1557
1558
1559 //=======================================================================
1560 //function : KnotDistribution
1561 //purpose  :
1562 //=======================================================================
1563
1564 GeomAbs_BSplKnotDistribution Law_BSpline::KnotDistribution () const
1565 {
1566   return knotSet;
1567 }
1568
1569
1570 //=======================================================================
1571 //function : Knots
1572 //purpose  :
1573 //=======================================================================
1574
1575 void Law_BSpline::Knots (TColStd_Array1OfReal& K) const
1576 {
1577   Standard_DimensionError_Raise_if
1578     (K.Length() != knots->Length(), "Law_BSpline::Knots");
1579   K = knots->Array1();
1580 }
1581
1582
1583 //=======================================================================
1584 //function : KnotSequence
1585 //purpose  :
1586 //=======================================================================
1587
1588 void Law_BSpline::KnotSequence (TColStd_Array1OfReal& K) const
1589 {
1590   Standard_DimensionError_Raise_if
1591     (K.Length() != flatknots->Length(), "Law_BSpline::KnotSequence");
1592   K = flatknots->Array1();
1593 }
1594
1595
1596
1597 //=======================================================================
1598 //function : LastUKnotIndex
1599 //purpose  :
1600 //=======================================================================
1601
1602 Standard_Integer Law_BSpline::LastUKnotIndex() const
1603 {
1604   if (periodic) return knots->Length();
1605   else return BSplCLib::LastUKnotIndex (deg, mults->Array1());
1606 }
1607
1608
1609 //=======================================================================
1610 //function : LastParameter
1611 //purpose  :
1612 //=======================================================================
1613
1614 Standard_Real Law_BSpline::LastParameter () const
1615 {
1616   return flatknots->Value (flatknots->Upper()-deg);
1617 }
1618
1619
1620 //=======================================================================
1621 //function : LocalValue
1622 //purpose  :
1623 //=======================================================================
1624
1625 Standard_Real Law_BSpline::LocalValue
1626 (const Standard_Real    U,
1627  const Standard_Integer FromK1,
1628  const Standard_Integer ToK2)   const
1629 {
1630   Standard_Real P;
1631   LocalD0(U,FromK1,ToK2,P);
1632   return P;
1633 }
1634
1635 //=======================================================================
1636 //function : LocalD0
1637 //purpose  :
1638 //=======================================================================
1639
1640 void  Law_BSpline::LocalD0
1641 (const Standard_Real    U,
1642  const Standard_Integer FromK1,
1643  const Standard_Integer ToK2,
1644  Standard_Real& P)   const
1645 {
1646   Standard_DomainError_Raise_if (FromK1 == ToK2,
1647                                  "Law_BSpline::LocalValue");
1648   Standard_Real u = U;
1649   Standard_Integer index = 0;
1650   BSplCLib::LocateParameter(deg, FKNOTS, U, periodic,FromK1,ToK2, index,u);
1651   index = BSplCLib::FlatIndex(deg,index,mults->Array1(),periodic);
1652   if (rational) {
1653     BSplCLib::D0(u,index,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,P);
1654   }
1655   else {
1656     BSplCLib::D0(u,index,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,P);
1657   }
1658 }
1659
1660 //=======================================================================
1661 //function : LocalD1
1662 //purpose  :
1663 //=======================================================================
1664
1665 void Law_BSpline::LocalD1 (const Standard_Real    U,
1666                            const Standard_Integer FromK1,
1667                            const Standard_Integer ToK2,
1668                            Standard_Real&    P,
1669                            Standard_Real&    V1)    const
1670 {
1671   Standard_DomainError_Raise_if (FromK1 == ToK2,
1672                                  "Law_BSpline::LocalD1");
1673   Standard_Real u = U;
1674   Standard_Integer index = 0;
1675   BSplCLib::LocateParameter(deg, FKNOTS, U, periodic, FromK1,ToK2, index, u);
1676   index = BSplCLib::FlatIndex(deg,index,mults->Array1(),periodic);
1677   if (rational) {
1678     BSplCLib::D1(u,index,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,P,V1);
1679   }
1680   else {
1681     BSplCLib::D1(u,index,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,P,V1);
1682   }
1683 }
1684
1685
1686 //=======================================================================
1687 //function : LocalD2
1688 //purpose  :
1689 //=======================================================================
1690
1691 void Law_BSpline::LocalD2
1692 (const Standard_Real    U,
1693  const Standard_Integer FromK1,
1694  const Standard_Integer ToK2,
1695  Standard_Real&    P,
1696  Standard_Real&    V1,
1697  Standard_Real&    V2) const
1698 {
1699   Standard_DomainError_Raise_if (FromK1 == ToK2,
1700                                  "Law_BSpline::LocalD2");
1701   Standard_Real u = U;
1702   Standard_Integer index = 0;
1703   BSplCLib::LocateParameter(deg, FKNOTS, U, periodic, FromK1,ToK2, index, u);
1704   index = BSplCLib::FlatIndex(deg,index,mults->Array1(),periodic);
1705   if (rational) {
1706     BSplCLib::D2(u,index,deg,periodic,POLES, weights->Array1(),FKNOTS,FMULTS,P,V1,V2);
1707   }
1708   else {
1709     BSplCLib::D2(u,index,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,P,V1,V2);
1710   }
1711 }
1712
1713
1714 //=======================================================================
1715 //function : LocalD3
1716 //purpose  :
1717 //=======================================================================
1718
1719 void Law_BSpline::LocalD3
1720 (const Standard_Real    U,
1721  const Standard_Integer FromK1,
1722  const Standard_Integer ToK2,
1723  Standard_Real&    P,
1724  Standard_Real&    V1,
1725  Standard_Real&    V2,
1726  Standard_Real&    V3) const
1727 {
1728   Standard_DomainError_Raise_if (FromK1 == ToK2,
1729                                  "Law_BSpline::LocalD3");
1730   Standard_Real u = U;
1731   Standard_Integer index = 0;
1732   BSplCLib::LocateParameter(deg, FKNOTS, U, periodic, FromK1,ToK2, index, u);
1733   index = BSplCLib::FlatIndex(deg,index,mults->Array1(),periodic);
1734   if (rational) {
1735     BSplCLib::D3(u,index,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,P,V1,V2,V3);
1736   }
1737   else {
1738     BSplCLib::D3(u,index,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,P,V1,V2,V3);
1739   }
1740 }
1741
1742
1743
1744 //=======================================================================
1745 //function : LocalDN
1746 //purpose  :
1747 //=======================================================================
1748
1749 Standard_Real Law_BSpline::LocalDN
1750 (const Standard_Real    U,
1751  const Standard_Integer FromK1,
1752  const Standard_Integer ToK2,
1753  const Standard_Integer N      ) const
1754 {
1755   Standard_DomainError_Raise_if (FromK1 == ToK2,
1756                                  "Law_BSpline::LocalD3");
1757   Standard_Real u = U;
1758   Standard_Integer index = 0;
1759   BSplCLib::LocateParameter(deg, FKNOTS, U, periodic, FromK1,ToK2, index, u);
1760   index = BSplCLib::FlatIndex(deg,index,mults->Array1(),periodic);
1761
1762   Standard_Real V;
1763   if (rational) {
1764     BSplCLib::DN(u,N,index,deg,periodic,POLES,weights->Array1(),FKNOTS,FMULTS,V);
1765   }
1766   else {
1767     BSplCLib::DN(u,N,index,deg,periodic,POLES,BSplCLib::NoWeights(),FKNOTS,FMULTS,V);
1768   }
1769   return V;
1770 }
1771
1772
1773
1774 //=======================================================================
1775 //function : Multiplicity
1776 //purpose  :
1777 //=======================================================================
1778
1779 Standard_Integer Law_BSpline::Multiplicity
1780 (const Standard_Integer Index) const
1781 {
1782   Standard_OutOfRange_Raise_if (Index < 1 || Index > mults->Length(),
1783                                 "Law_BSpline::Multiplicity");
1784   return mults->Value (Index);
1785 }
1786
1787
1788 //=======================================================================
1789 //function : Multiplicities
1790 //purpose  :
1791 //=======================================================================
1792
1793 void Law_BSpline::Multiplicities (TColStd_Array1OfInteger& M) const
1794 {
1795   Standard_DimensionError_Raise_if (M.Length() != mults->Length(),
1796                                     "Law_BSpline::Multiplicities");
1797   M = mults->Array1();
1798 }
1799
1800
1801 //=======================================================================
1802 //function : NbKnots
1803 //purpose  :
1804 //=======================================================================
1805
1806 Standard_Integer Law_BSpline::NbKnots () const
1807 { return knots->Length(); }
1808
1809 //=======================================================================
1810 //function : NbPoles
1811 //purpose  :
1812 //=======================================================================
1813
1814 Standard_Integer Law_BSpline::NbPoles () const
1815 { return poles->Length(); }
1816
1817
1818 //=======================================================================
1819 //function : Pole
1820 //purpose  :
1821 //=======================================================================
1822
1823 Standard_Real Law_BSpline::Pole (const Standard_Integer Index) const
1824 {
1825   Standard_OutOfRange_Raise_if (Index < 1 || Index > poles->Length(),
1826                                 "Law_BSpline::Pole");
1827   return poles->Value (Index);
1828 }
1829
1830
1831 //=======================================================================
1832 //function : Poles
1833 //purpose  :
1834 //=======================================================================
1835
1836 void Law_BSpline::Poles (TColStd_Array1OfReal& P) const
1837 {
1838   Standard_DimensionError_Raise_if (P.Length() != poles->Length(),
1839                                     "Law_BSpline::Poles");
1840   P = poles->Array1();
1841 }
1842
1843
1844 //=======================================================================
1845 //function : StartPoint
1846 //purpose  :
1847 //=======================================================================
1848
1849 Standard_Real Law_BSpline::StartPoint () const
1850 {
1851   if (mults->Value (1) == deg + 1)
1852     return poles->Value (1);
1853   else
1854     return Value(FirstParameter());
1855 }
1856
1857 //=======================================================================
1858 //function : Weight
1859 //purpose  :
1860 //=======================================================================
1861
1862 Standard_Real Law_BSpline::Weight
1863 (const Standard_Integer Index) const
1864 {
1865   Standard_OutOfRange_Raise_if (Index < 1 || Index > poles->Length(),
1866                                 "Law_BSpline::Weight");
1867   if (IsRational())
1868     return weights->Value (Index);
1869   else
1870     return 1.;
1871 }
1872
1873
1874
1875 //=======================================================================
1876 //function : Weights
1877 //purpose  :
1878 //=======================================================================
1879
1880 void Law_BSpline::Weights
1881 (TColStd_Array1OfReal& W) const
1882 {
1883   Standard_DimensionError_Raise_if (W.Length() != poles->Length(),
1884                                     "Law_BSpline::Weights");
1885   if (IsRational())
1886     W = weights->Array1();
1887   else {
1888     Standard_Integer i;
1889     for (i = W.Lower(); i <= W.Upper(); i++)
1890       W(i) = 1.;
1891   }
1892 }
1893
1894
1895
1896 //=======================================================================
1897 //function : IsRational
1898 //purpose  :
1899 //=======================================================================
1900
1901 Standard_Boolean Law_BSpline::IsRational () const
1902 {
1903   return !weights.IsNull();
1904 }
1905
1906
1907 //=======================================================================
1908 //function : LocateU
1909 //purpose  :
1910 //=======================================================================
1911
1912 void Law_BSpline::LocateU
1913 (const Standard_Real    U,
1914  const Standard_Real    ParametricTolerance,
1915  Standard_Integer&      I1,
1916  Standard_Integer&      I2,
1917  const Standard_Boolean WithKnotRepetition) const
1918 {
1919   Standard_Real NewU = U;
1920   Handle(TColStd_HArray1OfReal) TheKnots;
1921   if (WithKnotRepetition)  TheKnots = flatknots;
1922   else                     TheKnots = knots;
1923
1924   PeriodicNormalization(NewU); //Attention a la periode
1925
1926   const TColStd_Array1OfReal & CKnots = TheKnots->Array1();
1927   Standard_Real UFirst = CKnots (1);
1928   Standard_Real ULast  = CKnots (CKnots.Length());
1929   if (Abs (U - UFirst) <= Abs(ParametricTolerance)) { I1 = I2 = 1; }
1930   else if (Abs (U - ULast) <= Abs(ParametricTolerance)) {
1931     I1 = I2 = CKnots.Length();
1932   }
1933   else if (NewU < UFirst - Abs(ParametricTolerance)) {
1934     I2 = 1;
1935     I1 = 0;
1936   }
1937   else if (NewU > ULast + Abs(ParametricTolerance)) {
1938     I1 = CKnots.Length();
1939     I2 = I1 + 1;
1940   }
1941   else {
1942     I1 = 1;
1943     BSplCLib::Hunt (CKnots, NewU, I1);
1944     while ( Abs( CKnots(I1+1) - NewU) <= Abs(ParametricTolerance)) I1++;
1945     if ( Abs( CKnots(I1) - NewU) <= Abs(ParametricTolerance)) {
1946       I2 = I1;
1947     }
1948     else {
1949       I2 = I1 + 1;
1950     }
1951   }
1952 }
1953 //=======================================================================
1954 //function : MovePointAndTangent
1955 //purpose  :
1956 //=======================================================================
1957
1958 void Law_BSpline::
1959   MovePointAndTangent(const Standard_Real    U,
1960                       const Standard_Real    P,
1961                       const Standard_Real    Tangent,
1962                       const Standard_Real    Tolerance,
1963                       const Standard_Integer StartingCondition,
1964                       const Standard_Integer EndingCondition,
1965                       Standard_Integer&      ErrorStatus)
1966 {
1967   TColStd_Array1OfReal new_poles(1, poles->Length());
1968   Standard_Real delta,
1969   *poles_array,
1970   *new_poles_array,
1971   delta_derivative;
1972   const Standard_Integer dimension = 1 ;
1973   D1(U,
1974      delta,
1975      delta_derivative) ;
1976   delta = P - delta  ;
1977
1978   delta_derivative = Tangent - delta_derivative ;
1979   poles_array = (Standard_Real *)
1980     &poles->Array1()(1) ;
1981   new_poles_array = (Standard_Real *)
1982     &new_poles(1) ;
1983   BSplCLib::MovePointAndTangent(U,
1984                                 dimension,
1985                                 delta,
1986                                 delta_derivative,
1987                                 Tolerance,
1988                                 deg,
1989                                 rational,
1990                                 StartingCondition,
1991                                 EndingCondition,
1992                                 poles_array[0],
1993                                 weights->Array1(),
1994                                 flatknots->Array1(),
1995                                 new_poles_array[0],
1996                                 ErrorStatus) ;
1997   if (!ErrorStatus) {
1998     poles->ChangeArray1() = new_poles;
1999   }
2000 }
2001
2002
2003 //=======================================================================
2004 //function : Resolution
2005 //purpose  :
2006 //=======================================================================
2007
2008 void Law_BSpline::Resolution(const Standard_Real Tolerance3D,
2009                              Standard_Real &     UTolerance) const
2010 {
2011   void* bid = (void*)(&(poles->Value(1)));
2012   Standard_Real* bidr = (Standard_Real*)bid;
2013   if (rational) {
2014     BSplCLib::Resolution(*bidr,1,poles->Length(),
2015                          weights->Array1(),FKNOTS,deg,
2016                          Tolerance3D,
2017                          UTolerance) ;
2018   }
2019   else {
2020
2021     BSplCLib::Resolution(*bidr,1,poles->Length(),
2022                          BSplCLib::NoWeights(),FKNOTS,deg,
2023                          Tolerance3D,
2024                          UTolerance) ;
2025   }
2026 }