src/GEOMUtils/GEOMUtils.cxx

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