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