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(1,
 350                                  ContactOrderArray(ii),
 351                                  Order,
 352                                  FlatKnots,
 353                                  Parameters(ii),
 354
 355                                  FirstNonZeroBsplineIndex,
 356                                  BSplineBasis) ;
 357     if (ErrorCode != 0) {
 358       ReturnCode = 2 ;
 359       goto FINISH ;
 360     }
 361     Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
 362
 363     for (jj = 1 ; jj < Index ; jj++) {
 364       Matrix.Value(ii,jj) = 0.0e0 ;
 365     }
 366
 367     for (jj = 1 ; jj <= Order ; jj++) {
 368       Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
 369       Index += 1 ;
 370     }
 371
 372     for (jj = Index ; jj <= BandWidth ; jj++) {
 373       Matrix.Value(ii,jj) = 0.0e0 ;
 374     }
 375   }
 376   FINISH : ;
 377   return (ReturnCode) ;
 378 }
 379
 380 //=======================================================================
 381 //function : Makes LU decompositiomn without Pivoting
 382 //purpose  : Builds the Bspline Matrix
 383 //=======================================================================
 384
 385 Standard_Integer
 386 BSplCLib::FactorBandedMatrix(math_Matrix&   Matrix,
 387                              const Standard_Integer UpperBandWidth,
 388                              const Standard_Integer LowerBandWidth,
 389                              Standard_Integer&  PivotIndexProblem)
 390 {
 391   Standard_Integer ii,
 392   jj,
 393   kk,
 394   Index,
 395   MinIndex,
 396   MaxIndex,
 397   ReturnCode = 0,
 398   BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
 399
 400   Standard_Real Inverse ;
 401   PivotIndexProblem = 0 ;
 402
 403   for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
 404     MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
 405
 406     for (jj = MinIndex ; jj <= LowerBandWidth  ; jj++) {
 407       Index = ii - LowerBandWidth + jj - 1 ;
 408       Inverse = Matrix(Index ,LowerBandWidth + 1) ;
 409       if (Abs(Inverse) > RealSmall()) {
 410         Inverse = -1.0e0/Inverse ;
 411       }
 412       else {
 413         ReturnCode = 1 ;
 414         PivotIndexProblem = Index ;
 415         goto FINISH ;
 416       }
 417       Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
 418       MaxIndex = BandWidth + Index - ii ;
 419
 420       for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
 421         Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
 422       }
 423     }
 424   }
 425   FINISH :
 426     return (ReturnCode) ;
 427 }
 428
 429 //=======================================================================
 430 //function : Build BSpline Matrix
 431 //purpose  : Builds the Bspline Matrix
 432 //=======================================================================
 433
 434 Standard_Integer
 435 BSplCLib::EvalBsplineBasis
 436 //(const Standard_Integer              Side, // = 1 rigth side, -1 left side
 437 (const Standard_Integer              , // = 1 rigth side, -1 left side
 438  const  Standard_Integer              DerivativeRequest,
 439  const  Standard_Integer              Order,
 440  const  TColStd_Array1OfReal&         FlatKnots,
 441  const  Standard_Real                 Parameter,
 442  Standard_Integer&             FirstNonZeroBsplineIndex,
 443  math_Matrix&                  BsplineBasis,
 444  Standard_Boolean              isPeriodic)
 445 {
 446   // the matrix must have at least DerivativeRequest + 1
 447   //   row and Order columns
 448   // the result are stored in the following way in
 449   // the Bspline matrix
 450   // Let i be the FirstNonZeroBsplineIndex and
 451   // t be the parameter value, k the order of the
 452   // knot vector, r the DerivativeRequest :
 453   //
 454   //   B (t)   B (t)                     B (t)
 455   //    i       i+1                       i+k-1
 456   //
 457   //    (1)     (1)                       (1)
 458   //   B (t)   B (t)                     B (t)
 459   //    i       i+1                       i+k-1
 460   //
 461   //
 462   //
 463   //
 464   //    (r)     (r)                       (r)
 465   //   B (t)   B (t)                     B (t)
 466   //    i       i+1                       i+k-1
 467   //
 468   Standard_Integer
 469     ReturnCode,
 470   ii,
 471   pp,
 472   qq,
 473   ss,
 474   NumPoles,
 475   LocalRequest ;
 476 //  ,Index ;
 477
 478   Standard_Real NewParameter,
 479   Inverse,
 480   Factor,
 481   LocalInverse,
 482   Saved ;
 483 // , *FlatKnotsArray ;
 484
 485   ReturnCode = 0 ;
 486   FirstNonZeroBsplineIndex = 0 ;
 487   LocalRequest = DerivativeRequest ;
 488   if (DerivativeRequest >= Order) {
 489     LocalRequest = Order - 1 ;
 490   }
 491
 492   if (BsplineBasis.LowerCol() != 1    ||
 493       BsplineBasis.UpperCol() < Order ||
 494       BsplineBasis.LowerRow() != 1    ||
 495       BsplineBasis.UpperRow() <= LocalRequest) {
 496     ReturnCode = 1;
 497     goto FINISH ;
 498   }
 499   NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
 500   BSplCLib::LocateParameter(Order - 1,
 501                             FlatKnots,
 502                             Parameter,
 503                             isPeriodic,
 504                             Order,
 505                             NumPoles+1,
 506                             ii,
 507                             NewParameter) ;
 508
 509   FirstNonZeroBsplineIndex = ii - Order + 1 ;
 510
 511   BsplineBasis(1,1) = 1.0e0 ;
 512   LocalRequest = DerivativeRequest ;
 513   if (DerivativeRequest >= Order) {
 514     LocalRequest = Order - 1 ;
 515   }
 516
 517   for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
 518     BsplineBasis(1,qq) = 0.0e0 ;
 519
 520     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 521       //
 522       // this should be always invertible if ii is correctly computed
 523       //
 524       Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
 525         / (FlatKnots(ii + pp)   - FlatKnots(ii - qq + pp + 1)) ;
 526       Saved = Factor *    BsplineBasis(1,pp) ;
 527       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
 528       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
 529       BsplineBasis(1,qq) = Saved ;
 530     }
 531   }
 532
 533   for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
 534
 535     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 536       BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
 537     }
 538     BsplineBasis(1,qq) = 0.0e0 ;
 539
 540     for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
 541       BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
 542     }
 543
 544     for (pp = 1 ; pp <= qq - 1 ; pp++) {
 545       Inverse = 1.0e0 / (FlatKnots(ii + pp)  - FlatKnots(ii - qq + pp + 1)) ;
 546       Factor  =  (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
 547       Saved = Factor *                 BsplineBasis(1,pp) ;
 548       BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
 549       BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
 550       BsplineBasis(1,qq) = Saved ;
 551       LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
 552
 553       for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
 554         Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
 555         BsplineBasis(Order - ss + 2, pp) *= - LocalInverse  ;
 556         BsplineBasis(Order - ss + 2, pp) +=   BsplineBasis(Order - ss + 2,qq) ;
 557         BsplineBasis(Order - ss + 2,qq) = Saved ;
 558       }
 559     }
 560   }
 561   FINISH :
 562     return (ReturnCode) ;
 563 }
 564
 565 //=======================================================================
 566 //function : MovePointAndTangent
 567 //purpose  :
 568 //=======================================================================
 569
 570 void BSplCLib::MovePointAndTangent(const Standard_Real    U,
 571                                    const Standard_Integer ArrayDimension,
 572                                    Standard_Real    &Delta,
 573                                    Standard_Real    &DeltaDerivatives,
 574                                    const Standard_Real    Tolerance,
 575                                    const Standard_Integer Degree,
 576                                    const Standard_Boolean Rational,
 577                                    const Standard_Integer StartingCondition,
 578                                    const Standard_Integer EndingCondition,
 579                                    Standard_Real&         Poles,
 580                                    const TColStd_Array1OfReal&   Weights,
 581                                    const TColStd_Array1OfReal&   FlatKnots,
 582                                    Standard_Real&        NewPoles,
 583                                    Standard_Integer&     ErrorStatus)
 584 {
 585   Standard_Integer num_poles,
 586   num_knots,
 587   ii,
 588   jj,
 589   conditions,
 590   start_num_poles,
 591   end_num_poles,
 592   index,
 593   start_index,
 594   end_index,
 595   other_index,
 596   type,
 597   order ;
 598
 599   Standard_Real    new_parameter,
 600   value,
 601   divide,
 602   end_value,
 603   start_value,
 604   *poles_array,
 605   *new_poles_array,
 606   *delta_array,
 607   *derivatives_array,
 608   *weights_array ;
 609
 610   ErrorStatus = 0 ;
 611   weights_array = NULL ;
 612   if (Rational) {
 613     weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
 614   }
 615
 616   poles_array = &Poles ;
 617   new_poles_array = &NewPoles ;
 618   delta_array = &Delta ;
 619   derivatives_array = &DeltaDerivatives ;
 620   order = Degree + 1 ;
 621   num_knots = FlatKnots.Length() ;
 622   num_poles = num_knots - order ;
 623   conditions = StartingCondition + EndingCondition + 4 ;
 624   //
 625   // check validity of input data
 626   //
 627   if (StartingCondition >= -1 &&
 628       StartingCondition <= Degree &&
 629       EndingCondition >= -1 &&
 630       EndingCondition <= Degree &&
 631       conditions <= num_poles) {
 632     //
 633     // check the parameter is within bounds
 634     //
 635     start_index = FlatKnots.Lower() + Degree ;
 636     end_index = FlatKnots.Upper() - Degree ;
 637     //
 638     //  check if there is enough room to move the poles
 639     //
 640     conditions = 1 ;
 641     if (StartingCondition == -1) {
 642       conditions = conditions && (FlatKnots(start_index) <=  U) ;
 643     }
 644     else {
 645       conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
 646     }
 647     if (EndingCondition == -1) {
 648       conditions = conditions && (FlatKnots(end_index) >= U) ;
 649     }
 650     else {
 651       conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ;
 652     }
 653
 654     if (conditions) {
 655       //
 656       // build 2 auxialiary functions
 657       //
 658       TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
 659       TColStd_Array1OfReal first_function  (1,num_poles) ;
 660       TColStd_Array1OfReal second_function (1,num_poles) ;
 661
 662       BuildSchoenbergPoints(Degree,
 663                             FlatKnots,
 664                             schoenberg_points) ;
 665       start_index = StartingCondition + 2 ;
 666       end_index = num_poles - EndingCondition - 1 ;
 667       LocateParameter(schoenberg_points,
 668                       U,
 669                       Standard_False,
 670                       start_index,
 671                       end_index,
 672                       index,
 673                       new_parameter,
 674                       0, 1) ;
 675
 676       if (index == start_index) {
 677         other_index = index + 1 ;
 678       }
 679       else if (index == end_index) {
 680         other_index = index -1  ;
 681       }
 682       else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
 683         other_index = index - 1 ;
 684       }
 685       else {
 686         other_index = index + 1 ;
 687       }
 688       type = 3 ;
 689
 690       start_num_poles = StartingCondition + 2 ;
 691       end_num_poles   = num_poles - EndingCondition - 1 ;
 692       if (start_num_poles == 1) {
 693         start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
 694         start_value = schoenberg_points(1) - start_value ;
 695       }
 696       else {
 697         start_value = schoenberg_points(start_num_poles - 1) ;
 698       }
 699       if (end_num_poles == num_poles) {
 700         end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
 701         end_value = schoenberg_points(num_poles) + end_value ;
 702       }
 703       else {
 704         end_value = schoenberg_points(end_num_poles + 1) ;
 705       }
 706
 707       for (ii = 1 ; ii < start_num_poles ; ii++) {
 708         first_function(ii) = 0.e0 ;
 709         second_function(ii) = 0.0e0 ;
 710       }
 711
 712       for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
 713         first_function(ii) = 0.e0 ;
 714         second_function(ii) = 0.0e0 ;
 715       }
 716       divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
 717
 718       for (ii = start_num_poles  ; ii <= index ; ii++) {
 719         value = schoenberg_points(ii) - start_value ;
 720         value *= divide ;
 721         first_function(ii) = 1.0e0 ;
 722
 723         for (jj = 0 ; jj < type ; jj++) {
 724           first_function(ii) *= value ;
 725         }
 726       }
 727       divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
 728
 729       for (ii = index ; ii <= end_num_poles ; ii++) {
 730         value = end_value - schoenberg_points(ii) ;
 731         value *= divide ;
 732         first_function(ii) = 1.0e0 ;
 733
 734         for (jj = 0 ; jj < type ; jj++) {
 735           first_function(ii) *= value ;
 736         }
 737       }
 738
 739       divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
 740
 741       for (ii = start_num_poles  ; ii <= other_index ; ii++) {
 742         value = schoenberg_points(ii) - start_value ;
 743         value *= divide ;
 744         second_function(ii) = 1.0e0 ;
 745
 746         for (jj = 0 ; jj < type ; jj++) {
 747           second_function(ii) *= value ;
 748         }
 749       }
 750       divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
 751
 752       for (ii = other_index ; ii <= end_num_poles ; ii++) {
 753         value = end_value - schoenberg_points(ii) ;
 754         value *= divide ;
 755         second_function(ii) = 1.0e0 ;
 756
 757         for (jj = 0 ; jj < type ; jj++) {
 758           second_function(ii) *= value ;
 759         }
 760       }
 761
 762       //
 763       //  compute the point and derivatives of both functions
 764       //
 765       Standard_Real results[2][2],
 766       weights_results[2][2];
 767       Standard_Integer extrap_mode[2],
 768       derivative_request = 1,
 769       dimension = 1 ;
 770       Standard_Boolean periodic_flag = Standard_False ;
 771
 772       extrap_mode[0] = Degree ;
 773       extrap_mode[1] = Degree ;
 774       if (Rational) {
 775         //
 776         // evaluate in homogenised form
 777         //
 778         Eval(U,
 779              periodic_flag,
 780              derivative_request,
 781              extrap_mode[0],
 782              Degree,
 783              FlatKnots,
 784              dimension,
 785              first_function(1),
 786              weights_array[0],
 787              results[0][0],
 788              weights_results[0][0]) ;
 789
 790         Eval(U,
 791              periodic_flag,
 792              derivative_request,
 793              extrap_mode[0],
 794              Degree,
 795              FlatKnots,
 796              dimension,
 797              second_function(1),
 798              weights_array[0],
 799              results[1][0],
 800              weights_results[1][0]) ;
 801         //
 802         //  compute the rational derivatives values
 803         //
 804
 805         for (ii = 0 ; ii < 2 ; ii++) {
 806           PLib::RationalDerivatives(1,
 807                                     1,
 808                                     results[ii][0],
 809                                     weights_results[ii][0],
 810                                     results[ii][0]) ;
 811         }
 812       }
 813       else {
 814         Eval(U,
 815              Standard_False,
 816              1,
 817              extrap_mode[0],
 818              Degree,
 819              FlatKnots,
 820              1,
 821              first_function(1),
 822              results[0][0]) ;
 823
 824         Eval(U,
 825              Standard_False,
 826              1,
 827              extrap_mode[0],
 828              Degree,
 829              FlatKnots,
 830              1,
 831              second_function(1),
 832              results[1][0]) ;
 833       }
 834       gp_Mat2d  a_matrix ;
 835
 836       for (ii = 0 ; ii < 2 ; ii++) {
 837
 838         for (jj = 0 ; jj < 2 ; jj++) {
 839           a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
 840         }
 841       }
 842       a_matrix.Invert() ;
 843       TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
 844       TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
 845
 846       for( ii = 0 ; ii < ArrayDimension ; ii++) {
 847         the_a_vector(ii) =
 848           a_matrix.Value(1,1) * delta_array[ii] +
 849             a_matrix.Value(2,1) * derivatives_array[ii] ;
 850         the_b_vector(ii) =
 851           a_matrix.Value(1,2) * delta_array[ii] +
 852             a_matrix.Value(2,2) * derivatives_array[ii] ;
 853       }
 854       index = 0 ;
 855
 856       for (ii = 0 ; ii < num_poles ; ii++) {
 857
 858         for (jj = 0 ; jj < ArrayDimension ; jj++) {
 859           new_poles_array[index] = poles_array[index] ;
 860           new_poles_array[index] +=
 861             first_function(ii+1) * the_a_vector(jj) ;
 862           new_poles_array[index] +=
 863             second_function(ii+1) * the_b_vector(jj) ;
 864           index += 1 ;
 865         }
 866       }
 867     }
 868     else {
 869       ErrorStatus = 1 ;
 870     }
 871   }
 872   else {
 873     ErrorStatus = 2 ;
 874   }
 875 }
 876
 877 //=======================================================================
 878 //function : FunctionMultiply
 879 //purpose  :
 880 //=======================================================================
 881
 882 void BSplCLib::FunctionMultiply
 883 (const BSplCLib_EvaluatorFunction & FunctionPtr,
 884  const Standard_Integer             BSplineDegree,
 885  const TColStd_Array1OfReal &       BSplineFlatKnots,
 886  const Standard_Integer             PolesDimension,
 887  Standard_Real &                    Poles,
 888  const TColStd_Array1OfReal &       FlatKnots,
 889  const Standard_Integer             NewDegree,
 890  Standard_Real &                    NewPoles,
 891  Standard_Integer &                 Status)
 892 {
 893   Standard_Integer ii,
 894   jj,
 895   index ;
 896   Standard_Integer extrap_mode[2],
 897   error_code,
 898   num_new_poles,
 899   derivative_request = 0 ;
 900   Standard_Boolean  periodic_flag = Standard_False ;
 901   Standard_Real  result,
 902   start_end[2],
 903   *array_of_poles,
 904   *array_of_new_poles ;
 905
 906   array_of_poles = (Standard_Real *) &NewPoles ;
 907   extrap_mode[0] =
 908     extrap_mode[1] = BSplineDegree ;
 909   num_new_poles =
 910     FlatKnots.Length() - NewDegree - 1 ;
 911   start_end[0] = FlatKnots(NewDegree+1) ;
 912   start_end[1] = FlatKnots(num_new_poles+1) ;
 913   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
 914   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
 915   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
 916
 917   array_of_new_poles =
 918     (Standard_Real *) &new_poles_array(1) ;
 919   BuildSchoenbergPoints(NewDegree,
 920                         FlatKnots,
 921                         parameters) ;
 922   //
 923   // on recadre sur les bornes
 924   //
 925   if (parameters(1) < start_end[0]) {
 926     parameters(1) = start_end[0] ;
 927   }
 928   if (parameters(num_new_poles) > start_end[1]) {
 929     parameters(num_new_poles) = start_end[1] ;
 930   }
 931   index = 0 ;
 932
 933   for (ii = 1 ; ii <= num_new_poles ; ii++) {
 934     contact_order_array(ii) = 0 ;
 935     FunctionPtr.Evaluate (contact_order_array(ii),
 936                    start_end,
 937                    parameters(ii),
 938                    result,
 939                    error_code);
 940     if (error_code) {
 941       Status = 1 ;
 942       goto FINISH ;
 943     }
 944
 945     Eval(parameters(ii),
 946          periodic_flag,
 947          derivative_request,
 948          extrap_mode[0],
 949          BSplineDegree,
 950          BSplineFlatKnots,
 951          PolesDimension,
 952          Poles,
 953          array_of_new_poles[index]) ;
 954
 955     for (jj = 0 ; jj < PolesDimension ; jj++) {
 956       array_of_new_poles[index] *= result ;
 957       index += 1 ;
 958     }
 959   }
 960   Interpolate(NewDegree,
 961               FlatKnots,
 962               parameters,
 963               contact_order_array,
 964               PolesDimension,
 965               array_of_new_poles[0],
 966               Status) ;
 967
 968   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
 969     array_of_poles[ii] = array_of_new_poles[ii] ;
 970
 971   }
 972   FINISH :
 973     ;
 974 }
 975
 976 //=======================================================================
 977 // function : FunctionMultiply
 978 //purpose  :
 979 //=======================================================================
 980
 981 void BSplCLib::FunctionReparameterise
 982 (const BSplCLib_EvaluatorFunction & FunctionPtr,
 983  const Standard_Integer             BSplineDegree,
 984  const TColStd_Array1OfReal &       BSplineFlatKnots,
 985  const Standard_Integer             PolesDimension,
 986  Standard_Real &                    Poles,
 987  const TColStd_Array1OfReal &       FlatKnots,
 988  const Standard_Integer             NewDegree,
 989  Standard_Real &                    NewPoles,
 990  Standard_Integer &                 Status)
 991 {
 992   Standard_Integer ii,
 993 //  jj,
 994   index ;
 995   Standard_Integer extrap_mode[2],
 996   error_code,
 997   num_new_poles,
 998   derivative_request = 0 ;
 999   Standard_Boolean  periodic_flag = Standard_False ;
1000   Standard_Real  result,
1001   start_end[2],
1002   *array_of_poles,
1003   *array_of_new_poles ;
1004
1005   array_of_poles = (Standard_Real *) &NewPoles ;
1006   extrap_mode[0] =
1007     extrap_mode[1] = BSplineDegree ;
1008   num_new_poles =
1009     FlatKnots.Length() - NewDegree - 1 ;
1010   start_end[0] = FlatKnots(NewDegree+1) ;
1011   start_end[1] = FlatKnots(num_new_poles+1) ;
1012   TColStd_Array1OfReal  parameters(1,num_new_poles) ;
1013   TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1014   TColStd_Array1OfReal  new_poles_array(1, num_new_poles * PolesDimension) ;
1015
1016   array_of_new_poles =
1017     (Standard_Real *) &new_poles_array(1) ;
1018   BuildSchoenbergPoints(NewDegree,
1019                         FlatKnots,
1020                         parameters) ;
1021   index = 0 ;
1022
1023   for (ii = 1 ; ii <= num_new_poles ; ii++) {
1024     contact_order_array(ii) = 0 ;
1025     FunctionPtr.Evaluate (contact_order_array(ii),
1026                    start_end,
1027                    parameters(ii),
1028                    result,
1029                    error_code);
1030     if (error_code) {
1031       Status = 1 ;
1032       goto FINISH ;
1033     }
1034
1035     Eval(result,
1036          periodic_flag,
1037          derivative_request,
1038          extrap_mode[0],
1039          BSplineDegree,
1040          BSplineFlatKnots,
1041          PolesDimension,
1042          Poles,
1043          array_of_new_poles[index]) ;
1044     index += PolesDimension ;
1045   }
1046   Interpolate(NewDegree,
1047               FlatKnots,
1048               parameters,
1049               contact_order_array,
1050               PolesDimension,
1051               array_of_new_poles[0],
1052               Status) ;
1053
1054   for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1055     array_of_poles[ii] = array_of_new_poles[ii] ;
1056
1057   }
1058   FINISH :
1059     ;
1060 }
1061
1062 //=======================================================================
1063 //function : FunctionMultiply
1064 //purpose  :
1065 //=======================================================================
1066
1067 void BSplCLib::FunctionMultiply
1068 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1069  const Standard_Integer              BSplineDegree,
1070  const TColStd_Array1OfReal &        BSplineFlatKnots,
1071  const TColStd_Array1OfReal &        Poles,
1072  const TColStd_Array1OfReal &        FlatKnots,
1073  const Standard_Integer              NewDegree,
1074  TColStd_Array1OfReal &              NewPoles,
1075  Standard_Integer &                  Status)
1076 {
1077   Standard_Integer num_bspline_poles =
1078     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1079   Standard_Integer num_new_poles =
1080     FlatKnots.Length() - NewDegree - 1 ;
1081
1082   if (Poles.Length() != num_bspline_poles ||
1083       NewPoles.Length() != num_new_poles) {
1084     Standard_ConstructionError::Raise();
1085   }
1086   Standard_Real  * array_of_poles =
1087     (Standard_Real *) &Poles(Poles.Lower()) ;
1088   Standard_Real  * array_of_new_poles =
1089     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1090   BSplCLib::FunctionMultiply(FunctionPtr,
1091                              BSplineDegree,
1092                              BSplineFlatKnots,
1093                              1,
1094                              array_of_poles[0],
1095                              FlatKnots,
1096                              NewDegree,
1097                              array_of_new_poles[0],
1098                              Status) ;
1099 }
1100
1101 //=======================================================================
1102 //function : FunctionReparameterise
1103 //purpose  :
1104 //=======================================================================
1105
1106 void BSplCLib::FunctionReparameterise
1107 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1108  const Standard_Integer              BSplineDegree,
1109  const TColStd_Array1OfReal &        BSplineFlatKnots,
1110  const TColStd_Array1OfReal &        Poles,
1111  const TColStd_Array1OfReal &        FlatKnots,
1112  const Standard_Integer              NewDegree,
1113  TColStd_Array1OfReal &              NewPoles,
1114  Standard_Integer &                  Status)
1115 {
1116   Standard_Integer num_bspline_poles =
1117     BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1118   Standard_Integer num_new_poles =
1119     FlatKnots.Length() - NewDegree - 1 ;
1120
1121   if (Poles.Length() != num_bspline_poles ||
1122       NewPoles.Length() != num_new_poles) {
1123     Standard_ConstructionError::Raise();
1124   }
1125   Standard_Real  * array_of_poles =
1126     (Standard_Real *) &Poles(Poles.Lower()) ;
1127   Standard_Real  * array_of_new_poles =
1128     (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1129   BSplCLib::FunctionReparameterise(
1130                                    FunctionPtr,
1131                                    BSplineDegree,
1132                                    BSplineFlatKnots,
1133                                    1,
1134                                    array_of_poles[0],
1135                                    FlatKnots,
1136                                    NewDegree,
1137                                    array_of_new_poles[0],
1138                                    Status) ;
1139 }
1140
1141 //=======================================================================
1142 //function : MergeBSplineKnots
1143 //purpose  :
1144 //=======================================================================
1145 void BSplCLib::MergeBSplineKnots
1146 (const Standard_Real                 Tolerance,
1147  const Standard_Real                 StartValue,
1148  const Standard_Real                 EndValue,
1149  const Standard_Integer              Degree1,
1150  const TColStd_Array1OfReal          & Knots1,
1151  const TColStd_Array1OfInteger       & Mults1,
1152  const Standard_Integer              Degree2,
1153  const TColStd_Array1OfReal          & Knots2,
1154  const TColStd_Array1OfInteger       & Mults2,
1155  Standard_Integer &                  NumPoles,
1156  Handle(TColStd_HArray1OfReal) &     NewKnots,
1157  Handle(TColStd_HArray1OfInteger) &  NewMults)
1158 {
1159   Standard_Integer ii,
1160   jj,
1161   continuity,
1162   set_mults_flag,
1163   degree,
1164   index,
1165   num_knots ;
1166   if (StartValue < EndValue - Tolerance) {
1167     TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1168     TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1169     degree = Degree1 + Degree2 ;
1170     index = 1 ;
1171
1172     for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1173       knots1(index) = Knots1(ii) ;
1174       index += 1 ;
1175     }
1176     index = 1 ;
1177
1178     for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1179       knots2(index) = Knots2(ii) ;
1180       index += 1 ;
1181     }
1182     BSplCLib::Reparametrize(StartValue,
1183                             EndValue,
1184                             knots1) ;
1185
1186     BSplCLib::Reparametrize(StartValue,
1187                             EndValue,
1188                             knots2) ;
1189     num_knots = 0 ;
1190     jj =  1 ;
1191
1192     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1193
1194       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1195         jj += 1 ;
1196         num_knots += 1 ;
1197       }
1198
1199       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1200         jj += 1 ;
1201       }
1202       num_knots += 1 ;
1203     }
1204     NewKnots =
1205       new TColStd_HArray1OfReal(1,num_knots) ;
1206     NewMults =
1207       new TColStd_HArray1OfInteger(1,num_knots) ;
1208     num_knots = 1 ;
1209     jj = 1 ;
1210
1211     for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1212
1213       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1214         NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1215         NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1216         jj += 1 ;
1217         num_knots += 1 ;
1218       }
1219       set_mults_flag = 0 ;
1220
1221       while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1222         continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1223         set_mults_flag = 1 ;
1224         NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1225         jj += 1 ;
1226       }
1227
1228       NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1229       if (! set_mults_flag) {
1230         NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1231       }
1232       num_knots += 1 ;
1233     }
1234     num_knots -= 1 ;
1235     NewMults->ChangeArray1()(1) = degree + 1 ;
1236     NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1237     index = 0 ;
1238
1239     for (ii = 1 ; ii <= num_knots ; ii++) {
1240       index += NewMults->Value(ii) ;
1241     }
1242     NumPoles = index  - degree - 1  ;
1243   }
1244 }
1245