1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include <Standard_Stream.hxx>
25 #include <GEOMUtils.hxx>
27 #include <Basics_OCCTVersion.hxx>
29 #include <utilities.h>
31 #include <Utils_ExceptHandlers.hxx>
34 #include <BRepMesh_IncrementalMesh.hxx>
36 #include <BRepExtrema_DistShapeShape.hxx>
38 #include <BRep_Builder.hxx>
39 #include <BRep_Tool.hxx>
40 #include <BRepBndLib.hxx>
41 #include <BRepGProp.hxx>
42 #include <BRepTools.hxx>
44 #include <BRepClass3d_SolidClassifier.hxx>
46 #include <BRepBuilderAPI_MakeFace.hxx>
48 #include <Bnd_Box.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>
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>
72 #include <GeomLProp_CLProps.hxx>
73 #include <GeomLProp_SLProps.hxx>
75 #include <GProp_GProps.hxx>
76 #include <GProp_PrincipalProps.hxx>
78 #include <TColStd_Array1OfReal.hxx>
83 #include <ShapeAnalysis.hxx>
84 #include <ShapeFix_Shape.hxx>
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
92 #define STD_SORT_ALGO 1
94 //=======================================================================
95 //function : GetPosition
97 //=======================================================================
98 gp_Ax3 GEOMUtils::GetPosition (const TopoDS_Shape& theShape)
102 if (theShape.IsNull())
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)
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);
128 TopAbs_ShapeEnum aShType = theShape.ShapeType();
130 if (aShType == TopAbs_VERTEX) {
131 aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
134 if (aShType == TopAbs_COMPOUND) {
135 aShType = GetTypeOfSimplePart(theShape);
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);
144 BRepGProp::VolumeProperties(theShape, aSystem);
146 aPnt = aSystem.CentreOfMass();
149 aResult.SetLocation(aPnt);
154 //=======================================================================
155 //function : GetVector
157 //=======================================================================
158 gp_Vec GEOMUtils::GetVector (const TopoDS_Shape& theShape,
159 Standard_Boolean doConsiderOrientation)
161 if (theShape.IsNull())
162 Standard_NullObject::Raise("Null shape is given for a vector");
164 if (theShape.ShapeType() != TopAbs_EDGE)
165 Standard_TypeMismatch::Raise("Invalid shape is given, must be a vector or an edge");
167 TopoDS_Edge anE = TopoDS::Edge(theShape);
169 TopoDS_Vertex V1, V2;
170 TopExp::Vertices(anE, V1, V2, doConsiderOrientation);
172 if (V1.IsNull() || V2.IsNull())
173 Standard_NullObject::Raise("Invalid edge is given, it must have two points");
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");
183 //=======================================================================
184 //function : ShapeToDouble
185 //purpose : used by CompareShapes::operator()
186 //=======================================================================
187 std::pair<double, double> ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
189 // Computing of CentreOfMass
193 if (S.ShapeType() == TopAbs_VERTEX) {
194 GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
195 Len = (double)S.Orientation();
199 // BEGIN: fix for Mantis issue 0020842
201 BRepGProp::LinearProperties(S, GPr);
204 if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
205 BRepGProp::LinearProperties(S, GPr);
207 else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
208 BRepGProp::SurfaceProperties(S, GPr);
211 BRepGProp::VolumeProperties(S, GPr);
214 // END: fix for Mantis issue 0020842
215 GPoint = GPr.CentreOfMass();
219 double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
220 return std::make_pair(dMidXYZ, Len);
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)
230 if (!myMap.IsBound(theShape1)) {
231 myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
234 if (!myMap.IsBound(theShape2)) {
235 myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
238 std::pair<double, double> val1 = myMap.Find(theShape1);
239 std::pair<double, double> val2 = myMap.Find(theShape2);
241 double tol = Precision::Confusion();
242 bool exchange = Standard_False;
244 double dMidXYZ = val1.first - val2.first;
245 if (dMidXYZ >= tol) {
246 exchange = Standard_True;
248 else if (Abs(dMidXYZ) < tol) {
249 double dLength = val1.second - val2.second;
250 if (dLength >= tol) {
251 exchange = Standard_True;
253 else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
255 // equal values possible on shapes such as two halves of a sphere and
256 // a membrane inside the sphere
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;
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;
279 //return val1 < val2;
283 //=======================================================================
284 //function : SortShapes
286 //=======================================================================
287 void GEOMUtils::SortShapes (TopTools_ListOfShape& SL,
288 const Standard_Boolean isOldSorting)
291 std::vector<TopoDS_Shape> aShapesVec;
292 aShapesVec.reserve(SL.Extent());
294 TopTools_ListIteratorOfListOfShape it (SL);
295 for (; it.More(); it.Next()) {
296 aShapesVec.push_back(it.Value());
300 CompareShapes shComp (isOldSorting);
301 std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
302 //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
304 std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
305 for (; anIter != aShapesVec.end(); ++anIter) {
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;
316 // Computing of CentreOfMass
317 Standard_Integer Index;
320 TopTools_ListIteratorOfListOfShape it(SL);
321 for (Index=1; it.More(); Index++)
323 TopoDS_Shape S = it.Value();
324 SL.Remove( it ); // == it.Next()
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());
332 // BEGIN: fix for Mantis issue 0020842
334 BRepGProp::LinearProperties (S, GPr);
337 if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
338 BRepGProp::LinearProperties (S, GPr);
340 else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
341 BRepGProp::SurfaceProperties(S, GPr);
344 BRepGProp::VolumeProperties(S, GPr);
347 // END: fix for Mantis issue 0020842
348 GPoint = GPr.CentreOfMass();
349 Length.SetValue(Index, GPr.Mass());
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;
356 Standard_Integer aTemp;
357 Standard_Boolean exchange, Sort = Standard_True;
358 Standard_Real tol = Precision::Confusion();
361 Sort = Standard_False;
362 for (Index=1; Index < MaxShapes; Index++)
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;
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;
377 else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
378 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
380 // equal values possible on shapes such as two halves of a sphere and
381 // a membrane inside the sphere
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;
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;
401 //cout << "box: " << val1<<" > "<<val2 << endl;
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;
416 for (Index=1; Index <= MaxShapes; Index++)
417 SL.Append( aShapes( OrderInd(Index) ));
421 //=======================================================================
422 //function : CompsolidToCompound
424 //=======================================================================
425 TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
427 if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
431 TopoDS_Compound aCompound;
433 B.MakeCompound(aCompound);
435 TopTools_MapOfShape mapShape;
436 TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
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);
448 //=======================================================================
449 //function : AddSimpleShapes
451 //=======================================================================
452 void GEOMUtils::AddSimpleShapes (const TopoDS_Shape& theShape, TopTools_ListOfShape& theList)
454 if (theShape.ShapeType() != TopAbs_COMPOUND &&
455 theShape.ShapeType() != TopAbs_COMPSOLID) {
456 theList.Append(theShape);
460 TopTools_MapOfShape mapShape;
461 TopoDS_Iterator It (theShape, Standard_True, Standard_True);
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);
470 theList.Append(aShape_i);
476 //=======================================================================
477 //function : CheckTriangulation
479 //=======================================================================
480 bool GEOMUtils::CheckTriangulation (const TopoDS_Shape& aShape)
482 bool isTriangulation = true;
484 TopExp_Explorer exp (aShape, TopAbs_FACE);
487 TopLoc_Location aTopLoc;
488 Handle(Poly_Triangulation) aTRF;
489 aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
491 isTriangulation = false;
494 else // no faces, try edges
496 TopExp_Explorer expe (aShape, TopAbs_EDGE);
500 TopLoc_Location aLoc;
501 Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
503 isTriangulation = false;
507 if (!isTriangulation) {
508 // calculate deflection
509 Standard_Real aDeviationCoefficient = 0.001;
512 BRepBndLib::Add(aShape, B);
513 Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
514 B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
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;
520 BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
526 //=======================================================================
527 //function : GetTypeOfSimplePart
529 //=======================================================================
530 TopAbs_ShapeEnum GEOMUtils::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
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);
541 return GetTypeOfSimplePart(It.Value());
547 //=======================================================================
548 //function : GetEdgeNearPoint
550 //=======================================================================
551 TopoDS_Shape GEOMUtils::GetEdgeNearPoint (const TopoDS_Shape& theShape,
552 const TopoDS_Vertex& thePoint)
554 TopoDS_Shape aResult;
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())) {
567 Standard_NullObject::Raise("Given shape contains no edges");
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;
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");
583 aDistances(ind) = aDistTool.Value();
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) {
595 else if (aDistances(ind) < nearest) {
596 nearest = aDistances(ind);
597 aResult = anEdges(ind);
604 Standard_ConstructionError::Raise("Multiple edges near the given point are found");
606 else if (nbFound == 0) {
607 Standard_ConstructionError::Raise("There are no edges near the given point");
615 //=======================================================================
616 //function : PreciseBoundingBox
618 //=======================================================================
619 Standard_Boolean GEOMUtils::PreciseBoundingBox
620 (const TopoDS_Shape &theShape, Bnd_Box &theBox)
622 Standard_Real aBound[6];
624 theBox.Get(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
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] =
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
642 const gp_Dir aDir[3] = { gp::DX(), gp::DY(), gp::DZ() };
643 const Standard_Real aPlnSize[3] =
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
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]);
657 if (!aMkFace.IsDone()) {
658 return Standard_False;
661 TopoDS_Shape aFace = aMkFace.Shape();
663 // Get minimal distance between planar face and shape.
664 Standard_Real aMinDist =
665 GEOMUtils::GetMinDistance(aFace, theShape, aPMin[0], aPMin[1]);
668 return Standard_False;
671 aBound[i] = aPMin[1].Coord(iHalf + 1);
674 // Update Bounding box with the new values.
676 theBox.Update(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
678 return Standard_True;
681 //=======================================================================
682 //function : GetMinDistanceSingular
684 //=======================================================================
685 double GEOMUtils::GetMinDistanceSingular(const TopoDS_Shape& aSh1,
686 const TopoDS_Shape& aSh2,
687 gp_Pnt& Ptmp1, gp_Pnt& Ptmp2)
689 bool IsChange1 = false;
690 double AddDist1 = 0.0;
691 TopExp_Explorer anExp;
692 TopoDS_Shape tmpSh1, tmpSh2;
694 for ( anExp.Init( aSh1, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
696 tmpSh1 = anExp.Current();
699 TopoDS_Shape sh = aSh1;
700 while(sh.ShapeType()==TopAbs_COMPOUND) {
701 TopoDS_Iterator it(sh);
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 ) {
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);
720 TopoDS_Compound Comp;
721 B.MakeCompound(Comp);
722 B.MakeFace(F1,TrS1,1.e-7);
724 B.MakeFace(F2,TrS2,1.e-7);
726 Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
728 sfs->SetPrecision(1.e-6);
729 sfs->SetMaxTolerance(1.0);
731 tmpSh1 = sfs->Shape();
735 if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
736 Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
737 gp_Pnt PC = SS->Location();
740 B.MakeVertex(V,PC,1.e-7);
742 AddDist1 = SS->Radius();
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());
751 B.MakeEdge(E,C,1.e-7);
753 AddDist1 = TS->MinorRadius();
763 bool IsChange2 = false;
764 double AddDist2 = 0.0;
766 for ( anExp.Init( aSh2, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
768 tmpSh2 = anExp.Current();
771 TopoDS_Shape sh = aSh2;
772 while(sh.ShapeType()==TopAbs_COMPOUND) {
773 TopoDS_Iterator it(sh);
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 ) {
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);
790 TopoDS_Compound Comp;
791 B.MakeCompound(Comp);
792 B.MakeFace(F1,TrS1,1.e-7);
794 B.MakeFace(F2,TrS2,1.e-7);
796 Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
798 sfs->SetPrecision(1.e-6);
799 sfs->SetMaxTolerance(1.0);
801 tmpSh2 = sfs->Shape();
805 if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
806 Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
807 gp_Pnt PC = SS->Location();
810 B.MakeVertex(V,PC,1.e-7);
812 AddDist2 = SS->Radius();
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());
821 B.MakeEdge(E,C,1.e-7);
823 AddDist2 = TS->MinorRadius();
834 if( !IsChange1 && !IsChange2 )
837 BRepExtrema_DistShapeShape dst(tmpSh1,tmpSh2);
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) {
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);
868 Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
869 PMin1.Y() + aDir.Y()*AddDist1,
870 PMin1.Z() + aDir.Z()*AddDist1 );
874 Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
875 PMin2.Y() - aDir.Y()*AddDist2,
876 PMin2.Z() - aDir.Z()*AddDist2 );
881 double res = MinDist - AddDist1 - AddDist2;
882 if(res<0.) res = 0.0;
888 //=======================================================================
889 //function : GetMinDistance
891 //=======================================================================
892 Standard_Real GEOMUtils::GetMinDistance
893 (const TopoDS_Shape& theShape1,
894 const TopoDS_Shape& theShape2,
895 gp_Pnt& thePnt1, gp_Pnt& thePnt2)
897 Standard_Real aResult = 1.e9;
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
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;
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 ) {
927 // additional workaround for bugs 19899, 19908 and 19910 from Mantis
928 double dist = GEOMUtils::GetMinDistanceSingular
929 (theShape1, theShape2, thePnt1, thePnt2);
935 BRepExtrema_DistShapeShape dst (theShape1, theShape2);
939 for (int i = 1; i <= dst.NbSolution(); i++) {
940 P1 = dst.PointOnShape1(i);
941 P2 = dst.PointOnShape2(i);
943 Standard_Real Dist = P1.Distance(P2);
944 if (aResult > Dist) {