src/BndLib/BndLib_AddSurface.cxx

   1 // Created on: 1995-07-24
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503
  18
  19 #include <Adaptor3d_Surface.hxx>
  20 #include <Adaptor3d_Surface.hxx>
  21 #include <Bnd_Box.hxx>
  22 #include <BndLib.hxx>
  23 #include <BndLib_AddSurface.hxx>
  24 #include <ElSLib.hxx>
  25 #include <ElCLib.hxx>
  26 #include <Geom_BezierSurface.hxx>
  27 #include <Geom_BSplineSurface.hxx>
  28 #include <GeomAbs_SurfaceType.hxx>
  29 #include <gp_Pln.hxx>
  30 #include <gp_Pnt.hxx>
  31 #include <gp_Cone.hxx>
  32 #include <Precision.hxx>
  33 #include <TColgp_Array2OfPnt.hxx>
  34 #include <TColStd_Array1OfInteger.hxx>
  35 #include <TColStd_Array1OfReal.hxx>
  36 #include <math_PSO.hxx>
  37 #include <math_Powell.hxx>
  38 //
  39 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S,
  40                                    const Standard_Real Umin,
  41                                    const Standard_Real Umax);
  42 //
  43 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S,
  44                                    const Standard_Real Vmin,
  45                                    const Standard_Real Vmax);
  46 //
  47 static Standard_Real  AdjustExtr(const Adaptor3d_Surface& S,
  48                                  const Standard_Real UMin,
  49                                                    const Standard_Real UMax,
  50                                                    const Standard_Real VMin,
  51                                                    const Standard_Real VMax,
  52                                  const Standard_Real Extr0,
  53                                  const Standard_Integer CoordIndx,
  54                                  const Standard_Real Tol,
  55                                  const Standard_Boolean IsMin);
  56
  57
  58 static void ComputePolesIndexes(const TColStd_Array1OfReal &theKnots,
  59   const TColStd_Array1OfInteger &theMults,
  60   const Standard_Integer theDegree,
  61   const Standard_Real theMin,
  62   const Standard_Real theMax,
  63   const Standard_Integer theMaxPoleIdx,
  64   const Standard_Boolean theIsPeriodic,
  65   Standard_Integer &theOutMinIdx,
  66   Standard_Integer &theOutMaxIdx);
  67
  68 //=======================================================================
  69 //function : Add
  70 //purpose  :
  71 //=======================================================================
  72 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
  73                             const Standard_Real Tol,
  74                             Bnd_Box& B )
  75 {
  76
  77   BndLib_AddSurface::Add(S,
  78                          S.FirstUParameter(),
  79                          S.LastUParameter (),
  80                          S.FirstVParameter(),
  81                          S.LastVParameter (),Tol,B);
  82 }
  83 //=======================================================================
  84 //function : NbUSamples
  85 //purpose  :
  86 //=======================================================================
  87
  88 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S)
  89 {
  90   Standard_Integer N;
  91   GeomAbs_SurfaceType Type = S.GetType();
  92   switch (Type) {
  93   case GeomAbs_BezierSurface:
  94     {
  95       N = 2*S.NbUPoles();
  96       break;
  97     }
  98   case GeomAbs_BSplineSurface:
  99     {
 100       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
 101       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
 102       break;
 103     }
 104   default:
 105     N = 33;
 106   }
 107   return Min (50,N);
 108 }
 109
 110 //=======================================================================
 111 //function : NbVSamples
 112 //purpose  :
 113 //=======================================================================
 114
 115 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S)
 116 {
 117   Standard_Integer N;
 118   GeomAbs_SurfaceType Type = S.GetType();
 119   switch (Type) {
 120   case GeomAbs_BezierSurface:
 121     {
 122       N = 2*S.NbVPoles();
 123       break;
 124     }
 125   case GeomAbs_BSplineSurface:
 126     {
 127       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
 128       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
 129       break;
 130     }
 131   default:
 132     N = 33;
 133   }
 134   return Min(50,N);
 135 }
 136
 137 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
 138 static gp_Pnt BaryCenter(const gp_Pln        &aPlane,
 139                          const Standard_Real  aUMin,
 140                          const Standard_Real  aUMax,
 141                          const Standard_Real  aVMin,
 142                          const Standard_Real  aVMax)
 143 {
 144   Standard_Real aU, aV;
 145   Standard_Boolean isU1Inf = Precision::IsInfinite(aUMin);
 146   Standard_Boolean isU2Inf = Precision::IsInfinite(aUMax);
 147   Standard_Boolean isV1Inf = Precision::IsInfinite(aVMin);
 148   Standard_Boolean isV2Inf = Precision::IsInfinite(aVMax);
 149
 150   if (isU1Inf && isU2Inf)
 151     aU = 0;
 152   else if (isU1Inf)
 153     aU = aUMax - 10.;
 154   else if (isU2Inf)
 155     aU = aUMin + 10.;
 156   else
 157     aU = (aUMin + aUMax)/2.;
 158
 159   if (isV1Inf && isV2Inf)
 160     aV = 0;
 161   else if (isV1Inf)
 162     aV = aVMax - 10.;
 163   else if (isV2Inf)
 164     aV = aVMin + 10.;
 165   else
 166     aV = (aVMin + aVMax)/2.;
 167
 168   gp_Pnt aCenter = ElSLib::Value(aU, aV, aPlane);
 169
 170   return aCenter;
 171 }
 172
 173 static void TreatInfinitePlane(const gp_Pln        &aPlane,
 174                                const Standard_Real  aUMin,
 175                                const Standard_Real  aUMax,
 176                                const Standard_Real  aVMin,
 177                                const Standard_Real  aVMax,
 178                                const Standard_Real  aTol,
 179                                      Bnd_Box       &aB)
 180 {
 181   // Get 3 coordinate axes of the plane.
 182   const gp_Dir        &aNorm        = aPlane.Axis().Direction();
 183   const Standard_Real  anAngularTol = RealEpsilon();
 184
 185   // Get location of the plane as its barycenter
 186   gp_Pnt aLocation = BaryCenter(aPlane, aUMin, aUMax, aVMin, aVMax);
 187
 188   if (aNorm.IsParallel(gp::DX(), anAngularTol)) {
 189     aB.Add(aLocation);
 190     aB.OpenYmin();
 191     aB.OpenYmax();
 192     aB.OpenZmin();
 193     aB.OpenZmax();
 194   } else if (aNorm.IsParallel(gp::DY(), anAngularTol)) {
 195     aB.Add(aLocation);
 196     aB.OpenXmin();
 197     aB.OpenXmax();
 198     aB.OpenZmin();
 199     aB.OpenZmax();
 200   } else if (aNorm.IsParallel(gp::DZ(), anAngularTol)) {
 201     aB.Add(aLocation);
 202     aB.OpenXmin();
 203     aB.OpenXmax();
 204     aB.OpenYmin();
 205     aB.OpenYmax();
 206   } else {
 207     aB.SetWhole();
 208     return;
 209   }
 210
 211   aB.Enlarge(aTol);
 212 }
 213
 214 // Compute start and finish indexes used in convex hull.
 215 // theMinIdx - minimum poles index, that can be used.
 216 // theMaxIdx - maximum poles index, that can be used.
 217 // theShiftCoeff - shift between flatknots array and poles array.
 218 // This vaule should be equal to 1 in case of non periodic BSpline,
 219 // and (degree + 1) - mults(the lowest index).
 220
 221 void ComputePolesIndexes(const TColStd_Array1OfReal &theKnots,
 222   const TColStd_Array1OfInteger &theMults,
 223   const Standard_Integer theDegree,
 224   const Standard_Real theMin,
 225   const Standard_Real theMax,
 226   const Standard_Integer theMaxPoleIdx,
 227   const Standard_Boolean theIsPeriodic,
 228   Standard_Integer &theOutMinIdx,
 229   Standard_Integer &theOutMaxIdx)
 230 {
 231   BSplCLib::Hunt(theKnots, theMin, theOutMinIdx);
 232   theOutMinIdx = Max(theOutMinIdx, theKnots.Lower());
 233
 234   BSplCLib::Hunt(theKnots, theMax, theOutMaxIdx);
 235   theOutMaxIdx++;
 236   theOutMaxIdx = Min(theOutMaxIdx, theKnots.Upper());
 237   Standard_Integer mult = theMults(theOutMaxIdx);
 238
 239   theOutMinIdx = BSplCLib::PoleIndex(theDegree, theOutMinIdx, theIsPeriodic, theMults) + 1;
 240   theOutMinIdx = Max(theOutMinIdx, 1);
 241   theOutMaxIdx = BSplCLib::PoleIndex(theDegree, theOutMaxIdx, theIsPeriodic, theMults) + 1;
 242   theOutMaxIdx += theDegree - mult;
 243   if (!theIsPeriodic)
 244     theOutMaxIdx = Min(theOutMaxIdx, theMaxPoleIdx);
 245 }
 246
 247 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
 248 //=======================================================================
 249 //function : Add
 250 //purpose  :
 251 //=======================================================================
 252 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
 253                             const Standard_Real UMin,
 254                             const Standard_Real UMax,
 255                             const Standard_Real VMin,
 256                             const Standard_Real VMax,
 257                             const Standard_Real Tol,
 258                             Bnd_Box& B )
 259 {
 260   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
 261
 262   if (Precision::IsInfinite(VMin) ||
 263       Precision::IsInfinite(VMax) ||
 264       Precision::IsInfinite(UMin) ||
 265       Precision::IsInfinite(UMax)   ) {
 266 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
 267 //     B.SetWhole();
 268 //     return;
 269     switch (Type) {
 270     case GeomAbs_Plane:
 271       {
 272         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
 273         return;
 274       }
 275     default:
 276       {
 277         B.SetWhole();
 278         return;
 279       }
 280     }
 281 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
 282   }
 283
 284 //   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
 285
 286   switch (Type) {
 287
 288   case GeomAbs_Plane:
 289     {
 290       gp_Pln Plan = S.Plane();
 291       B.Add(ElSLib::Value(UMin,VMin,Plan));
 292       B.Add(ElSLib::Value(UMin,VMax,Plan));
 293       B.Add(ElSLib::Value(UMax,VMin,Plan));
 294       B.Add(ElSLib::Value(UMax,VMax,Plan));
 295       B.Enlarge(Tol);
 296       break;
 297     }
 298   case GeomAbs_Cylinder:
 299     {
 300       BndLib::Add(S.Cylinder(),UMin,UMax,VMin,VMax,Tol,B);
 301       break;
 302     }
 303   case GeomAbs_Cone:
 304     {
 305       BndLib::Add(S.Cone(),UMin,UMax,VMin,VMax,Tol,B);
 306       break;
 307     }
 308   case GeomAbs_Torus:
 309     {
 310       BndLib::Add(S.Torus(),UMin,UMax,VMin,VMax,Tol,B);
 311       break;
 312     }
 313   case GeomAbs_Sphere:
 314     {
 315       if (Abs(UMin) < Precision::Angular() &&
 316           Abs(UMax - 2.*M_PI) < Precision::Angular() &&
 317           Abs(VMin + M_PI/2.) < Precision::Angular() &&
 318           Abs(VMax - M_PI/2.) < Precision::Angular()) // a whole sphere
 319         BndLib::Add(S.Sphere(),Tol,B);
 320       else
 321         BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
 322       break;
 323     }
 324   case GeomAbs_OffsetSurface:
 325     {
 326       Handle(Adaptor3d_Surface) HS = S.BasisSurface();
 327       Add (*HS,UMin,UMax,VMin,VMax,Tol,B);
 328       B.Enlarge(S.OffsetValue());
 329       B.Enlarge(Tol);
 330       break;
 331     }
 332   case GeomAbs_BezierSurface:
 333   case GeomAbs_BSplineSurface:
 334     {
 335       Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
 336       Standard_Real PTol = Precision::Parametric(Precision::Confusion());
 337       // Borders of underlying geometry.
 338       Standard_Real anUMinParam = UMin, anUMaxParam = UMax,// BSpline case.
 339                      aVMinParam = VMin,  aVMaxParam = VMax;
 340       Handle(Geom_BSplineSurface) aBS;
 341       if (Type == GeomAbs_BezierSurface)
 342       {
 343         // Bezier surface:
 344         // All of poles used for any parameter,
 345         // that's why in case of trimmed parameters handled by grid algorithm.
 346
 347         if (Abs(UMin-S.FirstUParameter()) > PTol ||
 348             Abs(VMin-S.FirstVParameter()) > PTol ||
 349             Abs(UMax-S.LastUParameter ()) > PTol ||
 350             Abs(VMax-S.LastVParameter ()) > PTol )
 351         {
 352           // Borders not equal to topology borders.
 353           isUseConvexHullAlgorithm = Standard_False;
 354         }
 355       }
 356       else
 357       {
 358         // BSpline:
 359         // If Umin, Vmin, Umax, Vmax lies inside geometry bounds then:
 360         // use convex hull algorithm,
 361         // if Umin, VMin, Umax, Vmax lies outside then:
 362         // use grid algorithm on analytic continuation (default case).
 363         aBS = S.BSpline();
 364         aBS->Bounds(anUMinParam, anUMaxParam, aVMinParam, aVMaxParam);
 365         if ( (UMin - anUMinParam) < -PTol ||
 366              (VMin -  aVMinParam) < -PTol ||
 367              (UMax - anUMaxParam) >  PTol ||
 368              (VMax -  aVMaxParam) >  PTol )
 369         {
 370           // Out of geometry borders.
 371           isUseConvexHullAlgorithm = Standard_False;
 372         }
 373       }
 374
 375       if (isUseConvexHullAlgorithm)
 376       {
 377         Standard_Integer aNbUPoles = S.NbUPoles(), aNbVPoles = S.NbVPoles();
 378         TColgp_Array2OfPnt Tp(1, aNbUPoles, 1, aNbVPoles);
 379         Standard_Integer UMinIdx = 0, UMaxIdx = 0;
 380         Standard_Integer VMinIdx = 0, VMaxIdx = 0;
 381         Standard_Boolean isUPeriodic = S.IsUPeriodic(), isVPeriodic = S.IsVPeriodic();
 382         if (Type == GeomAbs_BezierSurface)
 383         {
 384           S.Bezier()->Poles(Tp);
 385           UMinIdx = 1; UMaxIdx = aNbUPoles;
 386           VMinIdx = 1; VMaxIdx = aNbVPoles;
 387         }
 388         else
 389         {
 390           aBS->Poles(Tp);
 391
 392           UMinIdx = 1;
 393           UMaxIdx = aNbUPoles;
 394           VMinIdx = 1;
 395           VMaxIdx = aNbVPoles;
 396
 397           if (UMin > anUMinParam ||
 398               UMax < anUMaxParam)
 399           {
 400             TColStd_Array1OfInteger aMults(1, aBS->NbUKnots());
 401             TColStd_Array1OfReal aKnots(1, aBS->NbUKnots());
 402             aBS->UKnots(aKnots);
 403             aBS->UMultiplicities(aMults);
 404
 405             ComputePolesIndexes(aKnots,
 406               aMults,
 407               aBS->UDegree(),
 408               UMin, UMax,
 409               aNbUPoles,
 410               isUPeriodic,
 411               UMinIdx, UMaxIdx); // the Output indexes
 412
 413           }
 414
 415           if (VMin > aVMinParam ||
 416             VMax < aVMaxParam)
 417           {
 418             TColStd_Array1OfInteger aMults(1, aBS->NbVKnots());
 419             TColStd_Array1OfReal aKnots(1, aBS->NbVKnots());
 420             aBS->VKnots(aKnots);
 421             aBS->VMultiplicities(aMults);
 422
 423             ComputePolesIndexes(aKnots,
 424               aMults,
 425               aBS->VDegree(),
 426               VMin, VMax,
 427               aNbVPoles,
 428               isVPeriodic,
 429               VMinIdx, VMaxIdx); // the Output indexes
 430           }
 431
 432         }
 433
 434         // Use poles to build convex hull.
 435         Standard_Integer ip, jp;
 436         for (Standard_Integer i = UMinIdx; i <= UMaxIdx; i++)
 437         {
 438           ip = i;
 439           if (isUPeriodic && ip > aNbUPoles)
 440           {
 441             ip = ip - aNbUPoles;
 442           }
 443           for (Standard_Integer j = VMinIdx; j <= VMaxIdx; j++)
 444           {
 445             jp = j;
 446             if (isVPeriodic && jp > aNbVPoles)
 447             {
 448               jp = jp - aNbVPoles;
 449             }
 450             B.Add(Tp(ip, jp));
 451           }
 452         }
 453
 454         B.Enlarge(Tol);
 455         break;
 456       }
 457   }
 458     Standard_FALLTHROUGH
 459   default:
 460     {
 461       Standard_Integer Nu = NbUSamples(S);
 462       Standard_Integer Nv = NbVSamples(S);
 463       gp_Pnt P;
 464       for (Standard_Integer i =1 ;i<=Nu;i++){
 465         Standard_Real U = UMin + ((UMax-UMin)*(i-1)/(Nu-1));
 466         for (Standard_Integer j=1 ;j<=Nv;j++){
 467           Standard_Real V = VMin + ((VMax-VMin)*(j-1)/(Nv-1));
 468           S.D0(U,V,P);
 469           B.Add(P);
 470         }
 471       }
 472       B.Enlarge(Tol);
 473     }
 474   }
 475 }
 476 //----- Methods for AddOptimal ---------------------------------------
 477
 478 //=======================================================================
 479 //function : AddOptimal
 480 //purpose  :
 481 //=======================================================================
 482 void BndLib_AddSurface::AddOptimal(const Adaptor3d_Surface& S,
 483                                                      const Standard_Real Tol,
 484                                                      Bnd_Box& B )
 485 {
 486
 487   BndLib_AddSurface::AddOptimal(S,
 488                                                   S.FirstUParameter(),
 489                                                   S.LastUParameter (),
 490                                                   S.FirstVParameter(),
 491                                                   S.LastVParameter (),Tol,B);
 492 }
 493 //=======================================================================
 494 //function : AddOptimal
 495 //purpose  :
 496 //=======================================================================
 497
 498 void BndLib_AddSurface::AddOptimal(const Adaptor3d_Surface& S,
 499                             const Standard_Real UMin,
 500                             const Standard_Real UMax,
 501                             const Standard_Real VMin,
 502                             const Standard_Real VMax,
 503                             const Standard_Real Tol,
 504                             Bnd_Box& B )
 505 {
 506   GeomAbs_SurfaceType Type = S.GetType();
 507
 508   if (Precision::IsInfinite(VMin) ||
 509       Precision::IsInfinite(VMax) ||
 510       Precision::IsInfinite(UMin) ||
 511       Precision::IsInfinite(UMax)   ) {
 512     switch (Type) {
 513       case GeomAbs_Plane:
 514       {
 515         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
 516         return;
 517       }
 518       default:
 519       {
 520               B.SetWhole();
 521               return;
 522       }
 523     }
 524   }
 525
 526   switch (Type) {
 527
 528     case GeomAbs_Plane:
 529     {
 530       gp_Pln Plan = S.Plane();
 531       B.Add(ElSLib::Value(UMin,VMin,Plan));
 532       B.Add(ElSLib::Value(UMin,VMax,Plan));
 533       B.Add(ElSLib::Value(UMax,VMin,Plan));
 534       B.Add(ElSLib::Value(UMax,VMax,Plan));
 535       B.Enlarge(Tol);
 536       break;
 537     }
 538     case GeomAbs_Cylinder:
 539     {
 540       BndLib::Add(S.Cylinder(), UMin, UMax, VMin, VMax, Tol, B);
 541       break;
 542     }
 543     case GeomAbs_Cone:
 544     {
 545       BndLib::Add(S.Cone(), UMin, UMax, VMin, VMax, Tol, B);
 546       break;
 547     }
 548     case GeomAbs_Sphere:
 549     {
 550       BndLib::Add(S.Sphere(), UMin, UMax, VMin, VMax, Tol, B);
 551       break;
 552     }
 553     default:
 554     {
 555       AddGenSurf(S, UMin, UMax, VMin, VMax, Tol, B);
 556     }
 557   }
 558 }
 559 //=======================================================================
 560 //function : AddGenSurf
 561 //purpose  :
 562 //=======================================================================
 563 void BndLib_AddSurface::AddGenSurf(const Adaptor3d_Surface& S,
 564                                    const Standard_Real UMin,
 565                                    const Standard_Real UMax,
 566                                    const Standard_Real VMin,
 567                                    const Standard_Real VMax,
 568                                    const Standard_Real Tol,
 569                                    Bnd_Box& B)
 570 {
 571   Standard_Integer Nu = NbUSamples(S, UMin, UMax);
 572   Standard_Integer Nv = NbVSamples(S, VMin, VMax);
 573   //
 574   Standard_Real CoordMin[3] = {RealLast(), RealLast(), RealLast()};
 575   Standard_Real CoordMax[3] = {-RealLast(), -RealLast(), -RealLast()};
 576   Standard_Real DeflMax[3] = {-RealLast(), -RealLast(), -RealLast()};
 577   //
 578   //
 579   Standard_Real du = (UMax-UMin)/(Nu-1), du2 = du / 2.;
 580   Standard_Real dv = (VMax-VMin)/(Nv-1), dv2 = dv / 2.;
 581   NCollection_Array2<gp_XYZ> aPnts(1, Nu, 1, Nv);
 582   Standard_Real u, v;
 583   Standard_Integer i, j, k;
 584   gp_Pnt P;
 585   for (i = 1, u = UMin; i <= Nu; i++, u += du){
 586     for (j = 1, v = VMin;j <= Nv; j++, v += dv){
 587       S.D0(u,v,P);
 588       aPnts(i, j) = P.XYZ();
 589       //
 590       for(k = 0; k < 3; ++k)
 591       {
 592         if(CoordMin[k] > P.Coord(k+1))
 593         {
 594           CoordMin[k] = P.Coord(k+1);
 595         }
 596         if(CoordMax[k] < P.Coord(k+1))
 597         {
 598           CoordMax[k] = P.Coord(k+1);
 599         }
 600       }
 601       //
 602       if(i > 1)
 603       {
 604         gp_XYZ aPm = 0.5 * (aPnts(i-1,j) + aPnts(i, j));
 605         S.D0(u - du2, v, P);
 606         gp_XYZ aD = (P.XYZ() - aPm);
 607         for(k = 0; k < 3; ++k)
 608         {
 609           if(CoordMin[k] > P.Coord(k+1))
 610           {
 611             CoordMin[k] = P.Coord(k+1);
 612           }
 613           if(CoordMax[k] < P.Coord(k+1))
 614           {
 615             CoordMax[k] = P.Coord(k+1);
 616           }
 617           Standard_Real d = Abs(aD.Coord(k+1));
 618           if(DeflMax[k] < d)
 619           {
 620             DeflMax[k] = d;
 621           }
 622         }
 623       }
 624       if(j > 1)
 625       {
 626         gp_XYZ aPm = 0.5 * (aPnts(i,j-1) + aPnts(i, j));
 627         S.D0(u , v - dv2, P);
 628         gp_XYZ aD = (P.XYZ() - aPm);
 629         for(k = 0; k < 3; ++k)
 630         {
 631           if(CoordMin[k] > P.Coord(k+1))
 632           {
 633             CoordMin[k] = P.Coord(k+1);
 634           }
 635           if(CoordMax[k] < P.Coord(k+1))
 636           {
 637             CoordMax[k] = P.Coord(k+1);
 638           }
 639           Standard_Real d = Abs(aD.Coord(k+1));
 640           if(DeflMax[k] < d)
 641           {
 642             DeflMax[k] = d;
 643           }
 644         }
 645       }
 646     }
 647   }
 648   //
 649   //Adjusting minmax
 650   Standard_Real eps = Max(Tol, Precision::Confusion());
 651   for(k = 0; k < 3; ++k)
 652   {
 653     Standard_Real d = DeflMax[k];
 654     if(d <= eps)
 655     {
 656       continue;
 657     }
 658
 659     Standard_Real CMin = CoordMin[k];
 660     Standard_Real CMax = CoordMax[k];
 661     for(i = 1; i <= Nu; ++i)
 662     {
 663       for(j = 1; j <= Nv; ++j)
 664       {
 665         if(aPnts(i,j).Coord(k+1) - CMin < d)
 666         {
 667           Standard_Real umin, umax, vmin, vmax;
 668           umin = UMin + Max(0, i-2) * du;
 669           umax = UMin + Min(Nu-1, i) * du;
 670           vmin = VMin + Max(0, j-2) * dv;
 671           vmax = VMin + Min(Nv-1, j) * dv;
 672           Standard_Real cmin = AdjustExtr(S, umin, umax, vmin, vmax,
 673                                           CMin, k + 1, eps, Standard_True);
 674           if(cmin < CMin)
 675           {
 676             CMin = cmin;
 677           }
 678         }
 679         else if(CMax - aPnts(i,j).Coord(k+1) < d)
 680         {
 681           Standard_Real umin, umax, vmin, vmax;
 682           umin = UMin + Max(0, i-2) * du;
 683           umax = UMin + Min(Nu-1, i) * du;
 684           vmin = VMin + Max(0, j-2) * dv;
 685           vmax = VMin + Min(Nv-1, j) * dv;
 686           Standard_Real cmax = AdjustExtr(S, umin, umax, vmin, vmax,
 687                                           CMax, k + 1, eps, Standard_False);
 688           if(cmax > CMax)
 689           {
 690             CMax = cmax;
 691           }
 692         }
 693       }
 694     }
 695     CoordMin[k] = CMin;
 696     CoordMax[k] = CMax;
 697
 698   }
 699
 700   B.Add(gp_Pnt(CoordMin[0], CoordMin[1], CoordMin[2]));
 701   B.Add(gp_Pnt(CoordMax[0], CoordMax[1], CoordMax[2]));
 702   B.Enlarge(eps);
 703 }
 704 //
 705
 706 //
 707 class SurfMaxMinCoord : public math_MultipleVarFunction
 708 {
 709 public:
 710   SurfMaxMinCoord(const Adaptor3d_Surface& theSurf,
 711               const Standard_Real UMin,
 712               const Standard_Real UMax,
 713               const Standard_Real VMin,
 714               const Standard_Real VMax,
 715               const Standard_Integer CoordIndx,
 716               const Standard_Real Sign)
 717 : mySurf(theSurf),
 718   myUMin(UMin),
 719   myUMax(UMax),
 720   myVMin(VMin),
 721   myVMax(VMax),
 722   myCoordIndx(CoordIndx),
 723   mySign(Sign),
 724   myPenalty(0.)
 725   {
 726     math_Vector X(1,2);
 727     X(1) = UMin;
 728     X(2) = (VMin + VMax) / 2.;
 729     Standard_Real F1, F2;
 730     Value(X, F1);
 731     X(1) = UMax;
 732     Value(X, F2);
 733     Standard_Real DU = Abs((F2 - F1) / (UMax - UMin));
 734     X(1) = (UMin + UMax) / 2.;
 735     X(2) = VMin;
 736     Value(X, F1);
 737     X(2) = VMax;
 738     Value(X, F2);
 739     Standard_Real DV = Abs((F2 - F1) / (VMax - VMin));
 740     myPenalty = 10. * Max(DU, DV);
 741     myPenalty = Max(myPenalty, 1.);
 742   }
 743
 744   Standard_Boolean Value (const math_Vector& X,
 745                                   Standard_Real& F)
 746   {
 747     if (CheckInputData(X))
 748     {
 749       gp_Pnt aP = mySurf.Value(X(1), X(2));
 750       F = mySign * aP.Coord(myCoordIndx);
 751     }
 752     else
 753     {
 754       Standard_Real UPen = 0., VPen = 0., u0, v0;
 755       if(X(1) < myUMin)
 756       {
 757         UPen = myPenalty * (myUMin - X(1));
 758         u0 = myUMin;
 759       }
 760       else if(X(1) > myUMax)
 761       {
 762         UPen = myPenalty * (X(1) - myUMax);
 763         u0 = myUMax;
 764       }
 765       else
 766       {
 767         u0 = X(1);
 768       }
 769       //
 770       if(X(2) < myVMin)
 771       {
 772         VPen = myPenalty * (myVMin - X(2));
 773         v0 = myVMin;
 774       }
 775       else if(X(2) > myVMax)
 776       {
 777         VPen = myPenalty * (X(2) - myVMax);
 778         v0 = myVMax;
 779       }
 780       else
 781       {
 782         v0 = X(2);
 783       }
 784       //
 785       gp_Pnt aP = mySurf.Value(u0, v0);
 786       F = mySign * aP.Coord(myCoordIndx) + UPen + VPen;
 787     }
 788
 789     return Standard_True;
 790   }
 791
 792
 793
 794   Standard_Integer NbVariables() const
 795   {
 796     return 2;
 797   }
 798
 799 private:
 800   SurfMaxMinCoord & operator = (const SurfMaxMinCoord & theOther);
 801
 802   Standard_Boolean CheckInputData(const math_Vector& theParams)
 803   {
 804     if (theParams(1) < myUMin ||
 805         theParams(1) > myUMax ||
 806         theParams(2) < myVMin ||
 807         theParams(2) > myVMax)
 808       return Standard_False;
 809     return Standard_True;
 810   }
 811
 812   const Adaptor3d_Surface& mySurf;
 813   Standard_Real myUMin;
 814   Standard_Real myUMax;
 815   Standard_Real myVMin;
 816   Standard_Real myVMax;
 817   Standard_Integer myCoordIndx;
 818   Standard_Real mySign;
 819   Standard_Real myPenalty;
 820 };
 821
 822 //=======================================================================
 823 //function : AdjustExtr
 824 //purpose  :
 825 //=======================================================================
 826
 827 Standard_Real AdjustExtr(const Adaptor3d_Surface& S,
 828                          const Standard_Real UMin,
 829                          const Standard_Real UMax,
 830                          const Standard_Real VMin,
 831                          const Standard_Real VMax,
 832                          const Standard_Real Extr0,
 833                          const Standard_Integer CoordIndx,
 834                          const Standard_Real Tol,
 835                          const Standard_Boolean IsMin)
 836 {
 837   Standard_Real aSign = IsMin ? 1.:-1.;
 838   Standard_Real extr = aSign * Extr0;
 839   Standard_Real relTol = 2.*Tol;
 840   if(Abs(extr) > Tol)
 841   {
 842     relTol /= Abs(extr);
 843   }
 844   Standard_Real Du = (S.LastUParameter() - S.FirstUParameter());
 845   Standard_Real Dv = (S.LastVParameter() - S.FirstVParameter());
 846   //
 847   math_Vector aT(1,2);
 848   math_Vector aLowBorder(1,2);
 849   math_Vector aUppBorder(1,2);
 850   math_Vector aSteps(1,2);
 851   aLowBorder(1) = UMin;
 852   aUppBorder(1) = UMax;
 853   aLowBorder(2) = VMin;
 854   aUppBorder(2) = VMax;
 855
 856   Standard_Integer aNbU = Max(8, RealToInt(32 * (UMax - UMin) / Du));
 857   Standard_Integer aNbV = Max(8, RealToInt(32 * (VMax - VMin) / Dv));
 858   Standard_Integer aNbParticles = aNbU * aNbV;
 859   Standard_Real aMaxUStep = (UMax - UMin) / (aNbU + 1);
 860   aSteps(1) = Min(0.1 * Du, aMaxUStep);
 861   Standard_Real aMaxVStep = (VMax - VMin) / (aNbV + 1);
 862   aSteps(2) = Min(0.1 * Dv, aMaxVStep);
 863
 864   SurfMaxMinCoord aFunc(S, UMin, UMax, VMin, VMax, CoordIndx, aSign);
 865   math_PSO aFinder(&aFunc, aLowBorder, aUppBorder, aSteps, aNbParticles);
 866   aFinder.Perform(aSteps, extr, aT);
 867
 868   //Refinement of extremal value
 869   math_Matrix aDir(1, 2, 1, 2, 0.0);
 870   aDir(1, 1) = 1.;
 871   aDir(2, 1) = 0.;
 872   aDir(1, 2) = 0.;
 873   aDir(2, 2) = 1.;
 874
 875   Standard_Integer aNbIter = 200;
 876   math_Powell powell(aFunc, relTol, aNbIter, Tol);
 877   powell.Perform(aFunc, aT, aDir);
 878
 879   if (powell.IsDone())
 880   {
 881     powell.Location(aT);
 882     extr = powell.Minimum();
 883   }
 884
 885   return aSign * extr;
 886 }
 887
 888 //=======================================================================
 889 //function : NbUSamples
 890 //purpose  :
 891 //=======================================================================
 892
 893 Standard_Integer NbUSamples(const Adaptor3d_Surface& S,
 894                             const Standard_Real Umin,
 895                             const Standard_Real Umax)
 896 {
 897   Standard_Integer N;
 898   GeomAbs_SurfaceType Type = S.GetType();
 899   switch (Type) {
 900   case GeomAbs_BezierSurface:
 901     {
 902       N = 2*S.NbUPoles();
 903       //By default parametric range of Bezier surf is [0, 1] [0, 1]
 904       Standard_Real du = Umax - Umin;
 905       if(du < .9)
 906       {
 907         N = RealToInt(du*N) + 1;
 908         N = Max(N, 5);
 909       }
 910       break;
 911     }
 912   case GeomAbs_BSplineSurface:
 913     {
 914       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
 915       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
 916       Standard_Real umin, umax, vmin, vmax;
 917       BS->Bounds(umin, umax, vmin, vmax);
 918       Standard_Real du = (Umax - Umin) / (umax - umin);
 919       if(du < .9)
 920       {
 921         N = RealToInt(du*N) + 1;
 922         N = Max(N, 5);
 923       }
 924       break;
 925     }
 926   default:
 927     N = 33;
 928   }
 929   return Min (50,N);
 930 }
 931
 932 //=======================================================================
 933 //function : NbVSamples
 934 //purpose  :
 935 //=======================================================================
 936
 937 Standard_Integer NbVSamples(const Adaptor3d_Surface& S,
 938                             const Standard_Real Vmin,
 939                             const Standard_Real Vmax)
 940 {
 941   Standard_Integer N;
 942   GeomAbs_SurfaceType Type = S.GetType();
 943   switch (Type) {
 944   case GeomAbs_BezierSurface:
 945     {
 946       N = 2*S.NbVPoles();
 947       //By default parametric range of Bezier surf is [0, 1] [0, 1]
 948       Standard_Real dv = Vmax - Vmin;
 949       if(dv < .9)
 950       {
 951         N = RealToInt(dv*N) + 1;
 952         N = Max(N, 5);
 953       }
 954       break;
 955     }
 956   case GeomAbs_BSplineSurface:
 957     {
 958       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
 959       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
 960       Standard_Real umin, umax, vmin, vmax;
 961       BS->Bounds(umin, umax, vmin, vmax);
 962       Standard_Real dv = (Vmax - Vmin) / (vmax - vmin);
 963       if(dv < .9)
 964       {
 965         N = RealToInt(dv*N) + 1;
 966         N = Max(N, 5);
 967       }
 968       break;
 969     }
 970   default:
 971     N = 33;
 972   }
 973   return Min(50,N);
 974 }
 975