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