src/BRepBndLib/BRepBndLib_1.cxx

   1 // Copyright (c) 1999-2017 OPEN CASCADE SAS
   2 //
   3 // This file is part of Open CASCADE Technology software library.
   4 //
   5 // This library is free software; you can redistribute it and/or modify it under
   6 // the terms of the GNU Lesser General Public License version 2.1 as published
   7 // by the Free Software Foundation, with special exception defined in the file
   8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
   9 // distribution for complete text of the license and disclaimer of any warranty.
  10 //
  11 // Alternatively, this file may be used under the terms of Open CASCADE
  12 // commercial license or contractual agreement.
  13
  14 #include <Adaptor3d_HCurve.hxx>
  15 #include <Adaptor3d_HSurface.hxx>
  16 #include <GeomAdaptor_HCurve.hxx>
  17 #include <BRepBndLib.hxx>
  18 #include <GProp_GProps.hxx>
  19 #include <TopoDS_Shape.hxx>
  20 #include <BRep_Tool.hxx>
  21 #include <TopoDS.hxx>
  22 #include <Bnd_OBB.hxx>
  23 #include <BRepGProp.hxx>
  24 #include <TopExp_Explorer.hxx>
  25 #include <GProp_PrincipalProps.hxx>
  26 #include <gp_Ax3.hxx>
  27 #include <BRepBuilderAPI_Transform.hxx>
  28 #include <Bnd_Box.hxx>
  29 #include <NCollection_List.hxx>
  30 #include <TColgp_Array1OfPnt.hxx>
  31 #include <TColStd_Array1OfReal.hxx>
  32 #include <Geom_Plane.hxx>
  33 #include <Geom_Line.hxx>
  34 #include <TColStd_Array1OfInteger.hxx>
  35 #include <BRepAdaptor_Curve.hxx>
  36 #include <BRepAdaptor_HSurface.hxx>
  37
  38 #include <Geom_OffsetCurve.hxx>
  39 #include <Geom_BSplineCurve.hxx>
  40 #include <Geom_BezierCurve.hxx>
  41 #include <Geom_BSplineSurface.hxx>
  42 #include <Geom_BezierSurface.hxx>
  43
  44 //=======================================================================
  45 // Function : IsLinear
  46 // purpose : Returns TRUE if theC is line-like.
  47 //=======================================================================
  48 static Standard_Boolean IsLinear(const Adaptor3d_Curve& theC)
  49 {
  50   const GeomAbs_CurveType aCT = theC.GetType();
  51   if(aCT == GeomAbs_OffsetCurve)
  52   {
  53     return IsLinear(GeomAdaptor_Curve(theC.OffsetCurve()->BasisCurve()));
  54   }
  55
  56   if((aCT == GeomAbs_BSplineCurve) || (aCT == GeomAbs_BezierCurve))
  57   {
  58     // Indeed, curves with C0-continuity and degree==1, may be
  59     // represented with set of points. It will be possible made
  60     // in the future.
  61
  62     return ((theC.Degree() == 1) &&
  63             (theC.Continuity() != GeomAbs_C0));
  64   }
  65
  66   if(aCT == GeomAbs_Line)
  67   {
  68     return Standard_True;
  69   }
  70
  71   return Standard_False;
  72 }
  73
  74 //=======================================================================
  75 // Function : IsPlanar
  76 // purpose : Returns TRUE if theS is plane-like.
  77 //=======================================================================
  78 static Standard_Boolean IsPlanar(const Adaptor3d_Surface& theS)
  79 {
  80   const GeomAbs_SurfaceType aST = theS.GetType();
  81   if(aST == GeomAbs_OffsetSurface)
  82   {
  83     return IsPlanar(theS.BasisSurface()->Surface());
  84   }
  85
  86   if(aST == GeomAbs_SurfaceOfExtrusion)
  87   {
  88     return IsLinear(theS.BasisCurve()->Curve());
  89   }
  90
  91   if((aST == GeomAbs_BSplineSurface) || (aST == GeomAbs_BezierSurface))
  92   {
  93     if((theS.UDegree() != 1) || (theS.VDegree() != 1))
  94       return Standard_False;
  95
  96     // Indeed, surfaces with C0-continuity and degree==1, may be
  97     // represented with set of points. It will be possible made
  98     // in the future.
  99
 100     return ((theS.UContinuity() != GeomAbs_C0) && (theS.VContinuity() != GeomAbs_C0));
 101   }
 102
 103   if(aST == GeomAbs_Plane)
 104   {
 105     return Standard_True;
 106   }
 107
 108   return Standard_False;
 109 }
 110
 111 //=======================================================================
 112 // Function : PointsForOBB
 113 // purpose : Returns number of points for array.
 114 //
 115 // Attention!!!
 116 //  1. Start index for thePts must be 0 strictly.
 117 //  2. Currently, infinite edges/faces (e.g. half-space) are not
 118 //      processed correctly because computation of UV-bounds is a costly operation.
 119 //=======================================================================
 120 static Standard_Integer PointsForOBB(const TopoDS_Shape& theS,
 121                                      const Standard_Boolean theIsTriangulationUsed,
 122                                      TColgp_Array1OfPnt* thePts = 0,
 123                                      TColStd_Array1OfReal* theArrOfToler = 0)
 124 {
 125   Standard_Integer aRetVal = 0;
 126   TopExp_Explorer anExpF, anExpE;
 127
 128   // get all vertices from the shape
 129   for(anExpF.Init(theS, TopAbs_VERTEX); anExpF.More(); anExpF.Next())
 130   {
 131     const TopoDS_Vertex &aVert = TopoDS::Vertex(anExpF.Current());
 132     if(thePts)
 133     {
 134       const gp_Pnt aP = BRep_Tool::Pnt(aVert);
 135       (*thePts)(aRetVal) = aP;
 136     }
 137
 138     if(theArrOfToler)
 139     {
 140       (*theArrOfToler) (aRetVal) = BRep_Tool::Tolerance(aVert);
 141     }
 142
 143     ++aRetVal;
 144   }
 145
 146   if(aRetVal == 0)
 147     return 0;
 148
 149   // analyze the faces of the shape on planarity and existence of triangulation
 150   TopLoc_Location aLoc;
 151   for(anExpF.Init(theS, TopAbs_FACE); anExpF.More(); anExpF.Next())
 152   {
 153     const TopoDS_Face &aF = TopoDS::Face(anExpF.Current());
 154     const BRepAdaptor_Surface anAS(aF, Standard_False);
 155
 156     if (!IsPlanar(anAS.Surface()))
 157     {
 158       if (!theIsTriangulationUsed)
 159         // not planar and triangulation usage disabled
 160         return 0;
 161     }
 162     else
 163     {
 164       // planar face
 165       for(anExpE.Init(aF, TopAbs_EDGE); anExpE.More(); anExpE.Next())
 166       {
 167         const TopoDS_Edge &anE = TopoDS::Edge(anExpE.Current());
 168         if (BRep_Tool::IsGeometric (anE))
 169         {
 170           const BRepAdaptor_Curve anAC(anE);
 171           if (!IsLinear(anAC))
 172           {
 173             if (!theIsTriangulationUsed)
 174               // not linear and triangulation usage disabled
 175               return 0;
 176
 177             break;
 178           }
 179         }
 180       }
 181
 182       if (!anExpE.More())
 183         // skip planar face with linear edges as its vertices have already been added
 184         continue;
 185     }
 186
 187     // Use triangulation of the face
 188     const Handle(Poly_Triangulation) &aTrng = BRep_Tool::Triangulation(aF, aLoc);
 189     if (aTrng.IsNull())
 190       // no triangulation on the face
 191       return 0;
 192
 193     const Standard_Integer aCNode = aTrng->NbNodes();
 194     const TColgp_Array1OfPnt& aNodesArr = aTrng->Nodes();
 195     for (Standard_Integer i = 1; i <= aCNode; i++)
 196     {
 197       if (thePts)
 198       {
 199         const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
 200           aNodesArr(i).Transformed(aLoc);
 201         (*thePts)(aRetVal) = aP;
 202       }
 203
 204       if (theArrOfToler)
 205       {
 206         (*theArrOfToler) (aRetVal) = aTrng->Deflection();
 207       }
 208
 209       ++aRetVal;
 210     }
 211   }
 212
 213   // Consider edges without faces
 214
 215   for(anExpE.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExpE.More(); anExpE.Next())
 216   {
 217     const TopoDS_Edge &anE = TopoDS::Edge(anExpE.Current());
 218     if (BRep_Tool::IsGeometric (anE))
 219     {
 220       const BRepAdaptor_Curve anAC(anE);
 221       if (IsLinear(anAC))
 222       {
 223         // skip linear edge as its vertices have already been added
 224         continue;
 225       }
 226     }
 227
 228     if (!theIsTriangulationUsed)
 229       // not linear and triangulation usage disabled
 230       return 0;
 231
 232     const Handle(Poly_Polygon3D) &aPolygon = BRep_Tool::Polygon3D(anE, aLoc);
 233     if (aPolygon.IsNull())
 234       return 0;
 235
 236     const Standard_Integer aCNode = aPolygon->NbNodes();
 237     const TColgp_Array1OfPnt& aNodesArr = aPolygon->Nodes();
 238     for (Standard_Integer i = 1; i <= aCNode; i++)
 239     {
 240       if (thePts)
 241       {
 242         const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
 243           aNodesArr(i).Transformed(aLoc);
 244         (*thePts)(aRetVal) = aP;
 245       }
 246
 247       if (theArrOfToler)
 248       {
 249         (*theArrOfToler) (aRetVal) = aPolygon->Deflection();
 250       }
 251
 252       ++aRetVal;
 253     }
 254   }
 255
 256   return aRetVal;
 257 }
 258
 259 //=======================================================================
 260 // Function : IsWCS
 261 // purpose : Returns 0 if the theDir does not match any axis of WCS.
 262 //            Otherwise, returns the index of correspond axis.
 263 //=======================================================================
 264 static Standard_Integer IsWCS(const gp_Dir& theDir)
 265 {
 266   const Standard_Real aToler = Precision::Angular()*Precision::Angular();
 267
 268   const Standard_Real aX = theDir.X(),
 269                       aY = theDir.Y(),
 270                       aZ = theDir.Z();
 271
 272   const Standard_Real aVx = aY*aY + aZ*aZ,
 273                       aVy = aX*aX + aZ*aZ,
 274                       aVz = aX*aX + aY*aY;
 275
 276   if(aVz < aToler)
 277     return 3; // Z-axis
 278
 279   if(aVy < aToler)
 280     return 2; // Y-axis
 281
 282   if(aVx < aToler)
 283     return 1; // X-axis
 284
 285   return 0;
 286 }
 287
 288 //=======================================================================
 289 // Function : CheckPoints
 290 // purpose : Collects points for DiTO algorithm for OBB construction on
 291 //            linear/planar shapes and shapes having triangulation
 292 //            (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf).
 293 //=======================================================================
 294 static Standard_Boolean CheckPoints(const TopoDS_Shape& theS,
 295                                     const Standard_Boolean theIsTriangulationUsed,
 296                                     const Standard_Boolean theIsOptimal,
 297                                     const Standard_Boolean theIsShapeToleranceUsed,
 298                                     Bnd_OBB& theOBB)
 299 {
 300   const Standard_Integer aNbPnts = PointsForOBB(theS, theIsTriangulationUsed);
 301
 302   if(aNbPnts < 1)
 303     return Standard_False;
 304
 305   TColgp_Array1OfPnt anArrPnts(0, theOBB.IsVoid() ? aNbPnts - 1 : aNbPnts + 7);
 306   TColStd_Array1OfReal anArrOfTolerances;
 307   if(theIsShapeToleranceUsed)
 308   {
 309     anArrOfTolerances.Resize(anArrPnts.Lower(), anArrPnts.Upper(), Standard_False);
 310     anArrOfTolerances.Init(0.0);
 311   }
 312
 313   TColStd_Array1OfReal *aPtrArrTol = theIsShapeToleranceUsed ? &anArrOfTolerances : 0;
 314
 315   PointsForOBB(theS, theIsTriangulationUsed, &anArrPnts, aPtrArrTol);
 316
 317   if(!theOBB.IsVoid())
 318   {
 319     // All points of old OBB have zero-tolerance
 320     theOBB.GetVertex(&anArrPnts(aNbPnts));
 321   }
 322
 323 #if 0
 324   for(Standard_Integer i = anArrPnts.Lower(); i <= anArrPnts.Upper(); i++)
 325   {
 326     const gp_Pnt &aP = anArrPnts(i);
 327     std::cout << "point p" << i << " " << aP.X() << ", " <<
 328                                           aP.Y() << ", " <<
 329                                           aP.Z() << ", "<< std::endl;
 330   }
 331 #endif
 332
 333   theOBB.ReBuild(anArrPnts, aPtrArrTol, theIsOptimal);
 334
 335   return (!theOBB.IsVoid());
 336 }
 337
 338 //=======================================================================
 339 // Function : ComputeProperties
 340 // purpose : Computes properties of theS.
 341 //=======================================================================
 342 static void ComputeProperties(const TopoDS_Shape& theS,
 343                               GProp_GProps& theGCommon)
 344 {
 345   TopExp_Explorer anExp;
 346   for(anExp.Init(theS, TopAbs_SOLID); anExp.More(); anExp.Next())
 347   {
 348     GProp_GProps aG;
 349     BRepGProp::VolumeProperties(anExp.Current(), aG, Standard_True);
 350     theGCommon.Add(aG);
 351   }
 352
 353   for(anExp.Init(theS, TopAbs_FACE, TopAbs_SOLID); anExp.More(); anExp.Next())
 354   {
 355     GProp_GProps aG;
 356     BRepGProp::SurfaceProperties(anExp.Current(), aG, Standard_True);
 357     theGCommon.Add(aG);
 358   }
 359
 360   for(anExp.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExp.More(); anExp.Next())
 361   {
 362     GProp_GProps aG;
 363     BRepGProp::LinearProperties(anExp.Current(), aG, Standard_True);
 364     theGCommon.Add(aG);
 365   }
 366
 367   for(anExp.Init(theS, TopAbs_VERTEX, TopAbs_EDGE); anExp.More(); anExp.Next())
 368   {
 369     GProp_GProps aG(BRep_Tool::Pnt(TopoDS::Vertex(anExp.Current())));
 370     theGCommon.Add(aG);
 371   }
 372 }
 373
 374 //=======================================================================
 375 // Function : ComputePCA
 376 // purpose : Creates OBB with axes of inertia.
 377 //=======================================================================
 378 static void ComputePCA(const TopoDS_Shape& theS,
 379                        Bnd_OBB& theOBB,
 380                        const Standard_Boolean theIsTriangulationUsed,
 381                        const Standard_Boolean theIsOptimal,
 382                        const Standard_Boolean theIsShapeToleranceUsed)
 383 {
 384   // Compute the transformation matrix to obtain more tight bounding box
 385   GProp_GProps aGCommon;
 386   ComputeProperties(theS, aGCommon);
 387
 388   // Transform the shape to the local coordinate system
 389   gp_Trsf aTrsf;
 390
 391   const Standard_Integer anIdx1 =
 392                   IsWCS(aGCommon.PrincipalProperties().FirstAxisOfInertia());
 393   const Standard_Integer anIdx2 =
 394                   IsWCS(aGCommon.PrincipalProperties().SecondAxisOfInertia());
 395
 396   if((anIdx1 == 0) || (anIdx2 == 0))
 397   {
 398     // Coordinate system in which the shape will have the optimal bounding box
 399     gp_Ax3 aLocCoordSys(aGCommon.CentreOfMass(),
 400                         aGCommon.PrincipalProperties().ThirdAxisOfInertia(),
 401                         aGCommon.PrincipalProperties().FirstAxisOfInertia());
 402     aTrsf.SetTransformation(aLocCoordSys);
 403   }
 404
 405   const TopoDS_Shape aST = (aTrsf.Form() == gp_Identity) ? theS :
 406                                               theS.Moved(TopLoc_Location(aTrsf));
 407
 408   // Initial axis-aligned BndBox
 409   Bnd_Box aShapeBox;
 410   if(theIsOptimal)
 411   {
 412     BRepBndLib::AddOptimal(aST, aShapeBox, theIsTriangulationUsed, theIsShapeToleranceUsed);
 413   }
 414   else
 415   {
 416     BRepBndLib::Add(aST, aShapeBox);
 417   }
 418   if (aShapeBox.IsVoid())
 419   {
 420     return;
 421   }
 422
 423   gp_Pnt aPMin = aShapeBox.CornerMin();
 424   gp_Pnt aPMax = aShapeBox.CornerMax();
 425
 426   gp_XYZ aXDir(1, 0, 0);
 427   gp_XYZ aYDir(0, 1, 0);
 428   gp_XYZ aZDir(0, 0, 1);
 429
 430   // Compute the center of the box
 431   gp_XYZ aCenter = (aPMin.XYZ() + aPMax.XYZ()) / 2.;
 432
 433   // Compute the half diagonal size of the box.
 434   // It takes into account the gap.
 435   gp_XYZ anOBBHSize = (aPMax.XYZ() - aPMin.XYZ()) / 2.;
 436
 437   // Apply transformation if necessary
 438   if(aTrsf.Form() != gp_Identity)
 439   {
 440     aTrsf.Invert();
 441     aTrsf.Transforms(aCenter);
 442
 443     // Make transformation
 444     const Standard_Real * aMat = &aTrsf.HVectorialPart().Value(1, 1);
 445     // Compute axes directions of the box
 446     aXDir = gp_XYZ(aMat[0], aMat[3], aMat[6]);
 447     aYDir = gp_XYZ(aMat[1], aMat[4], aMat[7]);
 448     aZDir = gp_XYZ(aMat[2], aMat[5], aMat[8]);
 449   }
 450
 451   if(theOBB.IsVoid())
 452   {
 453     // Create the OBB box
 454
 455     // Set parameters to the OBB
 456     theOBB.SetCenter(aCenter);
 457
 458     theOBB.SetXComponent(aXDir, anOBBHSize.X());
 459     theOBB.SetYComponent(aYDir, anOBBHSize.Y());
 460     theOBB.SetZComponent(aZDir, anOBBHSize.Z());
 461     theOBB.SetAABox(aTrsf.Form() == gp_Identity);
 462   }
 463   else
 464   {
 465     // Recreate the OBB box
 466
 467     TColgp_Array1OfPnt aListOfPnts(0, 15);
 468     theOBB.GetVertex(&aListOfPnts(0));
 469
 470     const Standard_Real aX = anOBBHSize.X();
 471     const Standard_Real aY = anOBBHSize.Y();
 472     const Standard_Real aZ = anOBBHSize.Z();
 473
 474     const gp_XYZ aXext = aX*aXDir,
 475                  aYext = aY*aYDir,
 476                  aZext = aZ*aZDir;
 477
 478     Standard_Integer aPntIdx = 8;
 479     aListOfPnts(aPntIdx++) = aCenter - aXext - aYext - aZext;
 480     aListOfPnts(aPntIdx++) = aCenter + aXext - aYext - aZext;
 481     aListOfPnts(aPntIdx++) = aCenter - aXext + aYext - aZext;
 482     aListOfPnts(aPntIdx++) = aCenter + aXext + aYext - aZext;
 483     aListOfPnts(aPntIdx++) = aCenter - aXext - aYext + aZext;
 484     aListOfPnts(aPntIdx++) = aCenter + aXext - aYext + aZext;
 485     aListOfPnts(aPntIdx++) = aCenter - aXext + aYext + aZext;
 486     aListOfPnts(aPntIdx++) = aCenter + aXext + aYext + aZext;
 487
 488     theOBB.ReBuild(aListOfPnts);
 489   }
 490 }
 491
 492 //=======================================================================
 493 // Function : AddOBB
 494 // purpose :
 495 //=======================================================================
 496 void BRepBndLib::AddOBB(const TopoDS_Shape& theS,
 497                         Bnd_OBB& theOBB,
 498                         const Standard_Boolean theIsTriangulationUsed,
 499                         const Standard_Boolean theIsOptimal,
 500                         const Standard_Boolean theIsShapeToleranceUsed)
 501 {
 502   if (CheckPoints(theS, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed, theOBB))
 503     return;
 504
 505   ComputePCA(theS, theOBB, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed);
 506 }