src/GEOMUtils/GEOMUtils.cxx

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