src/BSplCLib/BSplCLib_2.cxx

   1 // Created on: 1995-01-17
   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 // xab : modified 15-Mar-95 : added BuildBSplMatrix,FactorBandedMatrix,
  18 //                            EvalBsplineBasis,
  19 //                            EvalPolynomial : Horners method
  20
  21 #include <Standard_Stream.hxx>
  22
  23 #include <BSplCLib.hxx>
  24 #include <gp_Mat2d.hxx>
  25 #include <PLib.hxx>
  26 #include <Standard_ConstructionError.hxx>
  27 #include <TColStd_Array1OfReal.hxx>
  28 #include <TColStd_Array1OfInteger.hxx>
  29 #include <TColStd_HArray1OfReal.hxx>
  30 #include <TColStd_HArray1OfInteger.hxx>
  31
  32 #include <math_Matrix.hxx>
  33
  34 //=======================================================================
  35 //struct : BSplCLib_DataContainer
  36 //purpose: Auxiliary structure providing buffers for poles and knots used in
  37 //         evaluation of bspline (allocated in the stack)
  38 //=======================================================================
  39
  40 struct BSplCLib_DataContainer
  41 {
  42   BSplCLib_DataContainer(Standard_Integer Degree)
  43   {
  44     (void)Degree; // avoid compiler warning
  45     Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
  46         BSplCLib::MaxDegree() > 25,
  47         "BSplCLib: bspline degree is greater than maximum supported");
  48   }
  49
  50   Standard_Real poles[2*(25+1)];
  51   Standard_Real knots[2*25];
  52   Standard_Real ders[4];
  53 };
  54
  55 // methods for 1 dimensional BSplines
  56
  57 //=======================================================================
  58 //function : BuildEval
  59 //purpose  : builds the local array for evaluation
  60 //=======================================================================
  61
  62 void  BSplCLib::BuildEval(const Standard_Integer         Degree,
  63                           const Standard_Integer         Index,
  64                           const TColStd_Array1OfReal&    Poles,
  65                           const TColStd_Array1OfReal&    Weights,
  66                           Standard_Real&                 LP)
  67 {
  68   Standard_Integer PLower = Poles.Lower();
  69   Standard_Integer PUpper = Poles.Upper();
  70   Standard_Integer i;
  71   Standard_Integer ip = PLower + Index - 1;
  72   Standard_Real w, *pole = &LP;
  73   if (&Weights == NULL) {
  74
  75     for (i = 0; i <= Degree; i++) {
  76       ip++;
  77       if (ip > PUpper) ip = PLower;
  78       pole[0] = Poles(ip);
  79       pole += 1;
  80     }
  81   }
  82   else {
  83
  84     for (i = 0; i <= Degree; i++) {
  85       ip++;
  86       if (ip > PUpper) ip = PLower;
  87       pole[1] = w = Weights(ip);
  88       pole[0] = Poles(ip) * w;
  89       pole += 2;
  90     }
  91   }
  92 }
  93
  94 //=======================================================================
  95 //function : PrepareEval
  96 //purpose  : stores data for Eval in the local arrays
  97 //           dc.poles and dc.knots
  98 //=======================================================================
  99
 100 static void PrepareEval
 101 (Standard_Real&                 u,
 102  Standard_Integer&              index,
 103  Standard_Integer&              dim,
 104  Standard_Boolean&              rational,
 105  const Standard_Integer         Degree,
 106  const Standard_Boolean         Periodic,
 107  const TColStd_Array1OfReal&    Poles,
 108  const TColStd_Array1OfReal&    Weights,
 109  const TColStd_Array1OfReal&    Knots,
 110  const TColStd_Array1OfInteger& Mults,
 111  BSplCLib_DataContainer&        dc)
 112 {
 113   // Set the Index
 114   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
 115
 116   // make the knots
 117   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
 118   if (&Mults == NULL)
 119     index -= Knots.Lower() + Degree;
 120   else
 121     index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
 122
 123   // check truly rational
 124   rational = (&Weights != NULL);
 125   if (rational) {
 126     Standard_Integer WLower = Weights.Lower() + index;
 127     rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
 128   }
 129
 130   // make the poles
 131   if(rational) {
 132     dim = 2;
 133     BSplCLib::BuildEval(Degree, index, Poles, Weights              , *dc.poles);
 134   }
 135   else {
 136     dim = 1;
 137     BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
 138   }
 139 }
 140
 141 //=======================================================================
 142 //function : D0
 143 //purpose  :
 144 //=======================================================================
 145
 146 void BSplCLib::D0
 147 (const Standard_Real            U,
 148  const Standard_Integer         Index,
 149  const Standard_Integer         Degree,
 150  const Standard_Boolean         Periodic,
 151  const TColStd_Array1OfReal&    Poles,
 152  const TColStd_Array1OfReal&    Weights,
 153  const TColStd_Array1OfReal&    Knots,
 154  const TColStd_Array1OfInteger& Mults,
 155  Standard_Real&                 P)
 156 {
 157   Standard_Integer dim,index = Index;
 158   Standard_Real    u = U;
 159   Standard_Boolean rational;
 160   BSplCLib_DataContainer dc(Degree);
 161   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 162   BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
 163   if (rational) P = dc.poles[0] / dc.poles[1];
 164   else          P = dc.poles[0];
 165 }
 166
 167 //=======================================================================
 168 //function : D1
 169 //purpose  :
 170 //=======================================================================
 171
 172 void BSplCLib::D1
 173 (const Standard_Real            U,
 174  const Standard_Integer         Index,
 175  const Standard_Integer         Degree,
 176  const Standard_Boolean         Periodic,
 177  const TColStd_Array1OfReal&    Poles,
 178  const TColStd_Array1OfReal&    Weights,
 179  const TColStd_Array1OfReal&    Knots,
 180  const TColStd_Array1OfInteger& Mults,
 181  Standard_Real&                 P,
 182  Standard_Real&                 V)
 183 {
 184   Standard_Integer dim,index = Index;
 185   Standard_Real    u = U;
 186   Standard_Boolean rational;
 187   BSplCLib_DataContainer dc(Degree);
 188   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 189   BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
 190   Standard_Real *result = dc.poles;
 191   if (rational) {
 192     PLib::RationalDerivative(Degree,1,1,*dc.poles,*dc.ders);
 193     result = dc.ders;
 194   }
 195   P = result[0];
 196   V = result[1];
 197 }
 198
 199 //=======================================================================
 200 //function : D2
 201 //purpose  :
 202 //=======================================================================
 203
 204 void BSplCLib::D2
 205 (const Standard_Real            U,
 206  const Standard_Integer         Index,
 207  const Standard_Integer         Degree,
 208  const Standard_Boolean         Periodic,
 209  const TColStd_Array1OfReal&    Poles,
 210  const TColStd_Array1OfReal&    Weights,
 211  const TColStd_Array1OfReal&    Knots,
 212  const TColStd_Array1OfInteger& Mults,
 213  Standard_Real&                 P,
 214  Standard_Real&                 V1,
 215  Standard_Real&                 V2)
 216 {
 217   Standard_Integer dim,index = Index;
 218   Standard_Real    u = U;
 219   Standard_Boolean rational;
 220   BSplCLib_DataContainer dc(Degree);
 221   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 222   BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
 223   Standard_Real *result = dc.poles;
 224   if (rational) {
 225     PLib::RationalDerivative(Degree,2,1,*dc.poles,*dc.ders);
 226     result = dc.ders;
 227   }
 228   P  = result[0];
 229   V1 = result[1];
 230   if (!rational && (Degree < 2)) V2 = 0.;
 231   else                           V2 = result[2];
 232 }
 233
 234 //=======================================================================
 235 //function : D3
 236 //purpose  :
 237 //=======================================================================
 238
 239 void BSplCLib::D3
 240 (const Standard_Real            U,
 241  const Standard_Integer         Index,
 242  const Standard_Integer         Degree,
 243  const Standard_Boolean         Periodic,
 244  const TColStd_Array1OfReal&    Poles,
 245  const TColStd_Array1OfReal&    Weights,
 246  const TColStd_Array1OfReal&    Knots,
 247  const TColStd_Array1OfInteger& Mults,
 248  Standard_Real&                 P,
 249  Standard_Real&                 V1,
 250  Standard_Real&                 V2,
 251  Standard_Real&                 V3)
 252 {
 253   Standard_Integer dim,index = Index;
 254   Standard_Real    u = U;
 255   Standard_Boolean rational;
 256   BSplCLib_DataContainer dc(Degree);
 257   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 258   BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
 259   Standard_Real *result = dc.poles;
 260   if (rational) {
 261     PLib::RationalDerivative(Degree,3,1,*dc.poles,*dc.ders);
 262     result = dc.ders;
 263   }
 264   P  = result[0];
 265   V1 = result[1];
 266   if (!rational && (Degree < 2)) V2 = 0.;
 267   else                           V2 = result[2];
 268   if (!rational && (Degree < 3)) V3 = 0.;
 269   else                           V3 = result[3];
 270 }
 271
 272 //=======================================================================
 273 //function : DN
 274 //purpose  :
 275 //=======================================================================
 276
 277 void BSplCLib::DN
 278 (const Standard_Real            U,
 279  const Standard_Integer         N,
 280  const Standard_Integer         Index,
 281  const Standard_Integer         Degree,
 282  const Standard_Boolean         Periodic,
 283  const TColStd_Array1OfReal&    Poles,
 284  const TColStd_Array1OfReal&    Weights,
 285  const TColStd_Array1OfReal&    Knots,
 286  const TColStd_Array1OfInteger& Mults,
 287  Standard_Real&                 VN)
 288 {
 289   Standard_Integer dim,index = Index;
 290   Standard_Real    u = U;
 291   Standard_Boolean rational;
 292   BSplCLib_DataContainer dc(Degree);
 293   PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
 294   BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
 295   if (rational) {
 296     Standard_Real v;
 297     PLib::RationalDerivative(Degree,N,1,*dc.poles,v,Standard_False);
 298     VN = v;
 299   }
 300   else {
 301     if (N > Degree) VN = 0.;
 302     else            VN = dc.poles[N];
 303   }
 304 }
 305
 306 //=======================================================================
 307 //function : Build BSpline Matrix
 308 //purpose  : Builds the Bspline Matrix
 309 //=======================================================================
 310
 311 Standard_Integer
 312 BSplCLib::BuildBSpMatrix(const  TColStd_Array1OfReal&     Parameters,
 313                          const  TColStd_Array1OfInteger&  ContactOrderArray,
 314                          const  TColStd_Array1OfReal&     FlatKnots,
 315                          const  Standard_Integer          Degree,
 316                          math_Matrix&                     Matrix,
 317                          Standard_Integer&                UpperBandWidth,
 318                          Standard_Integer&                LowerBandWidth)
 319 {
 320   Standard_Integer ii,
 321   jj,
 322   Index,
 323   ErrorCode,
 324   ReturnCode = 0,
 325   FirstNonZeroBsplineIndex,
 326   BandWidth,
 327   MaxOrder = 21,
 328   Order ;
 329
 330   math_Matrix   BSplineBasis(1, MaxOrder,
 331                              1, MaxOrder) ;
 332
 333   Order = Degree + 1 ;
 334   UpperBandWidth = Degree ;
 335   LowerBandWidth = Degree ;
 336   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
 337   if (Matrix.LowerRow() != Parameters.Lower()  ||
 338       Matrix.UpperRow() != Parameters.Upper()  ||
 339       Matrix.LowerCol() != 1  ||
 340       Matrix.UpperCol() != BandWidth) {
 341     ReturnCode = 1;
 342     goto FINISH ;
 343   }
 344
 345   for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
 346     ErrorCode =
 347       BSplCLib::EvalBsplineBasis(1,
 348                                  ContactOrderArray(ii),
 349                                  Order,
 350                                  FlatKnots,
 351                                  Parameters(ii),
 352
 353                                  FirstNonZeroBsplineIndex,
 354                                  BSplineBasis) ;
 355     if (ErrorCode != 0) {
 356       ReturnCode = 2 ;
 357       goto FINISH ;
 358     }
 359     Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
 360
 361     for (jj = 1 ; jj < Index ; jj++) {
 362       Matrix.Value(ii,jj) = 0.0e0 ;
 363     }
 364
 365     for (jj = 1 ; jj <= Order ; jj++) {
 366       Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
 367       Index += 1 ;
 368     }
 369
 370     for (jj = Index ; jj <= BandWidth ; jj++) {
 371       Matrix.Value(ii,jj) = 0.0e0 ;
 372     }
 373   }
 374   FINISH : ;
 375   return (ReturnCode) ;
 376 }
 377
 378 //=======================================================================
 379 //function : Makes LU decompositiomn without Pivoting
 380 //purpose  : Builds the Bspline Matrix
 381 //=======================================================================
 382
 383 Standard_Integer
 384 BSplCLib::FactorBandedMatrix(math_Matrix&   Matrix,
 385                              const Standard_Integer UpperBandWidth,
 386                              const Standard_Integer LowerBandWidth,
 387                              Standard_Integer&  PivotIndexProblem)
 388 {
 389   Standard_Integer ii,
 390   jj,
 391   kk,
 392   Index,
 393   MinIndex,
 394   MaxIndex,
 395   ReturnCode = 0,
 396   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
 397
 398   Standard_Real Inverse ;
 399   PivotIndexProblem = 0 ;
 400
 401   for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
 402     MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
 403
 404     for (jj = MinIndex ; jj <= LowerBandWidth  ; jj++) {
 405       Index = ii - LowerBandWidth + jj - 1 ;
 406       Inverse = Matrix(Index ,LowerBandWidth + 1) ;
 407       if (Abs(Inverse) > RealSmall()) {
 408         Inverse = -1.0e0/Inverse ;
 409       }
 410       else {
 411         ReturnCode = 1 ;
 412         PivotIndexProblem = Index ;
 413         goto FINISH ;
 414       }
 415       Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
 416       MaxIndex = BandWidth + Index - ii ;
 417
 418       for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
 419         Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
 420       }
 421     }
 422   }
 423   FINISH :
 424     return (ReturnCode) ;
 425 }
 426
 427 //=======================================================================
 428 //function : Build BSpline Matrix
 429 //purpose  : Builds the Bspline Matrix
 430 //=======================================================================
 431
 432 Standard_Integer
 433 BSplCLib::EvalBsplineBasis
 434 //(const Standard_Integer              Side, // = 1 rigth side, -1 left side
 435 (const Standard_Integer              , // = 1 rigth side, -1 left side
 436  const  Standard_Integer              DerivativeRequest,
 437  const  Standard_Integer              Order,
 438  const  TColStd_Array1OfReal&         FlatKnots,
 439  const  Standard_Real                 Parameter,
 440  Standard_Integer&             FirstNonZeroBsplineIndex,
 441  math_Matrix&                  BsplineBasis,
 442  Standard_Boolean              isPeriodic)
 443 {
 444   // the matrix must have at least DerivativeRequest + 1
 445   //   row and Order columns
 446   // the result are stored in the following way in
 447   // the Bspline matrix
 448   // Let i be the FirstNonZeroBsplineIndex and
 449   // t be the parameter value, k the order of the
 450   // knot vector, r the DerivativeRequest :
 451   //
 452   //   B (t)   B (t)                     B (t)
 453   //    i       i+1                       i+k-1
 454   //
 455   //    (1)     (1)                       (1)
 456   //   B (t)   B (t)                     B (t)
 457   //    i       i+1                       i+k-1
 458   //
 459   //
 460   //
 461   //
 462   //    (r)     (r)                       (r)
 463   //   B (t)   B (t)                     B (t)
 464   //    i       i+1                       i+k-1
 465   //
 466   Standard_Integer
 467     ReturnCode,
 468   ii,
 469   pp,
 470   qq,
 471   ss,
 472   NumPoles,
 473   LocalRequest ;
 474 //  ,Index ;
 475
 476   Standard_Real NewParameter,
 477   Inverse,
 478   Factor,
 479   LocalInverse,
 480   Saved ;
 481 // , *FlatKnotsArray ;
 482
 483   ReturnCode = 0 ;
 484   FirstNonZeroBsplineIndex = 0 ;
 485   LocalRequest = DerivativeRequest ;
 486   if (DerivativeRequest >= Order) {
 487     LocalRequest = Order - 1 ;
 488   }
 489
 490   if (BsplineBasis.LowerCol() != 1    ||
 491       BsplineBasis.UpperCol() < Order ||
 492       BsplineBasis.LowerRow() != 1    ||
 493       BsplineBasis.UpperRow() <= LocalRequest) {
 494     ReturnCode = 1;
 495     goto FINISH ;
 496   }
 497   NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
 498   BSplCLib::LocateParameter(Order - 1,
 499                             FlatKnots,
 500                             Parameter,
 501                             isPeriodic,
 502                             Order,
 503                             NumPoles+1,
 504                             ii,
 505                             NewParameter) ;
 506
 507   FirstNonZeroBsplineIndex = ii - Order + 1 ;
 508
 509   BsplineBasis(1,1) = 1.0e0 ;
 510   LocalRequest = DerivativeRequest ;
 511   if (DerivativeRequest >= Order) {
 512     LocalRequest = Order - 1 ;
 513   }
 514
 515   for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
 516     BsplineBasis(1,qq) = 0.0e0 ;
 517
 518     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 519       //
 520       // this should be always invertible if ii is correctly computed
 521       //
 522       Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
 523         / (FlatKnots(ii + pp)   - FlatKnots(ii - qq + pp + 1)) ;
 524       Saved = Factor *    BsplineBasis(1,pp) ;
 525       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
 526       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
 527       BsplineBasis(1,qq) = Saved ;
 528     }
 529   }
 530
 531   for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
 532
 533     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 534       BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
 535     }
 536     BsplineBasis(1,qq) = 0.0e0 ;
 537
 538     for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
 539       BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
 540     }
 541
 542     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 543       Inverse = 1.0e0 / (FlatKnots(ii + pp)  - FlatKnots(ii - qq + pp + 1)) ;
 544       Factor  =  (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
 545       Saved = Factor *                 BsplineBasis(1,pp) ;
 546       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
 547       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
 548       BsplineBasis(1,qq) = Saved ;
 549       LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
 550
 551       for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
 552         Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
 553         BsplineBasis(Order - ss + 2, pp) *= - LocalInverse  ;
 554         BsplineBasis(Order - ss + 2, pp) +=   BsplineBasis(Order - ss + 2,qq) ;
 555         BsplineBasis(Order - ss + 2,qq) = Saved ;
 556       }
 557     }
 558   }
 559   FINISH :
 560     return (ReturnCode) ;
 561 }
 562
 563 //=======================================================================
 564 //function : MovePointAndTangent
 565 //purpose  :
 566 //=======================================================================
 567
 568 void BSplCLib::MovePointAndTangent(const Standard_Real    U,
 569                                    const Standard_Integer ArrayDimension,
 570                                    Standard_Real    &Delta,
 571                                    Standard_Real    &DeltaDerivatives,
 572                                    const Standard_Real    Tolerance,
 573                                    const Standard_Integer Degree,
 574                                    const Standard_Boolean Rational,
 575                                    const Standard_Integer StartingCondition,
 576                                    const Standard_Integer EndingCondition,
 577                                    Standard_Real&         Poles,
 578                                    const TColStd_Array1OfReal&   Weights,
 579                                    const TColStd_Array1OfReal&   FlatKnots,
 580                                    Standard_Real&        NewPoles,
 581                                    Standard_Integer&     ErrorStatus)
 582 {
 583   Standard_Integer num_poles,
 584   num_knots,
 585   ii,
 586   jj,
 587   conditions,
 588   start_num_poles,
 589   end_num_poles,
 590   index,
 591   start_index,
 592   end_index,
 593   other_index,
 594   type,
 595   order ;
 596
 597   Standard_Real    new_parameter,
 598   value,
 599   divide,
 600   end_value,
 601   start_value,
 602   *poles_array,
 603   *new_poles_array,
 604   *delta_array,
 605   *derivatives_array,
 606   *weights_array ;
 607
 608   ErrorStatus = 0 ;
 609   weights_array = NULL ;
 610   if (Rational) {
 611     weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
 612   }
 613
 614   poles_array = &Poles ;
 615   new_poles_array = &NewPoles ;
 616   delta_array = &Delta ;
 617   derivatives_array = &DeltaDerivatives ;
 618   order = Degree + 1 ;
 619   num_knots = FlatKnots.Length() ;
 620   num_poles = num_knots - order ;
 621   conditions = StartingCondition + EndingCondition + 4 ;
 622   //
 623   // check validity of input data
 624   //
 625   if (StartingCondition >= -1 &&
 626       StartingCondition <= Degree &&
 627       EndingCondition >= -1 &&
 628       EndingCondition <= Degree &&
 629       conditions <= num_poles) {
 630     //
 631     // check the parameter is within bounds
 632     //
 633     start_index = FlatKnots.Lower() + Degree ;
 634     end_index = FlatKnots.Upper() - Degree ;
 635     //
 636     //  check if there is enough room to move the poles
 637     //
 638     conditions = 1 ;
 639     if (StartingCondition == -1) {
 640       conditions = conditions && (FlatKnots(start_index) <=  U) ;
 641     }
 642     else {
 643       conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
 644     }
 645     if (EndingCondition == -1) {
 646       conditions = conditions && (FlatKnots(end_index) >= U) ;
 647     }
 648     else {
 649       conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ;
 650     }
 651
 652     if (conditions) {
 653       //
 654       // build 2 auxialiary functions
 655       //
 656       TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
 657       TColStd_Array1OfReal first_function  (1,num_poles) ;
 658       TColStd_Array1OfReal second_function (1,num_poles) ;
 659
 660       BuildSchoenbergPoints(Degree,
 661                             FlatKnots,
 662                             schoenberg_points) ;
 663       start_index = StartingCondition + 2 ;
 664       end_index = num_poles - EndingCondition - 1 ;
 665       LocateParameter(schoenberg_points,
 666                       U,
 667                       Standard_False,
 668                       start_index,
 669                       end_index,
 670                       index,
 671                       new_parameter,
 672                       0, 1) ;
 673
 674       if (index == start_index) {
 675         other_index = index + 1 ;
 676       }
 677       else if (index == end_index) {
 678         other_index = index -1  ;
 679       }
 680       else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
 681         other_index = index - 1 ;
 682       }
 683       else {
 684         other_index = index + 1 ;
 685       }
 686       type = 3 ;
 687
 688       start_num_poles = StartingCondition + 2 ;
 689       end_num_poles   = num_poles - EndingCondition - 1 ;
 690       if (start_num_poles == 1) {
 691         start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
 692         start_value = schoenberg_points(1) - start_value ;
 693       }
 694       else {
 695         start_value = schoenberg_points(start_num_poles - 1) ;
 696       }
 697       if (end_num_poles == num_poles) {
 698         end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
 699         end_value = schoenberg_points(num_poles) + end_value ;
 700       }
 701       else {
 702         end_value = schoenberg_points(end_num_poles + 1) ;
 703       }
 704
 705       for (ii = 1 ; ii < start_num_poles ; ii++) {
 706         first_function(ii) = 0.e0 ;
 707         second_function(ii) = 0.0e0 ;
 708       }
 709
 710       for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
 711         first_function(ii) = 0.e0 ;
 712         second_function(ii) = 0.0e0 ;
 713       }
 714       divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
 715
 716       for (ii = start_num_poles  ; ii <= index ; ii++) {
 717         value = schoenberg_points(ii) - start_value ;
 718         value *= divide ;
 719         first_function(ii) = 1.0e0 ;
 720
 721         for (jj = 0 ; jj < type ; jj++) {
 722           first_function(ii) *= value ;
 723         }
 724       }
 725       divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
 726
 727       for (ii = index ; ii <= end_num_poles ; ii++) {
 728         value = end_value - schoenberg_points(ii) ;
 729         value *= divide ;
 730         first_function(ii) = 1.0e0 ;
 731
 732         for (jj = 0 ; jj < type ; jj++) {
 733           first_function(ii) *= value ;
 734         }
 735       }
 736
 737       divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
 738
 739       for (ii = start_num_poles  ; ii <= other_index ; ii++) {
 740         value = schoenberg_points(ii) - start_value ;
 741         value *= divide ;
 742         second_function(ii) = 1.0e0 ;
 743
 744         for (jj = 0 ; jj < type ; jj++) {
 745           second_function(ii) *= value ;
 746         }
 747       }
 748       divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
 749
 750       for (ii = other_index ; ii <= end_num_poles ; ii++) {
 751         value = end_value - schoenberg_points(ii) ;
 752         value *= divide ;
 753         second_function(ii) = 1.0e0 ;
 754
 755         for (jj = 0 ; jj < type ; jj++) {
 756           second_function(ii) *= value ;
 757         }
 758       }
 759
 760       //
 761       //  compute the point and derivatives of both functions
 762       //
 763       Standard_Real results[2][2],
 764       weights_results[2][2];
 765       Standard_Integer extrap_mode[2],
 766       derivative_request = 1,
 767       dimension = 1 ;
 768       Standard_Boolean periodic_flag = Standard_False ;
 769
 770       extrap_mode[0] = Degree ;
 771       extrap_mode[1] = Degree ;
 772       if (Rational) {
 773         //
 774         // evaluate in homogenised form
 775         //
 776         Eval(U,
 777              periodic_flag,
 778              derivative_request,
 779              extrap_mode[0],
 780              Degree,
 781              FlatKnots,
 782              dimension,
 783              first_function(1),
 784              weights_array[0],
 785              results[0][0],
 786              weights_results[0][0]) ;
 787
 788         Eval(U,
 789              periodic_flag,
 790              derivative_request,
 791              extrap_mode[0],
 792              Degree,
 793              FlatKnots,
 794              dimension,
 795              second_function(1),
 796              weights_array[0],
 797              results[1][0],
 798              weights_results[1][0]) ;
 799         //
 800         //  compute the rational derivatives values
 801         //
 802
 803         for (ii = 0 ; ii < 2 ; ii++) {
 804           PLib::RationalDerivatives(1,
 805                                     1,
 806                                     results[ii][0],
 807                                     weights_results[ii][0],
 808                                     results[ii][0]) ;
 809         }
 810       }
 811       else {
 812         Eval(U,
 813              Standard_False,
 814              1,
 815              extrap_mode[0],
 816              Degree,
 817              FlatKnots,
 818              1,
 819              first_function(1),
 820              results[0][0]) ;
 821
 822         Eval(U,
 823              Standard_False,
 824              1,
 825              extrap_mode[0],
 826              Degree,
 827              FlatKnots,
 828              1,
 829              second_function(1),
 830              results[1][0]) ;
 831       }
 832       gp_Mat2d  a_matrix ;
 833
 834       for (ii = 0 ; ii < 2 ; ii++) {
 835
 836         for (jj = 0 ; jj < 2 ; jj++) {
 837           a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
 838         }
 839       }
 840       a_matrix.Invert() ;
 841       TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
 842       TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
 843
 844       for( ii = 0 ; ii < ArrayDimension ; ii++) {
 845         the_a_vector(ii) =
 846           a_matrix.Value(1,1) * delta_array[ii] +
 847             a_matrix.Value(2,1) * derivatives_array[ii] ;
 848         the_b_vector(ii) =
 849           a_matrix.Value(1,2) * delta_array[ii] +
 850             a_matrix.Value(2,2) * derivatives_array[ii] ;
 851       }
 852       index = 0 ;
 853
 854       for (ii = 0 ; ii < num_poles ; ii++) {
 855
 856         for (jj = 0 ; jj < ArrayDimension ; jj++) {
 857           new_poles_array[index] = poles_array[index] ;
 858           new_poles_array[index] +=
 859             first_function(ii+1) * the_a_vector(jj) ;
 860           new_poles_array[index] +=
 861             second_function(ii+1) * the_b_vector(jj) ;
 862           index += 1 ;
 863         }
 864       }
 865     }
 866     else {
 867       ErrorStatus = 1 ;
 868     }
 869   }
 870   else {
 871     ErrorStatus = 2 ;
 872   }
 873 }
 874
 875 //=======================================================================
 876 //function : FunctionMultiply
 877 //purpose  :
 878 //=======================================================================
 879
 880 void BSplCLib::FunctionMultiply
 881 (const BSplCLib_EvaluatorFunction & FunctionPtr,
 882  const Standard_Integer             BSplineDegree,
 883  const TColStd_Array1OfReal &       BSplineFlatKnots,
 884  const Standard_Integer             PolesDimension,
 885  Standard_Real &                    Poles,
 886  const TColStd_Array1OfReal &       FlatKnots,
 887  const Standard_Integer             NewDegree,
 888  Standard_Real &                    NewPoles,
 889  Standard_Integer &                 Status)
 890 {
 891   Standard_Integer ii,
 892   jj,
 893   index ;
 894   Standard_Integer extrap_mode[2],
 895   error_code,
 896   num_new_poles,
 897   derivative_request = 0 ;
 898   Standard_Boolean  periodic_flag = Standard_False ;
 899   Standard_Real  result,
 900   start_end[2],
 901   *array_of_poles,
 902   *array_of_new_poles ;
 903
 904   array_of_poles = (Standard_Real *) &NewPoles ;
 905   extrap_mode[0] =
 906     extrap_mode[1] = BSplineDegree ;
 907   num_new_poles =
 908     FlatKnots.Length() - NewDegree - 1 ;
 909   start_end[0] = FlatKnots(NewDegree+1) ;
 910   start_end[1] = FlatKnots(num_new_poles+1) ;
 911   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
 912   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
 913   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
 914
 915   array_of_new_poles =
 916     (Standard_Real *) &new_poles_array(1) ;
 917   BuildSchoenbergPoints(NewDegree,
 918                         FlatKnots,
 919                         parameters) ;
 920   //
 921   // on recadre sur les bornes
 922   //
 923   if (parameters(1) < start_end[0]) {
 924     parameters(1) = start_end[0] ;
 925   }
 926   if (parameters(num_new_poles) > start_end[1]) {
 927     parameters(num_new_poles) = start_end[1] ;
 928   }
 929   index = 0 ;
 930
 931   for (ii = 1 ; ii <= num_new_poles ; ii++) {
 932     contact_order_array(ii) = 0 ;
 933     FunctionPtr.Evaluate (contact_order_array(ii),
 934                    start_end,
 935                    parameters(ii),
 936                    result,
 937                    error_code);
 938     if (error_code) {
 939       Status = 1 ;
 940       goto FINISH ;
 941     }
 942
 943     Eval(parameters(ii),
 944          periodic_flag,
 945          derivative_request,
 946          extrap_mode[0],
 947          BSplineDegree,
 948          BSplineFlatKnots,
 949          PolesDimension,
 950          Poles,
 951          array_of_new_poles[index]) ;
 952
 953     for (jj = 0 ; jj < PolesDimension ; jj++) {
 954       array_of_new_poles[index] *= result ;
 955       index += 1 ;
 956     }
 957   }
 958   Interpolate(NewDegree,
 959               FlatKnots,
 960               parameters,
 961               contact_order_array,
 962               PolesDimension,
 963               array_of_new_poles[0],
 964               Status) ;
 965
 966   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
 967     array_of_poles[ii] = array_of_new_poles[ii] ;
 968
 969   }
 970   FINISH :
 971     ;
 972 }
 973
 974 //=======================================================================
 975 // function : FunctionMultiply
 976 //purpose  :
 977 //=======================================================================
 978
 979 void BSplCLib::FunctionReparameterise
 980 (const BSplCLib_EvaluatorFunction & FunctionPtr,
 981  const Standard_Integer             BSplineDegree,
 982  const TColStd_Array1OfReal &       BSplineFlatKnots,
 983  const Standard_Integer             PolesDimension,
 984  Standard_Real &                    Poles,
 985  const TColStd_Array1OfReal &       FlatKnots,
 986  const Standard_Integer             NewDegree,
 987  Standard_Real &                    NewPoles,
 988  Standard_Integer &                 Status)
 989 {
 990   Standard_Integer ii,
 991 //  jj,
 992   index ;
 993   Standard_Integer extrap_mode[2],
 994   error_code,
 995   num_new_poles,
 996   derivative_request = 0 ;
 997   Standard_Boolean  periodic_flag = Standard_False ;
 998   Standard_Real  result,
 999   start_end[2],
1000   *array_of_poles,
1001   *array_of_new_poles ;
1002
1003   array_of_poles = (Standard_Real *) &NewPoles ;
1004   extrap_mode[0] =
1005     extrap_mode[1] = BSplineDegree ;
1006   num_new_poles =
1007     FlatKnots.Length() - NewDegree - 1 ;
1008   start_end[0] = FlatKnots(NewDegree+1) ;
1009   start_end[1] = FlatKnots(num_new_poles+1) ;
1010   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
1011   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1012   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
1013
1014   array_of_new_poles =
1015     (Standard_Real *) &new_poles_array(1) ;
1016   BuildSchoenbergPoints(NewDegree,
1017                         FlatKnots,
1018                         parameters) ;
1019   index = 0 ;
1020
1021   for (ii = 1 ; ii <= num_new_poles ; ii++) {
1022     contact_order_array(ii) = 0 ;
1023     FunctionPtr.Evaluate (contact_order_array(ii),
1024                    start_end,
1025                    parameters(ii),
1026                    result,
1027                    error_code);
1028     if (error_code) {
1029       Status = 1 ;
1030       goto FINISH ;
1031     }
1032
1033     Eval(result,
1034          periodic_flag,
1035          derivative_request,
1036          extrap_mode[0],
1037          BSplineDegree,
1038          BSplineFlatKnots,
1039          PolesDimension,
1040          Poles,
1041          array_of_new_poles[index]) ;
1042     index += PolesDimension ;
1043   }
1044   Interpolate(NewDegree,
1045               FlatKnots,
1046               parameters,
1047               contact_order_array,
1048               PolesDimension,
1049               array_of_new_poles[0],
1050               Status) ;
1051
1052   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1053     array_of_poles[ii] = array_of_new_poles[ii] ;
1054
1055   }
1056   FINISH :
1057     ;
1058 }
1059
1060 //=======================================================================
1061 //function : FunctionMultiply
1062 //purpose  :
1063 //=======================================================================
1064
1065 void BSplCLib::FunctionMultiply
1066 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1067  const Standard_Integer              BSplineDegree,
1068  const TColStd_Array1OfReal &        BSplineFlatKnots,
1069  const TColStd_Array1OfReal &        Poles,
1070  const TColStd_Array1OfReal &        FlatKnots,
1071  const Standard_Integer              NewDegree,
1072  TColStd_Array1OfReal &              NewPoles,
1073  Standard_Integer &                  Status)
1074 {
1075   Standard_Integer num_bspline_poles =
1076     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1077   Standard_Integer num_new_poles =
1078     FlatKnots.Length() - NewDegree - 1 ;
1079
1080   if (Poles.Length() != num_bspline_poles ||
1081       NewPoles.Length() != num_new_poles) {
1082     Standard_ConstructionError::Raise();
1083   }
1084   Standard_Real  * array_of_poles =
1085     (Standard_Real *) &Poles(Poles.Lower()) ;
1086   Standard_Real  * array_of_new_poles =
1087     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1088   BSplCLib::FunctionMultiply(FunctionPtr,
1089                              BSplineDegree,
1090                              BSplineFlatKnots,
1091                              1,
1092                              array_of_poles[0],
1093                              FlatKnots,
1094                              NewDegree,
1095                              array_of_new_poles[0],
1096                              Status) ;
1097 }
1098
1099 //=======================================================================
1100 //function : FunctionReparameterise
1101 //purpose  :
1102 //=======================================================================
1103
1104 void BSplCLib::FunctionReparameterise
1105 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1106  const Standard_Integer              BSplineDegree,
1107  const TColStd_Array1OfReal &        BSplineFlatKnots,
1108  const TColStd_Array1OfReal &        Poles,
1109  const TColStd_Array1OfReal &        FlatKnots,
1110  const Standard_Integer              NewDegree,
1111  TColStd_Array1OfReal &              NewPoles,
1112  Standard_Integer &                  Status)
1113 {
1114   Standard_Integer num_bspline_poles =
1115     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1116   Standard_Integer num_new_poles =
1117     FlatKnots.Length() - NewDegree - 1 ;
1118
1119   if (Poles.Length() != num_bspline_poles ||
1120       NewPoles.Length() != num_new_poles) {
1121     Standard_ConstructionError::Raise();
1122   }
1123   Standard_Real  * array_of_poles =
1124     (Standard_Real *) &Poles(Poles.Lower()) ;
1125   Standard_Real  * array_of_new_poles =
1126     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1127   BSplCLib::FunctionReparameterise(
1128                                    FunctionPtr,
1129                                    BSplineDegree,
1130                                    BSplineFlatKnots,
1131                                    1,
1132                                    array_of_poles[0],
1133                                    FlatKnots,
1134                                    NewDegree,
1135                                    array_of_new_poles[0],
1136                                    Status) ;
1137 }
1138
1139 //=======================================================================
1140 //function : MergeBSplineKnots
1141 //purpose  :
1142 //=======================================================================
1143 void BSplCLib::MergeBSplineKnots
1144 (const Standard_Real                 Tolerance,
1145  const Standard_Real                 StartValue,
1146  const Standard_Real                 EndValue,
1147  const Standard_Integer              Degree1,
1148  const TColStd_Array1OfReal          & Knots1,
1149  const TColStd_Array1OfInteger       & Mults1,
1150  const Standard_Integer              Degree2,
1151  const TColStd_Array1OfReal          & Knots2,
1152  const TColStd_Array1OfInteger       & Mults2,
1153  Standard_Integer &                  NumPoles,
1154  Handle(TColStd_HArray1OfReal) &     NewKnots,
1155  Handle(TColStd_HArray1OfInteger) &  NewMults)
1156 {
1157   Standard_Integer ii,
1158   jj,
1159   continuity,
1160   set_mults_flag,
1161   degree,
1162   index,
1163   num_knots ;
1164   if (StartValue < EndValue - Tolerance) {
1165     TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1166     TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1167     degree = Degree1 + Degree2 ;
1168     index = 1 ;
1169
1170     for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1171       knots1(index) = Knots1(ii) ;
1172       index += 1 ;
1173     }
1174     index = 1 ;
1175
1176     for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1177       knots2(index) = Knots2(ii) ;
1178       index += 1 ;
1179     }
1180     BSplCLib::Reparametrize(StartValue,
1181                             EndValue,
1182                             knots1) ;
1183
1184     BSplCLib::Reparametrize(StartValue,
1185                             EndValue,
1186                             knots2) ;
1187     num_knots = 0 ;
1188     jj =  1 ;
1189
1190     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1191
1192       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1193         jj += 1 ;
1194         num_knots += 1 ;
1195       }
1196
1197       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1198         jj += 1 ;
1199       }
1200       num_knots += 1 ;
1201     }
1202     NewKnots =
1203       new TColStd_HArray1OfReal(1,num_knots) ;
1204     NewMults =
1205       new TColStd_HArray1OfInteger(1,num_knots) ;
1206     num_knots = 1 ;
1207     jj = 1 ;
1208
1209     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1210
1211       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1212         NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1213         NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1214         jj += 1 ;
1215         num_knots += 1 ;
1216       }
1217       set_mults_flag = 0 ;
1218
1219       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1220         continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1221         set_mults_flag = 1 ;
1222         NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1223         jj += 1 ;
1224       }
1225
1226       NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1227       if (! set_mults_flag) {
1228         NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1229       }
1230       num_knots += 1 ;
1231     }
1232     num_knots -= 1 ;
1233     NewMults->ChangeArray1()(1) = degree + 1 ;
1234     NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1235     index = 0 ;
1236
1237     for (ii = 1 ; ii <= num_knots ; ii++) {
1238       index += NewMults->Value(ii) ;
1239     }
1240     NumPoles = index  - degree - 1  ;
1241   }
1242 }
1243