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     const TColStd_Array2OfReal& refWeights = *WeightsArray;
2216     Standard_Real *
2217       WArray = (Standard_Real *)
2218         &refWeights(refWeights.LowerCol(),refWeights.LowerRow()) ;
2219     PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2220                          max_degree,
2221                          dimension,
2222                          max_degree*dimension,
2223                          WArray[0],
2224                          locpoles[0]) ;
2225
2226     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2227                          min_degree,
2228                          1,
2229                          min_degree,
2230                          locpoles[0],
2231                          inverse) ;
2232     inverse = 1.0e0/ inverse ;
2233
2234     myPoint[0] *= inverse ;
2235     myPoint[1] *= inverse ;
2236     myPoint[2] *= inverse ;
2237   }
2238 }
2239
2240 //=======================================================================
2241 //function : CacheD1
2242 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2243 //
2244 //=======================================================================
2245
2246 void  BSplSLib::CacheD1(const Standard_Real                  UParameter,
2247                         const Standard_Real                  VParameter,
2248                         const  Standard_Integer              UDegree,
2249                         const  Standard_Integer              VDegree,
2250                         const  Standard_Real                 UCacheParameter,
2251                         const  Standard_Real                 VCacheParameter,
2252                         const  Standard_Real                 USpanLenght,
2253                         const  Standard_Real                 VSpanLenght,
2254                         const  TColgp_Array2OfPnt&           PolesArray,
2255                         const  TColStd_Array2OfReal*         WeightsArray,
2256                         gp_Pnt&                              aPoint,
2257                         gp_Vec&                              aVecU,
2258                         gp_Vec&                              aVecV)
2259 {
2260   //
2261   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2262   // form
2263   // the SpanLenght     is the normalizing factor so that the polynomial is between
2264   // 0 and 1
2265   Standard_Integer
2266 //    ii,
2267 //  jj,
2268 //  kk,
2269   dimension,
2270   min_degree,
2271   max_degree  ;
2272
2273   Standard_Real
2274     inverse_min,
2275   inverse_max,
2276   new_parameter[2]  ;
2277
2278   Standard_Real *
2279     PArray = (Standard_Real *)
2280       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2281   Standard_Real local_poles_array[2][2][3],
2282   local_poles_and_weights_array[2][2][4],
2283   local_weights_array[2][2]  ;
2284   Standard_Real * my_vec_min,
2285   * my_vec_max,
2286   * my_point ;
2287   my_point = (Standard_Real *) &aPoint  ;
2288   //
2289   // initialize in case of rational evaluation
2290   // because RationalDerivative will use all
2291   // the coefficients
2292   //
2293   //
2294   if (WeightsArray != NULL) {
2295
2296     local_poles_array            [0][0][0] = 0.0e0 ;
2297     local_poles_array            [0][0][1] = 0.0e0 ;
2298     local_poles_array            [0][0][2] = 0.0e0 ;
2299     local_weights_array          [0][0]    = 0.0e0 ;
2300     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2301     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2302     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2303     local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2304
2305     local_poles_array            [0][1][0] = 0.0e0 ;
2306     local_poles_array            [0][1][1] = 0.0e0 ;
2307     local_poles_array            [0][1][2] = 0.0e0 ;
2308     local_weights_array          [0][1]    = 0.0e0 ;
2309     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2310     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2311     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2312     local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2313
2314     local_poles_array            [1][0][0] = 0.0e0 ;
2315     local_poles_array            [1][0][1] = 0.0e0 ;
2316     local_poles_array            [1][0][2] = 0.0e0 ;
2317     local_weights_array          [1][0]    = 0.0e0 ;
2318     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2319     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2320     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2321     local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2322
2323     local_poles_array            [1][1][0] = 0.0e0 ;
2324     local_poles_array            [1][1][1] = 0.0e0 ;
2325     local_poles_array            [1][1][2] = 0.0e0 ;
2326     local_weights_array          [1][1]    = 0.0e0 ;
2327     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2328     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2329     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2330     local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2331   }
2332
2333   if (UDegree <= VDegree) {
2334     min_degree = UDegree ;
2335     max_degree = VDegree ;
2336     inverse_min = 1.0e0/USpanLenght ;
2337     inverse_max = 1.0e0/VSpanLenght ;
2338     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2339     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2340
2341     dimension = 3 * (UDegree + 1) ;
2342     my_vec_min = (Standard_Real *) &aVecU ;
2343     my_vec_max = (Standard_Real *) &aVecV ;
2344   }
2345   else {
2346     min_degree = VDegree ;
2347     max_degree = UDegree ;
2348     inverse_min = 1.0e0/VSpanLenght ;
2349     inverse_max = 1.0e0/USpanLenght ;
2350     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2351     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2352     dimension = 3 * (VDegree + 1) ;
2353     my_vec_min = (Standard_Real *) &aVecV ;
2354     my_vec_max = (Standard_Real *) &aVecU ;
2355   }
2356
2357   NCollection_LocalArray<Standard_Real> locpoles (2 * dimension);
2358
2359   PLib::EvalPolynomial(new_parameter[0],
2360                        1,
2361                        max_degree,
2362                        dimension,
2363                        PArray[0],
2364                        locpoles[0]) ;
2365
2366   PLib::EvalPolynomial(new_parameter[1],
2367                        1,
2368                        min_degree,
2369                        3,
2370                        locpoles[0],
2371                        local_poles_array[0][0][0]) ;
2372   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2373                        min_degree,
2374                        3,
2375                        (min_degree << 1) + min_degree,
2376                        locpoles[dimension],
2377                        local_poles_array[1][0][0]) ;
2378
2379   if (WeightsArray != NULL) {
2380     dimension = min_degree + 1 ;
2381     const TColStd_Array2OfReal& refWeights = *WeightsArray;
2382     Standard_Real *
2383       WArray = (Standard_Real *)
2384         &refWeights(refWeights.LowerCol(),refWeights.LowerRow()) ;
2385     PLib::EvalPolynomial(new_parameter[0],
2386                          1,
2387                          max_degree,
2388                          dimension,
2389                          WArray[0],
2390                          locpoles[0]) ;
2391
2392     PLib::EvalPolynomial(new_parameter[1],
2393                          1,
2394                          min_degree,
2395                          1,
2396                          locpoles[0],
2397                          local_weights_array[0][0]) ;
2398     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2399                          min_degree,
2400                          1,
2401                          min_degree,
2402                          locpoles[dimension],
2403                          local_weights_array[1][0]) ;
2404
2405     local_poles_and_weights_array[0][0][0] = local_poles_array  [0][0][0] ;
2406     local_poles_and_weights_array[0][0][1] = local_poles_array  [0][0][1] ;
2407     local_poles_and_weights_array[0][0][2] = local_poles_array  [0][0][2] ;
2408     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0]    ;
2409
2410     local_poles_and_weights_array[0][1][0] = local_poles_array  [0][1][0] ;
2411     local_poles_and_weights_array[0][1][1] = local_poles_array  [0][1][1] ;
2412     local_poles_and_weights_array[0][1][2] = local_poles_array  [0][1][2] ;
2413     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1]    ;
2414
2415     local_poles_and_weights_array[1][0][0] = local_poles_array  [1][0][0] ;
2416     local_poles_and_weights_array[1][0][1] = local_poles_array  [1][0][1] ;
2417     local_poles_and_weights_array[1][0][2] = local_poles_array  [1][0][2] ;
2418     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0]    ;
2419
2420     local_poles_and_weights_array[1][1][0] = local_poles_array  [1][1][0] ;
2421     local_poles_and_weights_array[1][1][1] = local_poles_array  [1][1][1] ;
2422     local_poles_and_weights_array[1][1][2] = local_poles_array  [1][1][2] ;
2423     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1]    ;
2424
2425     BSplSLib::RationalDerivative(1,
2426                                  1,
2427                                  1,
2428                                  1,
2429                                  local_poles_and_weights_array[0][0][0],
2430                                  local_poles_array[0][0][0]) ;
2431   }
2432
2433   my_point  [0] = local_poles_array              [0][0][0] ;
2434   my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2435   my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2436
2437   my_point  [1] = local_poles_array              [0][0][1] ;
2438   my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2439   my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2440
2441   my_point  [2] = local_poles_array              [0][0][2] ;
2442   my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2443   my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2444 }
2445
2446 //=======================================================================
2447 //function : CacheD2
2448 //purpose  : Evaluates the polynomial cache of the Bspline Curve
2449 //
2450 //=======================================================================
2451
2452 void  BSplSLib::CacheD2(const Standard_Real                  UParameter,
2453                         const Standard_Real                  VParameter,
2454                         const  Standard_Integer              UDegree,
2455                         const  Standard_Integer              VDegree,
2456                         const  Standard_Real                 UCacheParameter,
2457                         const  Standard_Real                 VCacheParameter,
2458                         const  Standard_Real                 USpanLenght,
2459                         const  Standard_Real                 VSpanLenght,
2460                         const  TColgp_Array2OfPnt&           PolesArray,
2461                         const  TColStd_Array2OfReal*         WeightsArray,
2462                         gp_Pnt&                              aPoint,
2463                         gp_Vec&                              aVecU,
2464                         gp_Vec&                              aVecV,
2465                         gp_Vec&                              aVecUU,
2466                         gp_Vec&                              aVecUV,
2467                         gp_Vec&                              aVecVV)
2468 {
2469   //
2470   // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2471   // form
2472   // the SpanLenght     is the normalizing factor so that the polynomial is between
2473   // 0 and 1
2474   Standard_Integer
2475     ii,
2476 //  jj,
2477   kk,
2478   index,
2479   dimension,
2480   min_degree,
2481   max_degree  ;
2482
2483   Standard_Real
2484     inverse_min,
2485   inverse_max,
2486   new_parameter[2]  ;
2487
2488   Standard_Real *
2489     PArray = (Standard_Real *)
2490       &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2491   Standard_Real local_poles_array[3][3][3],
2492   local_poles_and_weights_array[3][3][4],
2493   local_weights_array[3][3]  ;
2494   Standard_Real * my_vec_min,
2495   * my_vec_max,
2496   * my_vec_min_min,
2497   * my_vec_max_max,
2498   * my_vec_min_max,
2499   * my_point ;
2500   my_point = (Standard_Real *) &aPoint  ;
2501
2502   //
2503   // initialize in case the min and max degree are less than 2
2504   //
2505   local_poles_array[0][0][0] = 0.0e0 ;
2506   local_poles_array[0][0][1] = 0.0e0 ;
2507   local_poles_array[0][0][2] = 0.0e0 ;
2508   local_poles_array[0][1][0] = 0.0e0 ;
2509   local_poles_array[0][1][1] = 0.0e0 ;
2510   local_poles_array[0][1][2] = 0.0e0 ;
2511   local_poles_array[0][2][0] = 0.0e0 ;
2512   local_poles_array[0][2][1] = 0.0e0 ;
2513   local_poles_array[0][2][2] = 0.0e0 ;
2514
2515   local_poles_array[1][0][0] = 0.0e0 ;
2516   local_poles_array[1][0][1] = 0.0e0 ;
2517   local_poles_array[1][0][2] = 0.0e0 ;
2518   local_poles_array[1][1][0] = 0.0e0 ;
2519   local_poles_array[1][1][1] = 0.0e0 ;
2520   local_poles_array[1][1][2] = 0.0e0 ;
2521   local_poles_array[1][2][0] = 0.0e0 ;
2522   local_poles_array[1][2][1] = 0.0e0 ;
2523   local_poles_array[1][2][2] = 0.0e0 ;
2524
2525   local_poles_array[2][0][0] = 0.0e0 ;
2526   local_poles_array[2][0][1] = 0.0e0 ;
2527   local_poles_array[2][0][2] = 0.0e0 ;
2528   local_poles_array[2][1][0] = 0.0e0 ;
2529   local_poles_array[2][1][1] = 0.0e0 ;
2530   local_poles_array[2][1][2] = 0.0e0 ;
2531   local_poles_array[2][2][0] = 0.0e0 ;
2532   local_poles_array[2][2][1] = 0.0e0 ;
2533   local_poles_array[2][2][2] = 0.0e0 ;
2534   //
2535   // initialize in case of rational evaluation
2536   // because RationalDerivative will use all
2537   // the coefficients
2538   //
2539   //
2540   if (WeightsArray != NULL) {
2541
2542     local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2543     local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2544     local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2545     local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2546     local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2547     local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2548     local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2549     local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2550     local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2551
2552     local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2553     local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2554     local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2555     local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2556     local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2557     local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2558     local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2559     local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2560     local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2561
2562     local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2563     local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2564     local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2565     local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2566     local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2567     local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2568     local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2569     local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2570     local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2571
2572     local_poles_and_weights_array[0][0][3] =
2573       local_weights_array[0][0] = 0.0e0 ;
2574     local_poles_and_weights_array[0][1][3] =
2575       local_weights_array[0][1] = 0.0e0 ;
2576     local_poles_and_weights_array[0][2][3] =
2577       local_weights_array[0][2] = 0.0e0 ;
2578     local_poles_and_weights_array[1][0][3] =
2579       local_weights_array[1][0] = 0.0e0 ;
2580     local_poles_and_weights_array[1][1][3] =
2581       local_weights_array[1][1] = 0.0e0 ;
2582     local_poles_and_weights_array[1][2][3] =
2583       local_weights_array[1][2] = 0.0e0 ;
2584     local_poles_and_weights_array[2][0][3] =
2585       local_weights_array[2][0] = 0.0e0 ;
2586     local_poles_and_weights_array[2][1][3] =
2587       local_weights_array[2][1] = 0.0e0 ;
2588     local_poles_and_weights_array[2][2][3] =
2589       local_weights_array[2][2] = 0.0e0 ;
2590   }
2591
2592   if (UDegree <= VDegree) {
2593     min_degree = UDegree ;
2594     max_degree = VDegree ;
2595     inverse_min = 1.0e0/USpanLenght ;
2596     inverse_max = 1.0e0/VSpanLenght ;
2597     new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2598     new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2599
2600     dimension = 3 * (UDegree + 1) ;
2601     my_vec_min     = (Standard_Real *) &aVecU ;
2602     my_vec_max     = (Standard_Real *) &aVecV ;
2603     my_vec_min_min = (Standard_Real *) &aVecUU ;
2604     my_vec_min_max = (Standard_Real *) &aVecUV ;
2605     my_vec_max_max = (Standard_Real *) &aVecVV ;
2606   }
2607   else {
2608     min_degree = VDegree ;
2609     max_degree = UDegree ;
2610     inverse_min = 1.0e0/VSpanLenght ;
2611     inverse_max = 1.0e0/USpanLenght ;
2612     new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2613     new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2614     dimension = 3 * (VDegree + 1) ;
2615     my_vec_min     = (Standard_Real *) &aVecV ;
2616     my_vec_max     = (Standard_Real *) &aVecU ;
2617     my_vec_min_min = (Standard_Real *) &aVecVV ;
2618     my_vec_min_max = (Standard_Real *) &aVecUV ;
2619     my_vec_max_max = (Standard_Real *) &aVecUU ;
2620   }
2621
2622   NCollection_LocalArray<Standard_Real> locpoles (3 * dimension);
2623
2624   //
2625   // initialize in case min or max degree are less than 2
2626   //
2627   Standard_Integer MinIndMax = 2;
2628   if ( max_degree < 2) MinIndMax = max_degree;
2629   Standard_Integer MinIndMin = 2;
2630   if ( min_degree < 2) MinIndMin = min_degree;
2631
2632   index =  MinIndMax * dimension ;
2633
2634   for (ii = MinIndMax ; ii <  3 ; ii++) {
2635
2636     for (kk = 0 ; kk < dimension ; kk++) {
2637       locpoles[index] = 0.0e0 ;
2638       index += 1 ;
2639     }
2640   }
2641
2642   PLib::EvalPolynomial(new_parameter[0],
2643                        MinIndMax,
2644                        max_degree,
2645                        dimension,
2646                        PArray[0],
2647                        locpoles[0]) ;
2648
2649   PLib::EvalPolynomial(new_parameter[1],
2650                        MinIndMin,
2651                        min_degree,
2652                        3,
2653                        locpoles[0],
2654                        local_poles_array[0][0][0]) ;
2655   PLib::EvalPolynomial(new_parameter[1],
2656                        1,
2657                        min_degree,
2658                        3,
2659                        locpoles[dimension],
2660                        local_poles_array[1][0][0]) ;
2661
2662   PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2663                        min_degree,
2664                        3,
2665                        (min_degree << 1) + min_degree,
2666                        locpoles[dimension + dimension],
2667                        local_poles_array[2][0][0]) ;
2668
2669   if (WeightsArray != NULL) {
2670     dimension = min_degree + 1 ;
2671     const TColStd_Array2OfReal& refWeights = *WeightsArray;
2672     Standard_Real *
2673       WArray = (Standard_Real *)
2674         &refWeights(refWeights.LowerCol(),refWeights.LowerRow()) ;
2675     PLib::EvalPolynomial(new_parameter[0],
2676                          MinIndMax,
2677                          max_degree,
2678                          dimension,
2679                          WArray[0],
2680                          locpoles[0]) ;
2681
2682     PLib::EvalPolynomial(new_parameter[1],
2683                          MinIndMin,
2684                          min_degree,
2685                          1,
2686                          locpoles[0],
2687                          local_weights_array[0][0]) ;
2688     PLib::EvalPolynomial(new_parameter[1],
2689                          1,
2690                          min_degree,
2691                          1,
2692                          locpoles[dimension],
2693                          local_weights_array[1][0]) ;
2694     PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2695                          min_degree,
2696                          1,
2697                          min_degree,
2698                          locpoles[dimension + dimension],
2699                          local_weights_array[2][0]) ;
2700
2701
2702     local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2703     local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2704     local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2705     local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2706     local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2707     local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2708     local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2709     local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2710     local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2711
2712     local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2713     local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2714     local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2715     local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2716     local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2717     local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2718     local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2719     local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2720     local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2721
2722     local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2723     local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2724     local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2725     local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2726     local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2727     local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2728     local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2729     local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2730     local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2731
2732
2733     local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2734     local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2735     local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2736     local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2737     local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2738     local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2739     local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2740     local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2741     local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2742
2743     BSplSLib::RationalDerivative(2,
2744                                  2,
2745                                  2,
2746                                  2,
2747                                  local_poles_and_weights_array[0][0][0],
2748                                  local_poles_array[0][0][0]) ;
2749   }
2750
2751
2752   Standard_Real minmin = inverse_min * inverse_min;
2753   Standard_Real minmax = inverse_min * inverse_max;
2754   Standard_Real maxmax = inverse_max * inverse_max;
2755
2756   my_point      [0] = local_poles_array              [0][0][0] ;
2757   my_vec_min    [0] = inverse_min * local_poles_array[0][1][0] ;
2758   my_vec_max    [0] = inverse_max * local_poles_array[1][0][0] ;
2759   my_vec_min_min[0] = minmin * local_poles_array     [0][2][0] ;
2760   my_vec_min_max[0] = minmax * local_poles_array     [1][1][0] ;
2761   my_vec_max_max[0] = maxmax * local_poles_array     [2][0][0] ;
2762
2763   my_point      [1] = local_poles_array              [0][0][1] ;
2764   my_vec_min    [1] = inverse_min * local_poles_array[0][1][1] ;
2765   my_vec_max    [1] = inverse_max * local_poles_array[1][0][1] ;
2766   my_vec_min_min[1] = minmin * local_poles_array     [0][2][1] ;
2767   my_vec_min_max[1] = minmax * local_poles_array     [1][1][1] ;
2768   my_vec_max_max[1] = maxmax * local_poles_array     [2][0][1] ;
2769
2770   my_point      [2] = local_poles_array              [0][0][2] ;
2771   my_vec_min    [2] = inverse_min * local_poles_array[0][1][2] ;
2772   my_vec_max    [2] = inverse_max * local_poles_array[1][0][2] ;
2773   my_vec_min_min[2] = minmin * local_poles_array     [0][2][2] ;
2774   my_vec_min_max[2] = minmax * local_poles_array     [1][1][2] ;
2775   my_vec_max_max[2] = maxmax * local_poles_array     [2][0][2] ;
2776 }
2777
2778 //=======================================================================
2779 //function : MovePoint
2780 //purpose  : Find the new poles which allows  an old point (with a
2781 //           given  u and v as parameters) to reach a new position
2782 //=======================================================================
2783
2784 void BSplSLib::MovePoint (const Standard_Real            U,
2785                           const Standard_Real            V,
2786                           const gp_Vec&                  Displ,
2787                           const Standard_Integer         UIndex1,
2788                           const Standard_Integer         UIndex2,
2789                           const Standard_Integer         VIndex1,
2790                           const Standard_Integer         VIndex2,
2791                           const Standard_Integer         UDegree,
2792                           const Standard_Integer         VDegree,
2793                           const Standard_Boolean         Rational,
2794                           const TColgp_Array2OfPnt&      Poles,
2795                           const TColStd_Array2OfReal&    Weights,
2796                           const TColStd_Array1OfReal&    UFlatKnots,
2797                           const TColStd_Array1OfReal&    VFlatKnots,
2798                           Standard_Integer&              UFirstIndex,
2799                           Standard_Integer&              ULastIndex,
2800                           Standard_Integer&              VFirstIndex,
2801                           Standard_Integer&              VLastIndex,
2802                           TColgp_Array2OfPnt&            NewPoles)
2803 {
2804   // calculate the UBSplineBasis in the parameter U
2805   Standard_Integer UFirstNonZeroBsplineIndex;
2806   math_Matrix UBSplineBasis(1, 1,
2807                             1, UDegree+1);
2808   Standard_Integer ErrorCod1 =  BSplCLib::EvalBsplineBasis(1,
2809                                                            0,
2810                                                            UDegree+1,
2811                                                            UFlatKnots,
2812                                                            U,
2813                                                            UFirstNonZeroBsplineIndex,
2814                                                            UBSplineBasis);
2815   // calculate the VBSplineBasis in the parameter V
2816   Standard_Integer VFirstNonZeroBsplineIndex;
2817   math_Matrix VBSplineBasis(1, 1,
2818                             1, VDegree+1);
2819   Standard_Integer ErrorCod2 =  BSplCLib::EvalBsplineBasis(1,
2820                                                            0,
2821                                                            VDegree+1,
2822                                                            VFlatKnots,
2823                                                            V,
2824                                                            VFirstNonZeroBsplineIndex,
2825                                                            VBSplineBasis);
2826   if (ErrorCod1 || ErrorCod2) {
2827     UFirstIndex = 0;
2828     ULastIndex = 0;
2829     VFirstIndex = 0;
2830     VLastIndex = 0;
2831     return;
2832   }
2833
2834   // find the span which is predominant for parameter U
2835   UFirstIndex = UFirstNonZeroBsplineIndex;
2836   ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2837   if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2838   if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2839
2840   Standard_Real maxValue = 0.0;
2841   Standard_Integer i, ukk1=0, ukk2;
2842
2843   for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2844     if (UBSplineBasis(1,i) > maxValue) {
2845       ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2846       maxValue = UBSplineBasis(1,i);
2847     }
2848   }
2849
2850   // find a ukk2 if symetriy
2851   ukk2 = ukk1;
2852   i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2853   if ((ukk1+1) <= ULastIndex) {
2854     if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2855       ukk2 = ukk1+1;
2856     }
2857   }
2858
2859   // find the span which is predominant for parameter V
2860   VFirstIndex = VFirstNonZeroBsplineIndex;
2861   VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2862
2863   if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2864   if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2865
2866   maxValue = 0.0;
2867   Standard_Integer j, vkk1=0, vkk2;
2868
2869   for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2870     if (VBSplineBasis(1,j) > maxValue) {
2871       vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2872       maxValue = VBSplineBasis(1,j);
2873     }
2874   }
2875
2876   // find a vkk2 if symetriy
2877   vkk2 = vkk1;
2878   j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2879   if ((vkk1+1) <= VLastIndex) {
2880     if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2881       vkk2 = vkk1+1;
2882     }
2883   }
2884
2885   // compute the vector of displacement
2886   Standard_Real D1 = 0.0;
2887   Standard_Real D2 = 0.0;
2888   Standard_Real hN, Coef, DvalU, DvalV;
2889
2890   Standard_Integer ii, jj;
2891
2892   for (i = 1; i <= UDegree+1; i++) {
2893     ii = i + UFirstNonZeroBsplineIndex - 1;
2894     if (ii < ukk1) {
2895       DvalU = ukk1-ii;
2896     }
2897     else if (ii > ukk2) {
2898       DvalU = ii - ukk2;
2899     }
2900     else {
2901       DvalU = 0.0;
2902     }
2903
2904     for (j = 1; j <= VDegree+1; j++) {
2905       jj = j + VFirstNonZeroBsplineIndex - 1;
2906       if (Rational) {
2907         hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2908         D2 += hN;
2909       }
2910       else {
2911         hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2912       }
2913       if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2914         if (jj < vkk1) {
2915           DvalV = vkk1-jj;
2916         }
2917         else if (jj > vkk2) {
2918           DvalV = jj - vkk2;
2919         }
2920         else {
2921           DvalV = 0.0;
2922         }
2923         D1 += 1./(DvalU + DvalV + 1.) * hN;
2924       }
2925     }
2926   }
2927
2928   if (Rational) {
2929     Coef = D2/D1;
2930   }
2931   else {
2932     Coef = 1./D1;
2933   }
2934
2935   // compute the new poles
2936
2937   for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2938     if (i < ukk1) {
2939       DvalU = ukk1-i;
2940     }
2941     else if (i > ukk2) {
2942       DvalU = i - ukk2;
2943     }
2944     else {
2945       DvalU = 0.0;
2946     }
2947
2948     for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2949       if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2950         if (j < vkk1) {
2951           DvalV = vkk1-j;
2952         }
2953         else if (j > vkk2) {
2954           DvalV = j - vkk2;
2955         }
2956         else {
2957           DvalV = 0.0;
2958         }
2959         NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2960       }
2961       else {
2962         NewPoles(i,j) = Poles(i,j);
2963       }
2964     }
2965   }
2966 }
2967
2968 //=======================================================================
2969 // function : Resolution
2970 // purpose  : this computes an estimate for the maximum of the
2971 // partial derivatives both in U and in V
2972 //
2973 //
2974 // The calculation resembles at the calculation of curves with
2975 // additional index for the control point. Let Si,j be the
2976 // control points for ls surface  and  Di,j  the weights.
2977 // The checking of upper bounds for the partial derivatives
2978 // will be omitted and Su is the next upper bound in the polynomial case :
2979 //
2980 //
2981 //
2982 //                        |  Si,j - Si-1,j  |
2983 //          d *   Max     |  -------------  |
2984 //                i = 2,n |     ti+d - ti   |
2985 //                i=1.m
2986 //
2987 //
2988 // and in the rational case :
2989 //
2990 //
2991 //
2992 //                         Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2993 //   Max   Max       d  *  -----------------------------------------------
2994 // k=1,n  i dans Rj                   ti+d  - ti
2995 // j=1,m
2996 //  ----------------------------------------------------------------------
2997 //
2998 //               Min    Di,j
2999 //              i=1,n
3000 //              j=1,m
3001 //
3002 //
3003 //
3004 // with Rj = {j-d, ....,  j+d+d+1}.
3005 //
3006 //
3007 //=======================================================================
3008
3009 void BSplSLib::Resolution(const TColgp_Array2OfPnt&      Poles,
3010                           const TColStd_Array2OfReal*    Weights,
3011                           const TColStd_Array1OfReal&    UKnots,
3012                           const TColStd_Array1OfReal&    VKnots,
3013                           const TColStd_Array1OfInteger& UMults,
3014                           const TColStd_Array1OfInteger& VMults,
3015                           const Standard_Integer         UDegree,
3016                           const Standard_Integer         VDegree,
3017                           const Standard_Boolean         URational,
3018                           const Standard_Boolean         VRational,
3019                           const Standard_Boolean         UPeriodic,
3020                           const Standard_Boolean         VPeriodic,
3021                           const Standard_Real            Tolerance3D,
3022                           Standard_Real&                 UTolerance,
3023                           Standard_Real&                 VTolerance)
3024 {
3025   Standard_Real Wij,Wmj,Wji,Wjm;
3026   Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
3027   Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
3028   Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
3029   Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
3030
3031   max_derivative[0] = max_derivative[1] = 0.0e0 ;
3032
3033   Standard_Integer PRowLength, PColLength;
3034   Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
3035   Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
3036   Standard_Integer num_poles[2],num_flat_knots[2];
3037
3038   num_flat_knots[0] =
3039     BSplCLib::KnotSequenceLength(UMults,
3040                                  UDegree,
3041                                  UPeriodic) ;
3042   num_flat_knots[1] =
3043     BSplCLib::KnotSequenceLength(VMults,
3044                                  VDegree,
3045                                  VPeriodic) ;
3046   TColStd_Array1OfReal  flat_knots_in_u(1,num_flat_knots[0]) ;
3047   TColStd_Array1OfReal  flat_knots_in_v(1,num_flat_knots[1]) ;
3048   BSplCLib::KnotSequence(UKnots,
3049                          UMults,
3050                          UDegree,
3051                          UPeriodic,
3052                          flat_knots_in_u) ;
3053   BSplCLib::KnotSequence(VKnots,
3054                          VMults,
3055                          VDegree,
3056                          VPeriodic,
3057                          flat_knots_in_v) ;
3058   PRowLength = Poles.RowLength();
3059   PColLength = Poles.ColLength();
3060   if (URational || VRational) {
3061     Standard_Integer Wsize = PRowLength * PColLength;
3062     const TColStd_Array2OfReal& refWights = *Weights;
3063     const Standard_Real * WG = &refWights(refWights.LowerRow(), refWights.LowerCol());
3064     min_weights = WG[0];
3065
3066     for (ii = 1 ; ii < Wsize ; ii++) {
3067       min = WG[ii];
3068       if (min_weights > min) min_weights = min;
3069     }
3070   }
3071   Standard_Integer UD1 = UDegree + 1;
3072   Standard_Integer VD1 = VDegree + 1;
3073   num_poles[0] = num_flat_knots[0] - UD1;
3074   num_poles[1] = num_flat_knots[1] - VD1;
3075   poles_length[0] = PColLength;
3076   poles_length[1] = PRowLength;
3077   if(URational) {
3078     Standard_Integer UD2 = UDegree << 1;
3079     Standard_Integer VD2 = VDegree << 1;
3080
3081     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3082       ii_index = (ii - 1) % poles_length[0] + 1 ;
3083       ii_minus = (ii - 2) % poles_length[0] + 1 ;
3084       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3085       inverse = 1.0e0 / inverse ;
3086       lower[0] = ii - UD1;
3087       if (lower[0] < 1) lower[0] = 1;
3088       upper[0] = ii + UD2 + 1;
3089       if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
3090
3091       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3092         jj_index = (jj - 1) % poles_length[1] + 1 ;
3093         lower[1] = jj - VD1;
3094         if (lower[1] < 1) lower[1] = 1;
3095         upper[1] = jj + VD2 + 1;
3096         if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
3097         const gp_Pnt& Pij = Poles  .Value(ii_index,jj_index);
3098         Wij               = Weights->Value(ii_index,jj_index);
3099         const gp_Pnt& Pmj = Poles  .Value(ii_minus,jj_index);
3100         Wmj               = Weights->Value(ii_minus,jj_index);
3101         Xij = Pij.X();
3102         Yij = Pij.Y();
3103         Zij = Pij.Z();
3104         Xmj = Pmj.X();
3105         Ymj = Pmj.Y();
3106         Zmj = Pmj.Z();
3107
3108         for (pp = lower[0] ; pp <= upper[0] ; pp++) {
3109           pp_index = (pp - 1) % poles_length[0] + 1 ;
3110
3111           for (qq = lower[1] ; qq <= upper[1] ; qq++) {
3112             value = 0.0e0 ;
3113             qq_index = (qq - 1) % poles_length[1] + 1 ;
3114             const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
3115             Xpq = Ppq.X();
3116             Ypq = Ppq.Y();
3117             Zpq = Ppq.Z();
3118             factor  = (Xpq - Xij) * Wij;
3119             factor -= (Xpq - Xmj) * Wmj;
3120             if (factor < 0) factor = - factor;
3121             value  += factor ;
3122             factor  = (Ypq - Yij) * Wij;
3123             factor -= (Ypq - Ymj) * Wmj;
3124             if (factor < 0) factor = - factor;
3125             value  += factor ;
3126             factor  = (Zpq - Zij) * Wij;
3127             factor -= (Zpq - Zmj) * Wmj;
3128             if (factor < 0) factor = - factor;
3129             value  += factor ;
3130             value *= inverse ;
3131             if (max_derivative[0] < value) max_derivative[0] = value ;
3132           }
3133         }
3134       }
3135     }
3136     max_derivative[0] /= min_weights ;
3137   }
3138   else {
3139
3140     for (ii = 2 ; ii <= num_poles[0] ; ii++) {
3141       ii_index = (ii - 1) % poles_length[0] + 1 ;
3142       ii_minus = (ii - 2) % poles_length[0] + 1 ;
3143       inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
3144       inverse = 1.0e0 / inverse ;
3145
3146       for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
3147         jj_index = (jj - 1) % poles_length[1] + 1 ;
3148         value = 0.0e0 ;
3149         const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
3150         const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
3151         factor = Pij.X() - Pmj.X();
3152         if (factor < 0) factor = - factor;
3153         value += factor;
3154         factor = Pij.Y() - Pmj.Y();
3155         if (factor < 0) factor = - factor;
3156         value += factor;
3157         factor = Pij.Z() - Pmj.Z();
3158         if (factor < 0) factor = - factor;
3159         value += factor;
3160         value *= inverse ;
3161         if (max_derivative[0] < value) max_derivative[0] = value ;
3162       }
3163     }
3164   }
3165   max_derivative[0] *= UDegree ;
3166   if(VRational) {
3167     Standard_Integer UD2 = UDegree << 1;
3168     Standard_Integer VD2 = VDegree << 1;
3169
3170     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3171       ii_index = (ii - 1) % poles_length[1] + 1 ;
3172       ii_minus = (ii - 2) % poles_length[1] + 1 ;
3173       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3174       inverse = 1.0e0 / inverse ;
3175       lower[0] = ii - VD1;
3176       if (lower[0] < 1) lower[0] = 1;
3177       upper[0] = ii + VD2 + 1;
3178       if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
3179
3180       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3181         jj_index = (jj - 1) % poles_length[0] + 1 ;
3182         lower[1] = jj - UD1;
3183         if (lower[1] < 1) lower[1] = 1;
3184         upper[1] = jj + UD2 + 1;
3185         if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
3186         const gp_Pnt& Pji = Poles  .Value(jj_index,ii_index);
3187         Wji               = Weights->Value(jj_index,ii_index);
3188         const gp_Pnt& Pjm = Poles  .Value(jj_index,ii_minus);
3189         Wjm               = Weights->Value(jj_index,ii_minus);
3190         Xji = Pji.X();
3191         Yji = Pji.Y();
3192         Zji = Pji.Z();
3193         Xjm = Pjm.X();
3194         Yjm = Pjm.Y();
3195         Zjm = Pjm.Z();
3196
3197         for (pp = lower[1] ; pp <= upper[1] ; pp++) {
3198           pp_index = (pp - 1) % poles_length[1] + 1 ;
3199
3200           for (qq = lower[0] ; qq <= upper[0] ; qq++) {
3201             value = 0.0e0 ;
3202             qq_index = (qq - 1) % poles_length[0] + 1 ;
3203             const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
3204             Xqp = Pqp.X();
3205             Yqp = Pqp.Y();
3206             Zqp = Pqp.Z();
3207             factor  = (Xqp - Xji) * Wji;
3208             factor -= (Xqp - Xjm) * Wjm;
3209             if (factor < 0) factor = - factor;
3210             value += factor ;
3211             factor  = (Yqp - Yji) * Wji;
3212             factor -= (Yqp - Yjm) * Wjm;
3213             if (factor < 0) factor = - factor;
3214             value += factor ;
3215             factor  = (Zqp - Zji) * Wji;
3216             factor -= (Zqp - Zjm) * Wjm;
3217             if (factor < 0) factor = - factor;
3218             value += factor ;
3219             value *= inverse ;
3220             if (max_derivative[1] < value) max_derivative[1] = value ;
3221           }
3222         }
3223       }
3224     }
3225     max_derivative[1] /= min_weights ;
3226   }
3227   else {
3228
3229     for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3230       ii_index = (ii - 1) % poles_length[1] + 1 ;
3231       ii_minus = (ii - 2) % poles_length[1] + 1 ;
3232       inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3233       inverse = 1.0e0 / inverse ;
3234
3235       for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3236         jj_index = (jj - 1) % poles_length[0] + 1 ;
3237         value = 0.0e0 ;
3238         const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3239         const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3240         factor = Pji.X() - Pjm.X() ;
3241         if (factor < 0) factor = - factor;
3242         value += factor;
3243         factor = Pji.Y() - Pjm.Y() ;
3244         if (factor < 0) factor = - factor;
3245         value += factor;
3246         factor = Pji.Z() - Pjm.Z() ;
3247         if (factor < 0) factor = - factor;
3248         value += factor;
3249         value *= inverse ;
3250         if (max_derivative[1] < value) max_derivative[1] = value ;
3251       }
3252     }
3253   }
3254   max_derivative[1] *= VDegree ;
3255   max_derivative[0] *= M_SQRT2 ;
3256   max_derivative[1] *= M_SQRT2 ;
3257   if(max_derivative[0] && max_derivative[1]) {
3258     UTolerance = Tolerance3D / max_derivative[0] ;
3259     VTolerance = Tolerance3D / max_derivative[1] ;
3260   }
3261   else {
3262     UTolerance=VTolerance=0.0;
3263 #ifdef OCCT_DEBUG
3264     cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3265 #endif
3266   }
3267 }
3268
3269 //=======================================================================
3270 //function : Interpolate
3271 //purpose  :
3272 //=======================================================================
3273
3274 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3275                            const Standard_Integer VDegree,
3276                            const TColStd_Array1OfReal& UFlatKnots,
3277                            const TColStd_Array1OfReal& VFlatKnots,
3278                            const TColStd_Array1OfReal& UParameters,
3279                            const TColStd_Array1OfReal& VParameters,
3280                            TColgp_Array2OfPnt& Poles,
3281                            TColStd_Array2OfReal& Weights,
3282                            Standard_Integer& InversionProblem)
3283 {
3284   Standard_Integer ii, jj, ll, kk, dimension;
3285   Standard_Integer ULength = UParameters.Length();
3286   Standard_Integer VLength = VParameters.Length();
3287   Standard_Real * poles_array;
3288
3289   // extraction of iso u
3290   dimension = 4*ULength;
3291   TColStd_Array2OfReal Points(1, VLength,
3292                               1, dimension);
3293
3294   Handle(TColStd_HArray1OfInteger) ContactOrder =
3295     new (TColStd_HArray1OfInteger)(1, VLength);
3296   ContactOrder->Init(0);
3297
3298   for (ii=1; ii <= VLength; ii++) {
3299
3300     for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3301       Points(ii,ll)   = Poles(jj, ii).X();
3302       Points(ii,ll+1) = Poles(jj, ii).Y();
3303       Points(ii,ll+2) = Poles(jj, ii).Z();
3304       Points(ii,ll+3) = Weights(jj,ii)  ;
3305     }
3306   }
3307
3308   // interpolation of iso u
3309   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3310   BSplCLib::Interpolate(VDegree,
3311                         VFlatKnots,
3312                         VParameters,
3313                         ContactOrder->Array1(),
3314                         dimension,
3315                         poles_array[0],
3316                         InversionProblem) ;
3317   if (InversionProblem != 0) return;
3318
3319   // extraction of iso v
3320
3321   dimension = VLength*4;
3322   TColStd_Array2OfReal IsoPoles(1, ULength,
3323                                 1, dimension);
3324
3325   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3326   ContactOrder->Init(0);
3327   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3328
3329   for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3330
3331     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3332       IsoPoles (ii,ll)   = Points(jj, kk);
3333       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3334       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3335       IsoPoles (ii,ll+3) = Points(jj, kk+3);
3336     }
3337   }
3338   // interpolation of iso v
3339   BSplCLib::Interpolate(UDegree,
3340                         UFlatKnots,
3341                         UParameters,
3342                         ContactOrder->Array1(),
3343                         dimension,
3344                         poles_array[0],
3345                         InversionProblem);
3346
3347   // return results
3348
3349   for (ii=1; ii <= ULength; ii++) {
3350
3351     for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3352       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3353       Poles.SetValue(ii, jj, Pnt);
3354       Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3355     }
3356   }
3357 }
3358
3359 //=======================================================================
3360 //function : Interpolate
3361 //purpose  :
3362 //=======================================================================
3363
3364 void BSplSLib::Interpolate(const Standard_Integer UDegree,
3365                            const Standard_Integer VDegree,
3366                            const TColStd_Array1OfReal& UFlatKnots,
3367                            const TColStd_Array1OfReal& VFlatKnots,
3368                            const TColStd_Array1OfReal& UParameters,
3369                            const TColStd_Array1OfReal& VParameters,
3370                            TColgp_Array2OfPnt& Poles,
3371                            Standard_Integer& InversionProblem)
3372 {
3373   Standard_Integer ii, jj, ll, kk, dimension;
3374   Standard_Integer ULength = UParameters.Length();
3375   Standard_Integer VLength = VParameters.Length();
3376   Standard_Real * poles_array;
3377
3378   // extraction of iso u
3379   dimension = 3*ULength;
3380   TColStd_Array2OfReal Points(1, VLength,
3381                               1, dimension);
3382
3383   Handle(TColStd_HArray1OfInteger) ContactOrder =
3384     new (TColStd_HArray1OfInteger)(1, VLength);
3385   ContactOrder->Init(0);
3386
3387   for (ii=1; ii <= VLength; ii++) {
3388
3389     for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3390       Points(ii,ll)   = Poles(jj, ii).X();
3391       Points(ii,ll+1) = Poles(jj, ii).Y();
3392       Points(ii,ll+2) = Poles(jj, ii).Z();
3393     }
3394   }
3395
3396   // interpolation of iso u
3397   poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3398   BSplCLib::Interpolate(VDegree,
3399                         VFlatKnots,
3400                         VParameters,
3401                         ContactOrder->Array1(),
3402                         dimension,
3403                         poles_array[0],
3404                         InversionProblem) ;
3405   if (InversionProblem != 0) return;
3406
3407   // extraction of iso v
3408
3409   dimension = VLength*3;
3410   TColStd_Array2OfReal IsoPoles(1, ULength,
3411                                 1, dimension);
3412
3413   ContactOrder =  new (TColStd_HArray1OfInteger)(1, ULength);
3414   ContactOrder->Init(0);
3415   poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3416
3417   for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3418
3419     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3420       IsoPoles (ii,ll)   = Points(jj, kk);
3421       IsoPoles (ii,ll+1) = Points(jj, kk+1);
3422       IsoPoles (ii,ll+2) = Points(jj, kk+2);
3423     }
3424   }
3425   // interpolation of iso v
3426   BSplCLib::Interpolate(UDegree,
3427                         UFlatKnots,
3428                         UParameters,
3429                         ContactOrder->Array1(),
3430                         dimension,
3431                         poles_array[0],
3432                         InversionProblem);
3433
3434   // return results
3435
3436   for (ii=1; ii <= ULength; ii++) {
3437
3438     for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3439       gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3440       Poles.SetValue(ii, jj, Pnt);
3441     }
3442   }
3443 }
3444
3445 //=======================================================================
3446 //function : FunctionMultiply
3447 //purpose  :
3448 //=======================================================================
3449
3450 void BSplSLib::FunctionMultiply
3451 (const BSplSLib_EvaluatorFunction&           Function,
3452  const Standard_Integer                      UBSplineDegree,
3453  const Standard_Integer                      VBSplineDegree,
3454  const TColStd_Array1OfReal&                 UBSplineKnots,
3455  const TColStd_Array1OfReal&                 VBSplineKnots,
3456  const TColStd_Array1OfInteger *             UMults,
3457  const TColStd_Array1OfInteger *             VMults,
3458  const TColgp_Array2OfPnt&                   Poles,
3459  const TColStd_Array2OfReal*                 Weights,
3460  const TColStd_Array1OfReal&                 UFlatKnots,
3461  const TColStd_Array1OfReal&                 VFlatKnots,
3462  const Standard_Integer                      UNewDegree,
3463  const Standard_Integer                      VNewDegree,
3464  TColgp_Array2OfPnt&                         NewNumerator,
3465  TColStd_Array2OfReal&                       NewDenominator,
3466  Standard_Integer&                           Status)
3467 {
3468   Standard_Integer num_uparameters,
3469 //  ii,jj,kk,
3470   ii,jj,
3471   error_code,
3472   num_vparameters ;
3473   Standard_Real    result ;
3474
3475   num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3476   num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3477   TColStd_Array1OfReal  UParameters(1,num_uparameters) ;
3478   TColStd_Array1OfReal  VParameters(1,num_vparameters) ;
3479
3480   if ((NewNumerator.ColLength() == num_uparameters) &&
3481       (NewNumerator.RowLength() == num_vparameters) &&
3482       (NewDenominator.ColLength() == num_uparameters) &&
3483       (NewDenominator.RowLength() == num_vparameters)) {
3484
3485
3486     BSplCLib::BuildSchoenbergPoints(UNewDegree,
3487                                     UFlatKnots,
3488                                     UParameters) ;
3489
3490     BSplCLib::BuildSchoenbergPoints(VNewDegree,
3491                                     VFlatKnots,
3492                                     VParameters) ;
3493
3494     for (ii = 1 ; ii <= num_uparameters ; ii++) {
3495
3496       for (jj = 1 ; jj <= num_vparameters ; jj++) {
3497         HomogeneousD0(UParameters(ii),
3498                       VParameters(jj),
3499                       0,
3500                       0,
3501                       Poles,
3502                       Weights,
3503                       UBSplineKnots,
3504                       VBSplineKnots,
3505                       UMults,
3506                       VMults,
3507                       UBSplineDegree,
3508                       VBSplineDegree,
3509                       Standard_True,
3510                       Standard_True,
3511                       Standard_False,
3512                       Standard_False,
3513                       NewDenominator(ii,jj),
3514                       NewNumerator(ii,jj)) ;
3515
3516         Function.Evaluate (0,
3517                  UParameters(ii),
3518                  VParameters(jj),
3519                  result,
3520                  error_code) ;
3521         if (error_code) {
3522           Standard_ConstructionError::Raise();
3523         }
3524         gp_Pnt& P = NewNumerator(ii,jj);
3525         P.SetX(P.X() * result);
3526         P.SetY(P.Y() * result);
3527         P.SetZ(P.Z() * result);
3528         NewDenominator(ii,jj) *= result ;
3529       }
3530     }
3531     Interpolate(UNewDegree,
3532                 VNewDegree,
3533                 UFlatKnots,
3534                 VFlatKnots,
3535                 UParameters,
3536                 VParameters,
3537                 NewNumerator,
3538                 NewDenominator,
3539                 Status) ;
3540   }
3541   else {
3542     Standard_ConstructionError::Raise();
3543   }
3544 }