src/GEOMUtils/GEOMUtils.cxx

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