src/GEOMUtils/GEOMUtils.cxx

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