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