src/GEOMUtils/GEOMUtils.cxx

   1 // Copyright (C) 2007-2015  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, or (at your option) any later version.
  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
  23 #include <Standard_Stream.hxx>
  24
  25 #include <GEOMUtils.hxx>
  26
  27 #include <Basics_OCCTVersion.hxx>
  28
  29 #include <OpUtil.hxx>
  30 #include <Utils_ExceptHandlers.hxx>
  31
  32 // OCCT Includes
  33 #include <BRepMesh_IncrementalMesh.hxx>
  34
  35 #include <BRepExtrema_DistShapeShape.hxx>
  36
  37 #include <BRep_Builder.hxx>
  38 #include <BRep_Tool.hxx>
  39 #include <BRepBndLib.hxx>
  40 #include <BRepGProp.hxx>
  41 #include <BRepTools.hxx>
  42
  43 #include <BRepClass3d_SolidClassifier.hxx>
  44
  45 #include <BRepBuilderAPI_MakeFace.hxx>
  46 #include <BRepBuilderAPI_Sewing.hxx>
  47
  48 #include <BRepCheck_Analyzer.hxx>
  49
  50 #include <Bnd_Box.hxx>
  51
  52 #include <BOPTools_AlgoTools.hxx>
  53
  54 #include <TopAbs.hxx>
  55 #include <TopExp.hxx>
  56 #include <TopoDS.hxx>
  57 #include <TopoDS_Edge.hxx>
  58 #include <TopoDS_Face.hxx>
  59 #include <TopoDS_Shape.hxx>
  60 #include <TopoDS_Vertex.hxx>
  61 #include <TopoDS_Compound.hxx>
  62 #include <TopoDS_Iterator.hxx>
  63 #include <TopExp_Explorer.hxx>
  64 #include <TopTools_MapOfShape.hxx>
  65 #include <TopTools_ListOfShape.hxx>
  66 #include <TopTools_ListIteratorOfListOfShape.hxx>
  67 #include <TopTools_Array1OfShape.hxx>
  68
  69 #include <Geom_Circle.hxx>
  70 #include <Geom_Surface.hxx>
  71 #include <Geom_Plane.hxx>
  72 #include <Geom_SphericalSurface.hxx>
  73 #include <Geom_ToroidalSurface.hxx>
  74 #include <Geom_RectangularTrimmedSurface.hxx>
  75
  76 #include <GeomLProp_CLProps.hxx>
  77 #include <GeomLProp_SLProps.hxx>
  78
  79 #include <GProp_GProps.hxx>
  80 #include <GProp_PrincipalProps.hxx>
  81
  82 #include <TColStd_Array1OfReal.hxx>
  83
  84 #include <gp_Pln.hxx>
  85 #include <gp_Lin.hxx>
  86
  87 #include <ShapeAnalysis.hxx>
  88 #include <ShapeFix_Shape.hxx>
  89 #include <ShapeFix_ShapeTolerance.hxx>
  90
  91 #include <ProjLib.hxx>
  92 #include <ElSLib.hxx>
  93
  94 #include <vector>
  95 #include <sstream>
  96 #include <algorithm>
  97
  98 #include <Standard_Failure.hxx>
  99 #include <Standard_NullObject.hxx>
 100 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 101
 102 #define MAX2(X, Y)    (Abs(X) > Abs(Y) ? Abs(X) : Abs(Y))
 103 #define MAX3(X, Y, Z) (MAX2(MAX2(X,Y), Z))
 104
 105 #define STD_SORT_ALGO 1
 106
 107 // When the following macro is defined, ShapeFix_ShapeTolerance function is used to set max tolerance of curve
 108 // in GEOMUtils::FixShapeCurves function; otherwise less restrictive BRep_Builder::UpdateEdge/UpdateVertex
 109 // approach is used
 110 // VSR (29/12/2014): macro disabled
 111 //#define USE_LIMIT_TOLERANCE
 112
 113 namespace
 114 {
 115   /**
 116    * This function constructs and returns modified shape from the original one
 117    * for singular cases. It is used for the method GetMinDistanceSingular.
 118    *
 119    * \param theShape the original shape
 120    * \param theModifiedShape output parameter. The modified shape.
 121    * \param theAddDist output parameter. The added distance for modified shape.
 122    * \retval true if the shape is modified; false otherwise.
 123    *
 124    * \internal
 125    */
 126   Standard_Boolean ModifyShape(const TopoDS_Shape  &theShape,
 127                                TopoDS_Shape  &theModifiedShape,
 128                                Standard_Real &theAddDist)
 129   {
 130     TopExp_Explorer anExp;
 131     int nbf = 0;
 132
 133     theAddDist = 0.;
 134     theModifiedShape.Nullify();
 135
 136     for ( anExp.Init( theShape, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
 137       nbf++;
 138       theModifiedShape = anExp.Current();
 139     }
 140     if(nbf==1) {
 141       TopoDS_Shape sh = theShape;
 142       while(sh.ShapeType()==TopAbs_COMPOUND) {
 143         TopoDS_Iterator it(sh);
 144         sh = it.Value();
 145       }
 146       Handle(Geom_Surface) S = BRep_Tool::Surface(TopoDS::Face(theModifiedShape));
 147       if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
 148           S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ||
 149           S->IsUPeriodic()) {
 150         const Standard_Boolean isShell =
 151           (sh.ShapeType()==TopAbs_SHELL || sh.ShapeType()==TopAbs_FACE);
 152
 153         if ( !isShell && S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
 154           Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
 155           gp_Pnt PC = SS->Location();
 156           BRep_Builder B;
 157           TopoDS_Vertex V;
 158           B.MakeVertex(V,PC,1.e-7);
 159           theModifiedShape = V;
 160           theAddDist = SS->Radius();
 161           return Standard_True;
 162         }
 163         if ( !isShell && S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
 164           Handle(Geom_ToroidalSurface) TS = Handle(Geom_ToroidalSurface)::DownCast(S);
 165           gp_Ax3 ax3 = TS->Position();
 166           Handle(Geom_Circle) C = new Geom_Circle(ax3.Ax2(),TS->MajorRadius());
 167           BRep_Builder B;
 168           TopoDS_Edge E;
 169           B.MakeEdge(E,C,1.e-7);
 170           theModifiedShape = E;
 171           theAddDist = TS->MinorRadius();
 172           return Standard_True;
 173         }
 174
 175         // non solid case or any periodic surface (Mantis 22454).
 176         double U1,U2,V1,V2;
 177         // changes for 0020677: EDF 1219 GEOM: MinDistance gives 0 instead of 20.88
 178         //S->Bounds(U1,U2,V1,V2); changed by
 179         ShapeAnalysis::GetFaceUVBounds(TopoDS::Face(theModifiedShape),U1,U2,V1,V2);
 180         // end of changes for 020677 (dmv)
 181         Handle(Geom_RectangularTrimmedSurface) TrS1 =
 182           new Geom_RectangularTrimmedSurface(S,U1,(U1+U2)/2.,V1,V2);
 183         Handle(Geom_RectangularTrimmedSurface) TrS2 =
 184           new Geom_RectangularTrimmedSurface(S,(U1+U2)/2.,U2,V1,V2);
 185         TopoDS_Shape aMShape;
 186
 187         TopoDS_Face F1 = BRepBuilderAPI_MakeFace(TrS1, Precision::Confusion());
 188         TopoDS_Face F2 = BRepBuilderAPI_MakeFace(TrS2, Precision::Confusion());
 189
 190         if (isShell) {
 191           BRep_Builder B;
 192           B.MakeCompound(TopoDS::Compound(aMShape));
 193           B.Add(aMShape, F1);
 194           B.Add(aMShape, F2);
 195         } else {
 196           // The original shape is a solid.
 197           BRepBuilderAPI_Sewing aSewing (Precision::Confusion()*10.0);
 198           aSewing.Add(F1);
 199           aSewing.Add(F2);
 200           aSewing.Perform();
 201           aMShape = aSewing.SewedShape();
 202           BRep_Builder B;
 203           TopoDS_Solid aSolid;
 204           B.MakeSolid(aSolid);
 205           B.Add(aSolid, aMShape);
 206           aMShape = aSolid;
 207         }
 208
 209         Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
 210         sfs->Init(aMShape);
 211         sfs->SetPrecision(1.e-6);
 212         sfs->SetMaxTolerance(1.0);
 213         sfs->Perform();
 214         theModifiedShape = sfs->Shape();
 215         return Standard_True;
 216       }
 217     }
 218
 219     theModifiedShape = theShape;
 220     return Standard_False;
 221   }
 222
 223   void parseWard( const GEOMUtils::LevelsList &theLevelList, std::string &treeStr )
 224   {
 225     treeStr.append( "{" );
 226     for( GEOMUtils::LevelsList::const_iterator j = theLevelList.begin();
 227          j != theLevelList.end(); ++j ) {
 228       if ( j != theLevelList.begin() ) {
 229         treeStr.append( ";" );
 230       }
 231       GEOMUtils::LevelInfo level = (*j);
 232       GEOMUtils::LevelInfo::iterator upIter;
 233       for ( upIter = level.begin(); upIter != level.end(); ++upIter ) {
 234         if ( upIter != level.begin() ) {
 235           treeStr.append( "," );
 236         }
 237         treeStr.append( upIter->first );
 238         for ( std::vector<std::string>::iterator k = upIter->second.begin(); k != upIter->second.end(); ++k ) {
 239           treeStr.append( "_" );
 240           treeStr.append( *k );
 241         }
 242       }
 243     }
 244     treeStr.append( "}" );
 245   }
 246
 247   GEOMUtils::LevelsList parseWard( const std::string& theData, std::size_t& theCursor )
 248   {
 249     std::size_t indexStart = theData.find( "{", theCursor ) + 1;
 250     std::size_t indexEnd = theData.find( "}", indexStart );
 251
 252     std::string ward = theData.substr( indexStart, indexEnd - indexStart );
 253     std::stringstream ss(ward);
 254     std::string substr;
 255     std::vector<std::string> levelsListStr;
 256     while ( std::getline( ss, substr, ';' ) ) {
 257       if ( !substr.empty() )
 258         levelsListStr.push_back( substr );
 259     }
 260     GEOMUtils::LevelsList levelsListData;
 261     for( int level = 0; level < levelsListStr.size(); level++ ) {
 262       std::vector<std::string> namesListStr;
 263       std::stringstream ss1( levelsListStr[level] );
 264       while ( std::getline( ss1, substr, ',' ) ) {
 265         if ( !substr.empty() )
 266           namesListStr.push_back( substr );
 267       }
 268       GEOMUtils::LevelInfo levelInfoData;
 269       for( int node = 0; node < namesListStr.size(); node++ ) {
 270         std::vector<std::string> linksListStr;
 271         std::stringstream ss2( namesListStr[node] );
 272         while ( std::getline( ss2, substr, '_' ) ) {
 273           if ( !substr.empty() )
 274             linksListStr.push_back( substr );
 275         }
 276         std::string nodeItem = linksListStr[0];
 277         if( !nodeItem.empty() ) {
 278           GEOMUtils::NodeLinks linksListData;
 279           for( int link = 1; link < linksListStr.size(); link++ ) {
 280             std::string linkItem = linksListStr[link];
 281             linksListData.push_back( linkItem );
 282           }// Links
 283           levelInfoData[nodeItem] = linksListData;
 284         }
 285       }// Level's objects
 286       levelsListData.push_back(levelInfoData);
 287     }// Levels
 288
 289     theCursor = indexEnd + 1;
 290     return levelsListData;
 291   }
 292
 293 }
 294
 295 //=======================================================================
 296 //function : ShapeToDouble
 297 //purpose  : used by CompareShapes::operator()
 298 //=======================================================================
 299 std::pair<double, double> GEOMUtils::ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
 300 {
 301   // Computing of CentreOfMass
 302   gp_Pnt GPoint;
 303   double Len;
 304
 305   if (S.ShapeType() == TopAbs_VERTEX) {
 306     GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
 307     Len = (double)S.Orientation();
 308   }
 309   else {
 310     GProp_GProps GPr;
 311     // BEGIN: fix for Mantis issue 0020842
 312     if (isOldSorting) {
 313       BRepGProp::LinearProperties(S, GPr);
 314     }
 315     else {
 316       if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 317         BRepGProp::LinearProperties(S, GPr);
 318       }
 319       else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 320         BRepGProp::SurfaceProperties(S, GPr);
 321       }
 322       else {
 323         BRepGProp::VolumeProperties(S, GPr);
 324       }
 325     }
 326     // END: fix for Mantis issue 0020842
 327     GPoint = GPr.CentreOfMass();
 328     Len = GPr.Mass();
 329   }
 330
 331   double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
 332   return std::make_pair(dMidXYZ, Len);
 333 }
 334
 335 //=======================================================================
 336 //function : GetPosition
 337 //purpose  :
 338 //=======================================================================
 339 gp_Ax3 GEOMUtils::GetPosition (const TopoDS_Shape& theShape)
 340 {
 341   gp_Ax3 aResult;
 342
 343   if (theShape.IsNull())
 344     return aResult;
 345
 346   // Axes
 347   aResult.Transform(theShape.Location().Transformation());
 348   if (theShape.ShapeType() == TopAbs_FACE) {
 349     Handle(Geom_Surface) aGS = BRep_Tool::Surface(TopoDS::Face(theShape));
 350     if (!aGS.IsNull() && aGS->IsKind(STANDARD_TYPE(Geom_Plane))) {
 351       Handle(Geom_Plane) aGPlane = Handle(Geom_Plane)::DownCast(aGS);
 352       gp_Pln aPln = aGPlane->Pln();
 353       aResult = aPln.Position();
 354       // In case of reverse orinetation of the face invert the plane normal
 355       // (the face's normal does not mathc the plane's normal in this case)
 356       if(theShape.Orientation() == TopAbs_REVERSED)
 357       {
 358         gp_Dir Vx =  aResult.XDirection();
 359         gp_Dir N  =  aResult.Direction().Mirrored(Vx);
 360         gp_Pnt P  =  aResult.Location();
 361         aResult = gp_Ax3(P, N, Vx);
 362       }
 363     }
 364   }
 365
 366   // Origin
 367   gp_Pnt aPnt;
 368
 369   TopAbs_ShapeEnum aShType = theShape.ShapeType();
 370
 371   if (aShType == TopAbs_VERTEX) {
 372     aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
 373   }
 374   else {
 375     if (aShType == TopAbs_COMPOUND) {
 376       aShType = GetTypeOfSimplePart(theShape);
 377     }
 378
 379     GProp_GProps aSystem;
 380     if (aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)
 381       BRepGProp::LinearProperties(theShape, aSystem);
 382     else if (aShType == TopAbs_FACE || aShType == TopAbs_SHELL)
 383       BRepGProp::SurfaceProperties(theShape, aSystem);
 384     else
 385       BRepGProp::VolumeProperties(theShape, aSystem);
 386
 387     aPnt = aSystem.CentreOfMass();
 388   }
 389
 390   aResult.SetLocation(aPnt);
 391
 392   return aResult;
 393 }
 394
 395 //=======================================================================
 396 //function : GetVector
 397 //purpose  :
 398 //=======================================================================
 399 gp_Vec GEOMUtils::GetVector (const TopoDS_Shape& theShape,
 400                              Standard_Boolean doConsiderOrientation)
 401 {
 402   if (theShape.IsNull())
 403     Standard_NullObject::Raise("Null shape is given for a vector");
 404
 405   if (theShape.ShapeType() != TopAbs_EDGE)
 406     Standard_TypeMismatch::Raise("Invalid shape is given, must be a vector or an edge");
 407
 408   TopoDS_Edge anE = TopoDS::Edge(theShape);
 409
 410   TopoDS_Vertex V1, V2;
 411   TopExp::Vertices(anE, V1, V2, doConsiderOrientation);
 412
 413   if (V1.IsNull() || V2.IsNull())
 414     Standard_NullObject::Raise("Invalid edge is given, it must have two points");
 415
 416   gp_Vec aV (BRep_Tool::Pnt(V1), BRep_Tool::Pnt(V2));
 417   if (aV.Magnitude() < gp::Resolution()) {
 418     Standard_ConstructionError::Raise("Vector of zero length is given");
 419   }
 420
 421   return aV;
 422 }
 423
 424 //=======================================================================
 425 //function : CompareShapes::operator()
 426 //purpose  : used by std::sort(), called from SortShapes()
 427 //=======================================================================
 428 bool GEOMUtils::CompareShapes::operator() (const TopoDS_Shape& theShape1,
 429                                            const TopoDS_Shape& theShape2)
 430 {
 431   if (!myMap.IsBound(theShape1)) {
 432     myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
 433   }
 434
 435   if (!myMap.IsBound(theShape2)) {
 436     myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
 437   }
 438
 439   std::pair<double, double> val1 = myMap.Find(theShape1);
 440   std::pair<double, double> val2 = myMap.Find(theShape2);
 441
 442   double tol = Precision::Confusion();
 443   bool exchange = Standard_False;
 444
 445   double dMidXYZ = val1.first - val2.first;
 446   if (dMidXYZ >= tol) {
 447     exchange = Standard_True;
 448   }
 449   else if (Abs(dMidXYZ) < tol) {
 450     double dLength = val1.second - val2.second;
 451     if (dLength >= tol) {
 452       exchange = Standard_True;
 453     }
 454     else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
 455       // PAL17233
 456       // equal values possible on shapes such as two halves of a sphere and
 457       // a membrane inside the sphere
 458       Bnd_Box box1,box2;
 459       BRepBndLib::Add(theShape1, box1);
 460       if (!box1.IsVoid()) {
 461         BRepBndLib::Add(theShape2, box2);
 462         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 463         if (dSquareExtent >= tol) {
 464           exchange = Standard_True;
 465         }
 466         else if (Abs(dSquareExtent) < tol) {
 467           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 468           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 469           val1 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 470           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 471           val2 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 472           if ((val1 - val2) >= tol) {
 473             exchange = Standard_True;
 474           }
 475         }
 476       }
 477     }
 478   }
 479
 480   //return val1 < val2;
 481   return !exchange;
 482 }
 483
 484 //=======================================================================
 485 //function : SortShapes
 486 //purpose  :
 487 //=======================================================================
 488 void GEOMUtils::SortShapes (TopTools_ListOfShape& SL,
 489                             const Standard_Boolean isOldSorting)
 490 {
 491 #ifdef STD_SORT_ALGO
 492   std::vector<TopoDS_Shape> aShapesVec;
 493   aShapesVec.reserve(SL.Extent());
 494
 495   TopTools_ListIteratorOfListOfShape it (SL);
 496   for (; it.More(); it.Next()) {
 497     aShapesVec.push_back(it.Value());
 498   }
 499   SL.Clear();
 500
 501   CompareShapes shComp (isOldSorting);
 502   std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 503   //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 504
 505   std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
 506   for (; anIter != aShapesVec.end(); ++anIter) {
 507     SL.Append(*anIter);
 508   }
 509 #else
 510   // old implementation
 511   Standard_Integer MaxShapes = SL.Extent();
 512   TopTools_Array1OfShape  aShapes (1,MaxShapes);
 513   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
 514   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
 515   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
 516
 517   // Computing of CentreOfMass
 518   Standard_Integer Index;
 519   GProp_GProps GPr;
 520   gp_Pnt GPoint;
 521   TopTools_ListIteratorOfListOfShape it(SL);
 522   for (Index=1;  it.More();  Index++)
 523   {
 524     TopoDS_Shape S = it.Value();
 525     SL.Remove( it ); // == it.Next()
 526     aShapes(Index) = S;
 527     OrderInd.SetValue (Index, Index);
 528     if (S.ShapeType() == TopAbs_VERTEX) {
 529       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
 530       Length.SetValue( Index, (Standard_Real) S.Orientation());
 531     }
 532     else {
 533       // BEGIN: fix for Mantis issue 0020842
 534       if (isOldSorting) {
 535         BRepGProp::LinearProperties (S, GPr);
 536       }
 537       else {
 538         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 539           BRepGProp::LinearProperties (S, GPr);
 540         }
 541         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 542           BRepGProp::SurfaceProperties(S, GPr);
 543         }
 544         else {
 545           BRepGProp::VolumeProperties(S, GPr);
 546         }
 547       }
 548       // END: fix for Mantis issue 0020842
 549       GPoint = GPr.CentreOfMass();
 550       Length.SetValue(Index, GPr.Mass());
 551     }
 552     MidXYZ.SetValue(Index, GPoint.X()*999.0 + GPoint.Y()*99.0 + GPoint.Z()*0.9);
 553     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
 554   }
 555
 556   // Sorting
 557   Standard_Integer aTemp;
 558   Standard_Boolean exchange, Sort = Standard_True;
 559   Standard_Real    tol = Precision::Confusion();
 560   while (Sort)
 561   {
 562     Sort = Standard_False;
 563     for (Index=1; Index < MaxShapes; Index++)
 564     {
 565       exchange = Standard_False;
 566       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
 567       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
 568       if ( dMidXYZ >= tol ) {
 569 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
 570 //              << " d: " << dMidXYZ << endl;
 571         exchange = Standard_True;
 572       }
 573       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
 574 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
 575 //              << " d: " << dLength << endl;
 576         exchange = Standard_True;
 577       }
 578       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
 579                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
 580         // PAL17233
 581         // equal values possible on shapes such as two halves of a sphere and
 582         // a membrane inside the sphere
 583         Bnd_Box box1,box2;
 584         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
 585         if ( box1.IsVoid() ) continue;
 586         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
 587         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 588         if ( dSquareExtent >= tol ) {
 589 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
 590           exchange = Standard_True;
 591         }
 592         else if ( Abs(dSquareExtent) < tol ) {
 593           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 594           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 595           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 596           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 597           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 598           //exchange = val1 > val2;
 599           if ((val1 - val2) >= tol) {
 600             exchange = Standard_True;
 601           }
 602           //cout << "box: " << val1<<" > "<<val2 << endl;
 603         }
 604       }
 605
 606       if (exchange)
 607       {
 608 //         cout << "exchange " << Index << " & " << Index+1 << endl;
 609         aTemp = OrderInd(Index);
 610         OrderInd(Index) = OrderInd(Index+1);
 611         OrderInd(Index+1) = aTemp;
 612         Sort = Standard_True;
 613       }
 614     }
 615   }
 616
 617   for (Index=1; Index <= MaxShapes; Index++)
 618     SL.Append( aShapes( OrderInd(Index) ));
 619 #endif
 620 }
 621
 622 //=======================================================================
 623 //function : CompsolidToCompound
 624 //purpose  :
 625 //=======================================================================
 626 TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
 627 {
 628   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
 629     return theCompsolid;
 630   }
 631
 632   TopoDS_Compound aCompound;
 633   BRep_Builder B;
 634   B.MakeCompound(aCompound);
 635
 636   TopTools_MapOfShape mapShape;
 637   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
 638
 639   for (; It.More(); It.Next()) {
 640     TopoDS_Shape aShape_i = It.Value();
 641     if (mapShape.Add(aShape_i)) {
 642       B.Add(aCompound, aShape_i);
 643     }
 644   }
 645
 646   return aCompound;
 647 }
 648
 649 //=======================================================================
 650 //function : AddSimpleShapes
 651 //purpose  :
 652 //=======================================================================
 653 void GEOMUtils::AddSimpleShapes (const TopoDS_Shape& theShape, TopTools_ListOfShape& theList)
 654 {
 655   if (theShape.ShapeType() != TopAbs_COMPOUND &&
 656       theShape.ShapeType() != TopAbs_COMPSOLID) {
 657     theList.Append(theShape);
 658     return;
 659   }
 660
 661   TopTools_MapOfShape mapShape;
 662   TopoDS_Iterator It (theShape, Standard_True, Standard_True);
 663
 664   for (; It.More(); It.Next()) {
 665     TopoDS_Shape aShape_i = It.Value();
 666     if (mapShape.Add(aShape_i)) {
 667       if (aShape_i.ShapeType() == TopAbs_COMPOUND ||
 668           aShape_i.ShapeType() == TopAbs_COMPSOLID) {
 669         AddSimpleShapes(aShape_i, theList);
 670       } else {
 671         theList.Append(aShape_i);
 672       }
 673     }
 674   }
 675 }
 676
 677 //=======================================================================
 678 //function : CheckTriangulation
 679 //purpose  :
 680 //=======================================================================
 681 bool GEOMUtils::CheckTriangulation (const TopoDS_Shape& aShape)
 682 {
 683   bool isTriangulation = true;
 684
 685   TopExp_Explorer exp (aShape, TopAbs_FACE);
 686   if (exp.More())
 687   {
 688     TopLoc_Location aTopLoc;
 689     Handle(Poly_Triangulation) aTRF;
 690     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
 691     if (aTRF.IsNull()) {
 692       isTriangulation = false;
 693     }
 694   }
 695   else // no faces, try edges
 696   {
 697     TopExp_Explorer expe (aShape, TopAbs_EDGE);
 698     if (!expe.More()) {
 699       return false;
 700     }
 701     TopLoc_Location aLoc;
 702     Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
 703     if (aPE.IsNull()) {
 704       isTriangulation = false;
 705     }
 706   }
 707
 708   if (!isTriangulation) {
 709     // calculate deflection
 710     Standard_Real aDeviationCoefficient = 0.001;
 711
 712     Bnd_Box B;
 713     BRepBndLib::Add(aShape, B);
 714     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 715     B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 716
 717     Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
 718     Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
 719     Standard_Real aHLRAngle = 0.349066;
 720
 721     BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
 722   }
 723
 724   return true;
 725 }
 726
 727 //=======================================================================
 728 //function : GetTypeOfSimplePart
 729 //purpose  :
 730 //=======================================================================
 731 TopAbs_ShapeEnum GEOMUtils::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
 732 {
 733   TopAbs_ShapeEnum aType = theShape.ShapeType();
 734   if      (aType == TopAbs_VERTEX)                             return TopAbs_VERTEX;
 735   else if (aType == TopAbs_EDGE  || aType == TopAbs_WIRE)      return TopAbs_EDGE;
 736   else if (aType == TopAbs_FACE  || aType == TopAbs_SHELL)     return TopAbs_FACE;
 737   else if (aType == TopAbs_SOLID || aType == TopAbs_COMPSOLID) return TopAbs_SOLID;
 738   else if (aType == TopAbs_COMPOUND) {
 739     // Only the iType of the first shape in the compound is taken into account
 740     TopoDS_Iterator It (theShape, Standard_False, Standard_False);
 741     if (It.More()) {
 742       return GetTypeOfSimplePart(It.Value());
 743     }
 744   }
 745   return TopAbs_SHAPE;
 746 }
 747
 748 //=======================================================================
 749 //function : GetEdgeNearPoint
 750 //purpose  :
 751 //=======================================================================
 752 TopoDS_Shape GEOMUtils::GetEdgeNearPoint (const TopoDS_Shape& theShape,
 753                                           const TopoDS_Vertex& thePoint)
 754 {
 755   TopoDS_Shape aResult;
 756
 757   // 1. Explode the shape on edges
 758   TopTools_MapOfShape mapShape;
 759   Standard_Integer nbEdges = 0;
 760   TopExp_Explorer exp (theShape, TopAbs_EDGE);
 761   for (; exp.More(); exp.Next()) {
 762     if (mapShape.Add(exp.Current())) {
 763       nbEdges++;
 764     }
 765   }
 766
 767   if (nbEdges == 0)
 768     Standard_NullObject::Raise("Given shape contains no edges");
 769
 770   mapShape.Clear();
 771   Standard_Integer ind = 1;
 772   TopTools_Array1OfShape anEdges (1, nbEdges);
 773   TColStd_Array1OfReal aDistances (1, nbEdges);
 774   for (exp.Init(theShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 775     if (mapShape.Add(exp.Current())) {
 776       TopoDS_Shape anEdge = exp.Current();
 777       anEdges(ind) = anEdge;
 778
 779       // 2. Classify the point relatively each edge
 780       BRepExtrema_DistShapeShape aDistTool (thePoint, anEdges(ind));
 781       if (!aDistTool.IsDone())
 782         Standard_ConstructionError::Raise("Cannot find a distance from the given point to one of edges");
 783
 784       aDistances(ind) = aDistTool.Value();
 785       ind++;
 786     }
 787   }
 788
 789   // 3. Define edge, having minimum distance to the point
 790   Standard_Real nearest = RealLast(), nbFound = 0;
 791   Standard_Real prec = Precision::Confusion();
 792   for (ind = 1; ind <= nbEdges; ind++) {
 793     if (Abs(aDistances(ind) - nearest) < prec) {
 794       nbFound++;
 795     }
 796     else if (aDistances(ind) < nearest) {
 797       nearest = aDistances(ind);
 798       aResult = anEdges(ind);
 799       nbFound = 1;
 800     }
 801     else {
 802     }
 803   }
 804   if (nbFound > 1) {
 805     Standard_ConstructionError::Raise("Multiple edges near the given point are found");
 806   }
 807   else if (nbFound == 0) {
 808     Standard_ConstructionError::Raise("There are no edges near the given point");
 809   }
 810   else {
 811   }
 812
 813   return aResult;
 814 }
 815
 816 //=======================================================================
 817 //function : PreciseBoundingBox
 818 //purpose  :
 819 //=======================================================================
 820 Standard_Boolean GEOMUtils::PreciseBoundingBox
 821                           (const TopoDS_Shape &theShape, Bnd_Box &theBox)
 822 {
 823   if ( theBox.IsVoid() ) BRepBndLib::Add( theShape, theBox );
 824   if ( theBox.IsVoid() ) return Standard_False;
 825
 826   Standard_Real aBound[6];
 827   theBox.Get(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 828
 829   Standard_Integer i;
 830   const gp_Pnt aMid(0.5*(aBound[1] + aBound[0]),  // XMid
 831                     0.5*(aBound[3] + aBound[2]),  // YMid
 832                     0.5*(aBound[5] + aBound[4])); // ZMid
 833   const gp_XYZ aSize(aBound[1] - aBound[0],       // DX
 834                      aBound[3] - aBound[2],       // DY
 835                      aBound[5] - aBound[4]);      // DZ
 836   const gp_Pnt aPnt[6] =
 837     {
 838       gp_Pnt(aBound[0] - (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMin
 839       gp_Pnt(aBound[1] + (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMax
 840       gp_Pnt(aMid.X(), aBound[2] - (aBound[3] - aBound[2]), aMid.Z()), // YMin
 841       gp_Pnt(aMid.X(), aBound[3] + (aBound[3] - aBound[2]), aMid.Z()), // YMax
 842       gp_Pnt(aMid.X(), aMid.Y(), aBound[4] - (aBound[5] - aBound[4])), // ZMin
 843       gp_Pnt(aMid.X(), aMid.Y(), aBound[5] + (aBound[5] - aBound[4]))  // ZMax
 844     };
 845   const gp_Dir aDir[3] = { gp::DX(), gp::DY(), gp::DZ() };
 846   const Standard_Real aPlnSize[3] =
 847     {
 848       0.5*Max(aSize.Y(), aSize.Z()), // XMin, XMax planes
 849       0.5*Max(aSize.X(), aSize.Z()), // YMin, YMax planes
 850       0.5*Max(aSize.X(), aSize.Y())  // ZMin, ZMax planes
 851     };
 852   gp_Pnt aPMin[2];
 853
 854   for (i = 0; i < 6; i++) {
 855     const Standard_Integer iHalf = i/2;
 856     const gp_Pln aPln(aPnt[i], aDir[iHalf]);
 857     BRepBuilderAPI_MakeFace aMkFace(aPln, -aPlnSize[iHalf], aPlnSize[iHalf],
 858                                     -aPlnSize[iHalf], aPlnSize[iHalf]);
 859
 860     if (!aMkFace.IsDone()) {
 861       return Standard_False;
 862     }
 863
 864     TopoDS_Shape aFace = aMkFace.Shape();
 865
 866     // Get minimal distance between planar face and shape.
 867     Standard_Real aMinDist =
 868       GEOMUtils::GetMinDistance(aFace, theShape, aPMin[0], aPMin[1]);
 869
 870     if (aMinDist < 0.) {
 871       return Standard_False;
 872     }
 873
 874     aBound[i] = aPMin[1].Coord(iHalf + 1);
 875   }
 876
 877   // Update Bounding box with the new values.
 878   theBox.SetVoid();
 879   theBox.Update(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 880
 881   return Standard_True;
 882 }
 883
 884 //=======================================================================
 885 //function : GetMinDistanceSingular
 886 //purpose  :
 887 //=======================================================================
 888 double GEOMUtils::GetMinDistanceSingular(const TopoDS_Shape& aSh1,
 889                                          const TopoDS_Shape& aSh2,
 890                                          gp_Pnt& Ptmp1, gp_Pnt& Ptmp2)
 891 {
 892   TopoDS_Shape     tmpSh1;
 893   TopoDS_Shape     tmpSh2;
 894   Standard_Real    AddDist1 = 0.;
 895   Standard_Real    AddDist2 = 0.;
 896   Standard_Boolean IsChange1 = ModifyShape(aSh1, tmpSh1, AddDist1);
 897   Standard_Boolean IsChange2 = ModifyShape(aSh2, tmpSh2, AddDist2);
 898
 899   if( !IsChange1 && !IsChange2 )
 900     return -2.0;
 901
 902   BRepExtrema_DistShapeShape dst(tmpSh1,tmpSh2);
 903   if (dst.IsDone()) {
 904     double MinDist = 1.e9;
 905     gp_Pnt PMin1, PMin2, P1, P2;
 906     for (int i = 1; i <= dst.NbSolution(); i++) {
 907       P1 = dst.PointOnShape1(i);
 908       P2 = dst.PointOnShape2(i);
 909       Standard_Real Dist = P1.Distance(P2);
 910       if (MinDist > Dist) {
 911         MinDist = Dist;
 912         PMin1 = P1;
 913         PMin2 = P2;
 914       }
 915     }
 916     if(MinDist<1.e-7) {
 917       Ptmp1 = PMin1;
 918       Ptmp2 = PMin2;
 919     }
 920     else {
 921       gp_Dir aDir(gp_Vec(PMin1,PMin2));
 922       if( MinDist > (AddDist1+AddDist2) ) {
 923         Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 924                         PMin1.Y() + aDir.Y()*AddDist1,
 925                         PMin1.Z() + aDir.Z()*AddDist1 );
 926         Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 927                         PMin2.Y() - aDir.Y()*AddDist2,
 928                         PMin2.Z() - aDir.Z()*AddDist2 );
 929         return (MinDist - AddDist1 - AddDist2);
 930       }
 931       else {
 932         if( AddDist1 > 0 ) {
 933           Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 934                           PMin1.Y() + aDir.Y()*AddDist1,
 935                           PMin1.Z() + aDir.Z()*AddDist1 );
 936           Ptmp2 = Ptmp1;
 937         }
 938         else {
 939           Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 940                           PMin2.Y() - aDir.Y()*AddDist2,
 941                           PMin2.Z() - aDir.Z()*AddDist2 );
 942           Ptmp1 = Ptmp2;
 943         }
 944       }
 945     }
 946     double res = MinDist - AddDist1 - AddDist2;
 947     if(res<0.) res = 0.0;
 948     return res;
 949   }
 950   return -2.0;
 951 }
 952
 953 //=======================================================================
 954 //function : GetMinDistance
 955 //purpose  :
 956 //=======================================================================
 957 Standard_Real GEOMUtils::GetMinDistance
 958                                (const TopoDS_Shape& theShape1,
 959                                 const TopoDS_Shape& theShape2,
 960                                 gp_Pnt& thePnt1, gp_Pnt& thePnt2)
 961 {
 962   Standard_Real aResult = 1.e9;
 963
 964   // Issue 0020231: A min distance bug with torus and vertex.
 965   // Make GetMinDistance() return zero if a sole VERTEX is inside any of SOLIDs
 966
 967   // which of shapes consists of only one vertex?
 968   TopExp_Explorer exp1(theShape1,TopAbs_VERTEX), exp2(theShape2,TopAbs_VERTEX);
 969   TopoDS_Shape V1 = exp1.More() ? exp1.Current() : TopoDS_Shape();
 970   TopoDS_Shape V2 = exp2.More() ? exp2.Current() : TopoDS_Shape();
 971   exp1.Next(); exp2.Next();
 972   if ( exp1.More() ) V1.Nullify();
 973   if ( exp2.More() ) V2.Nullify();
 974   // vertex and container of solids
 975   TopoDS_Shape V = V1.IsNull() ? V2 : V1;
 976   TopoDS_Shape S = V1.IsNull() ? theShape1 : theShape2;
 977   if ( !V.IsNull() ) {
 978     // classify vertex against solids
 979     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( V ) );
 980     for ( exp1.Init( S, TopAbs_SOLID ); exp1.More(); exp1.Next() ) {
 981       BRepClass3d_SolidClassifier classifier( exp1.Current(), p, 1e-6);
 982       if ( classifier.State() == TopAbs_IN ) {
 983         thePnt1 = p;
 984         thePnt2 = p;
 985         return 0.0;
 986       }
 987     }
 988   }
 989   // End Issue 0020231
 990
 991   // skl 30.06.2008
 992   // additional workaround for bugs 19899, 19908 and 19910 from Mantis
 993   double dist = GEOMUtils::GetMinDistanceSingular
 994       (theShape1, theShape2, thePnt1, thePnt2);
 995
 996   if (dist > -1.0) {
 997     return dist;
 998   }
 999
1000   BRepExtrema_DistShapeShape dst (theShape1, theShape2);
1001   if (dst.IsDone()) {
1002     gp_Pnt P1, P2;
1003
1004     for (int i = 1; i <= dst.NbSolution(); i++) {
1005       P1 = dst.PointOnShape1(i);
1006       P2 = dst.PointOnShape2(i);
1007
1008       Standard_Real Dist = P1.Distance(P2);
1009       if (aResult > Dist) {
1010         aResult = Dist;
1011         thePnt1 = P1;
1012         thePnt2 = P2;
1013       }
1014     }
1015   }
1016
1017   return aResult;
1018 }
1019
1020 //=======================================================================
1021 // function : ConvertClickToPoint()
1022 // purpose  : Returns the point clicked in 3D view
1023 //=======================================================================
1024 gp_Pnt GEOMUtils::ConvertClickToPoint( int x, int y, Handle(V3d_View) aView )
1025 {
1026   V3d_Coordinate XEye, YEye, ZEye, XAt, YAt, ZAt;
1027   aView->Eye( XEye, YEye, ZEye );
1028
1029   aView->At( XAt, YAt, ZAt );
1030   gp_Pnt EyePoint( XEye, YEye, ZEye );
1031   gp_Pnt AtPoint( XAt, YAt, ZAt );
1032
1033   gp_Vec EyeVector( EyePoint, AtPoint );
1034   gp_Dir EyeDir( EyeVector );
1035
1036   gp_Pln PlaneOfTheView = gp_Pln( AtPoint, EyeDir );
1037   Standard_Real X, Y, Z;
1038   aView->Convert( x, y, X, Y, Z );
1039   gp_Pnt ConvertedPoint( X, Y, Z );
1040
1041   gp_Pnt2d ConvertedPointOnPlane = ProjLib::Project( PlaneOfTheView, ConvertedPoint );
1042   gp_Pnt ResultPoint = ElSLib::Value( ConvertedPointOnPlane.X(), ConvertedPointOnPlane.Y(), PlaneOfTheView );
1043   return ResultPoint;
1044 }
1045
1046 //=======================================================================
1047 // function : ConvertTreeToString()
1048 // purpose  : Returns the string representation of dependency tree
1049 //=======================================================================
1050 void GEOMUtils::ConvertTreeToString( const TreeModel& tree,
1051                                      std::string& dependencyStr )
1052 {
1053   TreeModel::const_iterator i;
1054   for ( i = tree.begin(); i != tree.end(); ++i ) {
1055     dependencyStr.append( i->first );
1056     dependencyStr.append( "-" );
1057     std::vector<LevelInfo> upLevelList = i->second.first;
1058     dependencyStr.append( "upward" );
1059     parseWard( upLevelList, dependencyStr );
1060     std::vector<LevelInfo> downLevelList = i->second.second;
1061     dependencyStr.append( "downward" );
1062     parseWard( downLevelList, dependencyStr );
1063   }
1064 }
1065
1066 //=======================================================================
1067 // function : ConvertStringToTree()
1068 // purpose  : Returns the dependency tree
1069 //=======================================================================
1070 void GEOMUtils::ConvertStringToTree( const std::string& dependencyStr,
1071                                      TreeModel& tree )
1072 {
1073   std::size_t cursor = 0;
1074
1075   while( dependencyStr.find('-',cursor) != std::string::npos ) //find next selected object
1076   {
1077     std::size_t objectIndex = dependencyStr.find( '-', cursor );
1078     std::string objectEntry = dependencyStr.substr( cursor, objectIndex - cursor );
1079     cursor = objectIndex;
1080
1081     std::size_t upwardIndexBegin = dependencyStr.find("{",cursor) + 1;
1082     std::size_t upwardIndexFinish = dependencyStr.find("}",upwardIndexBegin);
1083     LevelsList upwardList = parseWard( dependencyStr, cursor );
1084
1085     LevelsList downwardList = parseWard( dependencyStr, cursor );
1086
1087     tree[objectEntry] = std::pair<LevelsList,LevelsList>( upwardList, downwardList );
1088   }
1089 }
1090
1091 bool GEOMUtils::CheckShape( TopoDS_Shape& shape,
1092                             bool checkGeometry )
1093 {
1094   BRepCheck_Analyzer analyzer( shape, checkGeometry );
1095   return analyzer.IsValid();
1096 }
1097
1098 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1099                                    TopAbs_ShapeEnum type,
1100                                    Standard_Real tolerance,
1101                                    bool checkGeometry )
1102 {
1103   ShapeFix_ShapeTolerance aSft;
1104   aSft.LimitTolerance( shape, tolerance, tolerance, type );
1105   Handle(ShapeFix_Shape) aSfs = new ShapeFix_Shape( shape );
1106   aSfs->Perform();
1107   shape = aSfs->Shape();
1108   return CheckShape( shape, checkGeometry );
1109 }
1110
1111 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1112                                    Standard_Real tolerance,
1113                                    bool checkGeometry )
1114 {
1115   return FixShapeTolerance( shape, TopAbs_SHAPE, tolerance, checkGeometry );
1116 }
1117
1118 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1119                                    bool checkGeometry )
1120 {
1121   return FixShapeTolerance( shape, Precision::Confusion(), checkGeometry );
1122 }
1123
1124 bool GEOMUtils::FixShapeCurves( TopoDS_Shape& shape )
1125 {
1126   Standard_Real aT, aTolE, aD, aDMax;
1127   TopExp_Explorer aExpF, aExpE;
1128   NCollection_DataMap<TopoDS_Edge, Standard_Real, TopTools_ShapeMapHasher> aDMETol;
1129   aExpF.Init(shape, TopAbs_FACE);
1130   for (; aExpF.More(); aExpF.Next()) {
1131     const TopoDS_Face& aF = *(TopoDS_Face*)&aExpF.Current();
1132     aExpE.Init(aF, TopAbs_EDGE);
1133     for (; aExpE.More(); aExpE.Next()) {
1134       const TopoDS_Edge& aE = *(TopoDS_Edge*)&aExpE.Current();
1135       try {
1136         if (!BOPTools_AlgoTools::ComputeTolerance(aF, aE, aDMax, aT)) {
1137           continue;
1138         }
1139       }
1140       catch(...) {
1141         continue;
1142       }
1143       aTolE = BRep_Tool::Tolerance(aE);
1144       if (aDMax < aTolE) continue;
1145       if (aDMETol.IsBound(aE)) {
1146         aD = aDMETol.Find(aE);
1147         if (aDMax > aD) {
1148           aDMETol.UnBind(aE);
1149           aDMETol.Bind(aE, aDMax);
1150         }
1151       }
1152       else {
1153         aDMETol.Bind(aE, aDMax);
1154       }
1155     }
1156   }
1157   NCollection_DataMap<TopoDS_Edge, Standard_Real, TopTools_ShapeMapHasher>::Iterator aDMETolIt(aDMETol);
1158 #ifdef USE_LIMIT_TOLERANCE
1159   ShapeFix_ShapeTolerance sat;
1160 #else
1161   BRep_Builder b;
1162 #endif
1163   for (; aDMETolIt.More(); aDMETolIt.Next()) {
1164 #ifdef USE_LIMIT_TOLERANCE
1165     sat.LimitTolerance(aDMETolIt.Key(), aDMETolIt.Value()*1.001);
1166 #else
1167     TopoDS_Iterator itv(aDMETolIt.Key());
1168     for (; itv.More(); itv.Next())
1169       b.UpdateVertex(TopoDS::Vertex(itv.Value()), aDMETolIt.Value()*1.001);
1170     b.UpdateEdge(aDMETolIt.Key(), aDMETolIt.Value()*1.001);
1171 #endif
1172   }
1173   return CheckShape( shape );
1174 }
1175
1176 bool GEOMUtils::Write( const TopoDS_Shape& shape, const char* fileName )
1177 {
1178   return BRepTools::Write( shape, fileName );
1179 }
1180
1181 TopoDS_Shape GEOMUtils::ReduceCompound( const TopoDS_Shape& shape )
1182 {
1183   TopoDS_Shape result = shape;
1184
1185   if ( shape.ShapeType() == TopAbs_COMPOUND ||
1186        shape.ShapeType() == TopAbs_COMPSOLID ) {
1187
1188     TopTools_ListOfShape l;
1189
1190     TopoDS_Iterator it ( shape );
1191     for ( ; it.More(); it.Next() )
1192       l.Append( it.Value() );
1193     if ( l.Extent() == 1 && l.First() != shape )
1194       result = ReduceCompound( l.First() );
1195   }
1196
1197   return result;
1198 }
1199
1200 void GEOMUtils::MeshShape( const TopoDS_Shape shape,
1201                            double deflection, bool theForced )
1202 {
1203   Standard_Real aDeflection = ( deflection <= 0 ) ? DefaultDeflection() : deflection;
1204
1205   // Is shape triangulated?
1206   Standard_Boolean alreadyMeshed = true;
1207   TopExp_Explorer ex;
1208   TopLoc_Location aLoc;
1209   for ( ex.Init( shape, TopAbs_FACE ); ex.More() && alreadyMeshed; ex.Next() ) {
1210     const TopoDS_Face& aFace = TopoDS::Face( ex.Current() );
1211     Handle(Poly_Triangulation) aPoly = BRep_Tool::Triangulation( aFace, aLoc );
1212     alreadyMeshed = !aPoly.IsNull();
1213   }
1214
1215   if ( !alreadyMeshed || theForced ) {
1216     // Compute bounding box
1217     Bnd_Box B;
1218     BRepBndLib::Add( shape, B );
1219     if ( B.IsVoid() )
1220       return; // NPAL15983 (Bug when displaying empty groups)
1221     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1222     B.Get( aXmin, aYmin, aZmin, aXmax, aYmax, aZmax );
1223
1224     // This magic line comes from Prs3d_ShadedShape.gxx in OCCT
1225     aDeflection = MAX3(aXmax-aXmin, aYmax-aYmin, aZmax-aZmin) * aDeflection * 4;
1226
1227     // Clean triangulation before compute incremental mesh
1228     BRepTools::Clean( shape );
1229
1230     // Compute triangulation
1231     BRepMesh_IncrementalMesh mesh( shape, aDeflection );
1232   }
1233 }
1234
1235 double GEOMUtils::DefaultDeflection()
1236 {
1237   return 0.001;
1238 }