src/Geom/Geom_BezierSurface.cxx

   1 // Created on: 1993-03-09
   2 // Created by: JCV
   3 // Copyright (c) 1993-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 //Passage en classe persistante - 23/01/91
  18 //Modif suite a la deuxieme revue de projet toolkit Geometry -23/01/91
  19 // pmn : 21/10/95 ; Correction de la methode segment (PRO5853)
  20 // pmn : 31-Dec-96; Bonne gestion des poids (bug PRO4622)
  21 // xab : 07-Jul-97; le cache est instable en degree 21
  22 //       a partir du degree 15 on ne l'utilise plus
  23 // RBD : 15/10/98 ; Le cache est desormais defini sur [-1,1] (pro15537).
  24 // pmn : 10/12/98 ; Update de la methode segment (suite a la modif de cache).
  25
  26 #define No_Standard_OutOfRange
  27 #define No_Standard_DimensionError
  28
  29
  30 #include <BSplCLib.hxx>
  31 #include <BSplSLib.hxx>
  32 #include <Geom_BezierCurve.hxx>
  33 #include <Geom_BezierSurface.hxx>
  34 #include <Geom_Curve.hxx>
  35 #include <Geom_Geometry.hxx>
  36 #include <gp.hxx>
  37 #include <gp_Pnt.hxx>
  38 #include <gp_Trsf.hxx>
  39 #include <gp_Vec.hxx>
  40 #include <gp_XYZ.hxx>
  41 #include <PLib.hxx>
  42 #include <Precision.hxx>
  43 #include <Standard_ConstructionError.hxx>
  44 #include <Standard_DimensionError.hxx>
  45 #include <Standard_OutOfRange.hxx>
  46 #include <Standard_RangeError.hxx>
  47 #include <Standard_Type.hxx>
  48 #include <TColStd_Array1OfInteger.hxx>
  49
  50 //=======================================================================
  51 //function : Rational
  52 //purpose  : check rationality of an array of weights
  53 //=======================================================================
  54 static void Rational(const TColStd_Array2OfReal& Weights,
  55                      Standard_Boolean& Urational,
  56                      Standard_Boolean& Vrational)
  57 {
  58   Standard_Integer I,J;
  59   J = Weights.LowerCol ();
  60   Vrational = Standard_False;
  61   while (!Vrational && J <= Weights.UpperCol()) {
  62     I = Weights.LowerRow();
  63     while (!Vrational && I <= Weights.UpperRow() - 1) {
  64       Vrational = (Abs(Weights (I, J) - Weights (I+1, J))
  65                    > Epsilon (Abs(Weights (I, J))));
  66       I++;
  67     }
  68     J++;
  69   }
  70
  71   I = Weights.LowerRow ();
  72   Urational = Standard_False;
  73   while (!Urational && I <= Weights.UpperRow()) {
  74     J = Weights.LowerCol();
  75     while (!Urational && J <= Weights.UpperCol() - 1) {
  76       Urational = (Abs(Weights (I, J) - Weights (I, J+1))
  77                    > Epsilon (Abs(Weights (I, J))));
  78       J++;
  79     }
  80     I++;
  81   }
  82 }
  83
  84 //=======================================================================
  85 //function : AddPoleCol
  86 //purpose  : Internal use only.
  87 //=======================================================================
  88
  89 static void AddPoleCol
  90   (const TColgp_Array2OfPnt& Poles,
  91    const TColgp_Array1OfPnt& PoleCol,
  92    const Standard_Integer    AfterIndex,
  93          TColgp_Array2OfPnt& NewPoles)
  94 {
  95   Standard_Integer InsertIndex = AfterIndex + NewPoles.LowerCol();
  96   Standard_Integer Offset = NewPoles.LowerRow() - PoleCol.Lower();
  97   Standard_Integer ColIndex = NewPoles.LowerCol();
  98   Standard_Integer RowIndex;
  99   while (ColIndex < InsertIndex) {
 100     RowIndex = NewPoles.LowerRow();
 101     while (RowIndex <= NewPoles.UpperRow()){
 102       NewPoles (RowIndex, ColIndex) =  Poles (RowIndex, ColIndex);
 103       RowIndex++;
 104     }
 105     ColIndex++;
 106   }
 107   RowIndex = NewPoles.LowerRow();
 108   while (RowIndex <= NewPoles.UpperRow()){
 109     NewPoles (RowIndex, ColIndex) = PoleCol (RowIndex - Offset);
 110     RowIndex++;
 111   }
 112   ColIndex++;
 113   while (ColIndex <= NewPoles.UpperCol()) {
 114     RowIndex = NewPoles.LowerRow();
 115     while (RowIndex <= NewPoles.UpperRow()){
 116       NewPoles (RowIndex, ColIndex) =  Poles (RowIndex, ColIndex - 1);
 117       RowIndex++;
 118     }
 119     ColIndex++;
 120   }
 121 }
 122
 123 //=======================================================================
 124 //function : AddRatPoleCol
 125 //purpose  : Internal use only.
 126 //=======================================================================
 127
 128 static void AddRatPoleCol
 129   (const TColgp_Array2OfPnt&   Poles,
 130    const TColStd_Array2OfReal& Weights,
 131    const TColgp_Array1OfPnt&   PoleCol,
 132    const TColStd_Array1OfReal& PoleWeightCol,
 133    const Standard_Integer      AfterIndex,
 134          TColgp_Array2OfPnt&   NewPoles,
 135          TColStd_Array2OfReal& NewWeights)
 136 {
 137   Standard_Integer InsertIndex = AfterIndex + NewPoles.LowerCol();
 138   Standard_Integer OffsetPol   = NewPoles.LowerRow() - PoleCol.Lower();
 139   Standard_Integer OffsetW     = NewWeights.LowerRow() - PoleWeightCol.Lower();
 140
 141   Standard_Integer ColIndex = NewPoles.LowerCol();
 142   Standard_Integer RowIndex;
 143   while (ColIndex < InsertIndex) {
 144     RowIndex = NewPoles.LowerRow();
 145     while (RowIndex <= NewPoles.UpperRow()){
 146       NewPoles (RowIndex, ColIndex) = Poles (RowIndex, ColIndex);
 147       NewWeights (RowIndex, ColIndex)  =  Weights (RowIndex, ColIndex);
 148       RowIndex++;
 149     }
 150     ColIndex++;
 151   }
 152   RowIndex = NewPoles.LowerRow();
 153   while (RowIndex <= NewPoles.UpperRow()){
 154     NewPoles (RowIndex, ColIndex) = PoleCol (RowIndex - OffsetPol);
 155     NewWeights (RowIndex, ColIndex) =  PoleWeightCol (RowIndex - OffsetW);
 156     RowIndex++;
 157   }
 158   ColIndex++;
 159   while (ColIndex <= NewPoles.UpperCol()) {
 160     RowIndex = NewPoles.LowerRow();
 161     while (RowIndex <= NewPoles.UpperRow()){
 162       NewPoles (RowIndex, ColIndex) = Poles (RowIndex, ColIndex - 1);
 163       RowIndex++;
 164       NewWeights (RowIndex, ColIndex)  = Weights (RowIndex, ColIndex - 1);
 165     }
 166     ColIndex++;
 167   }
 168 }
 169
 170 //=======================================================================
 171 //function : AddPoleRow
 172 //purpose  : Internal use only.
 173 //=======================================================================
 174
 175 static void AddPoleRow
 176   (const TColgp_Array2OfPnt& Poles,
 177    const TColgp_Array1OfPnt& PoleRow,
 178    const Standard_Integer    AfterIndex,
 179          TColgp_Array2OfPnt& NewPoles)
 180 {
 181   Standard_Integer InsertIndex = AfterIndex + NewPoles.LowerRow();
 182   Standard_Integer Offset = NewPoles.LowerCol() - PoleRow.Lower();
 183   Standard_Integer RowIndex = NewPoles.LowerRow();
 184   Standard_Integer ColIndex;
 185   while (RowIndex < InsertIndex) {
 186     ColIndex = NewPoles.LowerCol();
 187     while (ColIndex <= NewPoles.UpperCol()){
 188       NewPoles (RowIndex, ColIndex) = Poles (RowIndex, ColIndex);
 189       ColIndex++;
 190     }
 191     RowIndex++;
 192   }
 193   ColIndex = NewPoles.LowerCol();
 194   while (ColIndex <= NewPoles.UpperCol()){
 195     NewPoles (RowIndex, ColIndex) = PoleRow (ColIndex - Offset);
 196     ColIndex++;
 197   }
 198   RowIndex++;
 199   while (RowIndex <= NewPoles.UpperRow()) {
 200     ColIndex = NewPoles.LowerCol();
 201     while (ColIndex <= NewPoles.UpperCol()){
 202       NewPoles (RowIndex, ColIndex) = Poles (RowIndex - 1, ColIndex);
 203       ColIndex++;
 204     }
 205     RowIndex++;
 206   }
 207 }
 208
 209 //=======================================================================
 210 //function : AddRatPoleRow
 211 //purpose  :
 212 //=======================================================================
 213
 214 static void AddRatPoleRow
 215   (const TColgp_Array2OfPnt&   Poles,
 216    const TColStd_Array2OfReal& Weights,
 217    const TColgp_Array1OfPnt&   PoleRow,
 218    const TColStd_Array1OfReal& PoleWeightRow,
 219    const Standard_Integer      AfterIndex,
 220          TColgp_Array2OfPnt&   NewPoles,
 221          TColStd_Array2OfReal& NewWeights)
 222 {
 223   Standard_Integer InsertIndex = AfterIndex + NewPoles.LowerRow();
 224   Standard_Integer OffsetPol = NewPoles.LowerCol() - PoleRow.Lower();
 225   Standard_Integer OffsetW = NewWeights.LowerCol() - PoleWeightRow.Lower();
 226
 227   Standard_Integer ColIndex;
 228   Standard_Integer RowIndex = NewPoles.LowerRow();
 229   while (RowIndex < InsertIndex) {
 230     ColIndex = NewPoles.LowerCol();
 231     while (ColIndex <= NewPoles.UpperCol()){
 232       NewPoles (RowIndex, ColIndex) =  Poles (RowIndex, ColIndex);
 233       NewWeights (RowIndex, ColIndex)  =  Weights (RowIndex, ColIndex);
 234       ColIndex++;
 235     }
 236     RowIndex++;
 237   }
 238   ColIndex = NewPoles.LowerCol();
 239   while (ColIndex <= NewPoles.UpperCol()){
 240     NewPoles (RowIndex, ColIndex) = PoleRow (ColIndex - OffsetPol);
 241     NewWeights (RowIndex, ColIndex) = PoleWeightRow (ColIndex - OffsetW);
 242     ColIndex++;
 243   }
 244   RowIndex++;
 245   while (RowIndex <= NewPoles.UpperRow()) {
 246     ColIndex = NewPoles.LowerCol();
 247     while (ColIndex <= NewPoles.UpperCol()){
 248       NewPoles (RowIndex, ColIndex) = Poles (RowIndex - 1, ColIndex);
 249       NewWeights (RowIndex, ColIndex)  = Weights (RowIndex - 1, ColIndex);
 250       ColIndex++;
 251     }
 252     RowIndex++;
 253   }
 254 }
 255
 256 //=======================================================================
 257 //function : DeletePoleCol
 258 //purpose  :
 259 //=======================================================================
 260
 261 static void DeletePoleCol
 262   (const TColgp_Array2OfPnt& Poles,
 263    const Standard_Integer    Index,
 264          TColgp_Array2OfPnt& NewPoles)
 265 {
 266   Standard_Integer Offset = 0;
 267   Standard_Integer RowIndex;
 268   Standard_Integer ColIndex = NewPoles.LowerCol();
 269   while (ColIndex <= NewPoles.UpperCol()) {
 270     RowIndex = NewPoles.LowerRow();
 271     if (ColIndex == Index)  Offset = 1;
 272     while (RowIndex <= NewPoles.UpperRow()){
 273       NewPoles (RowIndex, ColIndex) = Poles (RowIndex, ColIndex + Offset);
 274       RowIndex++;
 275     }
 276     ColIndex++;
 277   }
 278 }
 279
 280 //=======================================================================
 281 //function : DeleteRatPoleCol
 282 //purpose  :
 283 //=======================================================================
 284
 285 static void DeleteRatPoleCol
 286   (const TColgp_Array2OfPnt&   Poles,
 287    const TColStd_Array2OfReal& Weights,
 288    const Standard_Integer      Index,
 289          TColgp_Array2OfPnt&   NewPoles,
 290          TColStd_Array2OfReal& NewWeights)
 291 {
 292   Standard_Integer Offset = 0;
 293   Standard_Integer RowIndex;
 294   Standard_Integer ColIndex = NewPoles.LowerCol();
 295   while (ColIndex <= NewPoles.UpperCol()) {
 296     RowIndex = NewPoles.LowerRow();
 297     if (ColIndex == Index) Offset = 1;
 298     while (RowIndex <= NewPoles.UpperRow()){
 299       NewPoles (RowIndex, ColIndex) = Poles (RowIndex, ColIndex + Offset);
 300       NewWeights (RowIndex, ColIndex) = Weights (RowIndex, ColIndex+Offset);
 301       RowIndex++;
 302     }
 303     ColIndex++;
 304   }
 305 }
 306
 307 //=======================================================================
 308 //function : DeletePoleRow
 309 //purpose  :
 310 //=======================================================================
 311
 312 static void DeletePoleRow
 313   (const TColgp_Array2OfPnt& Poles,
 314    const Standard_Integer    Index,
 315          TColgp_Array2OfPnt& NewPoles)
 316 {
 317   Standard_Integer Offset = 0;
 318   Standard_Integer ColIndex;
 319   Standard_Integer RowIndex = NewPoles.LowerRow();
 320   while (RowIndex <= NewPoles.UpperRow()) {
 321     ColIndex = NewPoles.LowerCol();
 322     if (RowIndex == Index)  Offset = 1;
 323     while (ColIndex <= NewPoles.UpperCol()){
 324       NewPoles (RowIndex, ColIndex) = Poles (RowIndex + Offset, ColIndex);
 325       ColIndex++;
 326     }
 327     RowIndex++;
 328   }
 329 }
 330
 331 //=======================================================================
 332 //function : DeleteRatPoleRow
 333 //purpose  :
 334 //=======================================================================
 335
 336 static void DeleteRatPoleRow
 337   (const TColgp_Array2OfPnt&   Poles,
 338    const TColStd_Array2OfReal& Weights,
 339    const Standard_Integer      Index,
 340          TColgp_Array2OfPnt&   NewPoles,
 341          TColStd_Array2OfReal& NewWeights)
 342 {
 343   Standard_Integer Offset = 0;
 344   Standard_Integer ColIndex;
 345   Standard_Integer RowIndex = NewPoles.LowerRow();
 346   while (RowIndex <= NewPoles.UpperRow()) {
 347     ColIndex = NewPoles.LowerCol();
 348     if (RowIndex == Index)  Offset = 1;
 349     while (ColIndex <= NewPoles.UpperCol()){
 350       NewPoles (RowIndex, ColIndex) =  Poles (RowIndex +  Offset, ColIndex);
 351       NewWeights (RowIndex, ColIndex) = Weights (RowIndex+Offset, ColIndex);
 352       ColIndex++;
 353     }
 354     RowIndex++;
 355   }
 356 }
 357
 358 //=======================================================================
 359 //function : Geom_BezierSurface
 360 //purpose  :
 361 //=======================================================================
 362
 363 Geom_BezierSurface::Geom_BezierSurface
 364 (const TColgp_Array2OfPnt& SurfacePoles):
 365  ucacheparameter(0.),
 366  vcacheparameter(0.),
 367  ucachespanlenght(1.),
 368  vcachespanlenght(1.),
 369  validcache(0),
 370  maxderivinvok(Standard_False)
 371 {
 372   Standard_Integer NbUPoles = SurfacePoles.ColLength();
 373   Standard_Integer NbVPoles = SurfacePoles.RowLength();
 374   if (NbUPoles < 2 || NbUPoles > MaxDegree()+1 ||
 375       NbVPoles < 2 || NbVPoles > MaxDegree()+1) {
 376     Standard_ConstructionError::Raise();
 377   }
 378
 379   Handle(TColgp_HArray2OfPnt) npoles =
 380     new TColgp_HArray2OfPnt   (1, NbUPoles, 1, NbVPoles);
 381
 382   urational = 0;
 383   vrational = 0;
 384
 385   npoles->ChangeArray2() = SurfacePoles;
 386
 387   // Init non rational
 388   Init(npoles,
 389        Handle(TColStd_HArray2OfReal)());
 390 }
 391
 392 //=======================================================================
 393 //function : Geom_BezierSurface
 394 //purpose  :
 395 //=======================================================================
 396
 397 Geom_BezierSurface::Geom_BezierSurface
 398   (const TColgp_Array2OfPnt&   SurfacePoles,
 399    const TColStd_Array2OfReal& PoleWeights  ):
 400  ucacheparameter(0.),
 401  vcacheparameter(0.),
 402  ucachespanlenght(1.),
 403  vcachespanlenght(1.),
 404  validcache(0),
 405  maxderivinvok(Standard_False)
 406 {
 407   Standard_Integer NbUPoles = SurfacePoles.ColLength();
 408   Standard_Integer NbVPoles = SurfacePoles.RowLength();
 409   if (NbUPoles < 2 || NbUPoles > MaxDegree()+1 ||
 410       NbVPoles < 2 || NbVPoles > MaxDegree()+1 ||
 411       NbVPoles != PoleWeights.RowLength()      ||
 412       NbUPoles != PoleWeights.ColLength()      ) {
 413     Standard_ConstructionError::Raise();
 414   }
 415
 416   Standard_Integer Row = PoleWeights.LowerRow();
 417   Standard_Integer Col = PoleWeights.LowerCol();
 418   while (Col <= PoleWeights.UpperCol()) {
 419     Row = PoleWeights.LowerRow();
 420     while (Row <= PoleWeights.UpperRow()) {
 421       if (PoleWeights(Row, Col) <= gp::Resolution()) {
 422         Standard_ConstructionError::Raise();
 423       }
 424       Row++;
 425     }
 426     Col++;
 427   }
 428
 429   Handle(TColgp_HArray2OfPnt)
 430     npoles = new TColgp_HArray2OfPnt   (1, NbUPoles, 1, NbVPoles);
 431   npoles->ChangeArray2() = SurfacePoles;
 432
 433   Standard_Integer I, J;
 434   urational = Standard_False;
 435   vrational = Standard_False;
 436   J = PoleWeights.LowerCol ();
 437   while (!vrational && J <= PoleWeights.UpperCol()) {
 438     I = PoleWeights.LowerRow();
 439     while (!vrational && I <= PoleWeights.UpperRow() - 1) {
 440       vrational = (Abs(PoleWeights (I, J) - PoleWeights (I+1, J))
 441                    > Epsilon (Abs(PoleWeights (I, J))));
 442       I++;
 443     }
 444     J++;
 445   }
 446   I = PoleWeights.LowerRow ();
 447   while (!urational && I <= PoleWeights.UpperRow()) {
 448     J = PoleWeights.LowerCol();
 449     while (!urational && J <= PoleWeights.UpperCol() - 1) {
 450       urational = (Abs(PoleWeights (I, J) - PoleWeights (I, J+1))
 451                    > Epsilon (Abs(PoleWeights (I, J))));
 452       J++;
 453     }
 454     I++;
 455   }
 456
 457
 458   Handle(TColStd_HArray2OfReal) nweights;
 459   if (urational || vrational) {
 460     nweights = new TColStd_HArray2OfReal (1, NbUPoles, 1, NbVPoles);
 461     nweights->ChangeArray2() = PoleWeights;
 462   }
 463
 464   // Init
 465   Init(npoles,nweights);
 466 }
 467
 468 //=======================================================================
 469 //function : Geom_BezierSurface
 470 //purpose  :
 471 //=======================================================================
 472
 473 Geom_BezierSurface::Geom_BezierSurface
 474   (const Handle(TColgp_HArray2OfPnt)&   SurfacePoles,
 475    const Handle(TColgp_HArray2OfPnt)&   SurfaceCoefs,
 476    const Handle(TColStd_HArray2OfReal)& PoleWeights,
 477    const Handle(TColStd_HArray2OfReal)& CoefWeights,
 478    const Standard_Boolean               IsURational,
 479    const Standard_Boolean               IsVRational)
 480 :maxderivinvok(Standard_False)
 481 {
 482   urational = IsURational;
 483   vrational = IsVRational;
 484   ucachespanlenght = 1.;
 485   vcachespanlenght = 1.;
 486   validcache = 1;
 487   ucacheparameter =0.;
 488   vcacheparameter = 0.;
 489   Standard_Integer NbUPoles = SurfacePoles->ColLength();
 490   Standard_Integer NbVPoles = SurfacePoles->RowLength();
 491
 492   poles    = new TColgp_HArray2OfPnt   (1,NbUPoles,
 493                                         1,NbVPoles) ;
 494   poles->ChangeArray2() = SurfacePoles->Array2();
 495
 496   coeffs   = new TColgp_HArray2OfPnt   (1,SurfaceCoefs->ColLength(),
 497                                         1,SurfaceCoefs->RowLength()) ;
 498   coeffs->ChangeArray2() = SurfaceCoefs->Array2();
 499
 500   if ( urational || vrational) {
 501     weights  = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
 502     weights->ChangeArray2() = PoleWeights->Array2();
 503
 504     wcoeffs  = new TColStd_HArray2OfReal (1,SurfaceCoefs->ColLength(),
 505                                           1,SurfaceCoefs->RowLength()) ;
 506     wcoeffs->ChangeArray2() = CoefWeights->Array2();
 507   }
 508 }
 509
 510 //=======================================================================
 511 //function : MaxDegree
 512 //purpose  :
 513 //=======================================================================
 514
 515 Standard_Integer Geom_BezierSurface::MaxDegree ()
 516 {
 517   return BSplCLib::MaxDegree();
 518 }
 519
 520 //=======================================================================
 521 //function : ExchangeUV
 522 //purpose  :
 523 //=======================================================================
 524
 525 void Geom_BezierSurface::ExchangeUV ()
 526 {
 527   Standard_Integer LC = poles->LowerCol();
 528   Standard_Integer UC = poles->UpperCol();
 529   Standard_Integer LR = poles->LowerRow();
 530   Standard_Integer UR = poles->UpperRow();
 531
 532   Handle(TColgp_HArray2OfPnt) npoles =
 533     new TColgp_HArray2OfPnt (LC, UC, LR, UR);
 534   Handle(TColStd_HArray2OfReal) nweights =
 535     new TColStd_HArray2OfReal (LC, UC, LR, UR);
 536
 537   const TColgp_Array2OfPnt   & spoles   = poles->Array2();
 538   const TColStd_Array2OfReal & sweights = weights->Array2();
 539
 540   TColgp_Array2OfPnt&   snpoles   = npoles->ChangeArray2();
 541   TColStd_Array2OfReal& snweights = nweights->ChangeArray2();
 542
 543   Standard_Integer i, j;
 544   for (i = LC; i <= UC; i++) {
 545     for (j = LR; j <= UR; j++) {
 546       snpoles   (i,j) = spoles   (j,i);
 547       snweights (i,j) = sweights (j,i);
 548     }
 549   }
 550
 551   poles   = npoles;
 552   weights = nweights;
 553
 554   Standard_Boolean temp = urational;
 555   urational = vrational;
 556   vrational = temp;
 557   coeffs  = new TColgp_HArray2OfPnt   (LC, UC, LR, UR);
 558   wcoeffs = new TColStd_HArray2OfReal (LC, UC, LR, UR);
 559
 560   UpdateCoefficients();
 561 }
 562
 563 //=======================================================================
 564 //function : Increase
 565 //purpose  :
 566 //=======================================================================
 567
 568 void Geom_BezierSurface::Increase (const Standard_Integer UDeg,
 569                                    const Standard_Integer VDeg)
 570 {
 571   if (UDeg < UDegree() || UDeg > Geom_BezierSurface::MaxDegree() ||
 572       VDeg < VDegree() || VDeg > Geom_BezierSurface::MaxDegree() )  {
 573     Standard_ConstructionError::Raise();
 574   }
 575
 576   Standard_Integer oldUDeg = UDegree();
 577   Standard_Integer oldVDeg = VDegree();
 578   Standard_Integer IncUDeg = UDeg - oldUDeg;
 579   Standard_Integer IncVDeg = VDeg - oldVDeg;
 580   if (IncUDeg == 0 && IncVDeg == 0) return;
 581
 582   TColStd_Array1OfReal biduknots(1,2); biduknots(1) = 0.; biduknots(2) = 1.;
 583   TColStd_Array1OfInteger bidumults(1,2); bidumults.Init(UDegree() + 1);
 584   TColStd_Array1OfReal bidvknots(1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
 585   TColStd_Array1OfInteger bidvmults(1,2); bidvmults.Init(VDegree() + 1);
 586   Handle(TColgp_HArray2OfPnt) npoles;
 587   Handle(TColStd_HArray2OfReal) nweights;
 588
 589   if(IncUDeg > 0){
 590     npoles = new TColgp_HArray2OfPnt( 1, UDeg + 1, 1, oldVDeg + 1);
 591
 592     if ( urational || vrational) {
 593       nweights = new TColStd_HArray2OfReal( 1, UDeg + 1, 1, VDegree() + 1);
 594
 595       BSplSLib::IncreaseDegree(1, oldUDeg, UDeg, 0,
 596                                poles->Array2(),
 597                                &weights->Array2(),
 598                                biduknots, bidumults,
 599                                npoles->ChangeArray2(),
 600                                &nweights->ChangeArray2(),
 601                                biduknots, bidumults);
 602       weights = nweights;
 603     }
 604     else {
 605       BSplSLib::IncreaseDegree(1, oldUDeg, UDeg, 0,
 606                                poles->Array2(),
 607                                BSplSLib::NoWeights(),
 608                                biduknots, bidumults,
 609                                npoles->ChangeArray2(),
 610                                BSplSLib::NoWeights(),
 611                                biduknots, bidumults);
 612     }
 613     poles   = npoles;
 614   }
 615   if(IncVDeg > 0){
 616     npoles = new TColgp_HArray2OfPnt( 1, UDeg + 1, 1, VDeg + 1);
 617
 618     if ( urational || vrational) {
 619       nweights = new TColStd_HArray2OfReal( 1, UDeg + 1, 1, VDeg + 1);
 620
 621       BSplSLib::IncreaseDegree(0, oldVDeg, VDeg, 0,
 622                                poles->Array2(),
 623                                &weights->Array2(),
 624                                bidvknots, bidvmults,
 625                                npoles->ChangeArray2(),
 626                                &nweights->ChangeArray2(),
 627                                bidvknots, bidvmults);
 628       weights = nweights;
 629     }
 630     else {
 631       BSplSLib::IncreaseDegree(0, oldVDeg, VDeg, 0,
 632                                poles->Array2(),
 633                                BSplSLib::NoWeights(),
 634                                bidvknots, bidvmults,
 635                                npoles->ChangeArray2(),
 636                                BSplSLib::NoWeights(),
 637                                bidvknots, bidvmults);
 638
 639     }
 640     poles   = npoles;
 641   }
 642   Init(npoles,nweights);
 643 }
 644
 645 //=======================================================================
 646 //function : InsertPoleColAfter
 647 //purpose  :
 648 //=======================================================================
 649
 650 void Geom_BezierSurface::InsertPoleColAfter
 651   (const Standard_Integer    VIndex,
 652    const TColgp_Array1OfPnt& CPoles)
 653 {
 654   const TColgp_Array2OfPnt & Poles = poles->Array2();
 655   if (VIndex < 1 || VIndex > Poles.RowLength())  Standard_OutOfRange::Raise();
 656   if (CPoles.Length() != Poles.ColLength()) {
 657     Standard_ConstructionError::Raise();
 658   }
 659
 660   Handle(TColgp_HArray2OfPnt) npoles =
 661     new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()+1);
 662
 663   Handle(TColStd_HArray2OfReal) nweights;
 664
 665   if (urational || vrational) {
 666     nweights =
 667       new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()+1);
 668
 669     TColStd_Array1OfReal CWeights(nweights->LowerRow(),nweights->UpperRow());
 670     CWeights.Init(1.);
 671
 672     AddRatPoleCol (poles->Array2(), weights->Array2(),
 673                    CPoles, CWeights, VIndex,
 674                    npoles->ChangeArray2(), nweights->ChangeArray2());
 675   }
 676   else {
 677     AddPoleCol (poles->Array2(),
 678                 CPoles, VIndex,
 679                 npoles->ChangeArray2());
 680   }
 681   poles   = npoles;
 682   weights = nweights;
 683   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 684                                     1,poles->RowLength());
 685   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 686                                       1,poles->RowLength());
 687   UpdateCoefficients();
 688 }
 689
 690 //=======================================================================
 691 //function : InsertPoleColAfter
 692 //purpose  :
 693 //=======================================================================
 694
 695 void Geom_BezierSurface::InsertPoleColAfter
 696   (const Standard_Integer      VIndex,
 697    const TColgp_Array1OfPnt&   CPoles,
 698    const TColStd_Array1OfReal& CPoleWeights)
 699 {
 700   const TColgp_Array2OfPnt & Poles = poles->Array2();
 701   if (VIndex < 1 || VIndex > Poles.RowLength())  Standard_OutOfRange::Raise();
 702   if (CPoles.Length() != Poles.ColLength() ||
 703       CPoleWeights.Length() != CPoles.Length()) {
 704     Standard_ConstructionError::Raise();
 705     }
 706   Standard_Integer Index = CPoleWeights.Lower();
 707   while (Index <= CPoleWeights.Upper()) {
 708     if (CPoleWeights (Index) <= gp::Resolution()) {
 709       Standard_ConstructionError::Raise();
 710       }
 711     Index++;
 712   }
 713
 714   Handle(TColgp_HArray2OfPnt) npoles =
 715     new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()+1);
 716
 717   Handle(TColStd_HArray2OfReal) nweights =
 718     new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()+1);
 719
 720   AddRatPoleCol (poles->Array2(), weights->Array2(),
 721                  CPoles, CPoleWeights, VIndex,
 722                  npoles->ChangeArray2(), nweights->ChangeArray2());
 723
 724   poles   = npoles;
 725   weights = nweights;
 726   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 727                                     1,poles->RowLength());
 728   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 729                                       1,poles->RowLength());
 730
 731   Rational(weights->Array2(), urational, vrational);
 732
 733   UpdateCoefficients();
 734 }
 735
 736 //=======================================================================
 737 //function : InsertPoleColBefore
 738 //purpose  :
 739 //=======================================================================
 740
 741 void Geom_BezierSurface::InsertPoleColBefore (const Standard_Integer VIndex,
 742                                               const TColgp_Array1OfPnt& CPoles)
 743 {
 744   InsertPoleColAfter(VIndex - 1, CPoles);
 745 }
 746
 747 //=======================================================================
 748 //function : InsertPoleColBefore
 749 //purpose  :
 750 //=======================================================================
 751
 752 void Geom_BezierSurface::InsertPoleColBefore
 753   (const Standard_Integer      VIndex,
 754    const TColgp_Array1OfPnt&   CPoles,
 755    const TColStd_Array1OfReal& CPoleWeights)
 756 {
 757   InsertPoleColAfter( VIndex - 1, CPoles, CPoleWeights);
 758 }
 759
 760 //=======================================================================
 761 //function : InsertPoleRowAfter
 762 //purpose  :
 763 //=======================================================================
 764
 765 void Geom_BezierSurface::InsertPoleRowAfter (const Standard_Integer UIndex,
 766                                              const TColgp_Array1OfPnt& CPoles)
 767 {
 768   const TColgp_Array2OfPnt & Poles = poles->Array2();
 769   if (UIndex < 1 || UIndex > Poles.ColLength())  Standard_OutOfRange::Raise();
 770   if (CPoles.Length() != Poles.RowLength()) {
 771     Standard_ConstructionError::Raise();
 772   }
 773
 774   Handle(TColgp_HArray2OfPnt) npoles =
 775     new TColgp_HArray2OfPnt(1,poles->ColLength()+1,1,poles->RowLength());
 776
 777   Handle(TColStd_HArray2OfReal) nweights;
 778
 779   if (urational || vrational) {
 780     nweights =
 781       new TColStd_HArray2OfReal(1,poles->ColLength()+1,1,poles->RowLength());
 782
 783 //    TColStd_Array1OfReal CWeights(nweights->LowerCol(),nweights->UpperCol(),
 784 //                                1.0); ???????????
 785     TColStd_Array1OfReal CWeights(1.0,
 786                                   nweights->LowerCol(),nweights->UpperCol());
 787
 788     AddRatPoleRow (poles->Array2(), weights->Array2(),
 789                    CPoles, CWeights, UIndex,
 790                    npoles->ChangeArray2(), nweights->ChangeArray2());
 791   }
 792   else {
 793     AddPoleRow (poles->Array2(),
 794                 CPoles, UIndex,
 795                 npoles->ChangeArray2());
 796   }
 797   poles   = npoles;
 798   weights = nweights;
 799   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 800                                     1,poles->RowLength());
 801   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 802                                       1,poles->RowLength());
 803
 804   UpdateCoefficients();
 805 }
 806
 807 //=======================================================================
 808 //function : InsertPoleRowAfter
 809 //purpose  :
 810 //=======================================================================
 811
 812 void Geom_BezierSurface::InsertPoleRowAfter
 813   (const Standard_Integer      UIndex,
 814    const TColgp_Array1OfPnt&   CPoles,
 815    const TColStd_Array1OfReal& CPoleWeights)
 816 {
 817   const TColgp_Array2OfPnt & Poles = poles->Array2();
 818   if (UIndex < 1 || UIndex > Poles.ColLength())  Standard_OutOfRange::Raise();
 819   if (CPoles.Length() != Poles.RowLength() ||
 820       CPoleWeights.Length() != CPoles.Length()) {
 821     Standard_ConstructionError::Raise();
 822   }
 823   Standard_Integer Index = CPoleWeights.Lower();
 824   while (Index <= CPoleWeights.Upper()) {
 825     if (CPoleWeights(Index) <= gp::Resolution()) {
 826       Standard_ConstructionError::Raise();
 827     }
 828     Index++;
 829   }
 830
 831   Handle(TColgp_HArray2OfPnt) npoles =
 832     new TColgp_HArray2OfPnt(1,poles->ColLength()+1,1,poles->RowLength());
 833
 834   Handle(TColStd_HArray2OfReal) nweights =
 835     new TColStd_HArray2OfReal(1,poles->ColLength()+1,1,poles->RowLength());
 836
 837   AddRatPoleCol (poles->Array2(), weights->Array2(),
 838                  CPoles, CPoleWeights, UIndex,
 839                  npoles->ChangeArray2(), nweights->ChangeArray2());
 840
 841   poles   = npoles;
 842   weights = nweights;
 843   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 844                                     1,poles->RowLength());
 845   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 846                                       1,poles->RowLength());
 847
 848   Rational(weights->Array2(), urational, vrational);
 849
 850   UpdateCoefficients();
 851 }
 852
 853 //=======================================================================
 854 //function : InsertPoleRowBefore
 855 //purpose  :
 856 //=======================================================================
 857
 858 void Geom_BezierSurface::InsertPoleRowBefore (const Standard_Integer UIndex,
 859                                               const TColgp_Array1OfPnt& CPoles)
 860 {
 861   InsertPoleRowAfter( UIndex - 1, CPoles);
 862 }
 863
 864 //=======================================================================
 865 //function : InsertPoleRowBefore
 866 //purpose  :
 867 //=======================================================================
 868
 869 void Geom_BezierSurface::InsertPoleRowBefore
 870   (const Standard_Integer      UIndex,
 871    const TColgp_Array1OfPnt&   CPoles,
 872    const TColStd_Array1OfReal& CPoleWeights)
 873 {
 874   InsertPoleRowAfter( UIndex - 1, CPoles, CPoleWeights);
 875 }
 876
 877 //=======================================================================
 878 //function : RemovePoleCol
 879 //purpose  :
 880 //=======================================================================
 881
 882 void Geom_BezierSurface::RemovePoleCol (const Standard_Integer VIndex)
 883 {
 884   const TColgp_Array2OfPnt & Poles = poles->Array2();
 885   if (VIndex < 1 || VIndex > Poles.RowLength())  Standard_OutOfRange::Raise();
 886   if (Poles.RowLength() <= 2)             Standard_ConstructionError::Raise();
 887
 888   Handle(TColgp_HArray2OfPnt) npoles =
 889     new TColgp_HArray2OfPnt(1,poles->ColLength(),1,poles->RowLength()-1);
 890
 891   Handle(TColStd_HArray2OfReal) nweights;
 892
 893   if (urational || vrational) {
 894     nweights =
 895       new TColStd_HArray2OfReal(1,poles->ColLength(),1,poles->RowLength()-1);
 896
 897     DeleteRatPoleCol (poles->Array2(), weights->Array2(),
 898                       VIndex,
 899                       npoles->ChangeArray2(), nweights->ChangeArray2());
 900     // Mise a jour de la rationalite
 901     Rational(nweights->Array2(), urational, vrational);
 902   }
 903   else {
 904     DeletePoleCol (poles->Array2(),
 905                    VIndex,
 906                    npoles->ChangeArray2());
 907   }
 908   poles   = npoles;
 909   weights = nweights;
 910   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 911                                     1,poles->RowLength());
 912   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 913                                       1,poles->RowLength());
 914   UpdateCoefficients();
 915 }
 916
 917 //=======================================================================
 918 //function : RemovePoleRow
 919 //purpose  :
 920 //=======================================================================
 921
 922 void Geom_BezierSurface::RemovePoleRow (const Standard_Integer UIndex)
 923 {
 924   const TColgp_Array2OfPnt & Poles = poles->Array2();
 925   if (UIndex < 1 || UIndex > Poles.ColLength()) Standard_OutOfRange::Raise();
 926   if (Poles.ColLength() <= 2)            Standard_ConstructionError::Raise();
 927
 928   Handle(TColgp_HArray2OfPnt) npoles =
 929     new TColgp_HArray2OfPnt(1,poles->ColLength()-1,1,poles->RowLength());
 930
 931   Handle(TColStd_HArray2OfReal) nweights;
 932
 933   if (urational || vrational) {
 934     nweights =
 935       new TColStd_HArray2OfReal(1,poles->ColLength()-1,1,poles->RowLength());
 936
 937     DeleteRatPoleRow (poles->Array2(), weights->Array2(),
 938                       UIndex,
 939                       npoles->ChangeArray2(), nweights->ChangeArray2());
 940
 941     // Mise a jour de la rationalite
 942     Rational(nweights->Array2(), urational, vrational);
 943   }
 944   else {
 945     DeletePoleRow (poles->Array2(),
 946                    UIndex,
 947                    npoles->ChangeArray2());
 948   }
 949   poles   = npoles;
 950   weights = nweights;
 951   coeffs  = new TColgp_HArray2OfPnt(1,poles->ColLength(),
 952                                     1,poles->RowLength());
 953   wcoeffs = new TColStd_HArray2OfReal(1,poles->ColLength(),
 954                                       1,poles->RowLength());
 955   UpdateCoefficients();
 956 }
 957
 958 //=======================================================================
 959 //function : Segment
 960 //purpose  :
 961 //=======================================================================
 962
 963 void Geom_BezierSurface::Segment
 964   (const Standard_Real U1,
 965    const Standard_Real U2,
 966    const Standard_Real V1,
 967    const Standard_Real V2)
 968 {
 969   Standard_Boolean rat = (urational || vrational);
 970   Handle(TColgp_HArray2OfPnt)  Coefs;
 971   Handle(TColStd_HArray2OfReal) WCoefs;
 972
 973   if (validcache == 0) UpdateCoefficients(0., 0.);
 974
 975   // Attention si udeg <= vdeg u et v sont intervertis
 976   // dans les coeffs, il faut donc tout transposer.
 977   if(UDegree() <= VDegree()) {
 978     Standard_Integer ii, jj;
 979     Coefs = new  (TColgp_HArray2OfPnt)(1,UDegree()+1,1,VDegree()+1);
 980     if (rat) {
 981       WCoefs = new  (TColStd_HArray2OfReal)(1,UDegree()+1,1,VDegree()+1);
 982     }
 983     for (ii=1; ii<=UDegree()+1; ii++)
 984       for (jj=1; jj<=VDegree()+1; jj++) {
 985         Coefs->SetValue(ii, jj, coeffs->Value(jj,ii));
 986         if (rat)  WCoefs->SetValue(ii, jj, wcoeffs->Value(jj,ii));
 987       }
 988   }
 989   else {
 990    Coefs = coeffs;
 991    if (rat)  {WCoefs =  wcoeffs;}
 992  }
 993
 994 // Trim dans la base cannonique et Update des Poles et Coeffs
 995
 996 // PMN : tranfo sur les parametres
 997   Standard_Real ufirst = 2*(U1 - 0.5),
 998                 ulast =  2*(U2 - 0.5),
 999                 vfirst = 2*(V1 - 0.5),
1000                 vlast  = 2*(V2 - 0.5);
1001   if (rat) {
1002     PLib::UTrimming (ufirst, ulast, Coefs->ChangeArray2(),
1003                      &WCoefs->ChangeArray2());
1004     PLib::VTrimming (vfirst, vlast, Coefs->ChangeArray2(),
1005                      &WCoefs->ChangeArray2());
1006     PLib::CoefficientsPoles(Coefs->Array2(),
1007                             &WCoefs->Array2(),
1008                             poles->ChangeArray2(),
1009                             &weights->ChangeArray2());
1010   }
1011   else {
1012     PLib::UTrimming (ufirst, ulast, Coefs->ChangeArray2(), PLib::NoWeights2());
1013     PLib::VTrimming (vfirst, vlast, Coefs->ChangeArray2(), PLib::NoWeights2());
1014     PLib::CoefficientsPoles (Coefs->Array2(), PLib::NoWeights2(),
1015                                                poles->ChangeArray2(), PLib::NoWeights2());
1016   }
1017   UpdateCoefficients();
1018 }
1019
1020 //=======================================================================
1021 //function : SetPole
1022 //purpose  :
1023 //=======================================================================
1024
1025 void Geom_BezierSurface::SetPole
1026   (const Standard_Integer UIndex,
1027    const Standard_Integer VIndex,
1028    const gp_Pnt&          P)
1029 {
1030   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1031   if (UIndex < 1                 ||
1032       UIndex > Poles.ColLength() ||
1033       VIndex < 1                 ||
1034       VIndex > Poles.RowLength() )    Standard_OutOfRange::Raise();
1035
1036   Poles (UIndex, VIndex) = P;
1037   UpdateCoefficients();
1038 }
1039
1040 //=======================================================================
1041 //function : SetPole
1042 //purpose  :
1043 //=======================================================================
1044
1045 void Geom_BezierSurface::SetPole
1046   (const Standard_Integer UIndex,
1047    const Standard_Integer VIndex,
1048    const gp_Pnt&          P,
1049    const Standard_Real    Weight)
1050 {
1051
1052   if (Weight <= gp::Resolution())
1053     Standard_ConstructionError::Raise("Geom_BezierSurface::SetPole");
1054   if (UIndex < 1                 ||
1055       UIndex > poles->ColLength() ||
1056       VIndex < 1                 ||
1057       VIndex > poles->RowLength())
1058     Standard_OutOfRange::Raise("Geom_BezierSurface::SetPole");
1059
1060   poles->SetValue(UIndex, VIndex, P);
1061
1062   SetWeight(UIndex, VIndex,  Weight); //L'update des coeff est fait la dedans
1063 }
1064
1065 //=======================================================================
1066 //function : SetPoleCol
1067 //purpose  :
1068 //=======================================================================
1069
1070 void Geom_BezierSurface::SetPoleCol
1071   (const Standard_Integer      VIndex,
1072    const TColgp_Array1OfPnt&   CPoles,
1073    const TColStd_Array1OfReal& CPoleWeights)
1074 {
1075   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1076   if (VIndex < 1 || VIndex > Poles.RowLength()) Standard_OutOfRange::Raise();
1077
1078   if (CPoles.Lower() < 1                     ||
1079       CPoles.Lower() > Poles.ColLength()     ||
1080       CPoles.Upper() < 1                     ||
1081       CPoles.Upper() > Poles.ColLength()     ||
1082       CPoleWeights.Lower() != CPoles.Lower() ||
1083       CPoleWeights.Upper() != CPoles.Upper()) {
1084     Standard_ConstructionError::Raise();
1085   }
1086
1087   Standard_Integer I;
1088   for (I = CPoles.Lower();  I <= CPoles.Upper(); I++) {
1089     Poles (I, VIndex) = CPoles (I);
1090   }
1091   SetWeightCol(VIndex, CPoleWeights); //Avec l'update
1092 }
1093
1094 //=======================================================================
1095 //function : SetPoleCol
1096 //purpose  :
1097 //=======================================================================
1098
1099 void Geom_BezierSurface::SetPoleCol (const Standard_Integer      VIndex,
1100                                      const TColgp_Array1OfPnt&   CPoles)
1101 {
1102   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1103   if (VIndex < 1 || VIndex > Poles.RowLength())  Standard_OutOfRange::Raise();
1104
1105   if (CPoles.Lower() < 1                 ||
1106       CPoles.Lower() > Poles.ColLength() ||
1107       CPoles.Upper() < 1                 ||
1108       CPoles.Upper() > Poles.ColLength()) {
1109     Standard_ConstructionError::Raise();
1110   }
1111   for (Standard_Integer I = CPoles.Lower(); I <= CPoles.Upper(); I++) {
1112     Poles (I, VIndex) = CPoles (I);
1113   }
1114
1115   UpdateCoefficients();
1116 }
1117
1118 //=======================================================================
1119 //function : SetPoleRow
1120 //purpose  :
1121 //=======================================================================
1122
1123 void Geom_BezierSurface::SetPoleRow (const Standard_Integer    UIndex,
1124                                      const TColgp_Array1OfPnt& CPoles)
1125 {
1126   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1127   if (UIndex < 1 || UIndex > Poles.ColLength())  Standard_OutOfRange::Raise();
1128
1129   if (CPoles.Lower() < 1                 ||
1130       CPoles.Lower() > Poles.RowLength() ||
1131       CPoles.Upper() < 1                 ||
1132       CPoles.Upper() > Poles.RowLength())  Standard_ConstructionError::Raise();
1133
1134   for (Standard_Integer I = CPoles.Lower(); I <= CPoles.Upper(); I++) {
1135     Poles (UIndex, I) = CPoles (I);
1136   }
1137   UpdateCoefficients();
1138 }
1139
1140 //=======================================================================
1141 //function : SetPoleRow
1142 //purpose  :
1143 //=======================================================================
1144
1145 void Geom_BezierSurface::SetPoleRow
1146   (const Standard_Integer      UIndex,
1147    const TColgp_Array1OfPnt&   CPoles,
1148    const TColStd_Array1OfReal& CPoleWeights)
1149 {
1150   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1151   if (UIndex < 1 || UIndex > Poles.ColLength())   Standard_OutOfRange::Raise();
1152
1153   if (CPoles.Lower() < 1                     ||
1154       CPoles.Lower() > Poles.RowLength()     ||
1155       CPoles.Upper() < 1                     ||
1156       CPoles.Upper() > Poles.RowLength()     ||
1157       CPoleWeights.Lower() != CPoles.Lower() ||
1158       CPoleWeights.Upper() != CPoles.Upper()) {
1159     Standard_ConstructionError::Raise();
1160   }
1161
1162   Standard_Integer I;
1163
1164   for (I = CPoles.Lower(); I <= CPoles.Upper(); I++) {
1165     Poles   (UIndex, I) = CPoles (I);
1166   }
1167
1168   SetWeightRow(UIndex, CPoleWeights); //Avec l'update
1169 }
1170
1171 //=======================================================================
1172 //function : SetWeight
1173 //purpose  :
1174 //=======================================================================
1175
1176 void Geom_BezierSurface::SetWeight (const Standard_Integer UIndex,
1177                                     const Standard_Integer VIndex,
1178                                     const Standard_Real    Weight)
1179 {
1180    // compute new rationality
1181   Standard_Boolean wasrat = (urational||vrational);
1182   if (!wasrat) {
1183     // a weight of 1. does not turn to rational
1184     if (Abs(Weight - 1.) <= gp::Resolution()) {
1185       UpdateCoefficients(); //Pour l'appel via SetPole
1186       return;
1187     }
1188
1189     // set weights of 1.
1190     weights = new TColStd_HArray2OfReal (1, poles->ColLength(),
1191                                          1, poles->RowLength(), 1.);
1192     wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(),
1193                                          1, poles->RowLength());
1194   }
1195
1196   TColStd_Array2OfReal & Weights = weights->ChangeArray2();
1197   if (Weight <= gp::Resolution())
1198     Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeight");
1199
1200   if (UIndex < 1                   ||
1201       UIndex > Weights.ColLength() ||
1202       VIndex < 1                   ||
1203       VIndex > Weights.RowLength())    Standard_OutOfRange::Raise();
1204
1205   if (Abs (Weight - Weights (UIndex, VIndex)) > gp::Resolution()) {
1206     Weights (UIndex, VIndex) = Weight;
1207     Rational(Weights, urational, vrational);
1208   }
1209
1210  // is it turning into non rational
1211   if (wasrat) {
1212     if (!(urational || vrational)) {
1213       weights.Nullify();
1214       wcoeffs.Nullify();
1215     }
1216   }
1217
1218   UpdateCoefficients(); //Dans tous cas :Attention a SetPoleCol !
1219 }
1220
1221 //=======================================================================
1222 //function : SetWeightCol
1223 //purpose  :
1224 //=======================================================================
1225
1226 void Geom_BezierSurface::SetWeightCol
1227   (const Standard_Integer      VIndex,
1228    const TColStd_Array1OfReal& CPoleWeights)
1229 {
1230   Standard_Integer I;
1231    // compute new rationality
1232   Standard_Boolean wasrat = (urational||vrational);
1233   if (!wasrat) {
1234     // set weights of 1.
1235     weights = new TColStd_HArray2OfReal (1, poles->ColLength(),
1236                                          1, poles->RowLength(), 1.);
1237     wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(),
1238                                          1, poles->RowLength());
1239   }
1240
1241   TColStd_Array2OfReal & Weights = weights->ChangeArray2();
1242   if (VIndex < 1 || VIndex > Weights.RowLength()) Standard_OutOfRange::Raise();
1243
1244   if (CPoleWeights.Length() !=  Weights.ColLength())  {
1245     Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeightCol");
1246   }
1247
1248   I = CPoleWeights.Lower();
1249   while (I <= CPoleWeights.Upper()) {
1250     if (CPoleWeights(I) <= gp::Resolution()) {
1251       Standard_ConstructionError::Raise();
1252     }
1253     Weights (I, VIndex) = CPoleWeights (I);
1254     I++;
1255   }
1256
1257  Rational(Weights, urational, vrational);
1258
1259  // is it turning into non rational
1260   if (wasrat) {
1261     if (!(urational || vrational)) {
1262       weights.Nullify();
1263       wcoeffs.Nullify();
1264     }
1265   }
1266
1267   UpdateCoefficients();
1268 }
1269
1270 //=======================================================================
1271 //function : SetWeightRow
1272 //purpose  :
1273 //=======================================================================
1274
1275 void Geom_BezierSurface::SetWeightRow
1276   (const Standard_Integer      UIndex,
1277    const TColStd_Array1OfReal& CPoleWeights)
1278 {
1279   Standard_Integer I;
1280    // compute new rationality
1281   Standard_Boolean wasrat = (urational||vrational);
1282   if (!wasrat) {
1283     // set weights of 1.
1284     weights = new TColStd_HArray2OfReal (1, poles->ColLength(),
1285                                          1, poles->RowLength(), 1.);
1286     wcoeffs = new TColStd_HArray2OfReal (1, poles->ColLength(),
1287                                          1, poles->RowLength());
1288   }
1289
1290   TColStd_Array2OfReal & Weights = weights->ChangeArray2();
1291   if (UIndex < 1 || UIndex > Weights.ColLength())
1292     Standard_OutOfRange::Raise("Geom_BezierSurface::SetWeightRow");
1293   if (CPoleWeights.Lower() < 1 ||
1294       CPoleWeights.Lower() > Weights.RowLength() ||
1295       CPoleWeights.Upper() < 1 ||
1296       CPoleWeights.Upper() > Weights.RowLength()  ) {
1297     Standard_ConstructionError::Raise("Geom_BezierSurface::SetWeightRow");
1298   }
1299
1300   I = CPoleWeights.Lower();
1301   while (I <= CPoleWeights.Upper()) {
1302     if (CPoleWeights(I) <= gp::Resolution())  {
1303       Standard_ConstructionError::Raise();
1304     }
1305     Weights (UIndex, I) = CPoleWeights (I);
1306     I++;
1307   }
1308
1309   Rational(Weights, urational, vrational);
1310
1311  // is it turning into non rational
1312   if (wasrat) {
1313     if (!(urational || vrational)) {
1314       weights.Nullify();
1315       wcoeffs.Nullify();
1316     }
1317   }
1318
1319   UpdateCoefficients();
1320 }
1321
1322 //=======================================================================
1323 //function : UReverse
1324 //purpose  :
1325 //=======================================================================
1326
1327 void Geom_BezierSurface::UReverse ()
1328 {
1329   gp_Pnt Pol;
1330   Standard_Integer Row,Col;
1331   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1332   if (urational || vrational) {
1333     TColStd_Array2OfReal & Weights = weights->ChangeArray2();
1334     Standard_Real W;
1335     for (Col = 1; Col <= Poles.RowLength(); Col++) {
1336       for (Row = 1; Row <= IntegerPart (Poles.ColLength() / 2); Row++) {
1337         W = Weights (Row, Col);
1338         Weights (Row, Col) = Weights (Poles.ColLength()- Row + 1, Col);
1339         Weights (Poles.ColLength() - Row + 1, Col) = W;
1340         Pol = Poles (Row, Col);
1341         Poles (Row, Col) = Poles (Poles.ColLength() - Row + 1, Col);
1342         Poles (Poles.ColLength() - Row + 1, Col) = Pol;
1343       }
1344     }
1345   }
1346   else {
1347     for (Col = 1; Col <= Poles.RowLength(); Col++) {
1348       for (Row = 1; Row <= IntegerPart (Poles.ColLength() / 2); Row++) {
1349         Pol = Poles (Row, Col);
1350         Poles (Row, Col) = Poles (Poles.ColLength() - Row + 1, Col);
1351         Poles (Poles.ColLength() - Row + 1, Col) = Pol;
1352       }
1353     }
1354   }
1355   UpdateCoefficients();
1356 }
1357
1358 //=======================================================================
1359 //function : UReversedParameter
1360 //purpose  :
1361 //=======================================================================
1362
1363 Standard_Real Geom_BezierSurface::UReversedParameter
1364   ( const Standard_Real U) const
1365 {
1366   return ( 1. - U);
1367 }
1368
1369 //=======================================================================
1370 //function : VReverse
1371 //purpose  :
1372 //=======================================================================
1373
1374 void Geom_BezierSurface::VReverse ()
1375 {
1376   gp_Pnt Pol;
1377   Standard_Integer Row,Col;
1378   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1379   if (urational || vrational) {
1380     TColStd_Array2OfReal & Weights = weights->ChangeArray2();
1381     Standard_Real W;
1382     for (Row = 1; Row <= Poles.ColLength(); Row++) {
1383       for (Col = 1; Col <= IntegerPart (Poles.RowLength()/2); Col++) {
1384         W = Weights (Row, Col);
1385         Weights (Row, Col) = Weights (Row, Poles.RowLength() - Col + 1);
1386         Weights (Row, Poles.RowLength() - Col + 1) = W;
1387         Pol = Poles (Row, Col);
1388         Poles (Row, Col) = Poles (Row, Poles.RowLength() - Col + 1);
1389         Poles (Row, Poles.RowLength() - Col + 1) = Pol;
1390       }
1391     }
1392   }
1393   else {
1394     for (Row = 1; Row <= Poles.ColLength(); Row++) {
1395       for (Col = 1; Col <= IntegerPart(Poles.RowLength()/2); Col++) {
1396         Pol = Poles (Row, Col);
1397         Poles (Row, Col)= Poles (Row, Poles.RowLength() - Col + 1);
1398         Poles (Row, Poles.RowLength() - Col + 1) = Pol;
1399       }
1400     }
1401   }
1402   UpdateCoefficients();
1403 }
1404
1405 //=======================================================================
1406 //function : VReversedParameter
1407 //purpose  :
1408 //=======================================================================
1409
1410 Standard_Real Geom_BezierSurface::VReversedParameter
1411   ( const Standard_Real V) const
1412 {
1413   return ( 1. - V);
1414 }
1415
1416 //=======================================================================
1417 //function : Bounds
1418 //purpose  :
1419 //=======================================================================
1420
1421 void Geom_BezierSurface::Bounds (Standard_Real& U1,
1422                                  Standard_Real& U2,
1423                                  Standard_Real& V1,
1424                                  Standard_Real& V2) const
1425 {
1426   U1 = 0.0;
1427   U2 = 1.0;
1428   V1 = 0.0;
1429   V2 = 1.0;
1430 }
1431
1432 //=======================================================================
1433 //function : Continuity
1434 //purpose  :
1435 //=======================================================================
1436
1437 GeomAbs_Shape Geom_BezierSurface::Continuity () const
1438 {
1439   return GeomAbs_CN;
1440 }
1441
1442 //=======================================================================
1443 //function : D0
1444 //purpose  :
1445 //=======================================================================
1446
1447 void Geom_BezierSurface::D0 (const Standard_Real U,
1448                              const Standard_Real V,
1449                                    gp_Pnt&       P ) const
1450 {
1451
1452   if (validcache == 1) {
1453     //
1454     // XAB : cet algorithme devient instable pour les hauts degres
1455     // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache
1456     //       sur [-1,1].
1457     //
1458     Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2,
1459                   uspanlenght_11 = ucachespanlenght/2,
1460                   vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2,
1461                   vspanlenght_11 = vcachespanlenght/2 ;
1462     if (urational || vrational) {
1463       BSplSLib::CacheD0(U, V, UDegree(), VDegree(),
1464                         uparameter_11, vparameter_11,
1465                         uspanlenght_11, vspanlenght_11,
1466                         coeffs->Array2(),
1467                         &wcoeffs->Array2(),
1468                         P);
1469     }
1470     else {
1471       BSplSLib::CacheD0(U, V, UDegree(), VDegree(),
1472                         uparameter_11, vparameter_11,
1473                         uspanlenght_11, vspanlenght_11,
1474                         coeffs->Array2(),
1475                         BSplSLib::NoWeights(),
1476                         P);
1477     }
1478   }
1479   else {
1480     Standard_Real array_u[2] ;
1481     Standard_Real array_v[2] ;
1482     Standard_Integer mult_u[2] ;
1483     Standard_Integer mult_v[2] ;
1484     TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1485     TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1);
1486     TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1487     TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1);
1488     if (urational || vrational) {
1489         BSplSLib::D0(U, V, 1,1,poles->Array2(),
1490                      &weights->Array2(),
1491                      biduknots,bidvknots,&bidumults,&bidvmults,
1492                      UDegree(),VDegree(),
1493                      urational,vrational,Standard_False,Standard_False,
1494                      P) ;
1495     }
1496     else {
1497
1498       BSplSLib::D0(U, V, 1,1,poles->Array2(),
1499                    BSplSLib::NoWeights(),
1500                    biduknots,bidvknots,&bidumults,&bidvmults,
1501                    UDegree(),VDegree(),
1502                    urational,vrational,Standard_False,Standard_False,
1503                    P) ;
1504     }
1505   }
1506 }
1507
1508 //=======================================================================
1509 //function : D1
1510 //purpose  :
1511 //=======================================================================
1512
1513 void Geom_BezierSurface::D1
1514   (const Standard_Real U,
1515    const Standard_Real V,
1516          gp_Pnt&       P,
1517          gp_Vec&       D1U,
1518          gp_Vec&       D1V ) const
1519 {
1520
1521   if (validcache == 1) {
1522     //
1523     // XAB : cet algorithme devient instable pour les hauts degres
1524     // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache
1525     //       sur [-1,1].
1526     //
1527     Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2,
1528                   uspanlenght_11 = ucachespanlenght/2,
1529                   vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2,
1530                   vspanlenght_11 = vcachespanlenght/2 ;
1531    if (urational || vrational) {
1532      BSplSLib::CacheD1(U, V, UDegree(), VDegree(),
1533                        uparameter_11, vparameter_11,
1534                        uspanlenght_11, vspanlenght_11,
1535                        coeffs->Array2(),
1536                        &wcoeffs->Array2(),
1537                        P, D1U, D1V);
1538    }
1539    else {
1540      BSplSLib::CacheD1(U, V, UDegree(), VDegree(),
1541                        uparameter_11, vparameter_11,
1542                        uspanlenght_11, vspanlenght_11,
1543                        coeffs->Array2(),
1544                        BSplSLib::NoWeights(),
1545                        P, D1U, D1V);
1546    }
1547   }
1548   else {
1549     Standard_Real array_u[2] ;
1550     Standard_Real array_v[2] ;
1551     Standard_Integer mult_u[2] ;
1552     Standard_Integer mult_v[2] ;
1553     TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1554     TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1);
1555     TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1556     TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1);
1557     if (urational || vrational) {
1558       BSplSLib::D1(U, V, 1,1,poles->Array2(),
1559                    &weights->Array2(),
1560                    biduknots,bidvknots,&bidumults,&bidvmults,
1561                    UDegree(),VDegree(),
1562                    urational,vrational,Standard_False,Standard_False,
1563                    P,D1U, D1V) ;
1564     }
1565     else {
1566       BSplSLib::D1(U, V, 1,1,poles->Array2(),
1567                    BSplSLib::NoWeights(),
1568                    biduknots,bidvknots,&bidumults,&bidvmults,
1569                    UDegree(),VDegree(),
1570                    urational,vrational,Standard_False,Standard_False,
1571                    P,D1U, D1V) ;
1572     }
1573   }
1574 }
1575
1576 //=======================================================================
1577 //function : D2
1578 //purpose  :
1579 //=======================================================================
1580
1581 void Geom_BezierSurface::D2
1582   (const Standard_Real U,
1583    const Standard_Real V,
1584          gp_Pnt&       P,
1585          gp_Vec&       D1U, gp_Vec& D1V,
1586          gp_Vec&       D2U, gp_Vec& D2V, gp_Vec& D2UV ) const
1587 {
1588
1589   if (validcache == 1) {
1590     //
1591     // XAB : cet algorithme devient instable pour les hauts degres
1592     // RBD : Beaucoup moins maintenant avec le calcul d'un nouveau cache
1593     //       sur [-1,1].
1594     //
1595     Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2,
1596                   uspanlenght_11 = ucachespanlenght/2,
1597                   vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2,
1598                   vspanlenght_11 = vcachespanlenght/2 ;
1599     if (urational || vrational) {
1600       //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB
1601       BSplSLib::CacheD2(U, V, UDegree(), VDegree(),
1602                         uparameter_11, vparameter_11,
1603                         uspanlenght_11, vspanlenght_11,
1604                         coeffs->Array2(),
1605                         &wcoeffs->Array2(),
1606                         P, D1U, D1V, D2U, D2UV , D2V);
1607     }
1608     else {
1609       //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB
1610       BSplSLib::CacheD2(U, V, UDegree(), VDegree(),
1611                         uparameter_11, vparameter_11,
1612                         uspanlenght_11, vspanlenght_11,
1613                         coeffs->Array2(),
1614                         BSplSLib::NoWeights(),
1615                         P, D1U, D1V, D2U, D2UV , D2V);
1616     }
1617   }
1618   else {
1619     Standard_Real array_u[2] ;
1620     Standard_Real array_v[2] ;
1621     Standard_Integer mult_u[2] ;
1622     Standard_Integer mult_v[2] ;
1623     TColStd_Array1OfReal biduknots(array_u[0],1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1624     TColStd_Array1OfInteger bidumults(mult_u[0],1,2); bidumults.Init(UDegree() + 1);
1625     TColStd_Array1OfReal bidvknots(array_v[0],1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1626     TColStd_Array1OfInteger bidvmults(mult_v[0],1,2); bidvmults.Init(VDegree() + 1);
1627     if (urational || vrational) {
1628       //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB
1629       BSplSLib::D2(U, V, 1,1,poles->Array2(),
1630                    &weights->Array2(),
1631                    biduknots,bidvknots,&bidumults,&bidvmults,
1632                    UDegree(),VDegree(),
1633                    urational,vrational,Standard_False,Standard_False,
1634                    P,D1U, D1V, D2U, D2V , D2UV) ;
1635     }
1636     else {
1637       //-- ATTENTION a l'ORDRE d'appel ds BSPLSLIB
1638       BSplSLib::D2(U, V, 1,1,poles->Array2(),
1639                    BSplSLib::NoWeights(),
1640                    biduknots,bidvknots,&bidumults,&bidvmults,
1641                    UDegree(),VDegree(),
1642                    urational,vrational,Standard_False,Standard_False,
1643                    P,D1U, D1V, D2U, D2V, D2UV ) ;
1644     }
1645   }
1646 }
1647
1648 //=======================================================================
1649 //function : D3
1650 //purpose  :
1651 //=======================================================================
1652
1653 void Geom_BezierSurface::D3
1654   (const Standard_Real U, const Standard_Real V,
1655    gp_Pnt& P,
1656    gp_Vec& D1U, gp_Vec& D1V,
1657    gp_Vec& D2U, gp_Vec& D2V, gp_Vec& D2UV,
1658    gp_Vec& D3U, gp_Vec& D3V, gp_Vec& D3UUV, gp_Vec& D3UVV) const
1659 {
1660   TColStd_Array1OfReal biduknots(1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1661   TColStd_Array1OfInteger bidumults(1,2); bidumults.Init(UDegree() + 1);
1662   TColStd_Array1OfReal bidvknots(1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1663   TColStd_Array1OfInteger bidvmults(1,2); bidvmults.Init(VDegree() + 1);
1664   if (urational || vrational) {
1665     BSplSLib::D3 (U, V, 0, 0, poles->Array2(),
1666                   &weights->Array2(),
1667                   biduknots, bidvknots, &bidumults, &bidvmults,
1668                   UDegree(), VDegree(), urational, vrational, 0, 0,
1669                   P,
1670                   D1U, D1V,
1671                   D2U, D2V, D2UV,
1672                   D3U, D3V, D3UUV, D3UVV);
1673   }
1674   else {
1675     BSplSLib::D3 (U, V, 0, 0, poles->Array2(),
1676                   BSplSLib::NoWeights(),
1677                   biduknots, bidvknots, &bidumults, &bidvmults,
1678                   UDegree(), VDegree(), urational, vrational, 0, 0,
1679                   P,
1680                   D1U, D1V,
1681                   D2U, D2V, D2UV,
1682                   D3U, D3V, D3UUV, D3UVV);
1683   }
1684 }
1685
1686 //=======================================================================
1687 //function : DN
1688 //purpose  :
1689 //=======================================================================
1690
1691 gp_Vec Geom_BezierSurface::DN
1692   (const Standard_Real    U,
1693    const Standard_Real    V,
1694    const Standard_Integer Nu,
1695    const Standard_Integer Nv) const
1696 {
1697   Standard_RangeError_Raise_if (Nu + Nv < 1 || Nv < 0 || Nu <0, " ");
1698   gp_Vec Derivative;
1699   TColStd_Array1OfReal biduknots(1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1700   TColStd_Array1OfInteger bidumults(1,2); bidumults.Init(UDegree() + 1);
1701   TColStd_Array1OfReal bidvknots(1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1702   TColStd_Array1OfInteger bidvmults(1,2); bidvmults.Init(VDegree() + 1);
1703   if (urational || vrational) {
1704     BSplSLib::DN (U, V, Nu, Nv, 0, 0, poles->Array2(),
1705                   &weights->Array2(),
1706                   biduknots, bidvknots, &bidumults, &bidvmults,
1707                   UDegree(), VDegree(), urational, vrational, 0, 0,
1708                   Derivative);
1709   }
1710   else {
1711     BSplSLib::DN (U, V, Nu, Nv, 0, 0, poles->Array2(),
1712                   BSplSLib::NoWeights(),
1713                   biduknots, bidvknots, &bidumults, &bidvmults,
1714                   UDegree(), VDegree(), urational, vrational, 0, 0,
1715                   Derivative);
1716   }
1717   return Derivative;
1718 }
1719
1720 //=======================================================================
1721 //function : NbUPoles
1722 //purpose  :
1723 //=======================================================================
1724
1725 Standard_Integer Geom_BezierSurface::NbUPoles () const
1726 {
1727   return poles->ColLength();
1728 }
1729
1730 //=======================================================================
1731 //function : NbVPoles
1732 //purpose  :
1733 //=======================================================================
1734
1735 Standard_Integer Geom_BezierSurface::NbVPoles () const
1736 {
1737   return poles->RowLength();
1738 }
1739
1740 //=======================================================================
1741 //function : Pole
1742 //purpose  :
1743 //=======================================================================
1744
1745 gp_Pnt Geom_BezierSurface::Pole (const Standard_Integer UIndex,
1746                                  const Standard_Integer VIndex) const
1747 {
1748   Standard_OutOfRange_Raise_if
1749     (UIndex < 1 || UIndex > poles->ColLength() ||
1750      VIndex < 1 || VIndex > poles->RowLength(), " ");
1751   return poles->Value (UIndex + poles->LowerRow() - 1,
1752                        VIndex + poles->LowerCol() - 1);
1753 }
1754
1755 //=======================================================================
1756 //function : Poles
1757 //purpose  :
1758 //=======================================================================
1759
1760 void Geom_BezierSurface::Poles (TColgp_Array2OfPnt& P) const
1761 {
1762   Standard_DimensionError_Raise_if
1763     (P.RowLength() != poles->RowLength() ||
1764      P.ColLength() != poles->ColLength(), " ");
1765   P = poles->Array2();
1766 }
1767
1768 //=======================================================================
1769 //function : UDegree
1770 //purpose  :
1771 //=======================================================================
1772
1773 Standard_Integer Geom_BezierSurface::UDegree () const
1774 {
1775   return poles->ColLength() - 1;
1776 }
1777
1778 //=======================================================================
1779 //function : UIso
1780 //purpose  :
1781 //=======================================================================
1782
1783 Handle(Geom_Curve) Geom_BezierSurface::UIso (const Standard_Real U) const
1784 {
1785   TColStd_Array1OfReal biduknots(1,2); biduknots(1) = 0.; biduknots(2) = 1.;
1786   TColStd_Array1OfInteger bidumults(1,2); bidumults.Init(UDegree() + 1);
1787
1788   Handle(Geom_BezierCurve) UIsoCurve;
1789   const TColgp_Array2OfPnt & Poles = poles->Array2();
1790   TColgp_Array1OfPnt VCurvePoles (Poles.LowerCol() , Poles.UpperCol());
1791   if (urational || vrational) {
1792     const TColStd_Array2OfReal & Weights = weights->Array2();
1793     TColStd_Array1OfReal VCurveWeights
1794       (Weights.LowerCol() , Weights.UpperCol());
1795     BSplSLib::Iso (U, 1, Poles,
1796                    &Weights,
1797                    biduknots, &bidumults,
1798                    UDegree(), 0, VCurvePoles, &VCurveWeights);
1799     if (urational)
1800       UIsoCurve = new Geom_BezierCurve (VCurvePoles, VCurveWeights);
1801     else
1802       UIsoCurve = new Geom_BezierCurve (VCurvePoles);
1803   }
1804   else {
1805     BSplSLib::Iso (U, 1, Poles,
1806                    BSplSLib::NoWeights(),
1807                    biduknots, &bidumults,
1808                    UDegree(), 0, VCurvePoles, PLib::NoWeights());
1809     UIsoCurve = new Geom_BezierCurve (VCurvePoles);
1810   }
1811   return UIsoCurve;
1812 }
1813
1814 //=======================================================================
1815 //function : VDegree
1816 //purpose  :
1817 //=======================================================================
1818
1819 Standard_Integer Geom_BezierSurface::VDegree () const
1820 {
1821   return poles->RowLength() - 1;
1822 }
1823
1824 //=======================================================================
1825 //function : VIso
1826 //purpose  :
1827 //=======================================================================
1828
1829 Handle(Geom_Curve) Geom_BezierSurface::VIso (const Standard_Real V) const
1830 {
1831   TColStd_Array1OfReal bidvknots(1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
1832   TColStd_Array1OfInteger bidvmults(1,2); bidvmults.Init(VDegree() + 1);
1833
1834   Handle(Geom_BezierCurve) VIsoCurve;
1835   const TColgp_Array2OfPnt & Poles = poles->Array2();
1836   TColgp_Array1OfPnt VCurvePoles (Poles.LowerRow() , Poles.UpperRow());
1837   if (vrational || urational) {
1838     const TColStd_Array2OfReal & Weights = weights->Array2();
1839     TColStd_Array1OfReal VCurveWeights
1840       (Weights.LowerRow() , Weights.UpperRow());
1841     BSplSLib::Iso (V, 0, Poles,
1842                    &Weights,
1843                    bidvknots, &bidvmults,
1844                    VDegree(), 0, VCurvePoles, &VCurveWeights);
1845     if (vrational)
1846       VIsoCurve = new Geom_BezierCurve (VCurvePoles, VCurveWeights);
1847     else
1848       VIsoCurve = new Geom_BezierCurve (VCurvePoles);
1849   }
1850   else {
1851     BSplSLib::Iso (V, 0, Poles,
1852                    BSplSLib::NoWeights(),
1853                    bidvknots, &bidvmults,
1854                    VDegree(), 0, VCurvePoles, PLib::NoWeights());
1855     VIsoCurve = new Geom_BezierCurve (VCurvePoles);
1856   }
1857   return VIsoCurve;
1858 }
1859
1860 //=======================================================================
1861 //function : Weight
1862 //purpose  :
1863 //=======================================================================
1864
1865 Standard_Real Geom_BezierSurface::Weight (const Standard_Integer UIndex,
1866                                           const Standard_Integer VIndex) const
1867 {
1868   Standard_OutOfRange_Raise_if (
1869            UIndex < 1 || UIndex > weights->ColLength() ||
1870            VIndex < 1 || VIndex > weights->RowLength(), " ");
1871
1872   if (urational || vrational)
1873     return weights->Value (UIndex, VIndex);
1874   else
1875     return 1;
1876 }
1877
1878 //=======================================================================
1879 //function : Weights
1880 //purpose  :
1881 //=======================================================================
1882
1883 void Geom_BezierSurface::Weights (TColStd_Array2OfReal& W ) const
1884 {
1885   Standard_DimensionError_Raise_if (
1886            W.RowLength() != weights->RowLength() ||
1887            W.ColLength() != weights->ColLength(), " " );
1888   if (urational || vrational)
1889     W = weights->Array2();
1890   else
1891     W.Init(1.);
1892 }
1893
1894 //=======================================================================
1895 //function : IsCNu
1896 //purpose  :
1897 //=======================================================================
1898
1899 Standard_Boolean Geom_BezierSurface::IsCNu (const Standard_Integer ) const
1900 {
1901   return Standard_True;
1902 }
1903
1904 //=======================================================================
1905 //function : IsCNv
1906 //purpose  :
1907 //=======================================================================
1908
1909 Standard_Boolean Geom_BezierSurface::IsCNv (const Standard_Integer ) const
1910 {
1911   return Standard_True;
1912 }
1913
1914 //=======================================================================
1915 //function : IsURational
1916 //purpose  :
1917 //=======================================================================
1918
1919 Standard_Boolean Geom_BezierSurface::IsURational () const
1920 {
1921   return urational;
1922 }
1923
1924 //=======================================================================
1925 //function : IsVRational
1926 //purpose  :
1927 //=======================================================================
1928
1929 Standard_Boolean Geom_BezierSurface::IsVRational () const
1930 {
1931   return vrational;
1932 }
1933
1934 //=======================================================================
1935 //function : Transform
1936 //purpose  :
1937 //=======================================================================
1938
1939 void Geom_BezierSurface::Transform (const gp_Trsf& T)
1940 {
1941   TColgp_Array2OfPnt & Poles = poles->ChangeArray2();
1942
1943   for (Standard_Integer I = 1; I <= Poles.ColLength(); I++) {
1944
1945     for (Standard_Integer J = 1; J <= Poles.RowLength(); J++) {
1946       Poles (I, J).Transform (T);
1947     }
1948   }
1949   UpdateCoefficients();
1950 }
1951
1952 //=======================================================================
1953 //function : IsUClosed
1954 //purpose  :
1955 //=======================================================================
1956
1957 Standard_Boolean Geom_BezierSurface::IsUClosed () const
1958 {
1959   const TColgp_Array2OfPnt & Poles = poles->Array2();
1960   Standard_Boolean Closed = Standard_True;
1961   Standard_Integer Lower  = Poles.LowerRow();
1962   Standard_Integer Upper  = Poles.UpperRow();
1963   Standard_Integer Length = Poles.RowLength();
1964   Standard_Integer j      = Poles.LowerCol();
1965
1966   while (Closed && j <= Length) {
1967     Closed = (Poles (Lower,j).Distance (Poles (Upper,j)) <= Precision::Confusion());
1968     j++;
1969   }
1970   return Closed;
1971 }
1972
1973 //=======================================================================
1974 //function : IsVClosed
1975 //purpose  :
1976 //=======================================================================
1977
1978 Standard_Boolean Geom_BezierSurface::IsVClosed () const
1979 {
1980   const TColgp_Array2OfPnt & Poles = poles->Array2();
1981   Standard_Boolean Closed = Standard_True;
1982   Standard_Integer Lower  = Poles.LowerCol();
1983   Standard_Integer Upper  = Poles.UpperCol();
1984   Standard_Integer Length = Poles.ColLength();
1985   Standard_Integer i      = Poles.LowerRow();
1986   while (Closed && i <= Length) {
1987     Closed = (Poles (i,Lower).Distance (Poles (i,Upper)) <= Precision::Confusion());
1988     i++;
1989   }
1990   return Closed;
1991 }
1992
1993 //=======================================================================
1994 //function : IsUPeriodic
1995 //purpose  :
1996 //=======================================================================
1997
1998 Standard_Boolean Geom_BezierSurface::IsUPeriodic () const
1999 {
2000   return Standard_False;
2001 }
2002
2003 //=======================================================================
2004 //function : IsVPeriodic
2005 //purpose  :
2006 //=======================================================================
2007
2008 Standard_Boolean Geom_BezierSurface::IsVPeriodic () const
2009 {
2010   return Standard_False;
2011 }
2012
2013 //=======================================================================
2014 //function : Resolution
2015 //purpose  :
2016 //=======================================================================
2017
2018 void Geom_BezierSurface::Resolution(const Standard_Real  Tolerance3D,
2019                                     Standard_Real&       UTolerance,
2020                                     Standard_Real&       VTolerance)
2021 {
2022   if(!maxderivinvok){
2023     TColStd_Array1OfReal biduknots(1,2); biduknots(1) = 0.; biduknots(2) = 1.;
2024     TColStd_Array1OfInteger bidumults(1,2); bidumults.Init(UDegree() + 1);
2025     TColStd_Array1OfReal bidvknots(1,2); bidvknots(1) = 0.; bidvknots(2) = 1.;
2026     TColStd_Array1OfInteger bidvmults(1,2); bidvmults.Init(VDegree() + 1);
2027     if(urational || vrational){
2028       BSplSLib::Resolution(poles->Array2(),
2029                            &weights->Array2(),
2030                            biduknots,
2031                            bidvknots,
2032                            bidumults,
2033                            bidvmults,
2034                            UDegree(),
2035                            VDegree(),
2036                            urational,
2037                            vrational,
2038                            0,
2039                            0,
2040                            1.,
2041                            umaxderivinv,
2042                            vmaxderivinv) ;
2043     }
2044     else{
2045       BSplSLib::Resolution(poles->Array2(),
2046                            BSplSLib::NoWeights(),
2047                            biduknots,
2048                            bidvknots,
2049                            bidumults,
2050                            bidvmults,
2051                            UDegree(),
2052                            VDegree(),
2053                            urational,
2054                            vrational,
2055                            0,
2056                            0,
2057                            1.,
2058                            umaxderivinv,
2059                            vmaxderivinv) ;
2060     }
2061     maxderivinvok = 1;
2062   }
2063   UTolerance = Tolerance3D * umaxderivinv;
2064   VTolerance = Tolerance3D * vmaxderivinv;
2065 }
2066
2067 //=======================================================================
2068 //function : Copy
2069 //purpose  :
2070 //=======================================================================
2071
2072 Handle(Geom_Geometry) Geom_BezierSurface::Copy() const
2073 {
2074   Handle(Geom_BezierSurface) S = new Geom_BezierSurface
2075     (poles, coeffs, weights, wcoeffs, urational, vrational);
2076   return S;
2077 }
2078
2079 //=======================================================================
2080 //function : Init
2081 //purpose  :
2082 //=======================================================================
2083
2084 void Geom_BezierSurface::Init
2085   (const Handle(TColgp_HArray2OfPnt)&   Poles,
2086    const Handle(TColStd_HArray2OfReal)& Weights)
2087 {
2088   Standard_Integer NbUPoles = Poles->ColLength();
2089   Standard_Integer NbVPoles = Poles->RowLength();
2090
2091   Standard_Integer maxcls = Max(NbUPoles, NbVPoles);
2092   Standard_Integer mincls = Min(NbUPoles, NbVPoles);
2093
2094   // set fields
2095   poles   = Poles;
2096   coeffs  = new TColgp_HArray2OfPnt  (1,maxcls,1,mincls);
2097
2098   if (urational || vrational) {
2099     weights = Weights;
2100     wcoeffs = new TColStd_HArray2OfReal (1,maxcls,1,mincls);
2101   }
2102   else {
2103     weights.Nullify();
2104     wcoeffs.Nullify();
2105   }
2106
2107   UpdateCoefficients();
2108 }
2109
2110
2111 //=======================================================================
2112 //function : UpdateCoefficients
2113 //purpose  :
2114 //=======================================================================
2115
2116 void  Geom_BezierSurface::UpdateCoefficients(const Standard_Real ,
2117                                              const Standard_Real )
2118 {
2119   maxderivinvok = Standard_False;
2120   ucacheparameter = 0.;
2121   TColStd_Array1OfReal biduflatknots(BSplCLib::FlatBezierKnots(UDegree()),
2122                                      1, 2*(UDegree()+1));
2123   vcacheparameter = 0.;
2124   TColStd_Array1OfReal bidvflatknots(BSplCLib::FlatBezierKnots(VDegree()),
2125                                      1, 2*(VDegree()+1));
2126
2127   Standard_Real uparameter_11 = (2*ucacheparameter + ucachespanlenght)/2,
2128                 uspanlenght_11 = ucachespanlenght/2,
2129                 vparameter_11 = (2*vcacheparameter + vcachespanlenght)/2,
2130                 vspanlenght_11 = vcachespanlenght/2 ;
2131
2132   if ( urational || vrational ) {
2133     BSplSLib::BuildCache(uparameter_11,vparameter_11,
2134                          uspanlenght_11,vspanlenght_11,0,0,
2135                          UDegree(),VDegree(),0,0,
2136                          biduflatknots,bidvflatknots,
2137                          poles->Array2(),
2138                          &weights->Array2(),
2139                          coeffs->ChangeArray2(),
2140                          &wcoeffs->ChangeArray2());
2141   }
2142   else {
2143     BSplSLib::BuildCache(uparameter_11,vparameter_11,
2144                          uspanlenght_11,vspanlenght_11,0,0,
2145                          UDegree(),VDegree(),0,0,
2146                          biduflatknots,bidvflatknots,
2147                          poles->Array2(),
2148                          BSplSLib::NoWeights(),
2149                          coeffs->ChangeArray2(),
2150                          BSplSLib::NoWeights());
2151   }
2152   validcache = 1;
2153 }
2154