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