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
  20 #include <BndLib_AddSurface.ixx>
  21
  22
  23 #include <Adaptor3d_HSurface.hxx>
  24 #include <GeomAbs_SurfaceType.hxx>
  25 #include <BndLib.hxx>
  26 #include <BSplCLib.hxx>
  27 #include <gp_Pnt.hxx>
  28 #include <gp_Pln.hxx>
  29 #include <ElSLib.hxx>
  30 #include <TColgp_Array2OfPnt.hxx>
  31 #include <TColStd_Array1OfReal.hxx>
  32 #include <TColStd_Array1OfInteger.hxx>
  33 #include <Geom_BSplineSurface.hxx>
  34 #include <Geom_BezierSurface.hxx>
  35
  36 #include <Precision.hxx>
  37
  38 //=======================================================================
  39 //function : Add
  40 //purpose  :
  41 //=======================================================================
  42
  43 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
  44                             const Standard_Real Tol,
  45                             Bnd_Box& B )
  46 {
  47
  48   BndLib_AddSurface::Add(S,
  49                          S.FirstUParameter(),
  50                          S.LastUParameter (),
  51                          S.FirstVParameter(),
  52                          S.LastVParameter (),Tol,B);
  53 }
  54 //=======================================================================
  55 //function : NbUSamples
  56 //purpose  :
  57 //=======================================================================
  58
  59 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S)
  60 {
  61   Standard_Integer N;
  62   GeomAbs_SurfaceType Type = S.GetType();
  63   switch (Type) {
  64   case GeomAbs_BezierSurface:
  65     {
  66       N = 2*S.NbUPoles();
  67       break;
  68     }
  69   case GeomAbs_BSplineSurface:
  70     {
  71       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
  72       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
  73       break;
  74     }
  75   default:
  76     N = 33;
  77   }
  78   return Min (50,N);
  79 }
  80
  81 //=======================================================================
  82 //function : NbVSamples
  83 //purpose  :
  84 //=======================================================================
  85
  86 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S)
  87 {
  88   Standard_Integer N;
  89   GeomAbs_SurfaceType Type = S.GetType();
  90   switch (Type) {
  91   case GeomAbs_BezierSurface:
  92     {
  93       N = 2*S.NbVPoles();
  94       break;
  95     }
  96   case GeomAbs_BSplineSurface:
  97     {
  98       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
  99       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
 100       break;
 101     }
 102   default:
 103     N = 33;
 104   }
 105   return Min(50,N);
 106 }
 107
 108 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
 109 static gp_Pnt BaryCenter(const gp_Pln        &aPlane,
 110                          const Standard_Real  aUMin,
 111                          const Standard_Real  aUMax,
 112                          const Standard_Real  aVMin,
 113                          const Standard_Real  aVMax)
 114 {
 115   Standard_Real aU, aV;
 116   Standard_Boolean isU1Inf = Precision::IsInfinite(aUMin);
 117   Standard_Boolean isU2Inf = Precision::IsInfinite(aUMax);
 118   Standard_Boolean isV1Inf = Precision::IsInfinite(aVMin);
 119   Standard_Boolean isV2Inf = Precision::IsInfinite(aVMax);
 120
 121   if (isU1Inf && isU2Inf)
 122     aU = 0;
 123   else if (isU1Inf)
 124     aU = aUMax - 10.;
 125   else if (isU2Inf)
 126     aU = aUMin + 10.;
 127   else
 128     aU = (aUMin + aUMax)/2.;
 129
 130   if (isV1Inf && isV2Inf)
 131     aV = 0;
 132   else if (isV1Inf)
 133     aV = aVMax - 10.;
 134   else if (isV2Inf)
 135     aV = aVMin + 10.;
 136   else
 137     aV = (aVMin + aVMax)/2.;
 138
 139   gp_Pnt aCenter = ElSLib::Value(aU, aV, aPlane);
 140
 141   return aCenter;
 142 }
 143
 144 static void TreatInfinitePlane(const gp_Pln        &aPlane,
 145                                const Standard_Real  aUMin,
 146                                const Standard_Real  aUMax,
 147                                const Standard_Real  aVMin,
 148                                const Standard_Real  aVMax,
 149                                const Standard_Real  aTol,
 150                                      Bnd_Box       &aB)
 151 {
 152   // Get 3 coordinate axes of the plane.
 153   const gp_Dir        &aNorm        = aPlane.Axis().Direction();
 154   const Standard_Real  anAngularTol = RealEpsilon();
 155
 156   // Get location of the plane as its barycenter
 157   gp_Pnt aLocation = BaryCenter(aPlane, aUMin, aUMax, aVMin, aVMax);
 158
 159   if (aNorm.IsParallel(gp::DX(), anAngularTol)) {
 160     aB.Add(aLocation);
 161     aB.OpenYmin();
 162     aB.OpenYmax();
 163     aB.OpenZmin();
 164     aB.OpenZmax();
 165   } else if (aNorm.IsParallel(gp::DY(), anAngularTol)) {
 166     aB.Add(aLocation);
 167     aB.OpenXmin();
 168     aB.OpenXmax();
 169     aB.OpenZmin();
 170     aB.OpenZmax();
 171   } else if (aNorm.IsParallel(gp::DZ(), anAngularTol)) {
 172     aB.Add(aLocation);
 173     aB.OpenXmin();
 174     aB.OpenXmax();
 175     aB.OpenYmin();
 176     aB.OpenYmax();
 177   } else {
 178     aB.SetWhole();
 179     return;
 180   }
 181
 182   aB.Enlarge(aTol);
 183 }
 184
 185 // Compute start and finish indexes used in convex hull.
 186 // theMinIdx - minimum poles index, that can be used.
 187 // theMaxIdx - maximum poles index, that can be used.
 188 // theShiftCoeff - shift between flatknots array and poles array.
 189 // This vaule should be equal to 1 in case of non periodic BSpline,
 190 // and (degree + 1) - mults(the lowest index).
 191 void ComputePolesIndexes(const TColStd_Array1OfReal &theFlatKnots,
 192                          const Standard_Integer theDegree,
 193                          const Standard_Real theMin,
 194                          const Standard_Real theMax,
 195                          const Standard_Integer theMinIdx,
 196                          const Standard_Integer theMaxIdx,
 197                          const Standard_Integer theShiftCoeff,
 198                          Standard_Integer &theOutMinIdx,
 199                          Standard_Integer &theOutMaxIdx)
 200 {
 201   // Compute first and last used flat knots.
 202   for(Standard_Integer aKnotIdx = theFlatKnots.Lower();
 203       aKnotIdx < theFlatKnots.Upper();
 204       aKnotIdx++)
 205   {
 206     if (theFlatKnots(aKnotIdx) <= theMin)
 207       theOutMinIdx = aKnotIdx;
 208
 209     if (theFlatKnots(theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower()) >= theMax)
 210       theOutMaxIdx = theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower();
 211   }
 212
 213   theOutMinIdx = Max(theOutMinIdx - 2 * theDegree + 2 - theShiftCoeff, theMinIdx);
 214   theOutMaxIdx = Min(theOutMaxIdx - 2 + theDegree + 1 - theShiftCoeff, theMaxIdx);
 215 }
 216
 217 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
 218 //=======================================================================
 219 //function : Add
 220 //purpose  :
 221 //=======================================================================
 222 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
 223                             const Standard_Real UMin,
 224                             const Standard_Real UMax,
 225                             const Standard_Real VMin,
 226                             const Standard_Real VMax,
 227                             const Standard_Real Tol,
 228                             Bnd_Box& B )
 229 {
 230   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
 231
 232   if (Precision::IsInfinite(VMin) ||
 233       Precision::IsInfinite(VMax) ||
 234       Precision::IsInfinite(UMin) ||
 235       Precision::IsInfinite(UMax)   ) {
 236 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
 237 //     B.SetWhole();
 238 //     return;
 239     switch (Type) {
 240     case GeomAbs_Plane:
 241       {
 242         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
 243         return;
 244       }
 245     default:
 246       {
 247         B.SetWhole();
 248         return;
 249       }
 250     }
 251 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
 252   }
 253
 254 //   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
 255
 256   switch (Type) {
 257
 258   case GeomAbs_Plane:
 259     {
 260       gp_Pln Plan = S.Plane();
 261       B.Add(ElSLib::Value(UMin,VMin,Plan));
 262       B.Add(ElSLib::Value(UMin,VMax,Plan));
 263       B.Add(ElSLib::Value(UMax,VMin,Plan));
 264       B.Add(ElSLib::Value(UMax,VMax,Plan));
 265       B.Enlarge(Tol);
 266       break;
 267     }
 268   case GeomAbs_Cylinder:
 269     {
 270       BndLib::Add(S.Cylinder(),UMin,UMax,VMin,VMax,Tol,B);
 271       break;
 272     }
 273   case GeomAbs_Cone:
 274     {
 275       BndLib::Add(S.Cone(),UMin,UMax,VMin,VMax,Tol,B);
 276       break;
 277     }
 278   case GeomAbs_Torus:
 279     {
 280       BndLib::Add(S.Torus(),UMin,UMax,VMin,VMax,Tol,B);
 281       break;
 282     }
 283   case GeomAbs_Sphere:
 284     {
 285       if (Abs(UMin) < Precision::Angular() &&
 286           Abs(UMax - 2.*M_PI) < Precision::Angular() &&
 287           Abs(VMin + M_PI/2.) < Precision::Angular() &&
 288           Abs(VMax - M_PI/2.) < Precision::Angular()) // a whole sphere
 289         BndLib::Add(S.Sphere(),Tol,B);
 290       else
 291         BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
 292       break;
 293     }
 294   case GeomAbs_OffsetSurface:
 295     {
 296       Handle(Adaptor3d_HSurface) HS = S.BasisSurface();
 297       Add (HS->Surface(),UMin,UMax,VMin,VMax,Tol,B);
 298       B.Enlarge(S.OffsetValue());
 299       B.Enlarge(Tol);
 300       break;
 301     }
 302   case GeomAbs_BezierSurface:
 303   case GeomAbs_BSplineSurface:
 304     {
 305       Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
 306       Standard_Real PTol = Precision::Parametric(Precision::Confusion());
 307       // Borders of underlying geometry.
 308       Standard_Real anUMinParam = UMin, anUMaxParam = UMax,// BSpline case.
 309                      aVMinParam = VMin,  aVMaxParam = VMax;
 310       if (Type == GeomAbs_BezierSurface)
 311       {
 312         // Bezier surface:
 313         // All of poles used for any parameter,
 314         // thats why in case of trimmed parameters handled by grid algorithm.
 315
 316         if (Abs(UMin-S.FirstUParameter()) > PTol ||
 317             Abs(VMin-S.FirstVParameter()) > PTol ||
 318             Abs(UMax-S.LastUParameter ()) > PTol ||
 319             Abs(VMax-S.LastVParameter ()) > PTol )
 320         {
 321           // Borders not equal to topology borders.
 322           isUseConvexHullAlgorithm = Standard_False;
 323         }
 324       }
 325       else
 326       {
 327         // BSpline:
 328         // If Umin, Vmin, Umax, Vmax lies inside geometry bounds then:
 329         // use convex hull algorithm,
 330         // if Umin, VMin, Umax, Vmax lies outside then:
 331         // use grid algorithm on analytic continuation (default case).
 332         S.BSpline()->Bounds(anUMinParam, anUMaxParam, aVMinParam, aVMaxParam);
 333
 334         if ( (UMin - anUMinParam) < -PTol ||
 335              (VMin -  aVMinParam) < -PTol ||
 336              (UMax - anUMaxParam) >  PTol ||
 337              (VMax -  aVMaxParam) >  PTol )
 338         {
 339           // Out of geometry borders.
 340           isUseConvexHullAlgorithm = Standard_False;
 341         }
 342       }
 343
 344       if (isUseConvexHullAlgorithm)
 345       {
 346           TColgp_Array2OfPnt Tp(1,S.NbUPoles(),1,S.NbVPoles());
 347           Standard_Integer UMinIdx = 0, UMaxIdx = 0;
 348           Standard_Integer VMinIdx = 0, VMaxIdx = 0;
 349           if (Type == GeomAbs_BezierSurface)
 350           {
 351             S.Bezier()->Poles(Tp);
 352
 353             UMinIdx = Tp.LowerRow();
 354             UMaxIdx = Tp.UpperRow();
 355             VMinIdx = Tp.LowerCol();
 356             VMaxIdx = Tp.UpperCol();
 357           }
 358           else
 359           {
 360             S.BSpline()->Poles(Tp);
 361
 362             UMinIdx = Tp.LowerRow();
 363             UMaxIdx = Tp.UpperRow();
 364             VMinIdx = Tp.LowerCol();
 365             VMaxIdx = Tp.UpperCol();
 366
 367             if (UMin > anUMinParam ||
 368                 UMax < anUMaxParam)
 369             {
 370               Standard_Integer anUFlatKnotsCount = S.BSpline()->NbUPoles() + S.BSpline()->UDegree() + 1;
 371               Standard_Integer aShift = 1;
 372
 373               if (S.BSpline()->IsUPeriodic())
 374               {
 375                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbUKnots());
 376                 S.BSpline()->UMultiplicities(aMults);
 377                 anUFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->UDegree(), Standard_True);
 378
 379                 aShift = S.BSpline()->UDegree() + 1 - S.BSpline()->UMultiplicity(1);
 380               }
 381
 382               TColStd_Array1OfReal anUFlatKnots(1, anUFlatKnotsCount);
 383               S.BSpline()->UKnotSequence(anUFlatKnots);
 384
 385               ComputePolesIndexes(anUFlatKnots,
 386                                   S.BSpline()->UDegree(),
 387                                   UMin, UMax,
 388                                   UMinIdx, UMaxIdx,  // Min and Max Indexes
 389                                   aShift,
 390                                   UMinIdx, UMaxIdx); // the Output indexes
 391             }
 392
 393             if (VMin > aVMinParam ||
 394                 VMax < aVMaxParam)
 395             {
 396               Standard_Integer anVFlatKnotsCount = S.BSpline()->NbVPoles() + S.BSpline()->VDegree() + 1;
 397               Standard_Integer aShift = 1;
 398
 399               if (S.BSpline()->IsVPeriodic())
 400               {
 401                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbVKnots());
 402                 S.BSpline()->VMultiplicities(aMults);
 403                 anVFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->VDegree(), Standard_True);
 404
 405                 aShift = S.BSpline()->VDegree() + 1 - S.BSpline()->VMultiplicity(1);
 406               }
 407
 408               TColStd_Array1OfReal anVFlatKnots(1, anVFlatKnotsCount);
 409               S.BSpline()->VKnotSequence(anVFlatKnots);
 410
 411               ComputePolesIndexes(anVFlatKnots,
 412                                   S.BSpline()->VDegree(),
 413                                   VMin, VMax,
 414                                   VMinIdx, VMaxIdx,  // Min and Max Indexes
 415                                   aShift,
 416                                   VMinIdx, VMaxIdx); // the Output indexes
 417             }
 418
 419           }
 420
 421           // Use poles to build convex hull.
 422           for (Standard_Integer i = UMinIdx; i <= UMaxIdx; i++)
 423           {
 424             for (Standard_Integer j = VMinIdx; j <= VMaxIdx; j++)
 425             {
 426               B.Add(Tp(i,j));
 427             }
 428           }
 429
 430           B.Enlarge(Tol);
 431           break;
 432       }
 433     }
 434   default:
 435     {
 436       Standard_Integer Nu = NbUSamples(S);
 437       Standard_Integer Nv = NbVSamples(S);
 438       gp_Pnt P;
 439       for (Standard_Integer i =1 ;i<=Nu;i++){
 440         Standard_Real U = UMin + ((UMax-UMin)*(i-1)/(Nu-1));
 441         for (Standard_Integer j=1 ;j<=Nv;j++){
 442           Standard_Real V = VMin + ((VMax-VMin)*(j-1)/(Nv-1));
 443           S.D0(U,V,P);
 444           B.Add(P);
 445         }
 446       }
 447       B.Enlarge(Tol);
 448     }
 449   }
 450 }