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