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