src/GEOMUtils/GEOMUtils.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21
  22 #include <Standard_Stream.hxx>
  23
  24 #include <GEOMUtils.hxx>
  25
  26 #include <Basics_OCCTVersion.hxx>
  27
  28 #include <utilities.h>
  29 #include <OpUtil.hxx>
  30 #include <Utils_ExceptHandlers.hxx>
  31
  32 // OCCT Includes
  33 #include <BRepMesh_IncrementalMesh.hxx>
  34
  35 #include <BRep_Builder.hxx>
  36 #include <BRep_Tool.hxx>
  37 #include <BRepBndLib.hxx>
  38 #include <BRepGProp.hxx>
  39 #include <BRepTools.hxx>
  40
  41 #include <Bnd_Box.hxx>
  42
  43 #include <TopAbs.hxx>
  44 #include <TopExp.hxx>
  45 #include <TopoDS.hxx>
  46 #include <TopoDS_Edge.hxx>
  47 #include <TopoDS_Face.hxx>
  48 #include <TopoDS_Shape.hxx>
  49 #include <TopoDS_Vertex.hxx>
  50 #include <TopoDS_Compound.hxx>
  51 #include <TopoDS_Iterator.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_MapOfShape.hxx>
  54 #include <TopTools_ListOfShape.hxx>
  55 #include <TopTools_ListIteratorOfListOfShape.hxx>
  56
  57 #include <Geom_Surface.hxx>
  58 #include <Geom_Plane.hxx>
  59
  60 #include <GeomLProp_CLProps.hxx>
  61 #include <GeomLProp_SLProps.hxx>
  62
  63 #include <GProp_GProps.hxx>
  64 #include <GProp_PrincipalProps.hxx>
  65
  66 #include <gp_Pln.hxx>
  67 #include <gp_Lin.hxx>
  68
  69 #include <vector>
  70
  71 #include <Standard_Failure.hxx>
  72 #include <Standard_NullObject.hxx>
  73 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
  74
  75 #define STD_SORT_ALGO 1
  76
  77 //=======================================================================
  78 //function : GetPosition
  79 //purpose  :
  80 //=======================================================================
  81 gp_Ax3 GEOMUtils::GetPosition (const TopoDS_Shape& theShape)
  82 {
  83   gp_Ax3 aResult;
  84
  85   if (theShape.IsNull())
  86     return aResult;
  87
  88   // Axes
  89   aResult.Transform(theShape.Location().Transformation());
  90   if (theShape.ShapeType() == TopAbs_FACE) {
  91     Handle(Geom_Surface) aGS = BRep_Tool::Surface(TopoDS::Face(theShape));
  92     if (!aGS.IsNull() && aGS->IsKind(STANDARD_TYPE(Geom_Plane))) {
  93       Handle(Geom_Plane) aGPlane = Handle(Geom_Plane)::DownCast(aGS);
  94       gp_Pln aPln = aGPlane->Pln();
  95       aResult = aPln.Position();
  96       // In case of reverse orinetation of the face invert the plane normal
  97       // (the face's normal does not mathc the plane's normal in this case)
  98       if(theShape.Orientation() == TopAbs_REVERSED)
  99       {
 100         gp_Dir Vx =  aResult.XDirection();
 101         gp_Dir N  =  aResult.Direction().Mirrored(Vx);
 102         gp_Pnt P  =  aResult.Location();
 103         aResult = gp_Ax3(P, N, Vx);
 104       }
 105     }
 106   }
 107
 108   // Origin
 109   gp_Pnt aPnt;
 110
 111   TopAbs_ShapeEnum aShType = theShape.ShapeType();
 112
 113   if (aShType == TopAbs_VERTEX) {
 114     aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
 115   }
 116   else {
 117     if (aShType == TopAbs_COMPOUND) {
 118       aShType = GetTypeOfSimplePart(theShape);
 119     }
 120
 121     GProp_GProps aSystem;
 122     if (aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)
 123       BRepGProp::LinearProperties(theShape, aSystem);
 124     else if (aShType == TopAbs_FACE || aShType == TopAbs_SHELL)
 125       BRepGProp::SurfaceProperties(theShape, aSystem);
 126     else
 127       BRepGProp::VolumeProperties(theShape, aSystem);
 128
 129     aPnt = aSystem.CentreOfMass();
 130   }
 131
 132   aResult.SetLocation(aPnt);
 133
 134   return aResult;
 135 }
 136
 137 //=======================================================================
 138 //function : GetVector
 139 //purpose  :
 140 //=======================================================================
 141 gp_Vec GEOMUtils::GetVector (const TopoDS_Shape& theShape)
 142 {
 143   if (theShape.IsNull())
 144     Standard_NullObject::Raise("Null shape is given for a vector");
 145
 146   if (theShape.ShapeType() != TopAbs_EDGE)
 147     Standard_TypeMismatch::Raise("Invalid shape is given, must be a vector or an edge");
 148
 149   TopoDS_Edge anE = TopoDS::Edge(theShape);
 150
 151   TopoDS_Vertex V1, V2;
 152   TopExp::Vertices(anE, V1, V2, Standard_True);
 153
 154   if (V1.IsNull() || V2.IsNull())
 155     Standard_NullObject::Raise("Invalid edge is given, it must have two points");
 156
 157   gp_Vec aV (BRep_Tool::Pnt(V1), BRep_Tool::Pnt(V2));
 158   if (aV.Magnitude() < gp::Resolution()) {
 159     Standard_ConstructionError::Raise("Vector of zero length is given");
 160   }
 161
 162   return aV;
 163 }
 164
 165 //=======================================================================
 166 //function : ShapeToDouble
 167 //purpose  : used by CompareShapes::operator()
 168 //=======================================================================
 169 std::pair<double, double> ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
 170 {
 171   // Computing of CentreOfMass
 172   gp_Pnt GPoint;
 173   double Len;
 174
 175   if (S.ShapeType() == TopAbs_VERTEX) {
 176     GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
 177     Len = (double)S.Orientation();
 178   }
 179   else {
 180     GProp_GProps GPr;
 181     // BEGIN: fix for Mantis issue 0020842
 182     if (isOldSorting) {
 183       BRepGProp::LinearProperties(S, GPr);
 184     }
 185     else {
 186       if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 187         BRepGProp::LinearProperties(S, GPr);
 188       }
 189       else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 190         BRepGProp::SurfaceProperties(S, GPr);
 191       }
 192       else {
 193         BRepGProp::VolumeProperties(S, GPr);
 194       }
 195     }
 196     // END: fix for Mantis issue 0020842
 197     GPoint = GPr.CentreOfMass();
 198     Len = GPr.Mass();
 199   }
 200
 201   double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
 202   return std::make_pair(dMidXYZ, Len);
 203 }
 204
 205 //=======================================================================
 206 //function : CompareShapes::operator()
 207 //purpose  : used by std::sort(), called from SortShapes()
 208 //=======================================================================
 209 bool GEOMUtils::CompareShapes::operator() (const TopoDS_Shape& theShape1,
 210                                            const TopoDS_Shape& theShape2)
 211 {
 212   if (!myMap.IsBound(theShape1)) {
 213     myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
 214   }
 215
 216   if (!myMap.IsBound(theShape2)) {
 217     myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
 218   }
 219
 220   std::pair<double, double> val1 = myMap.Find(theShape1);
 221   std::pair<double, double> val2 = myMap.Find(theShape2);
 222
 223   double tol = Precision::Confusion();
 224   bool exchange = Standard_False;
 225
 226   double dMidXYZ = val1.first - val2.first;
 227   if (dMidXYZ >= tol) {
 228     exchange = Standard_True;
 229   }
 230   else if (Abs(dMidXYZ) < tol) {
 231     double dLength = val1.second - val2.second;
 232     if (dLength >= tol) {
 233       exchange = Standard_True;
 234     }
 235     else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
 236       // PAL17233
 237       // equal values possible on shapes such as two halves of a sphere and
 238       // a membrane inside the sphere
 239       Bnd_Box box1,box2;
 240       BRepBndLib::Add(theShape1, box1);
 241       if (!box1.IsVoid()) {
 242         BRepBndLib::Add(theShape2, box2);
 243         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 244         if (dSquareExtent >= tol) {
 245           exchange = Standard_True;
 246         }
 247         else if (Abs(dSquareExtent) < tol) {
 248           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 249           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 250           val1 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 251           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 252           val2 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 253           if ((val1 - val2) >= tol) {
 254             exchange = Standard_True;
 255           }
 256         }
 257       }
 258     }
 259   }
 260
 261   //return val1 < val2;
 262   return !exchange;
 263 }
 264
 265 //=======================================================================
 266 //function : SortShapes
 267 //purpose  :
 268 //=======================================================================
 269 void GEOMUtils::SortShapes (TopTools_ListOfShape& SL,
 270                             const Standard_Boolean isOldSorting)
 271 {
 272 #ifdef STD_SORT_ALGO
 273   std::vector<TopoDS_Shape> aShapesVec;
 274   aShapesVec.reserve(SL.Extent());
 275
 276   TopTools_ListIteratorOfListOfShape it (SL);
 277   for (; it.More(); it.Next()) {
 278     aShapesVec.push_back(it.Value());
 279   }
 280   SL.Clear();
 281
 282   CompareShapes shComp (isOldSorting);
 283   std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 284   //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 285
 286   std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
 287   for (; anIter != aShapesVec.end(); ++anIter) {
 288     SL.Append(*anIter);
 289   }
 290 #else
 291   // old implementation
 292   Standard_Integer MaxShapes = SL.Extent();
 293   TopTools_Array1OfShape  aShapes (1,MaxShapes);
 294   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
 295   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
 296   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
 297
 298   // Computing of CentreOfMass
 299   Standard_Integer Index;
 300   GProp_GProps GPr;
 301   gp_Pnt GPoint;
 302   TopTools_ListIteratorOfListOfShape it(SL);
 303   for (Index=1;  it.More();  Index++)
 304   {
 305     TopoDS_Shape S = it.Value();
 306     SL.Remove( it ); // == it.Next()
 307     aShapes(Index) = S;
 308     OrderInd.SetValue (Index, Index);
 309     if (S.ShapeType() == TopAbs_VERTEX) {
 310       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
 311       Length.SetValue( Index, (Standard_Real) S.Orientation());
 312     }
 313     else {
 314       // BEGIN: fix for Mantis issue 0020842
 315       if (isOldSorting) {
 316         BRepGProp::LinearProperties (S, GPr);
 317       }
 318       else {
 319         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 320           BRepGProp::LinearProperties (S, GPr);
 321         }
 322         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 323           BRepGProp::SurfaceProperties(S, GPr);
 324         }
 325         else {
 326           BRepGProp::VolumeProperties(S, GPr);
 327         }
 328       }
 329       // END: fix for Mantis issue 0020842
 330       GPoint = GPr.CentreOfMass();
 331       Length.SetValue(Index, GPr.Mass());
 332     }
 333     MidXYZ.SetValue(Index, GPoint.X()*999.0 + GPoint.Y()*99.0 + GPoint.Z()*0.9);
 334     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
 335   }
 336
 337   // Sorting
 338   Standard_Integer aTemp;
 339   Standard_Boolean exchange, Sort = Standard_True;
 340   Standard_Real    tol = Precision::Confusion();
 341   while (Sort)
 342   {
 343     Sort = Standard_False;
 344     for (Index=1; Index < MaxShapes; Index++)
 345     {
 346       exchange = Standard_False;
 347       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
 348       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
 349       if ( dMidXYZ >= tol ) {
 350 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
 351 //              << " d: " << dMidXYZ << endl;
 352         exchange = Standard_True;
 353       }
 354       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
 355 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
 356 //              << " d: " << dLength << endl;
 357         exchange = Standard_True;
 358       }
 359       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
 360                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
 361         // PAL17233
 362         // equal values possible on shapes such as two halves of a sphere and
 363         // a membrane inside the sphere
 364         Bnd_Box box1,box2;
 365         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
 366         if ( box1.IsVoid() ) continue;
 367         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
 368         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 369         if ( dSquareExtent >= tol ) {
 370 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
 371           exchange = Standard_True;
 372         }
 373         else if ( Abs(dSquareExtent) < tol ) {
 374           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 375           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 376           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 377           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 378           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 379           //exchange = val1 > val2;
 380           if ((val1 - val2) >= tol) {
 381             exchange = Standard_True;
 382           }
 383           //cout << "box: " << val1<<" > "<<val2 << endl;
 384         }
 385       }
 386
 387       if (exchange)
 388       {
 389 //         cout << "exchange " << Index << " & " << Index+1 << endl;
 390         aTemp = OrderInd(Index);
 391         OrderInd(Index) = OrderInd(Index+1);
 392         OrderInd(Index+1) = aTemp;
 393         Sort = Standard_True;
 394       }
 395     }
 396   }
 397
 398   for (Index=1; Index <= MaxShapes; Index++)
 399     SL.Append( aShapes( OrderInd(Index) ));
 400 #endif
 401 }
 402
 403 //=======================================================================
 404 //function : CompsolidToCompound
 405 //purpose  :
 406 //=======================================================================
 407 TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
 408 {
 409   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
 410     return theCompsolid;
 411   }
 412
 413   TopoDS_Compound aCompound;
 414   BRep_Builder B;
 415   B.MakeCompound(aCompound);
 416
 417   TopTools_MapOfShape mapShape;
 418   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
 419
 420   for (; It.More(); It.Next()) {
 421     TopoDS_Shape aShape_i = It.Value();
 422     if (mapShape.Add(aShape_i)) {
 423       B.Add(aCompound, aShape_i);
 424     }
 425   }
 426
 427   return aCompound;
 428 }
 429
 430 //=======================================================================
 431 //function : CheckTriangulation
 432 //purpose  :
 433 //=======================================================================
 434 bool GEOMUtils::CheckTriangulation (const TopoDS_Shape& aShape)
 435 {
 436   bool isTriangulation = true;
 437
 438   TopExp_Explorer exp (aShape, TopAbs_FACE);
 439   if (exp.More())
 440   {
 441     TopLoc_Location aTopLoc;
 442     Handle(Poly_Triangulation) aTRF;
 443     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
 444     if (aTRF.IsNull()) {
 445       isTriangulation = false;
 446     }
 447   }
 448   else // no faces, try edges
 449   {
 450     TopExp_Explorer expe (aShape, TopAbs_EDGE);
 451     if (!expe.More()) {
 452       return false;
 453     }
 454     TopLoc_Location aLoc;
 455     Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
 456     if (aPE.IsNull()) {
 457       isTriangulation = false;
 458     }
 459   }
 460
 461   if (!isTriangulation) {
 462     // calculate deflection
 463     Standard_Real aDeviationCoefficient = 0.001;
 464
 465     Bnd_Box B;
 466     BRepBndLib::Add(aShape, B);
 467     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 468     B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 469
 470     Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
 471     Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
 472     Standard_Real aHLRAngle = 0.349066;
 473
 474     BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
 475   }
 476
 477   return true;
 478 }
 479
 480 //=======================================================================
 481 //function : GetTypeOfSimplePart
 482 //purpose  :
 483 //=======================================================================
 484 TopAbs_ShapeEnum GEOMUtils::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
 485 {
 486   TopAbs_ShapeEnum aType = theShape.ShapeType();
 487   if      (aType == TopAbs_VERTEX)                             return TopAbs_VERTEX;
 488   else if (aType == TopAbs_EDGE  || aType == TopAbs_WIRE)      return TopAbs_EDGE;
 489   else if (aType == TopAbs_FACE  || aType == TopAbs_SHELL)     return TopAbs_FACE;
 490   else if (aType == TopAbs_SOLID || aType == TopAbs_COMPSOLID) return TopAbs_SOLID;
 491   else if (aType == TopAbs_COMPOUND) {
 492     // Only the iType of the first shape in the compound is taken into account
 493     TopoDS_Iterator It (theShape, Standard_False, Standard_False);
 494     if (It.More()) {
 495       return GetTypeOfSimplePart(It.Value());
 496     }
 497   }
 498   return TopAbs_SHAPE;
 499 }