src/BSplCLib/BSplCLib_CurveComputation.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 // xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
  16 //                          Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
  17 //                          EvalPolynomial, RationalDerivatives.
  18 // xab : 22-Mar-94 : fixed rational problem in BuildCache
  19 // xab : 12-Mar-96 : added MovePointAndTangent
  20 #include <TColStd_Array1OfInteger.hxx>
  21 #include <TColStd_Array1OfReal.hxx>
  22 #include <TColStd_Array2OfReal.hxx>
  23 #include <gp_Vec2d.hxx>
  24 #include <Standard_ConstructionError.hxx>
  25 #include <PLib.hxx>
  26 #include <math_Matrix.hxx>
  27
  28 //=======================================================================
  29 //struct : BSplCLib_DataContainer
  30 //purpose: Auxiliary structure providing buffers for poles and knots used in
  31 //         evaluation of bspline (allocated in the stack)
  32 //=======================================================================
  33
  34 struct BSplCLib_DataContainer
  35 {
  36   BSplCLib_DataContainer(Standard_Integer Degree)
  37   {
  38     (void)Degree; // avoid compiler warning
  39     Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
  40         BSplCLib::MaxDegree() > 25,
  41         "BSplCLib: bspline degree is greater than maximum supported");
  42   }
  43
  44   Standard_Real poles[(25+1)*(Dimension_gen+1)];
  45   Standard_Real knots[2*25];
  46   Standard_Real ders[Dimension_gen*4];
  47 };
  48
  49 //=======================================================================
  50 //function : Reverse
  51 //purpose  :
  52 //=======================================================================
  53
  54 void  BSplCLib::Reverse(Array1OfPoints& Poles,
  55                         const Standard_Integer L)
  56 {
  57   Standard_Integer i, l = L;
  58   l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
  59
  60   Array1OfPoints temp(0,Poles.Length()-1);
  61
  62   for (i = Poles.Lower(); i <= l; i++)
  63     temp(l-i) = Poles(i);
  64
  65   for (i = l+1; i <= Poles.Upper(); i++)
  66     temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
  67
  68   for (i = Poles.Lower(); i <= Poles.Upper(); i++)
  69     Poles(i) = temp(i-Poles.Lower());
  70 }
  71
  72 //
  73 // CURVES MODIFICATIONS
  74 //
  75
  76 //=======================================================================
  77 //function : RemoveKnot
  78 //purpose  :
  79 //=======================================================================
  80
  81 Standard_Boolean BSplCLib::RemoveKnot
  82 (const Standard_Integer         Index,
  83  const Standard_Integer         Mult,
  84  const Standard_Integer         Degree,
  85  const Standard_Boolean         Periodic,
  86  const Array1OfPoints&          Poles,
  87  const TColStd_Array1OfReal*    Weights,
  88  const TColStd_Array1OfReal&    Knots,
  89  const TColStd_Array1OfInteger& Mults,
  90  Array1OfPoints&                NewPoles,
  91  TColStd_Array1OfReal*          NewWeights,
  92  TColStd_Array1OfReal&          NewKnots,
  93  TColStd_Array1OfInteger&       NewMults,
  94  const Standard_Real            Tolerance)
  95 {
  96   Standard_Boolean rational = Weights  != NULL;
  97   Standard_Integer dim;
  98   dim = Dimension_gen;
  99   if (rational) dim++;
 100
 101   TColStd_Array1OfReal poles(1,dim*Poles.Length());
 102   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
 103
 104   if (rational) PLib::SetPoles(Poles,*Weights,poles);
 105   else          PLib::SetPoles(Poles,poles);
 106
 107   if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
 108                   poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
 109     return Standard_False;
 110
 111   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
 112   else          PLib::GetPoles(newpoles,NewPoles);
 113   return Standard_True;
 114 }
 115
 116 //=======================================================================
 117 //function : InsertKnots
 118 //purpose  : insert an array of knots and multiplicities
 119 //=======================================================================
 120
 121 void BSplCLib::InsertKnots
 122 (const Standard_Integer         Degree,
 123  const Standard_Boolean         Periodic,
 124  const Array1OfPoints&          Poles,
 125  const TColStd_Array1OfReal*    Weights,
 126  const TColStd_Array1OfReal&    Knots,
 127  const TColStd_Array1OfInteger& Mults,
 128  const TColStd_Array1OfReal&    AddKnots,
 129  const TColStd_Array1OfInteger* AddMults,
 130  Array1OfPoints&                NewPoles,
 131  TColStd_Array1OfReal*          NewWeights,
 132  TColStd_Array1OfReal&          NewKnots,
 133  TColStd_Array1OfInteger&       NewMults,
 134  const Standard_Real            Epsilon,
 135  const Standard_Boolean         Add)
 136 {
 137   Standard_Boolean rational = Weights  != NULL;
 138   Standard_Integer dim;
 139   dim = Dimension_gen;
 140   if (rational) dim++;
 141
 142   TColStd_Array1OfReal poles(1,dim*Poles.Length());
 143   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
 144
 145   if (rational) PLib::SetPoles(Poles,*Weights,poles);
 146   else          PLib::SetPoles(Poles,poles);
 147
 148   InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
 149               AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
 150
 151   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
 152   else          PLib::GetPoles(newpoles,NewPoles);
 153 }
 154
 155 //=======================================================================
 156 //function : InsertKnot
 157 //purpose  :
 158 //=======================================================================
 159
 160 void  BSplCLib::InsertKnot(const Standard_Integer ,
 161                            const Standard_Real U,
 162                            const Standard_Integer UMult,
 163                            const Standard_Integer Degree,
 164                            const Standard_Boolean Periodic,
 165                            const Array1OfPoints& Poles,
 166                            const TColStd_Array1OfReal* Weights,
 167                            const TColStd_Array1OfReal& Knots,
 168                            const TColStd_Array1OfInteger& Mults,
 169                            Array1OfPoints& NewPoles,
 170                            TColStd_Array1OfReal* NewWeights)
 171 {
 172   TColStd_Array1OfReal k(1,1);
 173   k(1) = U;
 174   TColStd_Array1OfInteger m(1,1);
 175   m(1) = UMult;
 176   TColStd_Array1OfReal    nk(1,Knots.Length()+1);
 177   TColStd_Array1OfInteger nm(1,Knots.Length()+1);
 178   InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
 179               k,&m,NewPoles,NewWeights,nk,nm,Epsilon(U));
 180 }
 181
 182 //=======================================================================
 183 //function : RaiseMultiplicity
 184 //purpose  :
 185 //=======================================================================
 186
 187 void  BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
 188                                   const Standard_Integer Mult,
 189                                   const Standard_Integer Degree,
 190                                   const Standard_Boolean Periodic,
 191                                   const Array1OfPoints& Poles,
 192                                   const TColStd_Array1OfReal* Weights,
 193                                   const TColStd_Array1OfReal& Knots,
 194                                   const TColStd_Array1OfInteger& Mults,
 195                                   Array1OfPoints& NewPoles,
 196                                   TColStd_Array1OfReal* NewWeights)
 197 {
 198   TColStd_Array1OfReal k(1,1);
 199   k(1) = Knots(KnotIndex);
 200   TColStd_Array1OfInteger m(1,1);
 201   m(1) = Mult - Mults(KnotIndex);
 202   TColStd_Array1OfReal    nk(1,Knots.Length());
 203   TColStd_Array1OfInteger nm(1,Knots.Length());
 204   InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
 205               k,&m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
 206 }
 207
 208 //=======================================================================
 209 //function : IncreaseDegree
 210 //purpose  :
 211 //=======================================================================
 212
 213 void BSplCLib::IncreaseDegree
 214 (const Standard_Integer          Degree,
 215  const Standard_Integer          NewDegree,
 216  const Standard_Boolean          Periodic,
 217  const Array1OfPoints&           Poles,
 218  const TColStd_Array1OfReal*     Weights,
 219  const TColStd_Array1OfReal&     Knots,
 220  const TColStd_Array1OfInteger&  Mults,
 221  Array1OfPoints&                 NewPoles,
 222  TColStd_Array1OfReal*           NewWeights,
 223  TColStd_Array1OfReal&           NewKnots,
 224  TColStd_Array1OfInteger&        NewMults)
 225 {
 226   Standard_Boolean rational = Weights  != NULL;
 227   Standard_Integer dim;
 228   dim = Dimension_gen;
 229   if (rational) dim++;
 230
 231   TColStd_Array1OfReal poles(1,dim*Poles.Length());
 232   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
 233
 234   if (rational) PLib::SetPoles(Poles,*Weights,poles);
 235   else          PLib::SetPoles(Poles,poles);
 236
 237   IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
 238                  newpoles,NewKnots,NewMults);
 239
 240   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
 241   else          PLib::GetPoles(newpoles,NewPoles);
 242 }
 243
 244 //=======================================================================
 245 //function : Unperiodize
 246 //purpose  :
 247 //=======================================================================
 248
 249 void  BSplCLib::Unperiodize
 250 (const Standard_Integer         Degree,
 251  const TColStd_Array1OfInteger& Mults,
 252  const TColStd_Array1OfReal&    Knots,
 253  const Array1OfPoints&          Poles,
 254  const TColStd_Array1OfReal*    Weights,
 255  TColStd_Array1OfInteger& NewMults,
 256  TColStd_Array1OfReal&    NewKnots,
 257  Array1OfPoints&          NewPoles,
 258  TColStd_Array1OfReal*    NewWeights)
 259 {
 260   Standard_Boolean rational = Weights  != NULL;
 261   Standard_Integer dim;
 262   dim = Dimension_gen;
 263   if (rational) dim++;
 264
 265   TColStd_Array1OfReal poles(1,dim*Poles.Length());
 266   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
 267
 268   if (rational) PLib::SetPoles(Poles,*Weights,poles);
 269   else          PLib::SetPoles(Poles,poles);
 270
 271   Unperiodize(Degree, dim, Mults, Knots, poles,
 272               NewMults,NewKnots, newpoles);
 273
 274   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
 275   else          PLib::GetPoles(newpoles,NewPoles);
 276 }
 277
 278 //=======================================================================
 279 //function : Trimming
 280 //purpose  :
 281 //=======================================================================
 282
 283 void BSplCLib::Trimming(const Standard_Integer         Degree,
 284                         const Standard_Boolean         Periodic,
 285                         const TColStd_Array1OfReal&    Knots,
 286                         const TColStd_Array1OfInteger& Mults,
 287                         const Array1OfPoints&          Poles,
 288                         const TColStd_Array1OfReal*    Weights,
 289                         const Standard_Real            U1,
 290                         const Standard_Real            U2,
 291                         TColStd_Array1OfReal&    NewKnots,
 292                         TColStd_Array1OfInteger& NewMults,
 293                         Array1OfPoints&          NewPoles,
 294                         TColStd_Array1OfReal*    NewWeights)
 295 {
 296   Standard_Boolean rational = Weights  != NULL;
 297   Standard_Integer dim;
 298   dim = Dimension_gen;
 299   if (rational) dim++;
 300
 301   TColStd_Array1OfReal poles(1,dim*Poles.Length());
 302   TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
 303
 304   if (rational) PLib::SetPoles(Poles,*Weights,poles);
 305   else          PLib::SetPoles(Poles,poles);
 306
 307   Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
 308            NewKnots, NewMults, newpoles);
 309
 310   if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
 311   else          PLib::GetPoles(newpoles,NewPoles);
 312 }
 313
 314 //--------------------------------------------------------------------------
 315 //ELEMENTARY COMPUTATIONS
 316 //--------------------------------------------------------------------------
 317 //
 318 // All the following methods are methods of low level used in BSplCLib
 319 // but also in BSplSLib (but not in Geom)
 320 // At the creation time of this package there is no possibility
 321 // to declare private methods of package and to declare friend
 322 // methods of package.  It could interesting to declare that BSplSLib
 323 // is a package friend to BSplCLib because it uses all the basis methods
 324 // of BSplCLib.
 325 //--------------------------------------------------------------------------
 326
 327 //=======================================================================
 328 //function : BuildEval
 329 //purpose  : builds the local array for evaluation
 330 //=======================================================================
 331
 332 void  BSplCLib::BuildEval(const Standard_Integer         Degree,
 333                           const Standard_Integer         Index,
 334                           const Array1OfPoints&          Poles,
 335                           const TColStd_Array1OfReal*    Weights,
 336                           Standard_Real&                 LP)
 337 {
 338   Standard_Real w, *pole = &LP;
 339   Standard_Integer PLower = Poles.Lower();
 340   Standard_Integer PUpper = Poles.Upper();
 341   Standard_Integer i;
 342   Standard_Integer ip = PLower + Index - 1;
 343   if (Weights == NULL) {
 344     for (i = 0; i <= Degree; i++) {
 345       ip++;
 346       if (ip > PUpper) ip = PLower;
 347       const Point& P = Poles(ip);
 348       PointToCoords (pole,P,+0);
 349       pole += Dimension_gen;
 350     }
 351   }
 352   else {
 353     for (i = 0; i <= Degree; i++) {
 354       ip++;
 355       if (ip > PUpper) ip = PLower;
 356       const Point& P = Poles(ip);
 357       pole[Dimension_gen] = w = (*Weights)(ip);
 358       PointToCoords (pole, P, * w);
 359       pole += Dimension_gen + 1;
 360     }
 361   }
 362 }
 363
 364 //=======================================================================
 365 //function : PrepareEval
 366 //purpose  : stores data for Eval in the local arrays
 367 //           dc.poles and dc.knots
 368 //=======================================================================
 369
 370 static void PrepareEval
 371 (Standard_Real&                 u,
 372  Standard_Integer&              index,
 373  Standard_Integer&              dim,
 374  Standard_Boolean&              rational,
 375  const Standard_Integer         Degree,
 376  const Standard_Boolean         Periodic,
 377  const Array1OfPoints&          Poles,
 378  const TColStd_Array1OfReal*    Weights,
 379  const TColStd_Array1OfReal&    Knots,
 380  const TColStd_Array1OfInteger* Mults,
 381  BSplCLib_DataContainer& dc)
 382 {
 383   // Set the Index
 384   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
 385
 386   // make the knots
 387   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
 388   if (Mults == NULL)
 389     index -= Knots.Lower() + Degree;
 390   else
 391     index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults);
 392
 393   // check truly rational
 394   rational = (Weights != NULL);
 395   if (rational) {
 396     Standard_Integer WLower = Weights->Lower() + index;
 397     rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree);
 398   }
 399
 400   // make the poles
 401   if (rational) {
 402     dim = Dimension_gen + 1;
 403     BSplCLib::BuildEval(Degree, index, Poles, Weights              , *dc.poles);
 404   }
 405   else {
 406     dim = Dimension_gen;
 407     BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
 408   }
 409 }
 410
 411 //=======================================================================
 412 //function : D0
 413 //purpose  :
 414 //=======================================================================
 415
 416 void BSplCLib::D0
 417 (const Standard_Real            U,
 418  const Standard_Integer         Index,
 419  const Standard_Integer         Degree,
 420  const Standard_Boolean         Periodic,
 421  const Array1OfPoints&          Poles,
 422  const TColStd_Array1OfReal*    Weights,
 423  const TColStd_Array1OfReal&    Knots,
 424  const TColStd_Array1OfInteger* Mults,
 425  Point&                         P)
 426 {
 427 //  Standard_Integer k,dim,index = Index;
 428   Standard_Integer dim,index = Index;
 429   Standard_Real    u = U;
 430   Standard_Boolean rational;
 431   BSplCLib_DataContainer dc(Degree);
 432   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 433   BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
 434
 435   if (rational) {
 436     Standard_Real w = dc.poles[Dimension_gen];
 437     CoordsToPoint (P, dc.poles, / w);
 438   }
 439   else
 440     CoordsToPoint (P, dc.poles, );
 441 }
 442
 443 //=======================================================================
 444 //function : D1
 445 //purpose  :
 446 //=======================================================================
 447
 448 void BSplCLib::D1
 449 (const Standard_Real            U,
 450  const Standard_Integer         Index,
 451  const Standard_Integer         Degree,
 452  const Standard_Boolean         Periodic,
 453  const Array1OfPoints&          Poles,
 454  const TColStd_Array1OfReal*    Weights,
 455  const TColStd_Array1OfReal&    Knots,
 456  const TColStd_Array1OfInteger* Mults,
 457  Point&                         P,
 458  Vector&                        V)
 459 {
 460   Standard_Integer dim,index = Index;
 461   Standard_Real    u = U;
 462   Standard_Boolean rational;
 463   BSplCLib_DataContainer dc(Degree);
 464   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 465   BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
 466   Standard_Real *result = dc.poles;
 467   if (rational) {
 468     PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
 469     result = dc.ders;
 470   }
 471
 472   CoordsToPoint (P, result, );
 473   CoordsToPoint (V, result + Dimension_gen, );
 474 }
 475
 476 //=======================================================================
 477 //function : D2
 478 //purpose  :
 479 //=======================================================================
 480
 481 void BSplCLib::D2
 482 (const Standard_Real            U,
 483  const Standard_Integer         Index,
 484  const Standard_Integer         Degree,
 485  const Standard_Boolean         Periodic,
 486  const Array1OfPoints&          Poles,
 487  const TColStd_Array1OfReal*    Weights,
 488  const TColStd_Array1OfReal&    Knots,
 489  const TColStd_Array1OfInteger* Mults,
 490  Point&                         P,
 491  Vector&                        V1,
 492  Vector&                        V2)
 493 {
 494   Standard_Integer dim,index = Index;
 495   Standard_Real    u = U;
 496   Standard_Boolean rational;
 497   BSplCLib_DataContainer dc(Degree);
 498   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 499   BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
 500   Standard_Real *result = dc.poles;
 501   if (rational) {
 502     PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
 503     result = dc.ders;
 504   }
 505
 506   CoordsToPoint (P,  result, );
 507   CoordsToPoint (V1, result + Dimension_gen, );
 508   if (!rational && (Degree < 2))
 509     NullifyPoint (V2);
 510   else
 511     CoordsToPoint (V2, result + 2 * Dimension_gen, );
 512 }
 513
 514 //=======================================================================
 515 //function : D3
 516 //purpose  :
 517 //=======================================================================
 518
 519 void BSplCLib::D3
 520 (const Standard_Real            U,
 521  const Standard_Integer         Index,
 522  const Standard_Integer         Degree,
 523  const Standard_Boolean         Periodic,
 524  const Array1OfPoints&          Poles,
 525  const TColStd_Array1OfReal*    Weights,
 526  const TColStd_Array1OfReal&    Knots,
 527  const TColStd_Array1OfInteger* Mults,
 528  Point&                         P,
 529  Vector&                        V1,
 530  Vector&                        V2,
 531  Vector&                        V3)
 532 {
 533   Standard_Integer dim,index = Index;
 534   Standard_Real    u = U;
 535   Standard_Boolean rational;
 536   BSplCLib_DataContainer dc(Degree);
 537   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 538   BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
 539   Standard_Real *result = dc.poles;
 540   if (rational) {
 541     PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
 542     result = dc.ders;
 543   }
 544
 545   CoordsToPoint (P,  result, );
 546   CoordsToPoint (V1, result + Dimension_gen, );
 547   if (!rational && (Degree < 2))
 548     NullifyPoint (V2);
 549   else
 550     CoordsToPoint (V2, result + 2 * Dimension_gen, );
 551   if (!rational && (Degree < 3))
 552     NullifyPoint (V3);
 553   else
 554     CoordsToPoint (V3, result + 3 * Dimension_gen, );
 555 }
 556
 557 //=======================================================================
 558 //function : DN
 559 //purpose  :
 560 //=======================================================================
 561
 562 void BSplCLib::DN
 563 (const Standard_Real            U,
 564  const Standard_Integer         N,
 565  const Standard_Integer         Index,
 566  const Standard_Integer         Degree,
 567  const Standard_Boolean         Periodic,
 568  const Array1OfPoints&          Poles,
 569  const TColStd_Array1OfReal*    Weights,
 570  const TColStd_Array1OfReal&    Knots,
 571  const TColStd_Array1OfInteger* Mults,
 572  Vector&                        VN)
 573 {
 574   Standard_Integer dim,index = Index;
 575   Standard_Real    u = U;
 576   Standard_Boolean rational;
 577   BSplCLib_DataContainer dc(Degree);
 578   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 579   BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
 580
 581   if (rational) {
 582     Standard_Real v[Dimension_gen];
 583     PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
 584     CoordsToPoint (VN, v, );
 585   }
 586   else {
 587     if (N > Degree)
 588       NullifyPoint (VN);
 589     else {
 590       Standard_Real *DN = dc.poles + N * Dimension_gen;
 591       CoordsToPoint (VN, DN, );
 592     }
 593   }
 594 }
 595
 596 //=======================================================================
 597 //function : Solves a LU factored Matrix
 598 //purpose  :
 599 //=======================================================================
 600
 601 Standard_Integer
 602 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
 603                             const Standard_Integer UpperBandWidth,
 604                             const Standard_Integer LowerBandWidth,
 605                             Array1OfPoints&   PolesArray)
 606 {
 607   Standard_Real *PArray ;
 608   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
 609
 610   return BSplCLib::SolveBandedSystem(Matrix,
 611                                      UpperBandWidth,
 612                                      LowerBandWidth,
 613                                      Dimension_gen,
 614                                      PArray[0]) ;
 615 }
 616
 617 //=======================================================================
 618 //function : Solves a LU factored Matrix
 619 //purpose  : if HomogeneousFlag is 1 then the input and the output
 620 //           will be homogeneous that is no division or multiplication
 621 //           by weigths will happen. On the contrary if HomogeneousFlag
 622 //           is 0 then the poles will be multiplied first by the weights
 623 //           and after interpolation they will be devided by the weights
 624 //=======================================================================
 625
 626 Standard_Integer
 627 BSplCLib::SolveBandedSystem(const math_Matrix&  Matrix,
 628                             const Standard_Integer UpperBandWidth,
 629                             const Standard_Integer LowerBandWidth,
 630                             const Standard_Boolean HomogeneousFlag,
 631                             Array1OfPoints&        PolesArray,
 632                             TColStd_Array1OfReal&  WeightsArray)
 633 {
 634   Standard_Real *PArray,
 635   * WArray ;
 636   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
 637   WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
 638   return BSplCLib::SolveBandedSystem(Matrix,
 639                                      UpperBandWidth,
 640                                      LowerBandWidth,
 641                                      HomogeneousFlag,
 642                                      Dimension_gen,
 643                                      PArray[0],
 644                                      WArray[0]) ;
 645 }
 646
 647 //=======================================================================
 648 //function : Evaluates a Bspline function : uses the ExtrapMode
 649 //purpose  : the function is extrapolated using the Taylor expansion
 650 //           of degree ExtrapMode[0] to the left and the Taylor
 651 //           expansion of degree ExtrapMode[1] to the right
 652 //   if the HomogeneousFlag == True than the Poles are supposed
 653 //   to be stored homogeneously and the result will also be homogeneous
 654 //   Valid only if Weights
 655 //=======================================================================
 656 void  BSplCLib::Eval(const Standard_Real                   Parameter,
 657                      const Standard_Boolean                PeriodicFlag,
 658                      const Standard_Boolean                HomogeneousFlag,
 659                      Standard_Integer&                     ExtrapMode,
 660                      const Standard_Integer                Degree,
 661                      const  TColStd_Array1OfReal&          FlatKnots,
 662                      const  Array1OfPoints&                PolesArray,
 663                      const  TColStd_Array1OfReal&          WeightsArray,
 664                      Point&                                aPoint,
 665                      Standard_Real&                        aWeight)
 666 {
 667   Standard_Real  Inverse,
 668   P[Dimension_gen],
 669   *PArray,
 670   *WArray ;
 671   Standard_Integer kk ;
 672   PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
 673   WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
 674   if (HomogeneousFlag) {
 675     BSplCLib::Eval(Parameter,
 676                    PeriodicFlag,
 677                    0,
 678                    ExtrapMode,
 679                    Degree,
 680                    FlatKnots,
 681                    Dimension_gen,
 682                    PArray[0],
 683                    P[0]) ;
 684     BSplCLib::Eval(Parameter,
 685                    PeriodicFlag,
 686                    0,
 687                    ExtrapMode,
 688                    Degree,
 689                    FlatKnots,
 690                    1,
 691                    WArray[0],
 692                    aWeight) ;
 693   }
 694   else {
 695     BSplCLib::Eval(Parameter,
 696                    PeriodicFlag,
 697                    0,
 698                    ExtrapMode,
 699                    Degree,
 700                    FlatKnots,
 701                    Dimension_gen,
 702                    PArray[0],
 703                    WArray[0],
 704                    P[0],
 705                    aWeight) ;
 706     Inverse = 1.0e0 / aWeight ;
 707
 708     for (kk = 0 ; kk < Dimension_gen ; kk++) {
 709       P[kk] *= Inverse ;
 710     }
 711   }
 712
 713   for (kk = 0 ; kk < Dimension_gen ; kk++)
 714     aPoint.SetCoord(kk + 1, P[kk]) ;
 715 }
 716
 717 //=======================================================================
 718 //function : CacheD0
 719 //purpose  : Evaluates the polynomial cache of the Bspline Curve
 720 //
 721 //=======================================================================
 722 void  BSplCLib::CacheD0(const Standard_Real                  Parameter,
 723                         const  Standard_Integer               Degree,
 724                         const  Standard_Real                  CacheParameter,
 725                         const  Standard_Real                  SpanLenght,
 726                         const  Array1OfPoints&                PolesArray,
 727                         const  TColStd_Array1OfReal*          WeightsArray,
 728                         Point&                                aPoint)
 729 {
 730   //
 731   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
 732   // form
 733   // the SpanLenght     is the normalizing factor so that the polynomial is between
 734   // 0 and 1
 735   Standard_Real NewParameter, Inverse;
 736   Standard_Real * PArray  = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
 737   Standard_Real * myPoint = (Standard_Real *) &aPoint  ;
 738   NewParameter = (Parameter - CacheParameter) / SpanLenght ;
 739   PLib::NoDerivativeEvalPolynomial(NewParameter,
 740                        Degree,
 741                        Dimension_gen,
 742                        Degree * Dimension_gen,
 743                        PArray[0],
 744                        myPoint[0]) ;
 745   if (WeightsArray != NULL) {
 746     const TColStd_Array1OfReal& refWeights = *WeightsArray;
 747     Standard_Real *
 748       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
 749     PLib::NoDerivativeEvalPolynomial(NewParameter,
 750                          Degree,
 751                          1,
 752                          Degree,
 753                          WArray[0],
 754                          Inverse) ;
 755
 756     Inverse = 1.0e0 / Inverse;
 757     ModifyCoords (myPoint, *= Inverse);
 758   }
 759 }
 760
 761 //=======================================================================
 762 //function : CacheD1
 763 //purpose  : Evaluates the polynomial cache of the Bspline Curve
 764 //
 765 //=======================================================================
 766 void  BSplCLib::CacheD1(const Standard_Real                  Parameter,
 767                         const Standard_Integer                Degree,
 768                         const  Standard_Real                  CacheParameter,
 769                         const  Standard_Real                  SpanLenght,
 770                         const  Array1OfPoints&                PolesArray,
 771                         const  TColStd_Array1OfReal*          WeightsArray,
 772                         Point&                                aPoint,
 773                         Vector&                               aVector)
 774 {
 775   //
 776   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
 777   // form
 778   // the SpanLenght     is the normalizing factor so that the polynomial is between
 779   // 0 and 1
 780   Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
 781 //  Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
 782   Standard_Real LocalWDerivatives[2], NewParameter ;
 783
 784   Standard_Real *
 785     PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
 786   Standard_Real *
 787     myPoint = (Standard_Real *) &aPoint  ;
 788   Standard_Real *
 789     myVector = (Standard_Real *) &aVector ;
 790   NewParameter = (Parameter - CacheParameter) / SpanLenght ;
 791   PLib::EvalPolynomial(NewParameter,
 792                        1,
 793                        Degree,
 794                        Dimension_gen,
 795                        PArray[0],
 796                        LocalPDerivatives[0]) ;
 797   //
 798   // unormalize derivatives since those are computed normalized
 799   //
 800
 801   ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
 802
 803   if (WeightsArray != NULL) {
 804     const TColStd_Array1OfReal& refWeights = *WeightsArray;
 805     Standard_Real *
 806       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
 807     PLib::EvalPolynomial(NewParameter,
 808                          1,
 809                          Degree,
 810                          1,
 811                          WArray[0],
 812                          LocalWDerivatives[0]) ;
 813     //
 814     // unormalize the result since the polynomial stored in the cache
 815     // is normalized between 0 and 1
 816     //
 817     LocalWDerivatives[1] /= SpanLenght ;
 818
 819     PLib::RationalDerivatives(1,
 820                               Dimension_gen,
 821                               LocalPDerivatives[0],
 822                               LocalWDerivatives[0],
 823                               LocalPDerivatives[0]) ;
 824   }
 825
 826   CopyCoords (myPoint,  LocalPDerivatives);
 827   CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
 828 }
 829
 830 //=======================================================================
 831 //function : CacheD2
 832 //purpose  : Evaluates the polynomial cache of the Bspline Curve
 833 //
 834 //=======================================================================
 835 void  BSplCLib::CacheD2(const Standard_Real                  Parameter,
 836                         const Standard_Integer                Degree,
 837                         const  Standard_Real                  CacheParameter,
 838                         const  Standard_Real                  SpanLenght,
 839                         const  Array1OfPoints&                PolesArray,
 840                         const  TColStd_Array1OfReal*          WeightsArray,
 841                         Point&                                aPoint,
 842                         Vector&                               aVector1,
 843                         Vector&                               aVector2)
 844 {
 845   //
 846   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
 847   // form
 848   // the SpanLenght     is the normalizing factor so that the polynomial is between
 849   // 0 and 1
 850   Standard_Integer ii,Index,EndIndex;
 851   Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
 852 //  Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
 853   Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
 854   Standard_Real * PArray    = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
 855   Standard_Real * myPoint   = (Standard_Real *) &aPoint  ;
 856   Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
 857   Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
 858   NewParameter = (Parameter - CacheParameter) / SpanLenght ;
 859   PLib::EvalPolynomial(NewParameter,
 860                        2,
 861                        Degree,
 862                        Dimension_gen,
 863                        PArray[0],
 864                        LocalPDerivatives[0]) ;
 865   //
 866   // unormalize derivatives since those are computed normalized
 867   //
 868   Factor = 1.0e0/ SpanLenght ;
 869   Index = Dimension_gen ;
 870   EndIndex = Min (2, Degree) ;
 871
 872   for (ii = 1 ; ii <= EndIndex ; ii++) {
 873     ModifyCoords (LocalPDerivatives + Index, *= Factor);
 874     Factor /= SpanLenght ;
 875     Index  += Dimension_gen;
 876   }
 877
 878   Index = (Degree + 1) * Dimension_gen;
 879   for (ii = Degree ; ii < 2 ; ii++) {
 880     NullifyCoords (LocalPDerivatives + Index);
 881     Index += Dimension_gen;
 882   }
 883
 884   if (WeightsArray != NULL) {
 885     const TColStd_Array1OfReal& refWeights = *WeightsArray;
 886     Standard_Real *
 887       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
 888
 889     PLib::EvalPolynomial(NewParameter,
 890                          2,
 891                          Degree,
 892                          1,
 893                          WArray[0],
 894                          LocalWDerivatives[0]) ;
 895
 896     for (ii = Degree + 1  ; ii <= 2 ; ii++) {
 897       LocalWDerivatives[ii] = 0.0e0 ;
 898     }
 899     //
 900     // unormalize the result since the polynomial stored in the cache
 901     // is normalized between 0 and 1
 902     //
 903     Factor = 1.0e0 / SpanLenght ;
 904
 905     for (ii = 1 ; ii <= EndIndex ; ii++) {
 906       LocalWDerivatives[ii] *= Factor ;
 907       Factor /= SpanLenght ;
 908     }
 909     PLib::RationalDerivatives(2,
 910                               Dimension_gen,
 911                               LocalPDerivatives[0],
 912                               LocalWDerivatives[0],
 913                               LocalPDerivatives[0]) ;
 914   }
 915
 916   CopyCoords (myPoint,   LocalPDerivatives);
 917   CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
 918   CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
 919 }
 920
 921 //=======================================================================
 922 //function : CacheD3
 923 //purpose  : Evaluates the polynomial cache of the Bspline Curve
 924 //
 925 //=======================================================================
 926 void  BSplCLib::CacheD3(const Standard_Real                  Parameter,
 927                         const Standard_Integer                Degree,
 928                         const  Standard_Real                  CacheParameter,
 929                         const  Standard_Real                  SpanLenght,
 930                         const  Array1OfPoints&                PolesArray,
 931                         const  TColStd_Array1OfReal*          WeightsArray,
 932                         Point&                                aPoint,
 933                         Vector&                               aVector1,
 934                         Vector&                               aVector2,
 935                         Vector&                               aVector3)
 936 {
 937   //
 938   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
 939   // form
 940   // the SpanLenght     is the normalizing factor so that the polynomial is between
 941   // 0 and 1
 942   Standard_Integer ii, Index, EndIndex;
 943   Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
 944 //  Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
 945   Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
 946   Standard_Real * PArray    = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
 947   Standard_Real * myPoint   = (Standard_Real *) &aPoint  ;
 948   Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
 949   Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
 950   Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
 951   NewParameter = (Parameter - CacheParameter) / SpanLenght ;
 952   PLib::EvalPolynomial(NewParameter,
 953                        3,
 954                        Degree,
 955                        Dimension_gen,
 956                        PArray[0],
 957                        LocalPDerivatives[0]) ;
 958
 959   Index = (Degree + 1) * Dimension_gen;
 960   for (ii = Degree ; ii < 3 ; ii++) {
 961     NullifyCoords (LocalPDerivatives + Index);
 962     Index += Dimension_gen;
 963   }
 964
 965   //
 966   // unormalize derivatives since those are computed normalized
 967   //
 968   Factor = 1.0e0 / SpanLenght ;
 969   Index = Dimension_gen ;
 970   EndIndex = Min(3,Degree) ;
 971
 972   for (ii = 1 ; ii <= EndIndex ; ii++) {
 973     ModifyCoords (LocalPDerivatives + Index, *= Factor);
 974     Factor /= SpanLenght;
 975     Index  += Dimension_gen;
 976   }
 977
 978   if (WeightsArray != NULL) {
 979     const TColStd_Array1OfReal& refWeights = *WeightsArray;
 980     Standard_Real *
 981       WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
 982
 983     PLib::EvalPolynomial(NewParameter,
 984                          3,
 985                          Degree,
 986                          1,
 987                          WArray[0],
 988                          LocalWDerivatives[0]) ;
 989     //
 990     // unormalize the result since the polynomial stored in the cache
 991     // is normalized between 0 and 1
 992     //
 993     Factor = 1.0e0 / SpanLenght ;
 994
 995     for (ii = 1 ; ii <= EndIndex ; ii++) {
 996       LocalWDerivatives[ii] *= Factor ;
 997       Factor /= SpanLenght ;
 998     }
 999
1000     for (ii = (Degree + 1)  ; ii <= 3 ; ii++) {
1001       LocalWDerivatives[ii] = 0.0e0 ;
1002     }
1003     PLib::RationalDerivatives(3,
1004                               Dimension_gen,
1005                               LocalPDerivatives[0],
1006                               LocalWDerivatives[0],
1007                               LocalPDerivatives[0]) ;
1008   }
1009
1010   CopyCoords (myPoint,   LocalPDerivatives);
1011   CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1012   CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1013   CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1014 }
1015
1016 //=======================================================================
1017 //function : BuildCache
1018 //purpose  : Stores theTaylor expansion normalized between 0,1 in the
1019 //           Cache : in case of  a rational function the Poles are
1020 //           stored in homogeneous form
1021 //=======================================================================
1022
1023 void BSplCLib::BuildCache
1024 (const Standard_Real            U,
1025  const Standard_Real            SpanDomain,
1026  const Standard_Boolean         Periodic,
1027  const Standard_Integer         Degree,
1028  const TColStd_Array1OfReal&    FlatKnots,
1029  const Array1OfPoints&          Poles,
1030  const TColStd_Array1OfReal*    Weights,
1031  Array1OfPoints&                CachePoles,
1032  TColStd_Array1OfReal*          CacheWeights)
1033 {
1034   Standard_Integer ii,
1035   Dimension,
1036   LocalIndex,
1037   index = 0 ;
1038   Standard_Real    u = U,
1039   LocalValue ;
1040   Standard_Boolean rational;
1041
1042   BSplCLib_DataContainer dc(Degree);
1043   PrepareEval(u,
1044               index,
1045               Dimension,
1046               rational,
1047               Degree,
1048               Periodic,
1049               Poles,
1050               Weights,
1051               FlatKnots,
1052               (BSplCLib::NoMults()),
1053               dc);
1054   //
1055   // Watch out : PrepareEval checks if locally the Bspline is polynomial
1056   // therefore rational flag could be set to False even if there are
1057   // Weights. The Dimension is set accordingly.
1058   //
1059
1060   BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1061
1062   LocalValue = 1.0e0 ;
1063   LocalIndex = 0 ;
1064
1065   if (rational) {
1066
1067     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1068       CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1069       LocalIndex += Dimension_gen + 1;
1070       LocalValue *= SpanDomain / (Standard_Real) ii ;
1071     }
1072
1073     LocalIndex = Dimension_gen;
1074     LocalValue = 1.0e0 ;
1075     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1076       (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ;
1077       LocalIndex += Dimension_gen + 1;
1078       LocalValue *= SpanDomain / (Standard_Real) ii ;
1079     }
1080   }
1081   else {
1082
1083     for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1084       CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1085       LocalIndex += Dimension_gen;
1086       LocalValue *= SpanDomain / (Standard_Real) ii ;
1087     }
1088
1089     if (Weights != NULL) {
1090       for (ii = 1 ; ii <= Degree + 1 ; ii++)
1091         (*CacheWeights)(ii) = 0.0e0 ;
1092       (*CacheWeights)(1) = 1.0e0 ;
1093     }
1094   }
1095 }
1096
1097 void BSplCLib::BuildCache(const Standard_Real          theParameter,
1098                           const Standard_Real          theSpanDomain,
1099                           const Standard_Boolean       thePeriodicFlag,
1100                           const Standard_Integer       theDegree,
1101                           const Standard_Integer       theSpanIndex,
1102                           const TColStd_Array1OfReal&  theFlatKnots,
1103                           const Array1OfPoints&        thePoles,
1104                           const TColStd_Array1OfReal*  theWeights,
1105                                 TColStd_Array2OfReal&  theCacheArray)
1106 {
1107   Standard_Real    aParam = theParameter;
1108   Standard_Integer anIndex = theSpanIndex;
1109   Standard_Integer aDimension;
1110   Standard_Boolean isRational;
1111
1112   BSplCLib_DataContainer dc(theDegree);
1113   PrepareEval(aParam,
1114               anIndex,
1115               aDimension,
1116               isRational,
1117               theDegree,
1118               thePeriodicFlag,
1119               thePoles,
1120               theWeights,
1121               theFlatKnots,
1122               (BSplCLib::NoMults()),
1123               dc);
1124   //
1125   // Watch out : PrepareEval checks if locally the Bspline is polynomial
1126   // therefore rational flag could be set to False even if there are
1127   // Weights. The Dimension is set accordingly.
1128   //
1129
1130   Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
1131     (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
1132
1133   BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1134
1135   Standard_Real aCoeff = 1.0;
1136   Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1137   Standard_Real* aPolyCoeffs = dc.poles;
1138
1139   for (Standard_Integer i = 0; i <= theDegree; i++)
1140   {
1141     for (Standard_Integer j = 0; j< aDimension; j++)
1142       aCache[j] = aPolyCoeffs[j] * aCoeff;
1143     aCoeff *= theSpanDomain / (i+1);
1144     aPolyCoeffs += aDimension;
1145     aCache += aDimension;
1146     if (aCacheShift > aDimension)
1147     {
1148       aCache[0] = 0.0;
1149       aCache++;
1150     }
1151   }
1152
1153   if (aCacheShift > aDimension)
1154     theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1155 }
1156
1157 //=======================================================================
1158 //function : Interpolate
1159 //purpose  :
1160 //=======================================================================
1161
1162 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
1163                             const TColStd_Array1OfReal&    FlatKnots,
1164                             const TColStd_Array1OfReal&    Parameters,
1165                             const TColStd_Array1OfInteger& ContactOrderArray,
1166                             Array1OfPoints&                Poles,
1167                             Standard_Integer&              InversionProblem)
1168
1169 {
1170   Standard_Real *PArray ;
1171
1172   PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1173
1174   BSplCLib::Interpolate(Degree,
1175                         FlatKnots,
1176                         Parameters,
1177                         ContactOrderArray,
1178                         Dimension_gen,
1179                         PArray[0],
1180                         InversionProblem) ;
1181 }
1182
1183 //=======================================================================
1184 //function : Interpolate
1185 //purpose  :
1186 //=======================================================================
1187
1188 void  BSplCLib::Interpolate(const Standard_Integer         Degree,
1189                             const TColStd_Array1OfReal&    FlatKnots,
1190                             const TColStd_Array1OfReal&    Parameters,
1191                             const TColStd_Array1OfInteger& ContactOrderArray,
1192                             Array1OfPoints&                Poles,
1193                             TColStd_Array1OfReal&          Weights,
1194                             Standard_Integer&              InversionProblem)
1195 {
1196   Standard_Real *PArray,
1197   * WArray ;
1198   PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1199   WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1200   BSplCLib::Interpolate(Degree,
1201                         FlatKnots,
1202                         Parameters,
1203                         ContactOrderArray,
1204                         Dimension_gen,
1205                         PArray[0],
1206                         WArray[0],
1207                         InversionProblem) ;
1208 }
1209
1210 //=======================================================================
1211 //function : MovePoint
1212 //purpose  : Find the new poles which allows  an old point (with a
1213 //           given  u as parameter) to reach a new position
1214 //=======================================================================
1215
1216 void BSplCLib::MovePoint (const Standard_Real            U,
1217                           const Vector&                  Displ,
1218                           const Standard_Integer         Index1,
1219                           const Standard_Integer         Index2,
1220                           const Standard_Integer         Degree,
1221                           const Standard_Boolean         Rational,
1222                           const Array1OfPoints&          Poles,
1223                           const TColStd_Array1OfReal&    Weights,
1224                           const TColStd_Array1OfReal&    FlatKnots,
1225                           Standard_Integer&              FirstIndex,
1226                           Standard_Integer&              LastIndex,
1227                           Array1OfPoints&                NewPoles)
1228 {
1229   // calculate the BSplineBasis in the parameter U
1230   Standard_Integer FirstNonZeroBsplineIndex;
1231   math_Matrix BSplineBasis(1, 1,
1232                            1, Degree+1);
1233   Standard_Integer ErrorCode =
1234     BSplCLib::EvalBsplineBasis(0,
1235                                Degree+1,
1236                                FlatKnots,
1237                                U,
1238                                FirstNonZeroBsplineIndex,
1239                                BSplineBasis);
1240   if (ErrorCode != 0) {
1241     FirstIndex = 0;
1242     LastIndex = 0;
1243
1244     for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1245       NewPoles(i) = Poles(i);
1246     }
1247     return;
1248   }
1249
1250   // find the span which is predominant for parameter U
1251   FirstIndex = FirstNonZeroBsplineIndex;
1252   LastIndex = FirstNonZeroBsplineIndex + Degree ;
1253   if (FirstIndex < Index1) FirstIndex = Index1;
1254   if (LastIndex > Index2) LastIndex = Index2;
1255
1256   Standard_Real maxValue = 0.0;
1257   Standard_Integer i, kk1=0, kk2, ii;
1258
1259   for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1260        i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1261     if (BSplineBasis(1,i) > maxValue) {
1262       kk1 = i + FirstNonZeroBsplineIndex - 1;
1263       maxValue = BSplineBasis(1,i);
1264     }
1265   }
1266
1267   // find a kk2 if symetriy
1268   kk2 = kk1;
1269   i = kk1 - FirstNonZeroBsplineIndex + 2;
1270   if ((kk1+1) <= LastIndex) {
1271     if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1272       kk2 = kk1+1;
1273     }
1274   }
1275
1276   // compute the vector of displacement
1277   Standard_Real D1 = 0.0;
1278   Standard_Real D2 = 0.0;
1279   Standard_Real hN, Coef, Dval;
1280
1281   for (i = 1; i <= Degree+1; i++) {
1282     ii = i + FirstNonZeroBsplineIndex - 1;
1283     if (Rational) {
1284       hN = Weights(ii)*BSplineBasis(1, i);
1285       D2 += hN;
1286     }
1287     else {
1288       hN = BSplineBasis(1, i);
1289     }
1290     if (ii >= FirstIndex && ii <= LastIndex) {
1291       if (ii < kk1) {
1292         Dval = kk1-ii;
1293       }
1294       else if (ii > kk2) {
1295         Dval = ii - kk2;
1296       }
1297       else {
1298         Dval = 0.0;
1299       }
1300       D1 += 1./(Dval + 1.) * hN;
1301     }
1302   }
1303
1304   if (Rational) {
1305     Coef = D2/D1;
1306   }
1307   else {
1308     Coef = 1./D1;
1309   }
1310
1311   // compute the new poles
1312
1313   for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1314     if (i >= FirstIndex && i <= LastIndex) {
1315       if (i < kk1) {
1316         Dval = kk1-i;
1317       }
1318       else if (i > kk2) {
1319         Dval = i - kk2;
1320       }
1321       else {
1322         Dval = 0.0;
1323       }
1324       NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1325     }
1326     else {
1327       NewPoles(i) = Poles(i);
1328     }
1329   }
1330 }
1331
1332 //=======================================================================
1333 //function : MovePoint
1334 //purpose  : Find the new poles which allows  an old point (with a
1335 //           given  u as parameter) to reach a new position
1336 //=======================================================================
1337
1338 //=======================================================================
1339 //function : MovePointAndTangent
1340 //purpose  :
1341 //=======================================================================
1342
1343 void BSplCLib::MovePointAndTangent (const Standard_Real      U,
1344                                     const Vector&            Delta,
1345                                     const Vector&            DeltaDerivatives,
1346                                     const Standard_Real      Tolerance,
1347                                     const Standard_Integer   Degree,
1348                                     const Standard_Boolean   Rational,
1349                                     const Standard_Integer   StartingCondition,
1350                                     const Standard_Integer   EndingCondition,
1351                                     const Array1OfPoints&    Poles,
1352                                     const TColStd_Array1OfReal&    Weights,
1353                                     const TColStd_Array1OfReal&    FlatKnots,
1354                                     Array1OfPoints&                NewPoles,
1355                                     Standard_Integer &             ErrorStatus)
1356 {
1357   Standard_Real *delta_array,
1358   *delta_derivative_array,
1359   *poles_array,
1360   *new_poles_array ;
1361
1362   Standard_Integer num_poles ;
1363   num_poles = Poles.Length() ;
1364
1365   if (NewPoles.Length() != num_poles) {
1366     throw Standard_ConstructionError();
1367   }
1368   delta_array = (Standard_Real *)  &Delta ;
1369   delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1370   poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1371
1372   new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1373   MovePointAndTangent(U,
1374                       Dimension_gen,
1375                       delta_array[0],
1376                       delta_derivative_array[0],
1377                       Tolerance,
1378                       Degree,
1379                       Rational,
1380                       StartingCondition,
1381                       EndingCondition,
1382                       poles_array[0],
1383                       Weights,
1384                       FlatKnots,
1385                       new_poles_array[0],
1386                       ErrorStatus) ;
1387 }
1388
1389 //=======================================================================
1390 //function : Resolution
1391 //purpose  :
1392 //=======================================================================
1393
1394 void BSplCLib::Resolution(const Array1OfPoints&          Poles,
1395                           const TColStd_Array1OfReal*    Weights,
1396                           const Standard_Integer         NumPoles,
1397                           const TColStd_Array1OfReal&    FlatKnots,
1398                           const Standard_Integer         Degree,
1399                           const Standard_Real            Tolerance3D,
1400                           Standard_Real&                 UTolerance)
1401 {
1402   Standard_Real *PolesArray ;
1403   PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1404   BSplCLib::Resolution(PolesArray[0],
1405                        Dimension_gen,
1406                        NumPoles,
1407                        Weights,
1408                        FlatKnots,
1409                        Degree,
1410                        Tolerance3D,
1411                        UTolerance) ;
1412 }
1413
1414 //=======================================================================
1415 //function : FunctionMultiply
1416 //purpose  :
1417 //=======================================================================
1418
1419 void BSplCLib::FunctionMultiply
1420 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1421  const Standard_Integer              BSplineDegree,
1422  const TColStd_Array1OfReal &        BSplineFlatKnots,
1423  const Array1OfPoints &              Poles,
1424  const TColStd_Array1OfReal &        FlatKnots,
1425  const Standard_Integer              NewDegree,
1426  Array1OfPoints &                   NewPoles,
1427  Standard_Integer &                  Status)
1428 {
1429   Standard_Integer num_bspline_poles =
1430     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1431   Standard_Integer num_new_poles =
1432     FlatKnots.Length() - NewDegree - 1 ;
1433
1434   if (Poles.Length() != num_bspline_poles ||
1435       NewPoles.Length() != num_new_poles) {
1436     throw Standard_ConstructionError();
1437   }
1438   Standard_Real  * array_of_poles =
1439     (Standard_Real *) &Poles(Poles.Lower()) ;
1440   Standard_Real  * array_of_new_poles =
1441     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1442   BSplCLib::FunctionMultiply(FunctionPtr,
1443                              BSplineDegree,
1444                              BSplineFlatKnots,
1445                              Dimension_gen,
1446                              array_of_poles[0],
1447                              FlatKnots,
1448                              NewDegree,
1449                              array_of_new_poles[0],
1450                              Status) ;
1451 }
1452
1453 //=======================================================================
1454 //function : FunctionReparameterise
1455 //purpose  :
1456 //=======================================================================
1457
1458 void BSplCLib::FunctionReparameterise
1459 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1460  const Standard_Integer              BSplineDegree,
1461  const TColStd_Array1OfReal &        BSplineFlatKnots,
1462  const Array1OfPoints &              Poles,
1463  const TColStd_Array1OfReal &        FlatKnots,
1464  const Standard_Integer              NewDegree,
1465  Array1OfPoints &                   NewPoles,
1466  Standard_Integer &                  Status)
1467 {
1468   Standard_Integer num_bspline_poles =
1469     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1470   Standard_Integer num_new_poles =
1471     FlatKnots.Length() - NewDegree - 1 ;
1472
1473   if (Poles.Length() != num_bspline_poles ||
1474       NewPoles.Length() != num_new_poles) {
1475     throw Standard_ConstructionError();
1476   }
1477   Standard_Real  * array_of_poles =
1478     (Standard_Real *) &Poles(Poles.Lower()) ;
1479   Standard_Real  * array_of_new_poles =
1480     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1481   BSplCLib::FunctionReparameterise(FunctionPtr,
1482                                    BSplineDegree,
1483                                    BSplineFlatKnots,
1484                                    Dimension_gen,
1485                                    array_of_poles[0],
1486                                    FlatKnots,
1487                                    NewDegree,
1488                                    array_of_new_poles[0],
1489                                    Status) ;
1490 }