src/BSplSLib/BSplSLib.cxx

   1 // Created on: 1991-08-26
   2 // Created by: JCV
   3 // Copyright (c) 1991-1999 Matra Datavision
   4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 // Modifed RLE Aug 93 - Complete rewrite
  23 // xab  21-Mar-95  implemented cache mecanism
  24 // pmn  25-09-96   Interpolation
  25 // jct  25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
  26 // pmn  07-10-96 : Correction de DN dans le cas rationnal.
  27 // pmn  06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
  28
  29 #include <BSplSLib.ixx>
  30 #include <PLib.hxx>
  31 #include <NCollection_LocalArray.hxx>
  32 #include <BSplCLib.hxx>
  33 #include <TColgp_Array2OfXYZ.hxx>
  34 #include <TColgp_Array1OfXYZ.hxx>
  35 #include <TColStd_HArray1OfInteger.hxx>
  36 #include <Standard_NotImplemented.hxx>
  37 #include <Standard_ConstructionError.hxx>
  38 #include <math_Matrix.hxx>
  39
  40 // for null derivatives
  41 static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
  42 #ifdef WNT
  43 #define M_SQRT2 1.4142135623730950488016887
  44 #endif
  45
  46 //=======================================================================
  47 //struct : BSplCLib_DataContainer
  48 //purpose: Auxiliary structure providing buffers for poles and knots used in
  49 //         evaluation of bspline (allocated in the stack)
  50 //=======================================================================
  51
  52 struct BSplSLib_DataContainer
  53 {
  54   BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
  55   {
  56     (void)UDegree; (void)VDegree; // just to avoid compiler warning in Release mode
  57     Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
  58         VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
  59         "BSplSLib: bspline degree is greater than maximum supported");
  60   }
  61   Standard_Real poles[4*(25+1)*(25+1)];
  62   Standard_Real knots1[2*25];
  63   Standard_Real knots2[2*25];
  64   Standard_Real ders[48];
  65 };
  66
  67 //**************************************************************************
  68 //                     Evaluation methods
  69 //**************************************************************************
  70
  71 //=======================================================================
  72 //function : RationalDerivative
  73 //purpose  : computes the rational derivatives when whe have the
  74 //           the derivatives of the homogeneous numerator and the
  75 //           the derivatives of the denominator
  76 //=======================================================================
  77
  78 void  BSplSLib::RationalDerivative(const Standard_Integer UDeg,
  79                                    const Standard_Integer VDeg,
  80                                    const Standard_Integer N,
  81                                    const Standard_Integer M,
  82                                    Standard_Real& HDerivatives,
  83                                    Standard_Real& RDerivatives,
  84                                    const Standard_Boolean All)
  85 {
  86   //
  87   //  if All is True all derivatives are computed. if Not only
  88   //  the requested N, M is computed
  89   //
  90   //                                           Numerator(u,v)
  91   //   let f(u,v) be a rational function  = ------------------
  92   //                                         Denominator(u,v)
  93   //
  94   //
  95   //   Let (N,M) the order of the derivatives we want : then since
  96   //   we have :
  97   //
  98   //         Numerator = f * Denominator
  99   //
 100   //   we derive :
 101   //
 102   //   (N,M)         1           (         (N M)                  (p q)            (N -p M-q) )
 103   //  f       =  ------------    (  Numerator   -  SUM SUM  a   * f    * Denominator          )
 104   //                       (0,0) (                 p<N q<M   p q                              )
 105   //             Denominator
 106   //
 107   //   with :
 108   //
 109   //              ( N )  ( M )
 110   //    a      =  (   )  (   )
 111   //     p q      ( p )  ( q )
 112   //
 113   //
 114   //   HDerivatives is an array where derivatives are stored in the following form
 115   //   Numerator is assumee to have 3 functions that is a vector of dimension
 116   //   3
 117   //
 118   //             (0,0)           (0,0)                       (0, DegV)           (0, DegV)
 119   //     Numerator     Denominator      ...         Numerator          Denominator
 120   //
 121   //             (1,0)           (1,0)                       (1, DegV)           (1, DegV)
 122   //     Numerator     Denominator      ...         Numerator          Denominator
 123   //
 124   //         ...........................................................
 125   //
 126   //
 127   //            (DegU,0)        (DegU,0)                  (DegU, DegV)           (DegU, DegV)
 128   //     Numerator     Denominator      ...         Numerator          Denominator
 129   //
 130   //
 131   Standard_Integer ii,jj,pp,qq,index,index1,index2;
 132   Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
 133   Standard_Integer MinN,MinN1,MinM,MinM1;
 134   Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
 135
 136   M1 = M + 1;
 137   N1 = N + 1;
 138   ii = N1 * M1;
 139   M3 = (M1 << 1) + M1;
 140   M4 = (VDeg + 1) << 2;
 141
 142   NCollection_LocalArray<Standard_Real> StoreDerivatives (All ? 0 : ii * 3);
 143   Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
 144   NCollection_LocalArray<Standard_Real> StoreW (ii);
 145   Standard_Real *HomogeneousArray = &HDerivatives;
 146   Standard_Real denominator,Pii,Pip,Pjq;
 147
 148   denominator = 1.0e0 / HomogeneousArray[3];
 149   index_u  = 0;
 150   index_u1 = 0;
 151   if (UDeg < N) MinN = UDeg;
 152   else          MinN = N;
 153   if (VDeg < M) MinM = VDeg;
 154   else          MinM = M;
 155   MinN1 = MinN + 1;
 156   MinM1 = MinM + 1;
 157   iiM1 = - M1;
 158
 159   for (ii = 0 ; ii < MinN1 ; ii++) {
 160     iiM1    += M1;
 161     index_v  = index_u;
 162     index_v1 = index_u1;
 163     index_w  = iiM1;
 164
 165     for (jj = 0 ; jj < MinM1 ; jj++) {
 166       RArray[index_v++] = HomogeneousArray[index_v1++];
 167       RArray[index_v++] = HomogeneousArray[index_v1++];
 168       RArray[index_v++] = HomogeneousArray[index_v1++];
 169       StoreW[index_w++] = HomogeneousArray[index_v1++];
 170     }
 171
 172     for (jj = MinM1 ; jj < M1 ; jj++) {
 173       RArray[index_v++] = 0.;
 174       RArray[index_v++] = 0.;
 175       RArray[index_v++] = 0.;
 176       StoreW[index_w++] = 0.;
 177     }
 178     index_u1 += M4;
 179     index_u  += M3;
 180   }
 181   index_v = MinN1 * M3;
 182   index_w = MinN1 * M1;
 183
 184   for (ii = MinN1 ; ii < N1 ; ii++) {
 185
 186     for (jj = 0 ; jj < M1 ; jj++) {
 187       RArray[index_v++] = 0.0e0;
 188       RArray[index_v++] = 0.0e0;
 189       RArray[index_v++] = 0.0e0;
 190       StoreW[index_w++] = 0.0e0;
 191     }
 192   }
 193
 194   // ---------------  Calculation ----------------
 195
 196   iiM1 = - M1;
 197   iiM3 = - M3;
 198
 199   for (ii = 0 ; ii <= N  ; ii++) {
 200     iiM1  += M1;
 201     iiM3  += M3;
 202     index1 = iiM3 - 3;
 203     jjM1   = iiM1;
 204
 205     for (jj = 0 ; jj <= M ; jj++) {
 206       jjM1 ++;
 207       ppM1    = - M1;
 208       ppM3    = - M3;
 209       index1 += 3;
 210
 211       for (pp = 0 ; pp < ii ; pp++) {
 212         ppM1  += M1;
 213         ppM3  += M3;
 214         index  = ppM3;
 215         index2 = jjM1 - ppM1;
 216         Pip    = PLib::Bin(ii,pp);
 217
 218         for (qq = 0 ; qq <= jj ; qq++) {
 219           index2--;
 220           Pjq    = Pip * PLib::Bin(jj,qq) * StoreW[index2];
 221           RArray[index1] -= Pjq * RArray[index]; index++; index1++;
 222           RArray[index1] -= Pjq * RArray[index]; index++; index1++;
 223           RArray[index1] -= Pjq * RArray[index]; index++;
 224           index1 -= 2;
 225         }
 226       }
 227       index  = iiM3;
 228       index2 = jj + 1;
 229       Pii = PLib::Bin(ii,ii);
 230
 231       for (qq = 0 ; qq < jj ; qq++) {
 232         index2--;
 233         Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
 234         RArray[index1] -= Pjq * RArray[index]; index++; index1++;
 235         RArray[index1] -= Pjq * RArray[index]; index++; index1++;
 236         RArray[index1] -= Pjq * RArray[index]; index++;
 237         index1 -= 2;
 238       }
 239       RArray[index1] *= denominator; index1++;
 240       RArray[index1] *= denominator; index1++;
 241       RArray[index1] *= denominator;
 242       index1 -= 2;
 243     }
 244   }
 245   if (!All) {
 246     RArray = &RDerivatives;
 247     index = N * M1 + M;
 248     index = (index << 1) + index;
 249     RArray[0] = StoreDerivatives[index]; index++;
 250     RArray[1] = StoreDerivatives[index]; index++;
 251     RArray[2] = StoreDerivatives[index];
 252   }
 253 }
 254
 255 //=======================================================================
 256 //function : PrepareEval
 257 //purpose  :
 258 //=======================================================================
 259
 260 //
 261 // PrepareEval :
 262 //
 263 // Prepare all data for computing points :
 264 //  local arrays of knots
 265 //  local array  of poles (multiplied by the weights if rational)
 266 //
 267 //  The first direction to compute (smaller degree) is returned
 268 //  and the poles are stored according to this direction.
 269
 270 static Standard_Boolean  PrepareEval (const Standard_Real            U,
 271                                       const Standard_Real            V,
 272                                       const Standard_Integer         Uindex,
 273                                       const Standard_Integer         Vindex,
 274                                       const Standard_Integer         UDegree,
 275                                       const Standard_Integer         VDegree,
 276                                       const Standard_Boolean         URat,
 277                                       const Standard_Boolean         VRat,
 278                                       const Standard_Boolean         UPer,
 279                                       const Standard_Boolean         VPer,
 280                                       const TColgp_Array2OfPnt&      Poles,
 281                                       const TColStd_Array2OfReal&    Weights,
 282                                       const TColStd_Array1OfReal&    UKnots,
 283                                       const TColStd_Array1OfReal&    VKnots,
 284                                       const TColStd_Array1OfInteger& UMults,
 285                                       const TColStd_Array1OfInteger& VMults,
 286                                       Standard_Real& u1,     // first  parameter to use
 287                                       Standard_Real& u2,     // second parameter to use
 288                                       Standard_Integer& d1,  // first degree
 289                                       Standard_Integer& d2,  // second degree
 290                                       Standard_Boolean& rational,
 291                                       BSplSLib_DataContainer& dc)
 292   {
 293   rational = URat || VRat;
 294   Standard_Integer uindex = Uindex;
 295   Standard_Integer vindex = Vindex;
 296   Standard_Integer UKLower = UKnots.Lower();
 297   Standard_Integer UKUpper = UKnots.Upper();
 298   Standard_Integer VKLower = VKnots.Lower();
 299   Standard_Integer VKUpper = VKnots.Upper();
 300
 301   if (UDegree <= VDegree)
 302     {
 303     // compute the indices
 304     if (uindex < UKLower || uindex > UKUpper)
 305       BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
 306     else
 307       u1 = U;
 308
 309     if (vindex < VKLower || vindex > VKUpper)
 310       BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
 311     else
 312       u2 = V;
 313
 314     // get the knots
 315     d1 = UDegree;
 316     d2 = VDegree;
 317     BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
 318     BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
 319
 320     if (&UMults == NULL)
 321       uindex -= UKLower + UDegree;
 322     else
 323       uindex  = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
 324
 325     if (&VMults == NULL)
 326       vindex -= VKLower + VDegree;
 327     else
 328       vindex  = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
 329
 330     // get the poles
 331     Standard_Integer i,j,ip,jp;
 332     Standard_Real w, *pole = dc.poles;
 333     d1 = UDegree;
 334     d2 = VDegree;
 335     Standard_Integer PLowerRow = Poles.LowerRow();
 336     Standard_Integer PUpperRow = Poles.UpperRow();
 337     Standard_Integer PLowerCol = Poles.LowerCol();
 338     Standard_Integer PUpperCol = Poles.UpperCol();
 339
 340     // verify if locally non rational
 341     if (rational)
 342       {
 343       rational = Standard_False;
 344       ip = PLowerRow + uindex;
 345       jp = PLowerCol + vindex;
 346
 347       if(ip < PLowerRow)        ip = PUpperRow;
 348       if(jp < PLowerCol)        jp = PUpperCol;
 349
 350       w  = Weights.Value(ip,jp);
 351       Standard_Real eps = Epsilon(w);
 352       Standard_Real dw;
 353
 354       for (i = 0; i <= UDegree && !rational; i++)
 355         {
 356         jp = PLowerCol + vindex;
 357
 358         if(jp < PLowerCol)
 359           jp = PUpperCol;
 360
 361         for (j = 0; j <= VDegree && !rational; j++)
 362           {
 363           dw = Weights.Value(ip,jp) - w;
 364           if (dw < 0)
 365             dw = - dw;
 366
 367           rational = (dw > eps);
 368
 369           jp++;
 370
 371           if (jp > PUpperCol)
 372             jp = PLowerCol;
 373           }
 374
 375         ip++;
 376
 377         if (ip > PUpperRow)
 378           ip = PLowerRow;
 379
 380         }
 381       }
 382
 383     // copy the poles
 384     ip = PLowerRow + uindex;
 385
 386     if(ip < PLowerRow)
 387       ip = PUpperRow;
 388
 389     if (rational)
 390       {
 391       for (i = 0; i <= d1; i++)
 392         {
 393         jp = PLowerCol + vindex;
 394
 395         if(jp < PLowerCol)
 396           jp = PUpperCol;
 397
 398         for (j = 0; j <= d2; j++)
 399           {
 400           const gp_Pnt& P = Poles  .Value(ip,jp);
 401           pole[3] = w     = Weights.Value(ip,jp);
 402           pole[0] = P.X() * w;
 403           pole[1] = P.Y() * w;
 404           pole[2] = P.Z() * w;
 405           pole   += 4;
 406           jp++;
 407
 408           if (jp > PUpperCol)
 409             jp = PLowerCol;
 410           }
 411
 412         ip++;
 413
 414         if (ip > PUpperRow)
 415           ip = PLowerRow;
 416
 417         }
 418       }
 419     else
 420       {
 421       for (i = 0; i <= d1; i++)
 422         {
 423         jp = PLowerCol + vindex;
 424
 425         if(jp < PLowerCol)
 426           jp = PUpperCol;
 427
 428         for (j = 0; j <= d2; j++)
 429           {
 430           const gp_Pnt& P = Poles.Value(ip,jp);
 431           pole[0] = P.X();
 432           pole[1] = P.Y();
 433           pole[2] = P.Z();
 434           pole   += 3;
 435           jp++;
 436
 437           if (jp > PUpperCol)
 438             jp = PLowerCol;
 439           }
 440
 441         ip++;
 442
 443         if (ip > PUpperRow)
 444           ip = PLowerRow;
 445         }
 446       }
 447
 448     return Standard_True;
 449     }
 450   else
 451     {
 452     // compute the indices
 453     if (uindex < UKLower || uindex > UKUpper)
 454       BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
 455     else
 456       u2 = U;
 457
 458     if (vindex < VKLower || vindex > VKUpper)
 459       BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
 460     else
 461       u1 = V;
 462
 463     // get the knots
 464
 465     d2 = UDegree;
 466     d1 = VDegree;
 467
 468     BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
 469     BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
 470
 471     if (&UMults == NULL)
 472       uindex -= UKLower + UDegree;
 473     else
 474       uindex  = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults);
 475
 476     if (&VMults == NULL)
 477       vindex -= VKLower + VDegree;
 478     else
 479       vindex  = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults);
 480
 481     // get the poles
 482     Standard_Integer i,j,ip,jp;
 483     Standard_Real w, *pole = dc.poles;
 484     d1 = VDegree;
 485     d2 = UDegree;
 486     Standard_Integer PLowerRow = Poles.LowerRow();
 487     Standard_Integer PUpperRow = Poles.UpperRow();
 488     Standard_Integer PLowerCol = Poles.LowerCol();
 489     Standard_Integer PUpperCol = Poles.UpperCol();
 490
 491     // verify if locally non rational
 492     if (rational)
 493       {
 494       rational = Standard_False;
 495       ip = PLowerRow + uindex;
 496       jp = PLowerCol + vindex;
 497
 498       if(ip < PLowerRow)
 499         ip = PUpperRow;
 500
 501       if(jp < PLowerCol)
 502         jp = PUpperCol;
 503
 504       w  = Weights.Value(ip,jp);
 505       Standard_Real eps = Epsilon(w);
 506       Standard_Real dw;
 507
 508       for (i = 0; i <= UDegree && !rational; i++)
 509         {
 510         jp = PLowerCol + vindex;
 511
 512         if(jp < PLowerCol)
 513           jp = PUpperCol;
 514
 515         for (j = 0; j <= VDegree && !rational; j++)
 516           {
 517           dw = Weights.Value(ip,jp) - w;
 518           if (dw < 0) dw = - dw;
 519           rational = dw > eps;
 520
 521           jp++;
 522
 523           if (jp > PUpperCol)
 524             jp = PLowerCol;
 525           }
 526
 527         ip++;
 528
 529         if (ip > PUpperRow)
 530           ip = PLowerRow;
 531
 532         }
 533       }
 534
 535     // copy the poles
 536     jp = PLowerCol + vindex;
 537
 538     if(jp < PLowerCol)
 539       jp = PUpperCol;
 540
 541     if (rational)
 542       {
 543       for (i = 0; i <= d1; i++)
 544         {
 545         ip = PLowerRow + uindex;
 546
 547         if(ip < PLowerRow)
 548           ip = PUpperRow;
 549
 550         for (j = 0; j <= d2; j++)
 551           {
 552           const gp_Pnt& P = Poles.Value(ip,jp);
 553           pole[3] = w     = Weights.Value(ip,jp);
 554           pole[0] = P.X() * w;
 555           pole[1] = P.Y() * w;
 556           pole[2] = P.Z() * w;
 557           pole   += 4;
 558           ip++;
 559
 560           if (ip > PUpperRow)
 561             ip = PLowerRow;
 562
 563           }
 564
 565         jp++;
 566
 567         if (jp > PUpperCol)
 568           jp = PLowerCol;
 569
 570         }
 571       }
 572     else
 573       {
 574       for (i = 0; i <= d1; i++)
 575         {
 576         ip = PLowerRow + uindex;
 577
 578         if(ip < PLowerRow)
 579           ip = PUpperRow;
 580
 581         for (j = 0; j <= d2; j++)
 582           {
 583           const gp_Pnt& P = Poles.Value(ip,jp);
 584           pole[0] = P.X();
 585           pole[1] = P.Y();
 586           pole[2] = P.Z();
 587           pole   += 3;
 588           ip++;
 589
 590           if (ip > PUpperRow)
 591             ip = PLowerRow;
 592           }
 593
 594         jp++;
 595
 596         if (jp > PUpperCol)
 597           jp = PLowerCol;
 598
 599         }
 600       }
 601
 602     return Standard_False;
 603     }
 604   }
 605
 606 //=======================================================================
 607 //function : D0
 608 //purpose  :
 609 //=======================================================================
 610
 611 void  BSplSLib::D0
 612 (const Standard_Real U,
 613  const Standard_Real V,
 614  const Standard_Integer UIndex,
 615  const Standard_Integer VIndex,
 616  const TColgp_Array2OfPnt& Poles,
 617  const TColStd_Array2OfReal& Weights,
 618  const TColStd_Array1OfReal& UKnots,
 619  const TColStd_Array1OfReal& VKnots,
 620  const TColStd_Array1OfInteger& UMults,
 621  const TColStd_Array1OfInteger& VMults,
 622  const Standard_Integer UDegree,
 623  const Standard_Integer VDegree,
 624  const Standard_Boolean URat,
 625  const Standard_Boolean VRat,
 626  const Standard_Boolean UPer,
 627  const Standard_Boolean VPer,
 628  gp_Pnt& P)
 629 {
 630 //  Standard_Integer k ;
 631   Standard_Real W ;
 632   HomogeneousD0(U,
 633                 V,
 634                 UIndex,
 635                 VIndex,
 636                 Poles,
 637                 Weights,
 638                 UKnots,
 639                 VKnots,
 640                 UMults,
 641                 VMults,
 642                 UDegree,
 643                 VDegree,
 644                 URat,
 645                 VRat,
 646                 UPer,
 647                 VPer,
 648                 W,
 649                 P) ;
 650   P.SetX(P.X() / W);
 651   P.SetY(P.Y() / W);
 652   P.SetZ(P.Z() / W);
 653 }
 654
 655 //=======================================================================
 656 //function : D0
 657 //purpose  :
 658 //=======================================================================
 659
 660 void  BSplSLib::HomogeneousD0
 661 (const Standard_Real U,
 662  const Standard_Real V,
 663  const Standard_Integer UIndex,
 664  const Standard_Integer VIndex,
 665  const TColgp_Array2OfPnt& Poles,
 666  const TColStd_Array2OfReal& Weights,
 667  const TColStd_Array1OfReal& UKnots,
 668  const TColStd_Array1OfReal& VKnots,
 669  const TColStd_Array1OfInteger& UMults,
 670  const TColStd_Array1OfInteger& VMults,
 671  const Standard_Integer UDegree,
 672  const Standard_Integer VDegree,
 673  const Standard_Boolean URat,
 674  const Standard_Boolean VRat,
 675  const Standard_Boolean UPer,
 676  const Standard_Boolean VPer,
 677  Standard_Real & W,
 678  gp_Pnt& P)
 679 {
 680   Standard_Boolean rational;
 681 //  Standard_Integer k,dim;
 682   Standard_Integer dim;
 683   Standard_Real u1,u2;
 684   Standard_Integer d1,d2;
 685   W = 1.0e0 ;
 686
 687   BSplSLib_DataContainer dc (UDegree, VDegree);
 688   PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
 689               Poles,Weights,UKnots,VKnots,UMults,VMults,
 690               u1,u2,d1,d2,rational,dc);
 691   if (rational) {
 692     dim = 4;
 693     BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
 694     BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
 695     W = dc.poles[3];
 696     P.SetX(dc.poles[0]);
 697     P.SetY(dc.poles[1]);
 698     P.SetZ(dc.poles[2]);
 699   }
 700   else {
 701     dim = 3;
 702     BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
 703     BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
 704     P.SetX(dc.poles[0]);
 705     P.SetY(dc.poles[1]);
 706     P.SetZ(dc.poles[2]);
 707   }
 708 }
 709
 710 //=======================================================================
 711 //function : D1
 712 //purpose  :
 713 //=======================================================================
 714
 715 void  BSplSLib::D1
 716 (const Standard_Real U,
 717  const Standard_Real V,
 718  const Standard_Integer UIndex,
 719  const Standard_Integer VIndex,
 720  const TColgp_Array2OfPnt& Poles,
 721  const TColStd_Array2OfReal& Weights,
 722  const TColStd_Array1OfReal& UKnots,
 723  const TColStd_Array1OfReal& VKnots,
 724  const TColStd_Array1OfInteger& UMults,
 725  const TColStd_Array1OfInteger& VMults,
 726  const Standard_Integer UDegree,
 727  const Standard_Integer VDegree,
 728  const Standard_Boolean URat,
 729  const Standard_Boolean VRat,
 730  const Standard_Boolean UPer,
 731  const Standard_Boolean VPer,
 732  gp_Pnt& P,
 733  gp_Vec& Vu,
 734  gp_Vec& Vv)
 735 {
 736   Standard_Boolean rational;
 737 //  Standard_Integer k,dim,dim2;
 738   Standard_Integer dim,dim2;
 739   Standard_Real u1,u2;
 740   Standard_Integer d1,d2;
 741   Standard_Real *result, *resVu, *resVv;
 742   BSplSLib_DataContainer dc (UDegree, VDegree);
 743   if (PrepareEval
 744     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
 745      Poles,Weights,UKnots,VKnots,UMults,VMults,
 746      u1,u2,d1,d2,rational,dc)) {
 747     if (rational) {
 748       dim  = 4;
 749       dim2 = (d2 + 1) << 2;
 750       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
 751       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
 752       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
 753       BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
 754       result = dc.ders;
 755       resVu  = result + 6;
 756       resVv  = result + 3;
 757     }
 758     else {
 759       dim  = 3;
 760       dim2 = d2 + 1;
 761       dim2 = (dim2 << 1) + dim2;
 762       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
 763       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
 764       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
 765       result = dc.poles;
 766       resVu  = result + dim2;
 767       resVv  = result + 3;
 768     }
 769   }
 770   else {
 771     if (rational) {
 772       dim  = 4;
 773       dim2 = (d2 + 1) << 2;
 774       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
 775       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
 776       BSplCLib::Eval(u2,d2,  *dc.knots2,dim ,*(dc.poles + dim2));
 777       BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
 778       result = dc.ders;
 779       resVu  = result + 3;
 780       resVv  = result + 6;
 781     }
 782     else {
 783       dim  = 3;
 784       dim2 = d2 + 1;
 785       dim2 = (dim2 << 1) + dim2;
 786       BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
 787       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
 788       BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + dim2));
 789       result = dc.poles;
 790       resVu  = result + 3;
 791       resVv  = result + dim2;
 792     }
 793   }
 794
 795   P .SetX(result[0]);
 796   Vu.SetX(resVu [0]);
 797   Vv.SetX(resVv [0]);
 798
 799   P .SetY(result[1]);
 800   Vu.SetY(resVu [1]);
 801   Vv.SetY(resVv [1]);
 802
 803   P .SetZ(result[2]);
 804   Vu.SetZ(resVu [2]);
 805   Vv.SetZ(resVv [2]);
 806 }
 807
 808 //=======================================================================
 809 //function : D1
 810 //purpose  :
 811 //=======================================================================
 812
 813 void  BSplSLib::HomogeneousD1
 814 (const Standard_Real U,
 815  const Standard_Real V,
 816  const Standard_Integer UIndex,
 817  const Standard_Integer VIndex,
 818  const TColgp_Array2OfPnt& Poles,
 819  const TColStd_Array2OfReal& Weights,
 820  const TColStd_Array1OfReal& UKnots,
 821  const TColStd_Array1OfReal& VKnots,
 822  const TColStd_Array1OfInteger& UMults,
 823  const TColStd_Array1OfInteger& VMults,
 824  const Standard_Integer UDegree,
 825  const Standard_Integer VDegree,
 826  const Standard_Boolean URat,
 827  const Standard_Boolean VRat,
 828  const Standard_Boolean UPer,
 829  const Standard_Boolean VPer,
 830  gp_Pnt& N,
 831  gp_Vec& Nu,
 832  gp_Vec& Nv,
 833  Standard_Real& D,
 834  Standard_Real& Du,
 835  Standard_Real& Dv)
 836 {
 837   Standard_Boolean rational;
 838 //  Standard_Integer k,dim;
 839   Standard_Integer dim;
 840   Standard_Real u1,u2;
 841   Standard_Integer d1,d2;
 842
 843   D = 1.0e0 ;
 844   Du = 0.0e0 ;
 845   Dv = 0.0e0 ;
 846   BSplSLib_DataContainer dc (UDegree, VDegree);
 847   Standard_Boolean ufirst = PrepareEval
 848     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
 849      Poles,Weights,UKnots,VKnots,UMults,VMults,
 850      u1,u2,d1,d2,rational,dc);
 851   dim  = rational ? 4 : 3;
 852
 853   BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
 854   BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
 855   BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
 856
 857   Standard_Real *result, *resVu, *resVv;
 858   result = dc.poles;
 859   resVu  = result + (ufirst ? dim*(d2+1) : dim);
 860   resVv  = result + (ufirst ? dim : dim*(d2+1));
 861
 862   N .SetX(result[0]);
 863   Nu.SetX(resVu [0]);
 864   Nv.SetX(resVv [0]);
 865
 866   N .SetY(result[1]);
 867   Nu.SetY(resVu [1]);
 868   Nv.SetY(resVv [1]);
 869
 870   N .SetZ(result[2]);
 871   Nu.SetZ(resVu [2]);
 872   Nv.SetZ(resVv [2]);
 873
 874   if (rational) {
 875     D  = result[3];
 876     Du = resVu [3];
 877     Dv = resVv [3];
 878   }
 879 }
 880
 881 //=======================================================================
 882 //function : D2
 883 //purpose  :
 884 //=======================================================================
 885
 886 void  BSplSLib::D2
 887 (const Standard_Real U,
 888  const Standard_Real V,
 889  const Standard_Integer UIndex,
 890  const Standard_Integer VIndex,
 891  const TColgp_Array2OfPnt& Poles,
 892  const TColStd_Array2OfReal& Weights,
 893  const TColStd_Array1OfReal& UKnots,
 894  const TColStd_Array1OfReal& VKnots,
 895  const TColStd_Array1OfInteger& UMults,
 896  const TColStd_Array1OfInteger& VMults,
 897  const Standard_Integer UDegree,
 898  const Standard_Integer VDegree,
 899  const Standard_Boolean URat,
 900  const Standard_Boolean VRat,
 901  const Standard_Boolean UPer,
 902  const Standard_Boolean VPer,
 903  gp_Pnt& P,
 904  gp_Vec& Vu,
 905  gp_Vec& Vv,
 906  gp_Vec& Vuu,
 907  gp_Vec& Vvv,
 908  gp_Vec& Vuv)
 909 {
 910   Standard_Boolean rational;
 911 //  Standard_Integer k,dim,dim2;
 912   Standard_Integer dim,dim2;
 913   Standard_Real u1,u2;
 914   Standard_Integer d1,d2;
 915   Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
 916   BSplSLib_DataContainer dc (UDegree, VDegree);
 917   if (PrepareEval
 918       (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
 919        Poles,Weights,UKnots,VKnots,UMults,VMults,
 920        u1,u2,d1,d2,rational,dc)) {
 921     if (rational) {
 922       dim = 4;
 923       dim2 = (d2 + 1) << 2;
 924       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
 925       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
 926       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
 927       if (d1 > 1)
 928         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
 929       BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
 930       result = dc.ders;
 931       resVu  = result + 9;
 932       resVv  = result + 3;
 933       resVuu = result + 18;
 934       resVvv = result + 6;
 935       resVuv = result + 12;
 936     }
 937     else {
 938       dim = 3;
 939       dim2 = d2 + 1;
 940       dim2 = (dim2 << 1) + dim2;
 941       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
 942       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
 943       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
 944       if (d1 > 1)
 945         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
 946       result = dc.poles;
 947       resVu  = result + dim2;
 948       resVv  = result + 3;
 949       if (UDegree <= 1) resVuu = BSplSLib_zero;
 950       else              resVuu = result + (dim2 << 1);
 951       if (VDegree <= 1) resVvv = BSplSLib_zero;
 952       else              resVvv = result + 6;
 953       resVuv = result + (d2 << 1) + d2 + 6;
 954     }
 955   }
 956   else {
 957     if (rational) {
 958       dim = 4;
 959       dim2 = (d2 + 1) << 2;
 960       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
 961       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
 962       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
 963       if (d1 > 1)
 964         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
 965       BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
 966       result = dc.ders;
 967       resVu  = result + 3;
 968       resVv  = result + 9;
 969       resVuu = result + 6;
 970       resVvv = result + 18;
 971       resVuv = result + 12;
 972     }
 973     else {
 974       dim = 3;
 975       dim2 = d2 + 1;
 976       dim2 = (dim2 << 1) + dim2;
 977       BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
 978       BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
 979       BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
 980       if (d1 > 1)
 981         BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
 982       result = dc.poles;
 983       resVu  = result + 3;
 984       resVv  = result + dim2;
 985       if (UDegree <= 1) resVuu = BSplSLib_zero;
 986       else              resVuu = result + 6;
 987       if (VDegree <= 1) resVvv = BSplSLib_zero;
 988       else              resVvv = result + (dim2 << 1);
 989       resVuv = result + (d2 << 1) + d2 + 6;
 990     }
 991   }
 992
 993   P  .SetX(result[0]);
 994   Vu .SetX(resVu [0]);
 995   Vv .SetX(resVv [0]);
 996   Vuu.SetX(resVuu[0]);
 997   Vvv.SetX(resVvv[0]);
 998   Vuv.SetX(resVuv[0]);
 999
1000   P  .SetY(result[1]);
1001   Vu .SetY(resVu [1]);
1002   Vv .SetY(resVv [1]);
1003   Vuu.SetY(resVuu[1]);
1004   Vvv.SetY(resVvv[1]);
1005   Vuv.SetY(resVuv[1]);
1006
1007   P  .SetZ(result[2]);
1008   Vu .SetZ(resVu [2]);
1009   Vv .SetZ(resVv [2]);
1010   Vuu.SetZ(resVuu[2]);
1011   Vvv.SetZ(resVvv[2]);
1012   Vuv.SetZ(resVuv[2]);
1013 }
1014
1015 //=======================================================================
1016 //function : D3
1017 //purpose  :
1018 //=======================================================================
1019
1020 void  BSplSLib::D3
1021 (const Standard_Real U,
1022  const Standard_Real V,
1023  const Standard_Integer UIndex,
1024  const Standard_Integer VIndex,
1025  const TColgp_Array2OfPnt& Poles,
1026  const TColStd_Array2OfReal& Weights,
1027  const TColStd_Array1OfReal& UKnots,
1028  const TColStd_Array1OfReal& VKnots,
1029  const TColStd_Array1OfInteger& UMults,
1030  const TColStd_Array1OfInteger& VMults,
1031  const Standard_Integer UDegree,
1032  const Standard_Integer VDegree,
1033  const Standard_Boolean URat,
1034  const Standard_Boolean VRat,
1035  const Standard_Boolean UPer,
1036  const Standard_Boolean VPer,
1037  gp_Pnt& P,
1038  gp_Vec& Vu,
1039  gp_Vec& Vv,
1040  gp_Vec& Vuu,
1041  gp_Vec& Vvv,
1042  gp_Vec& Vuv,
1043  gp_Vec& Vuuu,
1044  gp_Vec& Vvvv,
1045  gp_Vec& Vuuv,
1046  gp_Vec& Vuvv)
1047 {
1048   Standard_Boolean rational;
1049 //  Standard_Integer k,dim,dim2;
1050   Standard_Integer dim,dim2;
1051   Standard_Real u1,u2;
1052   Standard_Integer d1,d2;
1053   Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
1054   *resVuuu, *resVvvv, *resVuuv, *resVuvv;
1055   BSplSLib_DataContainer dc (UDegree, VDegree);
1056   if (PrepareEval
1057     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1058      Poles,Weights,UKnots,VKnots,UMults,VMults,
1059      u1,u2,d1,d2,rational,dc)) {
1060     if (rational) {
1061       dim = 4;
1062       dim2 = (d2 + 1) << 2;
1063       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1064       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1065       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1066       if (d1 > 1)
1067         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1068       if (d1 > 2)
1069         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1070       BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1071       result  = dc.ders;
1072       resVu   = result + 12;
1073       resVv   = result + 3;
1074       resVuu  = result + 24;
1075       resVvv  = result + 6;
1076       resVuv  = result + 15;
1077       resVuuu = result + 36;
1078       resVvvv = result + 9;
1079       resVuuv = result + 27;
1080       resVuvv = result + 18;
1081     }
1082     else {
1083       dim = 3;
1084       dim2 = (d2 + 1);
1085       dim2 = (dim2 << 1) + dim2;
1086       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1087       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1088       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1089       if (d1 > 1)
1090         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1091       if (d1 > 2)
1092         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1093       result = dc.poles;
1094       resVu  = result + dim2;
1095       resVv  = result + 3;
1096       if (UDegree <= 1) {
1097         resVuu  = BSplSLib_zero;
1098         resVuuv = BSplSLib_zero;
1099       }
1100       else {
1101         resVuu  = result + (dim2 << 1);
1102         resVuuv = result + (dim2 << 1) + 3;
1103       }
1104       if (VDegree <= 1) {
1105         resVvv  = BSplSLib_zero;
1106         resVuvv = BSplSLib_zero;
1107       }
1108       else {
1109         resVvv  = result + 6;
1110         resVuvv = result + dim2 + 6;
1111       }
1112       resVuv = result + (d2 << 1) + d2 + 6;
1113       if (UDegree <= 2) resVuuu = BSplSLib_zero;
1114       else              resVuuu = result + (dim2 << 1) + dim2;
1115       if (VDegree <= 2) resVvvv = BSplSLib_zero;
1116       else              resVvvv = result + 9;
1117     }
1118   }
1119   else {
1120     if (rational) {
1121       dim = 4;
1122       dim2 = (d2 + 1) << 2;
1123       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1124       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1125       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1126       if (d1 > 1)
1127         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1128       if (d1 > 2)
1129         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1130       BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1131       result  = dc.ders;
1132       resVu   = result + 3;
1133       resVv   = result + 12;
1134       resVuu  = result + 6;
1135       resVvv  = result + 24;
1136       resVuv  = result + 15;
1137       resVuuu = result + 9;
1138       resVvvv = result + 36;
1139       resVuuv = result + 18;
1140       resVuvv = result + 27;
1141     }
1142     else {
1143       dim = 3;
1144       dim2 = (d2 + 1);
1145       dim2 = (dim2 << 1) + dim2;
1146       BSplCLib::Bohm  (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1147       BSplCLib::Bohm  (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1148       BSplCLib::Bohm  (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1149       if (d1 > 1)
1150         BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1151       if (d1 > 2)
1152         BSplCLib::Eval(u2,d2  ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1153       result = dc.poles;
1154       resVu  = result + 3;
1155       resVv  = result + dim2;
1156       if (UDegree <= 1) {
1157         resVuu  = BSplSLib_zero;
1158         resVuuv = BSplSLib_zero;
1159       }
1160       else {
1161         resVuu  = result + 6;
1162         resVuuv = result + dim2 + 6;
1163       }
1164       if (VDegree <= 1) {
1165         resVvv  = BSplSLib_zero;
1166         resVuvv = BSplSLib_zero;
1167       }
1168       else {
1169         resVvv  = result + (dim2 << 1);
1170         resVuvv = result + (dim2 << 1) + 3;
1171       }
1172       resVuv = result + (d2 << 1) + d2 + 6;
1173       if (UDegree <= 2) resVuuu = BSplSLib_zero;
1174       else              resVuuu = result + 9;
1175       if (VDegree <= 2) resVvvv = BSplSLib_zero;
1176       else              resVvvv = result + (dim2 << 1) + dim2;
1177     }
1178   }
1179
1180   P   .SetX(result [0]);
1181   Vu  .SetX(resVu  [0]);
1182   Vv  .SetX(resVv  [0]);
1183   Vuu .SetX(resVuu [0]);
1184   Vvv .SetX(resVvv [0]);
1185   Vuv .SetX(resVuv [0]);
1186   Vuuu.SetX(resVuuu[0]);
1187   Vvvv.SetX(resVvvv[0]);
1188   Vuuv.SetX(resVuuv[0]);
1189   Vuvv.SetX(resVuvv[0]);
1190
1191   P   .SetY(result [1]);
1192   Vu  .SetY(resVu  [1]);
1193   Vv  .SetY(resVv  [1]);
1194   Vuu .SetY(resVuu [1]);
1195   Vvv .SetY(resVvv [1]);
1196   Vuv .SetY(resVuv [1]);
1197   Vuuu.SetY(resVuuu[1]);
1198   Vvvv.SetY(resVvvv[1]);
1199   Vuuv.SetY(resVuuv[1]);
1200   Vuvv.SetY(resVuvv[1]);
1201
1202   P   .SetZ(result [2]);
1203   Vu  .SetZ(resVu  [2]);
1204   Vv  .SetZ(resVv  [2]);
1205   Vuu .SetZ(resVuu [2]);
1206   Vvv .SetZ(resVvv [2]);
1207   Vuv .SetZ(resVuv [2]);
1208   Vuuu.SetZ(resVuuu[2]);
1209   Vvvv.SetZ(resVvvv[2]);
1210   Vuuv.SetZ(resVuuv[2]);
1211   Vuvv.SetZ(resVuvv[2]);
1212 }
1213
1214 //=======================================================================
1215 //function : DN
1216 //purpose  :
1217 //=======================================================================
1218
1219 void  BSplSLib::DN
1220 (const Standard_Real U,
1221  const Standard_Real V,
1222  const Standard_Integer Nu,
1223  const Standard_Integer Nv,
1224  const Standard_Integer UIndex,
1225  const Standard_Integer VIndex,
1226  const TColgp_Array2OfPnt& Poles,
1227  const TColStd_Array2OfReal& Weights,
1228  const TColStd_Array1OfReal& UKnots,
1229  const TColStd_Array1OfReal& VKnots,
1230  const TColStd_Array1OfInteger& UMults,
1231  const TColStd_Array1OfInteger& VMults,
1232  const Standard_Integer UDegree,
1233  const Standard_Integer VDegree,
1234  const Standard_Boolean URat,
1235  const Standard_Boolean VRat,
1236  const Standard_Boolean UPer,
1237  const Standard_Boolean VPer,
1238  gp_Vec& Vn)
1239 {
1240   Standard_Boolean rational;
1241   Standard_Integer k,dim;
1242   Standard_Real u1,u2;
1243   Standard_Integer d1,d2;
1244
1245   BSplSLib_DataContainer dc (UDegree, VDegree);
1246   Standard_Boolean ufirst = PrepareEval
1247     (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1248      Poles,Weights,UKnots,VKnots,UMults,VMults,
1249      u1,u2,d1,d2,rational,dc);
1250   dim  = rational ? 4 : 3;
1251
1252   if (!rational) {
1253     if ((Nu > UDegree) || (Nv > VDegree)) {
1254       Vn.SetX(0.);
1255       Vn.SetY(0.);
1256       Vn.SetZ(0.);
1257       return;
1258     }
1259   }
1260
1261   Standard_Integer n1 = ufirst ? Nu : Nv;
1262   Standard_Integer n2 = ufirst ? Nv : Nu;
1263
1264   BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1265
1266   for (k = 0; k <= Min(n1,d1); k++)
1267     BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1268
1269   Standard_Real *result;
1270   if (rational) {
1271     BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1272     result = dc.ders; // because Standard_False ci-dessus.
1273
1274   }
1275   else {
1276     result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1277   }
1278
1279   Vn.SetX(result[0]);
1280   Vn.SetY(result[1]);
1281   Vn.SetZ(result[2]);
1282 }
1283
1284 //
1285 // Surface modifications
1286 //
1287 // a surface is processed as a curve of curves.
1288 // i.e : a curve of parameter u where the current point is the set of poles
1289 //       of the iso.
1290 //
1291
1292 //=======================================================================
1293 //function : Iso
1294 //purpose  :
1295 //=======================================================================
1296
1297 void  BSplSLib::Iso(const Standard_Real            Param,
1298                     const Standard_Boolean         IsU,
1299                     const TColgp_Array2OfPnt&      Poles,
1300                     const TColStd_Array2OfReal&    Weights,
1301                     const TColStd_Array1OfReal&    Knots,
1302                     const TColStd_Array1OfInteger& Mults,
1303                     const Standard_Integer         Degree,
1304                     const Standard_Boolean         Periodic,
1305                     TColgp_Array1OfPnt&            CPoles,
1306                     TColStd_Array1OfReal&          CWeights)
1307 {
1308   Standard_Integer index = 0;
1309   Standard_Real    u = Param;
1310   Standard_Boolean rational = &Weights != NULL;
1311   Standard_Integer dim = rational ? 4 : 3;
1312
1313   // compute local knots
1314
1315   NCollection_LocalArray<Standard_Real> locknots1 (2*Degree);
1316   BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1317   BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1318   if (&Mults == NULL)
1319     index -= Knots.Lower() + Degree;
1320   else
1321     index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1322
1323
1324   // copy the local poles
1325
1326 //  Standard_Integer f1,l1,f2,l2,i,j,k;
1327   Standard_Integer f1,l1,f2,l2,i,j;
1328
1329   if (IsU) {
1330     f1 = Poles.LowerRow();
1331     l1 = Poles.UpperRow();
1332     f2 = Poles.LowerCol();
1333     l2 = Poles.UpperCol();
1334   }
1335   else {
1336     f1 = Poles.LowerCol();
1337     l1 = Poles.UpperCol();
1338     f2 = Poles.LowerRow();
1339     l2 = Poles.UpperRow();
1340   }
1341
1342   NCollection_LocalArray<Standard_Real> locpoles ((Degree+1) * (l2-f2+1) * dim);
1343
1344   Standard_Real w, *pole = locpoles;
1345   index += f1;
1346
1347   for (i = 0; i <= Degree; i++) {
1348
1349     for (j = f2; j <= l2; j++) {
1350
1351       const gp_Pnt& P  = IsU ? Poles(index,j)   : Poles(j,index);
1352       if (rational) {
1353         pole[3] = w      = IsU ? Weights(index,j) : Weights(j,index);
1354         pole[0] = P.X() * w;
1355         pole[1] = P.Y() * w;
1356         pole[2] = P.Z() * w;
1357       }
1358       else {
1359         pole[0] = P.X();
1360         pole[1] = P.Y();
1361         pole[2] = P.Z();
1362       }
1363       pole += dim;
1364     }
1365     index++;
1366     if (index > l1) index = f1;
1367   }
1368
1369   // compute the iso
1370   BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1371
1372   // get the result
1373   pole = locpoles;
1374
1375   for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1376     gp_Pnt& P = CPoles(i);
1377     if (rational) {
1378       CWeights(i) = w = pole[3];
1379       P.SetX( pole[0] / w);
1380       P.SetY( pole[1] / w);
1381       P.SetZ( pole[2] / w);
1382     }
1383     else {
1384       P.SetX( pole[0]);
1385       P.SetY( pole[1]);
1386       P.SetZ( pole[2]);
1387     }
1388     pole += dim;
1389   }
1390
1391   // if the input is not rational but weights are wanted
1392   if (!rational && (&CWeights != NULL)) {
1393
1394     for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1395       CWeights(i) = 1.;
1396   }
1397 }
1398
1399 //=======================================================================
1400 //function : Reverse
1401 //purpose  :
1402 //=======================================================================
1403
1404 void  BSplSLib::Reverse(      TColgp_Array2OfPnt& Poles,
1405                         const Standard_Integer    Last,
1406                         const Standard_Boolean    UDirection)
1407 {
1408   Standard_Integer i,j, l = Last;
1409   if ( UDirection) {
1410     l = Poles.LowerRow() +
1411       (l - Poles.LowerRow())%(Poles.ColLength());
1412     TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1413                             Poles.LowerCol(), Poles.UpperCol());
1414
1415     for (i = Poles.LowerRow(); i <= l; i++) {
1416
1417       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1418         temp(l-i,j) = Poles(i,j);
1419       }
1420     }
1421
1422     for (i = l+1; i <= Poles.UpperRow(); i++) {
1423
1424       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1425         temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1426       }
1427     }
1428
1429     for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1430
1431       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1432         Poles(i,j) = temp (i-Poles.LowerRow(),j);
1433       }
1434     }
1435   }
1436   else {
1437     l = Poles.LowerCol() +
1438       (l - Poles.LowerCol())%(Poles.RowLength());
1439     TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1440                             0, Poles.RowLength()-1);
1441
1442     for (j = Poles.LowerCol(); j <= l; j++) {
1443
1444       for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1445         temp(i,l-j) = Poles(i,j);
1446       }
1447     }
1448
1449     for (j = l+1; j <= Poles.UpperCol(); j++) {
1450
1451       for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1452         temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1453       }
1454     }
1455
1456     for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1457
1458       for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1459         Poles(i,j) = temp (i,j-Poles.LowerCol());
1460       }
1461     }
1462   }
1463 }
1464
1465 //=======================================================================
1466 //function : Reverse
1467 //purpose  :
1468 //=======================================================================
1469
1470 void  BSplSLib::Reverse(      TColStd_Array2OfReal& Weights,
1471                         const Standard_Integer      Last,
1472                         const Standard_Boolean      UDirection)
1473 {
1474   Standard_Integer i,j, l = Last;
1475   if ( UDirection) {
1476     l = Weights.LowerRow() +
1477       (l - Weights.LowerRow())%(Weights.ColLength());
1478     TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1479                               Weights.LowerCol(), Weights.UpperCol());
1480
1481     for (i = Weights.LowerRow(); i <= l; i++) {
1482
1483       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1484         temp(l-i,j) = Weights(i,j);
1485       }
1486     }
1487
1488     for (i = l+1; i <= Weights.UpperRow(); i++) {
1489
1490       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1491         temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1492       }
1493     }
1494
1495     for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1496
1497       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1498         Weights(i,j) = temp (i-Weights.LowerRow(),j);
1499       }
1500     }
1501   }
1502   else {
1503     l = Weights.LowerCol() +
1504       (l - Weights.LowerCol())%(Weights.RowLength());
1505     TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1506                               0, Weights.RowLength()-1);
1507
1508     for (j = Weights.LowerCol(); j <= l; j++) {
1509
1510       for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1511         temp(i,l-j) = Weights(i,j);
1512       }
1513     }
1514
1515     for (j = l+1; j <= Weights.UpperCol(); j++) {
1516
1517       for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1518         temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1519       }
1520     }
1521
1522     for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1523
1524       for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1525         Weights(i,j) = temp (i,j-Weights.LowerCol());
1526       }
1527     }
1528   }
1529 }
1530
1531 //=======================================================================
1532 //function : IsRational
1533 //purpose  :
1534 //=======================================================================
1535
1536 Standard_Boolean BSplSLib::IsRational
1537 (const TColStd_Array2OfReal& Weights,
1538  const Standard_Integer      I1,
1539  const Standard_Integer      I2,
1540  const Standard_Integer      J1,
1541  const Standard_Integer      J2,
1542  const Standard_Real         Epsi)
1543 {
1544   Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1545   Standard_Integer i,j;
1546   Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1547   Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1548
1549   for (i = I1 - fi; i < I2 - fi; i++) {
1550
1551     for (j = J1 - fj; j < J2 - fj; j++) {
1552       if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1553         return Standard_True;
1554     }
1555   }
1556   return Standard_False;
1557 }
1558
1559 //=======================================================================
1560 //function : SetPoles
1561 //purpose  :
1562 //=======================================================================
1563
1564 void  BSplSLib::SetPoles(const TColgp_Array2OfPnt&   Poles,
1565                          TColStd_Array1OfReal& FP,
1566                          const Standard_Boolean      UDirection)
1567 {
1568   Standard_Integer i, j, l = FP.Lower();
1569   Standard_Integer PLowerRow = Poles.LowerRow();
1570   Standard_Integer PUpperRow = Poles.UpperRow();
1571   Standard_Integer PLowerCol = Poles.LowerCol();
1572   Standard_Integer PUpperCol = Poles.UpperCol();
1573   if (UDirection) {
1574
1575     for ( i = PLowerRow; i <= PUpperRow; i++) {
1576
1577       for ( j = PLowerCol; j <= PUpperCol; j++) {
1578         const gp_Pnt& P = Poles.Value(i,j);
1579         FP(l) = P.X(); l++;
1580         FP(l) = P.Y(); l++;
1581         FP(l) = P.Z(); l++;
1582       }
1583     }
1584   }
1585   else {
1586
1587     for ( j = PLowerCol; j <= PUpperCol; j++) {
1588
1589       for ( i = PLowerRow; i <= PUpperRow; i++) {
1590         const gp_Pnt& P = Poles.Value(i,j);
1591         FP(l) = P.X(); l++;
1592         FP(l) = P.Y(); l++;
1593         FP(l) = P.Z(); l++;
1594       }
1595     }
1596   }
1597 }
1598
1599 //=======================================================================
1600 //function : SetPoles
1601 //purpose  :
1602 //=======================================================================
1603
1604 void  BSplSLib::SetPoles(const TColgp_Array2OfPnt&   Poles,
1605                          const TColStd_Array2OfReal& Weights,
1606                          TColStd_Array1OfReal& FP,
1607                          const Standard_Boolean      UDirection)
1608 {
1609   Standard_Integer i, j, l = FP.Lower();
1610   Standard_Integer PLowerRow = Poles.LowerRow();
1611   Standard_Integer PUpperRow = Poles.UpperRow();
1612   Standard_Integer PLowerCol = Poles.LowerCol();
1613   Standard_Integer PUpperCol = Poles.UpperCol();
1614   if (UDirection) {
1615
1616     for ( i = PLowerRow; i <= PUpperRow; i++) {
1617
1618       for ( j = PLowerCol; j <= PUpperCol; j++) {
1619         const gp_Pnt& P = Poles  .Value(i,j);
1620         Standard_Real w = Weights.Value(i,j);
1621         FP(l) = P.X() * w; l++;
1622         FP(l) = P.Y() * w; l++;
1623         FP(l) = P.Z() * w; l++;
1624         FP(l) = w; l++;
1625       }
1626     }
1627   }
1628   else {
1629
1630     for ( j = PLowerCol; j <= PUpperCol; j++) {
1631
1632       for ( i = PLowerRow; i <= PUpperRow; i++) {
1633         const gp_Pnt& P = Poles  .Value(i,j);
1634         Standard_Real w = Weights.Value(i,j);
1635         FP(l) = P.X() * w; l++;
1636         FP(l) = P.Y() * w; l++;
1637         FP(l) = P.Z() * w; l++;
1638         FP(l) = w; l++;
1639       }
1640     }
1641   }
1642 }
1643
1644 //=======================================================================
1645 //function : GetPoles
1646 //purpose  :
1647 //=======================================================================
1648
1649 void  BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1650                          TColgp_Array2OfPnt&   Poles,
1651                          const Standard_Boolean      UDirection)
1652 {
1653   Standard_Integer i, j, l = FP.Lower();
1654   Standard_Integer PLowerRow = Poles.LowerRow();
1655   Standard_Integer PUpperRow = Poles.UpperRow();
1656   Standard_Integer PLowerCol = Poles.LowerCol();
1657   Standard_Integer PUpperCol = Poles.UpperCol();
1658   if (UDirection) {
1659
1660     for ( i = PLowerRow; i <= PUpperRow; i++) {
1661
1662       for ( j = PLowerCol; j <= PUpperCol; j++) {
1663         gp_Pnt& P = Poles.ChangeValue(i,j);
1664         P.SetX(FP(l)); l++;
1665         P.SetY(FP(l)); l++;
1666         P.SetZ(FP(l)); l++;
1667       }
1668     }
1669   }
1670   else {
1671
1672     for ( j = PLowerCol; j <= PUpperCol; j++) {
1673
1674       for ( i = PLowerRow; i <= PUpperRow; i++) {
1675         gp_Pnt& P = Poles.ChangeValue(i,j);
1676         P.SetX(FP(l)); l++;
1677         P.SetY(FP(l)); l++;
1678         P.SetZ(FP(l)); l++;
1679       }
1680     }
1681   }
1682 }
1683
1684 //=======================================================================
1685 //function : GetPoles
1686 //purpose  :
1687 //=======================================================================
1688
1689 void  BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1690                          TColgp_Array2OfPnt&   Poles,
1691                          TColStd_Array2OfReal& Weights,
1692                          const Standard_Boolean      UDirection)
1693 {
1694   Standard_Integer i, j, l = FP.Lower();
1695   Standard_Integer PLowerRow = Poles.LowerRow();
1696   Standard_Integer PUpperRow = Poles.UpperRow();
1697   Standard_Integer PLowerCol = Poles.LowerCol();
1698   Standard_Integer PUpperCol = Poles.UpperCol();
1699   if (UDirection) {
1700
1701     for ( i = PLowerRow; i <= PUpperRow; i++) {
1702
1703       for ( j = PLowerCol; j <= PUpperCol; j++) {
1704         Standard_Real w = FP( l + 3);
1705         Weights(i,j) = w;
1706         gp_Pnt& P = Poles.ChangeValue(i,j);
1707         P.SetX(FP(l) / w); l++;
1708         P.SetY(FP(l) / w); l++;
1709         P.SetZ(FP(l) / w); l++;
1710         l++;
1711       }
1712     }
1713   }
1714   else {
1715
1716     for ( j = PLowerCol; j <= PUpperCol; j++) {
1717
1718       for ( i = PLowerRow; i <= PUpperRow; i++) {
1719         Standard_Real w = FP( l + 3);
1720         Weights(i,j) = w;
1721         gp_Pnt& P = Poles.ChangeValue(i,j);
1722         P.SetX(FP(l) / w); l++;
1723         P.SetY(FP(l) / w); l++;
1724         P.SetZ(FP(l) / w); l++;
1725         l++;
1726       }
1727     }
1728   }
1729 }
1730
1731 //=======================================================================
1732 //function : InsertKnots
1733 //purpose  :
1734 //=======================================================================
1735
1736 void  BSplSLib::InsertKnots(const Standard_Boolean         UDirection,
1737                             const Standard_Integer         Degree,
1738                             const Standard_Boolean         Periodic,
1739                             const TColgp_Array2OfPnt&      Poles,
1740                             const TColStd_Array2OfReal&    Weights,
1741                             const TColStd_Array1OfReal&    Knots,
1742                             const TColStd_Array1OfInteger& Mults,
1743                             const TColStd_Array1OfReal&    AddKnots,
1744                             const TColStd_Array1OfInteger& AddMults,
1745                             TColgp_Array2OfPnt&      NewPoles,
1746                             TColStd_Array2OfReal&    NewWeights,
1747                             TColStd_Array1OfReal&    NewKnots,
1748                             TColStd_Array1OfInteger& NewMults,
1749                             const Standard_Real            Epsilon,
1750                             const Standard_Boolean         Add )
1751 {
1752   Standard_Boolean rational = &Weights != NULL;
1753   Standard_Integer dim = 3;
1754   if (rational) dim++;
1755
1756   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1757   TColStd_Array1OfReal
1758     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1759
1760   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1761   else          SetPoles(Poles,poles,UDirection);
1762
1763   if (UDirection) {
1764     dim *= Poles.RowLength();
1765   }
1766   else {
1767     dim *= Poles.ColLength();
1768   }
1769   BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1770                         AddKnots,AddMults,newpoles,NewKnots,NewMults,
1771                         Epsilon,Add);
1772
1773   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1774   else          GetPoles(newpoles,NewPoles,UDirection);
1775 }
1776
1777 //=======================================================================
1778 //function : RemoveKnot
1779 //purpose  :
1780 //=======================================================================
1781
1782 Standard_Boolean  BSplSLib::RemoveKnot
1783 (const Standard_Boolean         UDirection,
1784  const Standard_Integer         Index,
1785  const Standard_Integer         Mult,
1786  const Standard_Integer         Degree,
1787  const Standard_Boolean         Periodic,
1788  const TColgp_Array2OfPnt&      Poles,
1789  const TColStd_Array2OfReal&    Weights,
1790  const TColStd_Array1OfReal&    Knots,
1791  const TColStd_Array1OfInteger& Mults,
1792  TColgp_Array2OfPnt&      NewPoles,
1793  TColStd_Array2OfReal&    NewWeights,
1794  TColStd_Array1OfReal&    NewKnots,
1795  TColStd_Array1OfInteger& NewMults,
1796  const Standard_Real            Tolerance)
1797 {
1798   Standard_Boolean rational = &Weights != NULL;
1799   Standard_Integer dim = 3;
1800   if (rational) dim++;
1801
1802   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1803   TColStd_Array1OfReal
1804     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1805
1806   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1807   else          SetPoles(Poles,poles,UDirection);
1808
1809   if (UDirection) {
1810     dim *= Poles.RowLength();
1811   }
1812   else {
1813     dim *= Poles.ColLength();
1814   }
1815
1816   if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1817                              poles,Knots,Mults,newpoles,NewKnots,NewMults,
1818                              Tolerance))
1819     return Standard_False;
1820
1821   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1822   else          GetPoles(newpoles,NewPoles,UDirection);
1823   return Standard_True;
1824 }
1825
1826 //=======================================================================
1827 //function : IncreaseDegree
1828 //purpose  :
1829 //=======================================================================
1830
1831 void  BSplSLib::IncreaseDegree
1832 (const Standard_Boolean         UDirection,
1833  const Standard_Integer         Degree,
1834  const Standard_Integer         NewDegree,
1835  const Standard_Boolean         Periodic,
1836  const TColgp_Array2OfPnt&      Poles,
1837  const TColStd_Array2OfReal&    Weights,
1838  const TColStd_Array1OfReal&    Knots,
1839  const TColStd_Array1OfInteger& Mults,
1840  TColgp_Array2OfPnt&      NewPoles,
1841  TColStd_Array2OfReal&    NewWeights,
1842  TColStd_Array1OfReal&    NewKnots,
1843  TColStd_Array1OfInteger& NewMults)
1844 {
1845   Standard_Boolean rational = &Weights != NULL;
1846   Standard_Integer dim = 3;
1847   if (rational) dim++;
1848
1849   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1850   TColStd_Array1OfReal
1851     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1852
1853   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1854   else          SetPoles(Poles,poles,UDirection);
1855
1856   if (UDirection) {
1857     dim *= Poles.RowLength();
1858   }
1859   else {
1860     dim *= Poles.ColLength();
1861   }
1862
1863   BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1864                            newpoles,NewKnots,NewMults);
1865
1866   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1867   else          GetPoles(newpoles,NewPoles,UDirection);
1868 }
1869
1870 //=======================================================================
1871 //function : Unperiodize
1872 //purpose  :
1873 //=======================================================================
1874
1875 void  BSplSLib::Unperiodize
1876 (const Standard_Boolean         UDirection,
1877  const Standard_Integer         Degree,
1878  const TColStd_Array1OfInteger& Mults,
1879  const TColStd_Array1OfReal&    Knots,
1880  const TColgp_Array2OfPnt&      Poles,
1881  const TColStd_Array2OfReal&    Weights,
1882  TColStd_Array1OfInteger& NewMults,
1883  TColStd_Array1OfReal&    NewKnots,
1884  TColgp_Array2OfPnt&      NewPoles,
1885  TColStd_Array2OfReal&    NewWeights)
1886 {
1887   Standard_Boolean rational = &Weights != NULL;
1888   Standard_Integer dim = 3;
1889   if (rational) dim++;
1890
1891   TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1892   TColStd_Array1OfReal
1893     newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1894
1895   if (rational) SetPoles(Poles,Weights,poles,UDirection);
1896   else          SetPoles(Poles,poles,UDirection);
1897
1898   if (UDirection) {
1899     dim *= Poles.RowLength();
1900   }
1901   else {
1902     dim *= Poles.ColLength();
1903   }
1904
1905   BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1906                         NewMults, NewKnots, newpoles);
1907
1908   if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1909   else          GetPoles(newpoles,NewPoles,UDirection);
1910 }
1911
1912 //=======================================================================
1913 //function : BuildCache
1914 //purpose  : Stores theTaylor expansion normalized between 0,1 in the
1915 //           Cache : in case of  a rational function the Poles are
1916 //           stored in homogeneous form
1917 //=======================================================================
1918
1919 void BSplSLib::BuildCache
1920 (const Standard_Real            U,
1921  const Standard_Real            V,
1922  const Standard_Real            USpanDomain,
1923  const Standard_Real            VSpanDomain,
1924  const Standard_Boolean         UPeriodic,
1925  const Standard_Boolean         VPeriodic,
1926  const Standard_Integer         UDegree,
1927  const Standard_Integer         VDegree,
1928  const Standard_Integer         UIndex,
1929  const Standard_Integer         VIndex,
1930  const TColStd_Array1OfReal&    UFlatKnots,
1931  const TColStd_Array1OfReal&    VFlatKnots,
1932  const TColgp_Array2OfPnt&      Poles,
1933  const TColStd_Array2OfReal&    Weights,
1934  TColgp_Array2OfPnt&            CachePoles,
1935  TColStd_Array2OfReal&          CacheWeights)
1936 {
1937   Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1938   Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1939   Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1940   if (&Weights != NULL)
1941     rational_u = rational_v = Standard_True;
1942   else
1943     rational_u = rational_v = Standard_False;
1944   BSplSLib_DataContainer dc (UDegree, VDegree);
1945   flag_u_or_v =
1946     PrepareEval  (U,
1947                   V,
1948                   UIndex,
1949                   VIndex,
1950                   UDegree,
1951                   VDegree,
1952                   rational_u,
1953                   rational_v,
1954                   UPeriodic,
1955                   VPeriodic,
1956                   Poles,
1957                   Weights,
1958                   UFlatKnots,
1959                   VFlatKnots,
1960                   (BSplCLib::NoMults()),
1961                   (BSplCLib::NoMults()),
1962                   u1,
1963                   u2,
1964                   d1,
1965                   d2,
1966                   rational,
1967           dc);
1968   d1p1 = d1 + 1;
1969   d2p1 = d2 + 1;
1970   if (rational) {
1971     BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1972
1973     for (kk = 0; kk <= d1 ; kk++)
1974       BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1975     if (flag_u_or_v) {
1976       min_degree_domain = USpanDomain ;
1977       max_degree_domain = VSpanDomain ;
1978     }
1979     else {
1980       min_degree_domain = VSpanDomain ;
1981       max_degree_domain = USpanDomain ;
1982     }
1983     factor[0] = 1.0e0 ;
1984
1985     for (ii = 0 ; ii <= d2 ; ii++) {
1986       iii = ii + 1;
1987       factor[1] = 1.0e0 ;
1988
1989       for (jj = 0 ; jj <= d1 ; jj++) {
1990         jjj = jj + 1;
1991         Index = jj * d2p1 + ii ;
1992         Index = Index << 2;
1993         gp_Pnt& P = CachePoles(iii,jjj);
1994         f = factor[0] * factor[1];
1995         P.SetX( f * dc.poles[Index]); Index++;
1996         P.SetY( f * dc.poles[Index]); Index++;
1997         P.SetZ( f * dc.poles[Index]); Index++;
1998         CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1999         factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
2000       }
2001       factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
2002     }
2003   }
2004   else {
2005     BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
2006
2007     for (kk = 0; kk <= d1 ; kk++)
2008       BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
2009     if (flag_u_or_v) {
2010       min_degree_domain = USpanDomain ;
2011       max_degree_domain = VSpanDomain ;
2012     }
2013     else {
2014       min_degree_domain = VSpanDomain ;
2015       max_degree_domain = USpanDomain ;
2016     }
2017     factor[0] = 1.0e0 ;
2018
2019     for (ii = 0 ; ii <= d2 ; ii++) {
2020       iii = ii + 1;
2021       factor[1] = 1.0e0 ;
2022
2023       for (jj = 0 ; jj <= d1 ; jj++) {
2024         jjj = jj + 1;
2025         Index = jj * d2p1 + ii ;
2026         Index = (Index << 1) + Index;
2027         gp_Pnt& P = CachePoles(iii,jjj);
2028         f = factor[0] * factor[1];
2029         P.SetX( f * dc.poles[Index]); Index++;
2030         P.SetY( f * dc.poles[Index]); Index++;
2031         P.SetZ( f * dc.poles[Index]);
2032         factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
2033       }
2034       factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
2035     }
2036     if (&Weights != NULL) {
2037       //
2038       // means that PrepareEval did found out that the surface was
2039       // locally polynomial but since the surface is constructed
2040       // with some weights we need to set the weight polynomial to constant
2041       //
2042
2043       for (ii = 1 ; ii <= d2p1 ; ii++) {
2044
2045         for (jj = 1 ; jj <= d1p1 ; jj++) {
2046           CacheWeights(ii,jj) = 0.0e0 ;
2047         }
2048       }
2049       CacheWeights(1,1) = 1.0e0 ;
2050     }
2051   }
2052 }
2053
2054 //=======================================================================
2055 //function : CacheD0
2056 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2057 //
2058 //=======================================================================
2059
2060 void  BSplSLib::CacheD0(const Standard_Real                  UParameter,
2061                         const Standard_Real                  VParameter,
2062                         const  Standard_Integer              UDegree,
2063                         const  Standard_Integer              VDegree,
2064                         const  Standard_Real                 UCacheParameter,
2065                         const  Standard_Real                 VCacheParameter,
2066                         const  Standard_Real                 USpanLenght,
2067                         const  Standard_Real                 VSpanLenght,
2068                         const  TColgp_Array2OfPnt&           PolesArray,
2069                         const  TColStd_Array2OfReal&         WeightsArray,
2070                         gp_Pnt&                              aPoint)
2071 {
2072   //
2073   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2074   // form
2075   // the SpanLenght     is the normalizing factor so that the polynomial is between
2076   // 0 and 1
2077   Standard_Integer
2078 //    ii,
2079   dimension,
2080   min_degree,
2081   max_degree  ;
2082
2083   Standard_Real
2084     new_parameter[2] ,
2085   inverse ;
2086
2087   Standard_Real *
2088     PArray = (Standard_Real *)
2089       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2090   Standard_Real *
2091     myPoint = (Standard_Real *) &aPoint  ;
2092   if (UDegree <= VDegree) {
2093     min_degree = UDegree ;
2094     max_degree = VDegree ;
2095     new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
2096     new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
2097     dimension = 3 * (UDegree + 1) ;
2098   }
2099   else {
2100     min_degree = VDegree ;
2101     max_degree = UDegree ;
2102     new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
2103     new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
2104     dimension = 3 * (VDegree + 1) ;
2105   }
2106   NCollection_LocalArray<Standard_Real> locpoles(dimension);
2107
2108   PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2109                        max_degree,
2110                        dimension,
2111                        max_degree*dimension,
2112                        PArray[0],
2113                        locpoles[0]) ;
2114
2115   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2116                        min_degree,
2117                        3,
2118                        (min_degree << 1) + min_degree,
2119                        locpoles[0],
2120                        myPoint[0]) ;
2121   if (&WeightsArray != NULL) {
2122     dimension = min_degree + 1 ;
2123     Standard_Real *
2124       WArray = (Standard_Real *)
2125         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2126     PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2127                          max_degree,
2128                          dimension,
2129                          max_degree*dimension,
2130                          WArray[0],
2131                          locpoles[0]) ;
2132
2133     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2134                          min_degree,
2135                          1,
2136                          min_degree,
2137                          locpoles[0],
2138                          inverse) ;
2139     inverse = 1.0e0/ inverse ;
2140
2141     myPoint[0] *= inverse ;
2142     myPoint[1] *= inverse ;
2143     myPoint[2] *= inverse ;
2144   }
2145 }
2146
2147 //=======================================================================
2148 //function : CacheD1
2149 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2150 //
2151 //=======================================================================
2152
2153 void  BSplSLib::CacheD1(const Standard_Real                  UParameter,
2154                         const Standard_Real                  VParameter,
2155                         const  Standard_Integer              UDegree,
2156                         const  Standard_Integer              VDegree,
2157                         const  Standard_Real                 UCacheParameter,
2158                         const  Standard_Real                 VCacheParameter,
2159                         const  Standard_Real                 USpanLenght,
2160                         const  Standard_Real                 VSpanLenght,
2161                         const  TColgp_Array2OfPnt&           PolesArray,
2162                         const  TColStd_Array2OfReal&         WeightsArray,
2163                         gp_Pnt&                              aPoint,
2164                         gp_Vec&                              aVecU,
2165                         gp_Vec&                              aVecV)
2166 {
2167   //
2168   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2169   // form
2170   // the SpanLenght     is the normalizing factor so that the polynomial is between
2171   // 0 and 1
2172   Standard_Integer
2173 //    ii,
2174 //  jj,
2175 //  kk,
2176   dimension,
2177   min_degree,
2178   max_degree  ;
2179
2180   Standard_Real
2181     inverse_min,
2182   inverse_max,
2183   new_parameter[2]  ;
2184
2185   Standard_Real *
2186     PArray = (Standard_Real *)
2187       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2188   Standard_Real local_poles_array[2][2][3],
2189   local_poles_and_weights_array[2][2][4],
2190   local_weights_array[2][2]  ;
2191   Standard_Real * my_vec_min,
2192   * my_vec_max,
2193   * my_point ;
2194   my_point = (Standard_Real *) &aPoint  ;
2195   //
2196   // initialize in case of rational evaluation
2197   // because RationalDerivative will use all
2198   // the coefficients
2199   //
2200   //
2201   if (&WeightsArray != NULL) {
2202
2203     local_poles_array            [0][0][0] = 0.0e0 ;
2204     local_poles_array            [0][0][1] = 0.0e0 ;
2205     local_poles_array            [0][0][2] = 0.0e0 ;
2206     local_weights_array          [0][0]    = 0.0e0 ;
2207     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2208     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2209     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2210     local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2211
2212     local_poles_array            [0][1][0] = 0.0e0 ;
2213     local_poles_array            [0][1][1] = 0.0e0 ;
2214     local_poles_array            [0][1][2] = 0.0e0 ;
2215     local_weights_array          [0][1]    = 0.0e0 ;
2216     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2217     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2218     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2219     local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2220
2221     local_poles_array            [1][0][0] = 0.0e0 ;
2222     local_poles_array            [1][0][1] = 0.0e0 ;
2223     local_poles_array            [1][0][2] = 0.0e0 ;
2224     local_weights_array          [1][0]    = 0.0e0 ;
2225     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2226     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2227     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2228     local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2229
2230     local_poles_array            [1][1][0] = 0.0e0 ;
2231     local_poles_array            [1][1][1] = 0.0e0 ;
2232     local_poles_array            [1][1][2] = 0.0e0 ;
2233     local_weights_array          [1][1]    = 0.0e0 ;
2234     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2235     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2236     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2237     local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2238   }
2239
2240   if (UDegree <= VDegree) {
2241     min_degree = UDegree ;
2242     max_degree = VDegree ;
2243     inverse_min = 1.0e0/USpanLenght ;
2244     inverse_max = 1.0e0/VSpanLenght ;
2245     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2246     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2247
2248     dimension = 3 * (UDegree + 1) ;
2249     my_vec_min = (Standard_Real *) &aVecU ;
2250     my_vec_max = (Standard_Real *) &aVecV ;
2251   }
2252   else {
2253     min_degree = VDegree ;
2254     max_degree = UDegree ;
2255     inverse_min = 1.0e0/VSpanLenght ;
2256     inverse_max = 1.0e0/USpanLenght ;
2257     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2258     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2259     dimension = 3 * (VDegree + 1) ;
2260     my_vec_min = (Standard_Real *) &aVecV ;
2261     my_vec_max = (Standard_Real *) &aVecU ;
2262   }
2263
2264   NCollection_LocalArray<Standard_Real> locpoles (2 * dimension);
2265
2266   PLib::EvalPolynomial(new_parameter[0],
2267                        1,
2268                        max_degree,
2269                        dimension,
2270                        PArray[0],
2271                        locpoles[0]) ;
2272
2273   PLib::EvalPolynomial(new_parameter[1],
2274                        1,
2275                        min_degree,
2276                        3,
2277                        locpoles[0],
2278                        local_poles_array[0][0][0]) ;
2279   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2280                        min_degree,
2281                        3,
2282                        (min_degree << 1) + min_degree,
2283                        locpoles[dimension],
2284                        local_poles_array[1][0][0]) ;
2285
2286   if (&WeightsArray != NULL) {
2287     dimension = min_degree + 1 ;
2288     Standard_Real *
2289       WArray = (Standard_Real *)
2290         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2291     PLib::EvalPolynomial(new_parameter[0],
2292                          1,
2293                          max_degree,
2294                          dimension,
2295                          WArray[0],
2296                          locpoles[0]) ;
2297
2298     PLib::EvalPolynomial(new_parameter[1],
2299                          1,
2300                          min_degree,
2301                          1,
2302                          locpoles[0],
2303                          local_weights_array[0][0]) ;
2304     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2305                          min_degree,
2306                          1,
2307                          min_degree,
2308                          locpoles[dimension],
2309                          local_weights_array[1][0]) ;
2310
2311     local_poles_and_weights_array[0][0][0] = local_poles_array  [0][0][0] ;
2312     local_poles_and_weights_array[0][0][1] = local_poles_array  [0][0][1] ;
2313     local_poles_and_weights_array[0][0][2] = local_poles_array  [0][0][2] ;
2314     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0]    ;
2315
2316     local_poles_and_weights_array[0][1][0] = local_poles_array  [0][1][0] ;
2317     local_poles_and_weights_array[0][1][1] = local_poles_array  [0][1][1] ;
2318     local_poles_and_weights_array[0][1][2] = local_poles_array  [0][1][2] ;
2319     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1]    ;
2320
2321     local_poles_and_weights_array[1][0][0] = local_poles_array  [1][0][0] ;
2322     local_poles_and_weights_array[1][0][1] = local_poles_array  [1][0][1] ;
2323     local_poles_and_weights_array[1][0][2] = local_poles_array  [1][0][2] ;
2324     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0]    ;
2325
2326     local_poles_and_weights_array[1][1][0] = local_poles_array  [1][1][0] ;
2327     local_poles_and_weights_array[1][1][1] = local_poles_array  [1][1][1] ;
2328     local_poles_and_weights_array[1][1][2] = local_poles_array  [1][1][2] ;
2329     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1]    ;
2330
2331     BSplSLib::RationalDerivative(1,
2332                                  1,
2333                                  1,
2334                                  1,
2335                                  local_poles_and_weights_array[0][0][0],
2336                                  local_poles_array[0][0][0]) ;
2337   }
2338
2339   my_point  [0] = local_poles_array              [0][0][0] ;
2340   my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2341   my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2342
2343   my_point  [1] = local_poles_array              [0][0][1] ;
2344   my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2345   my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2346
2347   my_point  [2] = local_poles_array              [0][0][2] ;
2348   my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2349   my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2350 }
2351
2352 //=======================================================================
2353 //function : CacheD2
2354 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2355 //
2356 //=======================================================================
2357
2358 void  BSplSLib::CacheD2(const Standard_Real                  UParameter,
2359                         const Standard_Real                  VParameter,
2360                         const  Standard_Integer              UDegree,
2361                         const  Standard_Integer              VDegree,
2362                         const  Standard_Real                 UCacheParameter,
2363                         const  Standard_Real                 VCacheParameter,
2364                         const  Standard_Real                 USpanLenght,
2365                         const  Standard_Real                 VSpanLenght,
2366                         const  TColgp_Array2OfPnt&           PolesArray,
2367                         const  TColStd_Array2OfReal&         WeightsArray,
2368                         gp_Pnt&                              aPoint,
2369                         gp_Vec&                              aVecU,
2370                         gp_Vec&                              aVecV,
2371                         gp_Vec&                              aVecUU,
2372                         gp_Vec&                              aVecUV,
2373                         gp_Vec&                              aVecVV)
2374 {
2375   //
2376   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2377   // form
2378   // the SpanLenght     is the normalizing factor so that the polynomial is between
2379   // 0 and 1
2380   Standard_Integer
2381     ii,
2382 //  jj,
2383   kk,
2384   index,
2385   dimension,
2386   min_degree,
2387   max_degree  ;
2388
2389   Standard_Real
2390     inverse_min,
2391   inverse_max,
2392   new_parameter[2]  ;
2393
2394   Standard_Real *
2395     PArray = (Standard_Real *)
2396       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2397   Standard_Real local_poles_array[3][3][3],
2398   local_poles_and_weights_array[3][3][4],
2399   local_weights_array[3][3]  ;
2400   Standard_Real * my_vec_min,
2401   * my_vec_max,
2402   * my_vec_min_min,
2403   * my_vec_max_max,
2404   * my_vec_min_max,
2405   * my_point ;
2406   my_point = (Standard_Real *) &aPoint  ;
2407
2408   //
2409   // initialize in case the min and max degree are less than 2
2410   //
2411   local_poles_array[0][0][0] = 0.0e0 ;
2412   local_poles_array[0][0][1] = 0.0e0 ;
2413   local_poles_array[0][0][2] = 0.0e0 ;
2414   local_poles_array[0][1][0] = 0.0e0 ;
2415   local_poles_array[0][1][1] = 0.0e0 ;
2416   local_poles_array[0][1][2] = 0.0e0 ;
2417   local_poles_array[0][2][0] = 0.0e0 ;
2418   local_poles_array[0][2][1] = 0.0e0 ;
2419   local_poles_array[0][2][2] = 0.0e0 ;
2420
2421   local_poles_array[1][0][0] = 0.0e0 ;
2422   local_poles_array[1][0][1] = 0.0e0 ;
2423   local_poles_array[1][0][2] = 0.0e0 ;
2424   local_poles_array[1][1][0] = 0.0e0 ;
2425   local_poles_array[1][1][1] = 0.0e0 ;
2426   local_poles_array[1][1][2] = 0.0e0 ;
2427   local_poles_array[1][2][0] = 0.0e0 ;
2428   local_poles_array[1][2][1] = 0.0e0 ;
2429   local_poles_array[1][2][2] = 0.0e0 ;
2430
2431   local_poles_array[2][0][0] = 0.0e0 ;
2432   local_poles_array[2][0][1] = 0.0e0 ;
2433   local_poles_array[2][0][2] = 0.0e0 ;
2434   local_poles_array[2][1][0] = 0.0e0 ;
2435   local_poles_array[2][1][1] = 0.0e0 ;
2436   local_poles_array[2][1][2] = 0.0e0 ;
2437   local_poles_array[2][2][0] = 0.0e0 ;
2438   local_poles_array[2][2][1] = 0.0e0 ;
2439   local_poles_array[2][2][2] = 0.0e0 ;
2440   //
2441   // initialize in case of rational evaluation
2442   // because RationalDerivative will use all
2443   // the coefficients
2444   //
2445   //
2446   if (&WeightsArray != NULL) {
2447
2448     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2449     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2450     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2451     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2452     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2453     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2454     local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2455     local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2456     local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2457
2458     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2459     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2460     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2461     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2462     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2463     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2464     local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2465     local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2466     local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2467
2468     local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2469     local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2470     local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2471     local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2472     local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2473     local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2474     local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2475     local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2476     local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2477
2478     local_poles_and_weights_array[0][0][3] =
2479       local_weights_array[0][0] = 0.0e0 ;
2480     local_poles_and_weights_array[0][1][3] =
2481       local_weights_array[0][1] = 0.0e0 ;
2482     local_poles_and_weights_array[0][2][3] =
2483       local_weights_array[0][2] = 0.0e0 ;
2484     local_poles_and_weights_array[1][0][3] =
2485       local_weights_array[1][0] = 0.0e0 ;
2486     local_poles_and_weights_array[1][1][3] =
2487       local_weights_array[1][1] = 0.0e0 ;
2488     local_poles_and_weights_array[1][2][3] =
2489       local_weights_array[1][2] = 0.0e0 ;
2490     local_poles_and_weights_array[2][0][3] =
2491       local_weights_array[2][0] = 0.0e0 ;
2492     local_poles_and_weights_array[2][1][3] =
2493       local_weights_array[2][1] = 0.0e0 ;
2494     local_poles_and_weights_array[2][2][3] =
2495       local_weights_array[2][2] = 0.0e0 ;
2496   }
2497
2498   if (UDegree <= VDegree) {
2499     min_degree = UDegree ;
2500     max_degree = VDegree ;
2501     inverse_min = 1.0e0/USpanLenght ;
2502     inverse_max = 1.0e0/VSpanLenght ;
2503     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2504     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2505
2506     dimension = 3 * (UDegree + 1) ;
2507     my_vec_min     = (Standard_Real *) &aVecU ;
2508     my_vec_max     = (Standard_Real *) &aVecV ;
2509     my_vec_min_min = (Standard_Real *) &aVecUU ;
2510     my_vec_min_max = (Standard_Real *) &aVecUV ;
2511     my_vec_max_max = (Standard_Real *) &aVecVV ;
2512   }
2513   else {
2514     min_degree = VDegree ;
2515     max_degree = UDegree ;
2516     inverse_min = 1.0e0/VSpanLenght ;
2517     inverse_max = 1.0e0/USpanLenght ;
2518     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2519     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2520     dimension = 3 * (VDegree + 1) ;
2521     my_vec_min     = (Standard_Real *) &aVecV ;
2522     my_vec_max     = (Standard_Real *) &aVecU ;
2523     my_vec_min_min = (Standard_Real *) &aVecVV ;
2524     my_vec_min_max = (Standard_Real *) &aVecUV ;
2525     my_vec_max_max = (Standard_Real *) &aVecUU ;
2526   }
2527
2528   NCollection_LocalArray<Standard_Real> locpoles (3 * dimension);
2529
2530   //
2531   // initialize in case min or max degree are less than 2
2532   //
2533   Standard_Integer MinIndMax = 2;
2534   if ( max_degree < 2) MinIndMax = max_degree;
2535   Standard_Integer MinIndMin = 2;
2536   if ( min_degree < 2) MinIndMin = min_degree;
2537
2538   index =  MinIndMax * dimension ;
2539
2540   for (ii = MinIndMax ; ii <  3 ; ii++) {
2541
2542     for (kk = 0 ; kk < dimension ; kk++) {
2543       locpoles[index] = 0.0e0 ;
2544       index += 1 ;
2545     }
2546   }
2547
2548   PLib::EvalPolynomial(new_parameter[0],
2549                        MinIndMax,
2550                        max_degree,
2551                        dimension,
2552                        PArray[0],
2553                        locpoles[0]) ;
2554
2555   PLib::EvalPolynomial(new_parameter[1],
2556                        MinIndMin,
2557                        min_degree,
2558                        3,
2559                        locpoles[0],
2560                        local_poles_array[0][0][0]) ;
2561   PLib::EvalPolynomial(new_parameter[1],
2562                        1,
2563                        min_degree,
2564                        3,
2565                        locpoles[dimension],
2566                        local_poles_array[1][0][0]) ;
2567
2568   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2569                        min_degree,
2570                        3,
2571                        (min_degree << 1) + min_degree,
2572                        locpoles[dimension + dimension],
2573                        local_poles_array[2][0][0]) ;
2574
2575   if (&WeightsArray != NULL) {
2576     dimension = min_degree + 1 ;
2577     Standard_Real *
2578       WArray = (Standard_Real *)
2579         &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2580     PLib::EvalPolynomial(new_parameter[0],
2581                          MinIndMax,
2582                          max_degree,
2583                          dimension,
2584                          WArray[0],
2585                          locpoles[0]) ;
2586
2587     PLib::EvalPolynomial(new_parameter[1],
2588                          MinIndMin,
2589                          min_degree,
2590                          1,
2591                          locpoles[0],
2592                          local_weights_array[0][0]) ;
2593     PLib::EvalPolynomial(new_parameter[1],
2594                          1,
2595                          min_degree,
2596                          1,
2597                          locpoles[dimension],
2598                          local_weights_array[1][0]) ;
2599     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2600                          min_degree,
2601                          1,
2602                          min_degree,
2603                          locpoles[dimension + dimension],
2604                          local_weights_array[2][0]) ;
2605
2606
2607     local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2608     local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2609     local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2610     local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2611     local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2612     local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2613     local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2614     local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2615     local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2616
2617     local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2618     local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2619     local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2620     local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2621     local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2622     local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2623     local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2624     local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2625     local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2626
2627     local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2628     local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2629     local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2630     local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2631     local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2632     local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2633     local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2634     local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2635     local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2636
2637
2638     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2639     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2640     local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2641     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2642     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2643     local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2644     local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2645     local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2646     local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2647
2648     BSplSLib::RationalDerivative(2,
2649                                  2,
2650                                  2,
2651                                  2,
2652                                  local_poles_and_weights_array[0][0][0],
2653                                  local_poles_array[0][0][0]) ;
2654   }
2655
2656
2657   Standard_Real minmin = inverse_min * inverse_min;
2658   Standard_Real minmax = inverse_min * inverse_max;
2659   Standard_Real maxmax = inverse_max * inverse_max;
2660
2661   my_point      [0] = local_poles_array              [0][0][0] ;
2662   my_vec_min    [0] = inverse_min * local_poles_array[0][1][0] ;
2663   my_vec_max    [0] = inverse_max * local_poles_array[1][0][0] ;
2664   my_vec_min_min[0] = minmin * local_poles_array     [0][2][0] ;
2665   my_vec_min_max[0] = minmax * local_poles_array     [1][1][0] ;
2666   my_vec_max_max[0] = maxmax * local_poles_array     [2][0][0] ;
2667
2668   my_point      [1] = local_poles_array              [0][0][1] ;
2669   my_vec_min    [1] = inverse_min * local_poles_array[0][1][1] ;
2670   my_vec_max    [1] = inverse_max * local_poles_array[1][0][1] ;
2671   my_vec_min_min[1] = minmin * local_poles_array     [0][2][1] ;
2672   my_vec_min_max[1] = minmax * local_poles_array     [1][1][1] ;
2673   my_vec_max_max[1] = maxmax * local_poles_array     [2][0][1] ;
2674
2675   my_point      [2] = local_poles_array              [0][0][2] ;
2676   my_vec_min    [2] = inverse_min * local_poles_array[0][1][2] ;
2677   my_vec_max    [2] = inverse_max * local_poles_array[1][0][2] ;
2678   my_vec_min_min[2] = minmin * local_poles_array     [0][2][2] ;
2679   my_vec_min_max[2] = minmax * local_poles_array     [1][1][2] ;
2680   my_vec_max_max[2] = maxmax * local_poles_array     [2][0][2] ;
2681 }
2682
2683 //=======================================================================
2684 //function : MovePoint
2685 //purpose  : Find the new poles which allows  an old point (with a
2686 //           given  u and v as parameters) to reach a new position
2687 //=======================================================================
2688
2689 void BSplSLib::MovePoint (const Standard_Real            U,
2690                           const Standard_Real            V,
2691                           const gp_Vec&                  Displ,
2692                           const Standard_Integer         UIndex1,
2693                           const Standard_Integer         UIndex2,
2694                           const Standard_Integer         VIndex1,
2695                           const Standard_Integer         VIndex2,
2696                           const Standard_Integer         UDegree,
2697                           const Standard_Integer         VDegree,
2698                           const Standard_Boolean         Rational,
2699                           const TColgp_Array2OfPnt&      Poles,
2700                           const TColStd_Array2OfReal&    Weights,
2701                           const TColStd_Array1OfReal&    UFlatKnots,
2702                           const TColStd_Array1OfReal&    VFlatKnots,
2703                           Standard_Integer&              UFirstIndex,
2704                           Standard_Integer&              ULastIndex,
2705                           Standard_Integer&              VFirstIndex,
2706                           Standard_Integer&              VLastIndex,
2707                           TColgp_Array2OfPnt&            NewPoles)
2708 {
2709   // calculate the UBSplineBasis in the parameter U
2710   Standard_Integer UFirstNonZeroBsplineIndex;
2711   math_Matrix UBSplineBasis(1, 1,
2712                             1, UDegree+1);
2713   Standard_Integer ErrorCod1 =  BSplCLib::EvalBsplineBasis(1,
2714                                                            0,
2715                                                            UDegree+1,
2716                                                            UFlatKnots,
2717                                                            U,
2718                                                            UFirstNonZeroBsplineIndex,
2719                                                            UBSplineBasis);
2720   // calculate the VBSplineBasis in the parameter V
2721   Standard_Integer VFirstNonZeroBsplineIndex;
2722   math_Matrix VBSplineBasis(1, 1,
2723                             1, VDegree+1);
2724   Standard_Integer ErrorCod2 =  BSplCLib::EvalBsplineBasis(1,
2725                                                            0,
2726                                                            VDegree+1,
2727                                                            VFlatKnots,
2728                                                            V,
2729                                                            VFirstNonZeroBsplineIndex,
2730                                                            VBSplineBasis);
2731   if (ErrorCod1 || ErrorCod2) {
2732     UFirstIndex = 0;
2733     ULastIndex = 0;
2734     VFirstIndex = 0;
2735     VLastIndex = 0;
2736     return;
2737   }
2738
2739   // find the span which is predominant for parameter U
2740   UFirstIndex = UFirstNonZeroBsplineIndex;
2741   ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2742   if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2743   if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2744
2745   Standard_Real maxValue = 0.0;
2746   Standard_Integer i, ukk1=0, ukk2;
2747
2748   for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2749     if (UBSplineBasis(1,i) > maxValue) {
2750       ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2751       maxValue = UBSplineBasis(1,i);
2752     }
2753   }
2754
2755   // find a ukk2 if symetriy
2756   ukk2 = ukk1;
2757   i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2758   if ((ukk1+1) <= ULastIndex) {
2759     if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2760       ukk2 = ukk1+1;
2761     }
2762   }
2763
2764   // find the span which is predominant for parameter V
2765   VFirstIndex = VFirstNonZeroBsplineIndex;
2766   VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2767
2768   if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2769   if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2770
2771   maxValue = 0.0;
2772   Standard_Integer j, vkk1=0, vkk2;
2773
2774   for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2775     if (VBSplineBasis(1,j) > maxValue) {
2776       vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2777       maxValue = VBSplineBasis(1,j);
2778     }
2779   }
2780
2781   // find a vkk2 if symetriy
2782   vkk2 = vkk1;
2783   j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2784   if ((vkk1+1) <= VLastIndex) {
2785     if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2786       vkk2 = vkk1+1;
2787     }
2788   }
2789
2790   // compute the vector of displacement
2791   Standard_Real D1 = 0.0;
2792   Standard_Real D2 = 0.0;
2793   Standard_Real hN, Coef, DvalU, DvalV;
2794
2795   Standard_Integer ii, jj;
2796
2797   for (i = 1; i <= UDegree+1; i++) {
2798     ii = i + UFirstNonZeroBsplineIndex - 1;
2799     if (ii < ukk1) {
2800       DvalU = ukk1-ii;
2801     }
2802     else if (ii > ukk2) {
2803       DvalU = ii - ukk2;
2804     }
2805     else {
2806       DvalU = 0.0;
2807     }
2808
2809     for (j = 1; j <= VDegree+1; j++) {
2810       jj = j + VFirstNonZeroBsplineIndex - 1;
2811       if (Rational) {
2812         hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2813         D2 += hN;
2814       }
2815       else {
2816         hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2817       }
2818       if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2819         if (jj < vkk1) {
2820           DvalV = vkk1-jj;
2821         }
2822         else if (jj > vkk2) {
2823           DvalV = jj - vkk2;
2824         }
2825         else {
2826           DvalV = 0.0;
2827         }
2828         D1 += 1./(DvalU + DvalV + 1.) * hN;
2829       }
2830     }
2831   }
2832
2833   if (Rational) {
2834     Coef = D2/D1;
2835   }
2836   else {
2837     Coef = 1./D1;
2838   }
2839
2840   // compute the new poles
2841
2842   for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2843     if (i < ukk1) {
2844       DvalU = ukk1-i;
2845     }
2846     else if (i > ukk2) {
2847       DvalU = i - ukk2;
2848     }
2849     else {
2850       DvalU = 0.0;
2851     }
2852
2853     for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2854       if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2855         if (j < vkk1) {
2856           DvalV = vkk1-j;
2857         }
2858         else if (j > vkk2) {
2859           DvalV = j - vkk2;
2860         }
2861         else {
2862           DvalV = 0.0;
2863         }
2864         NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2865       }
2866       else {
2867         NewPoles(i,j) = Poles(i,j);
2868       }
2869     }
2870   }
2871 }
2872
2873 //=======================================================================
2874 // function : Resolution
2875 // purpose  : this computes an estimate for the maximum of the
2876 // partial derivatives both in U and in V
2877 //
2878 //
2879 // The calculation resembles at the calculation of curves with
2880 // additional index for the control point. Let Si,j be the
2881 // control points for ls surface  and  Di,j  the weights.
2882 // The checking of upper bounds for the partial derivatives
2883 // will be omitted and Su is the next upper bound in the polynomial case :
2884 //
2885 //
2886 //
2887 //                        |  Si,j - Si-1,j  |
2888 //          d *   Max     |  -------------  |
2889 //                i = 2,n |     ti+d - ti   |
2890 //                i=1.m
2891 //
2892 //
2893 // and in the rational case :
2894 //
2895 //
2896 //
2897 //                         Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2898 //   Max   Max       d  *  -----------------------------------------------
2899 // k=1,n  i dans Rj                   ti+d  - ti
2900 // j=1,m
2901 //  ----------------------------------------------------------------------
2902 //
2903 //               Min    Di,j
2904 //              i=1,n
2905 //              j=1,m
2906 //
2907 //
2908 //
2909 // with Rj = {j-d, ....,  j+d+d+1}.
2910 //
2911 //
2912 //=======================================================================
2913
2914 void BSplSLib::Resolution(const TColgp_Array2OfPnt&      Poles,
2915                           const TColStd_Array2OfReal&    Weights,
2916                           const TColStd_Array1OfReal&    UKnots,
2917                           const TColStd_Array1OfReal&    VKnots,
2918                           const TColStd_Array1OfInteger& UMults,
2919                           const TColStd_Array1OfInteger& VMults,
2920                           const Standard_Integer         UDegree,
2921                           const Standard_Integer         VDegree,
2922                           const Standard_Boolean         URational,
2923                           const Standard_Boolean         VRational,
2924                           const Standard_Boolean         UPeriodic,
2925                           const Standard_Boolean         VPeriodic,
2926                           const Standard_Real            Tolerance3D,
2927                           Standard_Real&                 UTolerance,
2928                           Standard_Real&                 VTolerance)
2929 {
2930   Standard_Real Wij,Wmj,Wji,Wjm;
2931   Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2932   Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2933   Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2934   Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2935
2936   max_derivative[0] = max_derivative[1] = 0.0e0 ;
2937
2938   Standard_Integer PRowLength, PColLength;
2939   Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2940   Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2941   Standard_Integer num_poles[2],num_flat_knots[2];
2942
2943   num_flat_knots[0] =
2944     BSplCLib::KnotSequenceLength(UMults,
2945                                  UDegree,
2946                                  UPeriodic) ;
2947   num_flat_knots[1] =
2948     BSplCLib::KnotSequenceLength(VMults,
2949                                  VDegree,
2950                                  VPeriodic) ;
2951   TColStd_Array1OfReal  flat_knots_in_u(1,num_flat_knots[0]) ;
2952   TColStd_Array1OfReal  flat_knots_in_v(1,num_flat_knots[1]) ;
2953   BSplCLib::KnotSequence(UKnots,
2954                          UMults,
2955                          UDegree,
2956                          UPeriodic,
2957                          flat_knots_in_u) ;
2958   BSplCLib::KnotSequence(VKnots,
2959                          VMults,
2960                          VDegree,
2961                          VPeriodic,
2962                          flat_knots_in_v) ;
2963   PRowLength = Poles.RowLength();
2964   PColLength = Poles.ColLength();
2965   if (URational || VRational) {
2966     Standard_Integer Wsize = PRowLength * PColLength;
2967     const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2968     min_weights = WG[0];
2969
2970     for (ii = 1 ; ii < Wsize ; ii++) {
2971       min = WG[ii];
2972       if (min_weights > min) min_weights = min;
2973     }
2974   }
2975   Standard_Integer UD1 = UDegree + 1;
2976   Standard_Integer VD1 = VDegree + 1;
2977   num_poles[0] = num_flat_knots[0] - UD1;
2978   num_poles[1] = num_flat_knots[1] - VD1;
2979   poles_length[0] = PColLength;
2980   poles_length[1] = PRowLength;
2981   if(URational) {
2982     Standard_Integer UD2 = UDegree << 1;
2983     Standard_Integer VD2 = VDegree << 1;
2984
2985     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2986       ii_index = (ii - 1) % poles_length[0] + 1 ;
2987       ii_minus = (ii - 2) % poles_length[0] + 1 ;
2988       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2989       inverse = 1.0e0 / inverse ;
2990       lower[0] = ii - UD1;
2991       if (lower[0] < 1) lower[0] = 1;
2992       upper[0] = ii + UD2 + 1;
2993       if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2994
2995       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2996         jj_index = (jj - 1) % poles_length[1] + 1 ;
2997         lower[1] = jj - VD1;
2998         if (lower[1] < 1) lower[1] = 1;
2999         upper[1] = jj + VD2 + 1;
3000         if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
3001         const gp_Pnt& Pij = Poles  .Value(ii_index,jj_index);
3002         Wij               = Weights.Value(ii_index,jj_index);
3003         const gp_Pnt& Pmj = Poles  .Value(ii_minus,jj_index);
3004         Wmj               = Weights.Value(ii_minus,jj_index);
3005         Xij = Pij.X();
3006         Yij = Pij.Y();
3007         Zij = Pij.Z();
3008         Xmj = Pmj.X();
3009         Ymj = Pmj.Y();
3010         Zmj = Pmj.Z();
3011
3012         for (pp = lower[0] ; pp <= upper[0] ; pp++) {
3013           pp_index = (pp - 1) % poles_length[0] + 1 ;
3014
3015           for (qq = lower[1] ; qq <= upper[1] ; qq++) {
3016             value = 0.0e0 ;
3017             qq_index = (qq - 1) % poles_length[1] + 1 ;
3018             const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
3019             Xpq = Ppq.X();
3020             Ypq = Ppq.Y();
3021             Zpq = Ppq.Z();
3022             factor  = (Xpq - Xij) * Wij;
3023             factor -= (Xpq - Xmj) * Wmj;
3024             if (factor < 0) factor = - factor;
3025             value  += factor ;
3026             factor  = (Ypq - Yij) * Wij;
3027             factor -= (Ypq - Ymj) * Wmj;
3028             if (factor < 0) factor = - factor;
3029             value  += factor ;
3030             factor  = (Zpq - Zij) * Wij;
3031             factor -= (Zpq - Zmj) * Wmj;
3032             if (factor < 0) factor = - factor;
3033             value  += factor ;
3034             value *= inverse ;
3035             if (max_derivative[0] < value) max_derivative[0] = value ;
3036           }
3037         }
3038       }
3039     }
3040     max_derivative[0] /= min_weights ;
3041   }
3042   else {
3043
3044     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3045       ii_index = (ii - 1) % poles_length[0] + 1 ;
3046       ii_minus = (ii - 2) % poles_length[0] + 1 ;
3047       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3048       inverse = 1.0e0 / inverse ;
3049
3050       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3051         jj_index = (jj - 1) % poles_length[1] + 1 ;
3052         value = 0.0e0 ;
3053         const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
3054         const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
3055         factor = Pij.X() - Pmj.X();
3056         if (factor < 0) factor = - factor;
3057         value += factor;
3058         factor = Pij.Y() - Pmj.Y();
3059         if (factor < 0) factor = - factor;
3060         value += factor;
3061         factor = Pij.Z() - Pmj.Z();
3062         if (factor < 0) factor = - factor;
3063         value += factor;
3064         value *= inverse ;
3065         if (max_derivative[0] < value) max_derivative[0] = value ;
3066       }
3067     }
3068   }
3069   max_derivative[0] *= UDegree ;
3070   if(VRational) {
3071     Standard_Integer UD2 = UDegree << 1;
3072     Standard_Integer VD2 = VDegree << 1;
3073
3074     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3075       ii_index = (ii - 1) % poles_length[1] + 1 ;
3076       ii_minus = (ii - 2) % poles_length[1] + 1 ;
3077       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3078       inverse = 1.0e0 / inverse ;
3079       lower[0] = ii - VD1;
3080       if (lower[0] < 1) lower[0] = 1;
3081       upper[0] = ii + VD2 + 1;
3082       if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
3083
3084       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3085         jj_index = (jj - 1) % poles_length[0] + 1 ;
3086         lower[1] = jj - UD1;
3087         if (lower[1] < 1) lower[1] = 1;
3088         upper[1] = jj + UD2 + 1;
3089         if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
3090         const gp_Pnt& Pji = Poles  .Value(jj_index,ii_index);
3091         Wji               = Weights.Value(jj_index,ii_index);
3092         const gp_Pnt& Pjm = Poles  .Value(jj_index,ii_minus);
3093         Wjm               = Weights.Value(jj_index,ii_minus);
3094         Xji = Pji.X();
3095         Yji = Pji.Y();
3096         Zji = Pji.Z();
3097         Xjm = Pjm.X();
3098         Yjm = Pjm.Y();
3099         Zjm = Pjm.Z();
3100
3101         for (pp = lower[1] ; pp <= upper[1] ; pp++) {
3102           pp_index = (pp - 1) % poles_length[1] + 1 ;
3103
3104           for (qq = lower[0] ; qq <= upper[0] ; qq++) {
3105             value = 0.0e0 ;
3106             qq_index = (qq - 1) % poles_length[0] + 1 ;
3107             const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
3108             Xqp = Pqp.X();
3109             Yqp = Pqp.Y();
3110             Zqp = Pqp.Z();
3111             factor  = (Xqp - Xji) * Wji;
3112             factor -= (Xqp - Xjm) * Wjm;
3113             if (factor < 0) factor = - factor;
3114             value += factor ;
3115             factor  = (Yqp - Yji) * Wji;
3116             factor -= (Yqp - Yjm) * Wjm;
3117             if (factor < 0) factor = - factor;
3118             value += factor ;
3119             factor  = (Zqp - Zji) * Wji;
3120             factor -= (Zqp - Zjm) * Wjm;
3121             if (factor < 0) factor = - factor;
3122             value += factor ;
3123             value *= inverse ;
3124             if (max_derivative[1] < value) max_derivative[1] = value ;
3125           }
3126         }
3127       }
3128     }
3129     max_derivative[1] /= min_weights ;
3130   }
3131   else {
3132
3133     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3134       ii_index = (ii - 1) % poles_length[1] + 1 ;
3135       ii_minus = (ii - 2) % poles_length[1] + 1 ;
3136       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3137       inverse = 1.0e0 / inverse ;
3138
3139       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3140         jj_index = (jj - 1) % poles_length[0] + 1 ;
3141         value = 0.0e0 ;
3142         const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3143         const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3144         factor = Pji.X() - Pjm.X() ;
3145         if (factor < 0) factor = - factor;
3146         value += factor;
3147         factor = Pji.Y() - Pjm.Y() ;
3148         if (factor < 0) factor = - factor;
3149         value += factor;
3150         factor = Pji.Z() - Pjm.Z() ;
3151         if (factor < 0) factor = - factor;
3152         value += factor;
3153         value *= inverse ;
3154         if (max_derivative[1] < value) max_derivative[1] = value ;
3155       }
3156     }
3157   }
3158   max_derivative[1] *= VDegree ;
3159   max_derivative[0] *= M_SQRT2 ;
3160   max_derivative[1] *= M_SQRT2 ;
3161   if(max_derivative[0] && max_derivative[1]) {
3162     UTolerance = Tolerance3D / max_derivative[0] ;
3163     VTolerance = Tolerance3D / max_derivative[1] ;
3164   }
3165   else {
3166     UTolerance=VTolerance=0.0;
3167 #ifdef DEB
3168     cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3169 #endif
3170   }
3171 }
3172
3173 //=======================================================================
3174 //function : Interpolate
3175 //purpose  :
3176 //=======================================================================
3177
3178 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3179                            const Standard_Integer VDegree,
3180                            const TColStd_Array1OfReal& UFlatKnots,
3181                            const TColStd_Array1OfReal& VFlatKnots,
3182                            const TColStd_Array1OfReal& UParameters,
3183                            const TColStd_Array1OfReal& VParameters,
3184                            TColgp_Array2OfPnt& Poles,
3185                            TColStd_Array2OfReal& Weights,
3186                            Standard_Integer& InversionProblem)
3187 {
3188   Standard_Integer ii, jj, ll, kk, dimension;
3189   Standard_Integer ULength = UParameters.Length();
3190   Standard_Integer VLength = VParameters.Length();
3191   Standard_Real * poles_array;
3192
3193   // extraction of iso u
3194   dimension = 4*ULength;
3195   TColStd_Array2OfReal Points(1, VLength,
3196                               1, dimension);
3197
3198   Handle(TColStd_HArray1OfInteger) ContactOrder =
3199     new (TColStd_HArray1OfInteger)(1, VLength);
3200   ContactOrder->Init(0);
3201
3202   for (ii=1; ii <= VLength; ii++) {
3203
3204     for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3205       Points(ii,ll)   = Poles(jj, ii).X();
3206       Points(ii,ll+1) = Poles(jj, ii).Y();
3207       Points(ii,ll+2) = Poles(jj, ii).Z();
3208       Points(ii,ll+3) = Weights(jj,ii)  ;
3209     }
3210   }
3211
3212   // interpolation of iso u
3213   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3214   BSplCLib::Interpolate(VDegree,
3215                         VFlatKnots,
3216                         VParameters,
3217                         ContactOrder->Array1(),
3218                         dimension,
3219                         poles_array[0],
3220                         InversionProblem) ;
3221   if (InversionProblem != 0) return;
3222
3223   // extraction of iso v
3224
3225   dimension = VLength*4;
3226   TColStd_Array2OfReal IsoPoles(1, ULength,
3227                                 1, dimension);
3228
3229   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3230   ContactOrder->Init(0);
3231   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3232
3233   for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3234
3235     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3236       IsoPoles (ii,ll)   = Points(jj, kk);
3237       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3238       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3239       IsoPoles (ii,ll+3) = Points(jj, kk+3);
3240     }
3241   }
3242   // interpolation of iso v
3243   BSplCLib::Interpolate(UDegree,
3244                         UFlatKnots,
3245                         UParameters,
3246                         ContactOrder->Array1(),
3247                         dimension,
3248                         poles_array[0],
3249                         InversionProblem);
3250
3251   // return results
3252
3253   for (ii=1; ii <= ULength; ii++) {
3254
3255     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3256       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3257       Poles.SetValue(ii, jj, Pnt);
3258       Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3259     }
3260   }
3261 }
3262
3263 //=======================================================================
3264 //function : Interpolate
3265 //purpose  :
3266 //=======================================================================
3267
3268 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3269                            const Standard_Integer VDegree,
3270                            const TColStd_Array1OfReal& UFlatKnots,
3271                            const TColStd_Array1OfReal& VFlatKnots,
3272                            const TColStd_Array1OfReal& UParameters,
3273                            const TColStd_Array1OfReal& VParameters,
3274                            TColgp_Array2OfPnt& Poles,
3275                            Standard_Integer& InversionProblem)
3276 {
3277   Standard_Integer ii, jj, ll, kk, dimension;
3278   Standard_Integer ULength = UParameters.Length();
3279   Standard_Integer VLength = VParameters.Length();
3280   Standard_Real * poles_array;
3281
3282   // extraction of iso u
3283   dimension = 3*ULength;
3284   TColStd_Array2OfReal Points(1, VLength,
3285                               1, dimension);
3286
3287   Handle(TColStd_HArray1OfInteger) ContactOrder =
3288     new (TColStd_HArray1OfInteger)(1, VLength);
3289   ContactOrder->Init(0);
3290
3291   for (ii=1; ii <= VLength; ii++) {
3292
3293     for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3294       Points(ii,ll)   = Poles(jj, ii).X();
3295       Points(ii,ll+1) = Poles(jj, ii).Y();
3296       Points(ii,ll+2) = Poles(jj, ii).Z();
3297     }
3298   }
3299
3300   // interpolation of iso u
3301   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3302   BSplCLib::Interpolate(VDegree,
3303                         VFlatKnots,
3304                         VParameters,
3305                         ContactOrder->Array1(),
3306                         dimension,
3307                         poles_array[0],
3308                         InversionProblem) ;
3309   if (InversionProblem != 0) return;
3310
3311   // extraction of iso v
3312
3313   dimension = VLength*3;
3314   TColStd_Array2OfReal IsoPoles(1, ULength,
3315                                 1, dimension);
3316
3317   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3318   ContactOrder->Init(0);
3319   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3320
3321   for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3322
3323     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3324       IsoPoles (ii,ll)   = Points(jj, kk);
3325       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3326       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3327     }
3328   }
3329   // interpolation of iso v
3330   BSplCLib::Interpolate(UDegree,
3331                         UFlatKnots,
3332                         UParameters,
3333                         ContactOrder->Array1(),
3334                         dimension,
3335                         poles_array[0],
3336                         InversionProblem);
3337
3338   // return results
3339
3340   for (ii=1; ii <= ULength; ii++) {
3341
3342     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3343       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3344       Poles.SetValue(ii, jj, Pnt);
3345     }
3346   }
3347 }
3348
3349 //=======================================================================
3350 //function : FunctionMultiply
3351 //purpose  :
3352 //=======================================================================
3353
3354 void BSplSLib::FunctionMultiply
3355 (const BSplSLib_EvaluatorFunction&           Function,
3356  const Standard_Integer                      UBSplineDegree,
3357  const Standard_Integer                      VBSplineDegree,
3358  const TColStd_Array1OfReal&                 UBSplineKnots,
3359  const TColStd_Array1OfReal&                 VBSplineKnots,
3360  const TColStd_Array1OfInteger &             UMults,
3361  const TColStd_Array1OfInteger &             VMults,
3362  const TColgp_Array2OfPnt&                   Poles,
3363  const TColStd_Array2OfReal&                 Weights,
3364  const TColStd_Array1OfReal&                 UFlatKnots,
3365  const TColStd_Array1OfReal&                 VFlatKnots,
3366  const Standard_Integer                      UNewDegree,
3367  const Standard_Integer                      VNewDegree,
3368  TColgp_Array2OfPnt&                         NewNumerator,
3369  TColStd_Array2OfReal&                       NewDenominator,
3370  Standard_Integer&                           Status)
3371 {
3372   Standard_Integer num_uparameters,
3373 //  ii,jj,kk,
3374   ii,jj,
3375   error_code,
3376   num_vparameters ;
3377   Standard_Real    result ;
3378
3379   num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3380   num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3381   TColStd_Array1OfReal  UParameters(1,num_uparameters) ;
3382   TColStd_Array1OfReal  VParameters(1,num_vparameters) ;
3383
3384   if ((NewNumerator.ColLength() == num_uparameters) &&
3385       (NewNumerator.RowLength() == num_vparameters) &&
3386       (NewDenominator.ColLength() == num_uparameters) &&
3387       (NewDenominator.RowLength() == num_vparameters)) {
3388
3389
3390     BSplCLib::BuildSchoenbergPoints(UNewDegree,
3391                                     UFlatKnots,
3392                                     UParameters) ;
3393
3394     BSplCLib::BuildSchoenbergPoints(VNewDegree,
3395                                     VFlatKnots,
3396                                     VParameters) ;
3397
3398     for (ii = 1 ; ii <= num_uparameters ; ii++) {
3399
3400       for (jj = 1 ; jj <= num_vparameters ; jj++) {
3401         HomogeneousD0(UParameters(ii),
3402                       VParameters(jj),
3403                       0,
3404                       0,
3405                       Poles,
3406                       Weights,
3407                       UBSplineKnots,
3408                       VBSplineKnots,
3409                       UMults,
3410                       VMults,
3411                       UBSplineDegree,
3412                       VBSplineDegree,
3413                       Standard_True,
3414                       Standard_True,
3415                       Standard_False,
3416                       Standard_False,
3417                       NewDenominator(ii,jj),
3418                       NewNumerator(ii,jj)) ;
3419
3420         Function.Evaluate (0,
3421                  UParameters(ii),
3422                  VParameters(jj),
3423                  result,
3424                  error_code) ;
3425         if (error_code) {
3426           Standard_ConstructionError::Raise();
3427         }
3428         gp_Pnt& P = NewNumerator(ii,jj);
3429         P.SetX(P.X() * result);
3430         P.SetY(P.Y() * result);
3431         P.SetZ(P.Z() * result);
3432         NewDenominator(ii,jj) *= result ;
3433       }
3434     }
3435     Interpolate(UNewDegree,
3436                 VNewDegree,
3437                 UFlatKnots,
3438                 VFlatKnots,
3439                 UParameters,
3440                 VParameters,
3441                 NewNumerator,
3442                 NewDenominator,
3443                 Status) ;
3444   }
3445   else {
3446     Standard_ConstructionError::Raise();
3447   }
3448 }