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