Salome HOME
21854: Add persistent dimensions
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IMeasureOperations.cxx
index 98abba29cb334a09d5a9ccd8300692767b77a087..aaa0b4ed86df74b97ec38973f0acc4b00ae37a9c 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 #include <GEOMImpl_IMeasure.hxx>
 #include <GEOMImpl_IShapesOperations.hxx>
 
+#include <GEOMUtils.hxx>
+
 #include <GEOMAlgo_ShapeInfo.hxx>
 #include <GEOMAlgo_ShapeInfoFiller.hxx>
 
 #include <GEOM_Function.hxx>
 #include <GEOM_PythonDump.hxx>
 
-#include <NMTTools_CheckerSI.hxx>
-
-#include <NMTDS_Tools.hxx>
-#include <NMTDS_InterfPool.hxx>
-#include <NMTDS_PInterfPool.hxx>
-//#include <NMTDS_PassKeyBoolean.hxx>
-#include <NMTDS_PairBoolean.hxx>
-#include <NMTDS_ShapesDataStructure.hxx>
-//#include <NMTDS_ListIteratorOfListOfPassKeyBoolean.hxx>
-#include <NMTDS_ListIteratorOfListOfPairBoolean.hxx>
-
 #include <Basics_OCCTVersion.hxx>
 
 #include <utilities.h>
 #include <TFunction_Logbook.hxx>
 #include <TDF_Tool.hxx>
 
+#include <BRep_Builder.hxx>
+#include <BRep_TFace.hxx>
 #include <BRep_Tool.hxx>
 #include <BRepAdaptor_Surface.hxx>
+#include <BRepBuilderAPI_Copy.hxx>
 #include <BRepBndLib.hxx>
 #include <BRepCheck.hxx>
 #include <BRepCheck_ListIteratorOfListOfStatus.hxx>
 #include <BRepCheck_Result.hxx>
 #include <BRepCheck_Shell.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
+#include <BRepClass_FaceClassifier.hxx>
 #include <BRepExtrema_DistShapeShape.hxx>
 #include <BRepGProp.hxx>
 #include <BRepTools.hxx>
 
 #include <Bnd_Box.hxx>
 
-#include <GProp_GProps.hxx>
-#include <GProp_PrincipalProps.hxx>
-
 #include <TopAbs.hxx>
 #include <TopExp.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Shape.hxx>
 #include <TopoDS_Vertex.hxx>
+#include <TopoDS_Compound.hxx>
 #include <TopoDS_Iterator.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 
+#include <ShapeAnalysis.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+
+#include <GeomAPI_IntSS.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
 #include <GeomAbs_SurfaceType.hxx>
-#include <Geom_Surface.hxx>
-#include <Geom_Plane.hxx>
-#include <Geom_SphericalSurface.hxx>
-#include <Geom_CylindricalSurface.hxx>
-#include <Geom_ToroidalSurface.hxx>
-#include <Geom_ConicalSurface.hxx>
-#include <Geom_SurfaceOfLinearExtrusion.hxx>
-#include <Geom_SurfaceOfRevolution.hxx>
-#include <Geom_BezierSurface.hxx>
-#include <Geom_BSplineSurface.hxx>
-#include <Geom_RectangularTrimmedSurface.hxx>
-#include <Geom_OffsetSurface.hxx>
+
 #include <Geom_Line.hxx>
+#include <Geom_Surface.hxx>
+
+#include <GeomLProp_CLProps.hxx>
+#include <GeomLProp_SLProps.hxx>
+
+#include <GProp_GProps.hxx>
+#include <GProp_PrincipalProps.hxx>
 
 #include <gp_Pln.hxx>
 #include <gp_Lin.hxx>
 
-#include <GeomAPI_ProjectPointOnCurve.hxx>
-#include <GeomLProp_CLProps.hxx>
-#include <GeomLProp_SLProps.hxx>
-#include <ShapeAnalysis.hxx>
-#include <ShapeAnalysis_Surface.hxx>
+#include <BOPCol_ListOfShape.hxx>
+#include <BOPDS_DS.hxx>
+#include <BOPDS_MapOfPassKey.hxx>
+#include <BOPDS_PassKey.hxx>
+#include <GEOMAlgo_AlgoTools.hxx>
+#include <BOPAlgo_CheckerSI.hxx>
 
 #include <Standard_Failure.hxx>
 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 
-#include <BRepClass3d_SolidClassifier.hxx>
-#include <BRep_Builder.hxx>
-#include <GeomAPI_IntSS.hxx>
-#include <Geom_Circle.hxx>
-#include <Geom_SphericalSurface.hxx>
-#include <Geom_ToroidalSurface.hxx>
-#include <ShapeFix_Shape.hxx>
-#include <TopoDS_Compound.hxx>
-
-
 //=============================================================================
 /*!
  *  Constructor
@@ -767,69 +755,6 @@ GEOMImpl_IMeasureOperations::ShapeKind GEOMImpl_IMeasureOperations::KindOfShape
   return aKind;
 }
 
-//=============================================================================
-/*! Get LCS, corresponding to the given shape.
- *  Origin of the LCS is situated at the shape's center of mass.
- *  Axes of the LCS are obtained from shape's location or,
- *  if the shape is a planar face, from position of its plane.
- */
-//=============================================================================
-gp_Ax3 GEOMImpl_IMeasureOperations::GetPosition (const TopoDS_Shape& theShape)
-{
-  gp_Ax3 aResult;
-
-  if (theShape.IsNull())
-    return aResult;
-
-  // Axes
-  aResult.Transform(theShape.Location().Transformation());
-  if (theShape.ShapeType() == TopAbs_FACE) {
-    Handle(Geom_Surface) aGS = BRep_Tool::Surface(TopoDS::Face(theShape));
-    if (!aGS.IsNull() && aGS->IsKind(STANDARD_TYPE(Geom_Plane))) {
-      Handle(Geom_Plane) aGPlane = Handle(Geom_Plane)::DownCast(aGS);
-      gp_Pln aPln = aGPlane->Pln();
-      aResult = aPln.Position();
-      // In case of reverse orinetation of the face invert the plane normal
-      // (the face's normal does not mathc the plane's normal in this case)
-      if(theShape.Orientation() == TopAbs_REVERSED)
-      {
-        gp_Dir Vx =  aResult.XDirection();
-        gp_Dir N  =  aResult.Direction().Mirrored(Vx); 
-        gp_Pnt P  =  aResult.Location();
-        aResult = gp_Ax3(P, N, Vx);
-      }       
-    }
-  }
-
-  // Origin
-  gp_Pnt aPnt;
-
-  TopAbs_ShapeEnum aShType = theShape.ShapeType();
-
-  if (aShType == TopAbs_VERTEX) {
-    aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
-  }
-  else {
-    if (aShType == TopAbs_COMPOUND) {
-      aShType = GEOMImpl_IShapesOperations::GetTypeOfSimplePart(theShape);
-    }
-
-    GProp_GProps aSystem;
-    if (aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)
-      BRepGProp::LinearProperties(theShape, aSystem);
-    else if (aShType == TopAbs_FACE || aShType == TopAbs_SHELL)
-      BRepGProp::SurfaceProperties(theShape, aSystem);
-    else
-      BRepGProp::VolumeProperties(theShape, aSystem);
-
-    aPnt = aSystem.CentreOfMass();
-  }
-
-  aResult.SetLocation(aPnt);
-
-  return aResult;
-}
-
 //=============================================================================
 /*!
  *  GetPosition
@@ -863,7 +788,7 @@ void GEOMImpl_IMeasureOperations::GetPosition
     OCC_CATCH_SIGNALS;
 #endif
 
-    gp_Ax3 anAx3 = GetPosition(aShape);
+    gp_Ax3 anAx3 = GEOMUtils::GetPosition(aShape);
 
     gp_Pnt anOri = anAx3.Location();
     gp_Dir aDirZ = anAx3.Direction();
@@ -1182,6 +1107,7 @@ void GEOMImpl_IMeasureOperations::GetInertia
 //=============================================================================
 void GEOMImpl_IMeasureOperations::GetBoundingBox
                                      (Handle(GEOM_Object) theShape,
+                                      const Standard_Boolean precise,
                                       Standard_Real& Xmin, Standard_Real& Xmax,
                                       Standard_Real& Ymin, Standard_Real& Ymax,
                                       Standard_Real& Zmin, Standard_Real& Zmax)
@@ -1206,7 +1132,26 @@ void GEOMImpl_IMeasureOperations::GetBoundingBox
 #if OCC_VERSION_LARGE > 0x06010000
     OCC_CATCH_SIGNALS;
 #endif
+    BRepBuilderAPI_Copy aCopyTool (aShape);
+    if (!aCopyTool.IsDone()) {
+      SetErrorCode("GetBoundingBox Error: Bad shape detected");
+      return;
+    }
+
+    aShape = aCopyTool.Shape();
+
+    // remove triangulation to obtain more exact boundaries
+    BRepTools::Clean(aShape);
+
     BRepBndLib::Add(aShape, B);
+
+    if (precise) {
+      if (!GEOMUtils::PreciseBoundingBox(aShape, B)) {
+        SetErrorCode("GetBoundingBox Error: Bounding box cannot be precised");
+        return;
+      }
+    }
+
     B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
   }
   catch (Standard_Failure) {
@@ -1218,6 +1163,69 @@ void GEOMImpl_IMeasureOperations::GetBoundingBox
   SetErrorCode(OK);
 }
 
+//=============================================================================
+/*!
+ *  GetBoundingBox
+ */
+//=============================================================================
+Handle(GEOM_Object) GEOMImpl_IMeasureOperations::GetBoundingBox
+                                                (Handle(GEOM_Object) theShape,
+                                                 const Standard_Boolean precise)
+{
+  SetErrorCode(KO);
+
+  if (theShape.IsNull()) return NULL;
+
+  //Add a new BoundingBox object
+  Handle(GEOM_Object) aBnd = GetEngine()->AddObject(GetDocID(), GEOM_BOX);
+
+  //Add a new BoundingBox function
+  const int aType = (precise ? BND_BOX_MEASURE_PRECISE : BND_BOX_MEASURE);
+  Handle(GEOM_Function) aFunction =
+    aBnd->AddFunction(GEOMImpl_MeasureDriver::GetID(), aType);
+  if (aFunction.IsNull()) return NULL;
+
+  //Check if the function is set correctly
+  if (aFunction->GetDriverGUID() != GEOMImpl_MeasureDriver::GetID()) return NULL;
+
+  GEOMImpl_IMeasure aCI (aFunction);
+
+  Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
+  if (aRefShape.IsNull()) return NULL;
+
+  aCI.SetBase(aRefShape);
+
+  //Compute the BoundingBox value
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
+    if (!GetSolver()->ComputeFunction(aFunction)) {
+      SetErrorCode("Measure driver failed to compute a bounding box");
+      return NULL;
+    }
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return NULL;
+  }
+
+  //Make a Python command
+  GEOM::TPythonDump aPd(aFunction);
+  
+  aPd << aBnd << " = geompy.MakeBoundingBox(" << theShape;
+
+  if (precise) {
+    aPd << ", True";
+  }
+
+  aPd << ")";
+
+  SetErrorCode(OK);
+  return aBnd;
+}
+
 //=============================================================================
 /*!
  *  GetTolerance
@@ -1361,8 +1369,9 @@ bool GEOMImpl_IMeasureOperations::CheckSelfIntersections
   BRep_Builder aBB;
   TopoDS_Compound aCS;
   TopoDS_Shape aScopy;
-  NMTDS_Tools::CopyShape(aShape, aScopy);
-
+  //
+  GEOMAlgo_AlgoTools::CopyShape(aShape, aScopy);
+  //
   // Map sub-shapes and their indices
   TopTools_IndexedMapOfShape anIndices;
   TopExp::MapShapes(aScopy, anIndices);
@@ -1370,40 +1379,38 @@ bool GEOMImpl_IMeasureOperations::CheckSelfIntersections
   aBB.MakeCompound(aCS);
   aBB.Add(aCS, aScopy);
 
-  NMTTools_CheckerSI aCSI; // checker of self-interferences
-  aCSI.SetCompositeShape(aCS);
+  BOPCol_ListOfShape aLCS;
+  aLCS.Append(aScopy);
+  //
+  BOPAlgo_CheckerSI aCSI; // checker of self-interferences
+  aCSI.SetArguments(aLCS);
 
   // 1. Launch the checker
   aCSI.Perform();
-  Standard_Integer iErr = aCSI.StopStatus();
+  Standard_Integer iErr = aCSI.ErrorStatus();
   if (iErr) {
     return false; // Error
   }
 
   isGood = true;
-
+  //
+  Standard_Integer aNbS, n1, n2;
+  BOPDS_MapIteratorMapOfPassKey aItMPK;
+  //
   // 2. Take the shapes from DS
-  const NMTDS_ShapesDataStructure& aDS = *(aCSI.DS());
-  Standard_Integer aNbS = aDS.NumberOfShapesOfTheObject();
-
-  // 3. Get the pairs of interfered shapes 
-  NMTDS_PInterfPool pIP = aCSI.IP();
-  //const NMTDS_ListOfPassKeyBoolean& aLPKB = pIP->Get();
-  const NMTDS_ListOfPairBoolean& aLPKB = pIP->Get();
-
-  Standard_Integer n1, n2;
-  //NMTDS_ListIteratorOfListOfPassKeyBoolean aIt;
-  NMTDS_ListIteratorOfListOfPairBoolean aIt;
-
-  aIt.Initialize(aLPKB);
-  for (; aIt.More(); aIt.Next()) {
-    //const NMTDS_PassKeyBoolean& aPKB = aIt.Value();
-    const NMTDS_PairBoolean& aPKB = aIt.Value();
-    aPKB.Ids(n1, n2);
-
-    if (n1 > aNbS || n2 > aNbS)
+  const BOPDS_DS& aDS = aCSI.DS();
+  aNbS=aDS.NbShapes();
+  //
+  // 3. Get the pairs of interfered shapes
+  const BOPDS_MapOfPassKey& aMPK=aDS.Interferences();
+  aItMPK.Initialize(aMPK);
+  for (; aItMPK.More(); aItMPK.Next()) {
+    const BOPDS_PassKey& aPK=aItMPK.Value();
+    aPK.Ids(n1, n2);
+    //
+    if (n1 > aNbS || n2 > aNbS){
       return false; // Error
-
+    }
     const TopoDS_Shape& aS1 = aDS.Shape(n1);
     const TopoDS_Shape& aS2 = aDS.Shape(n2);
 
@@ -1541,339 +1548,121 @@ TCollection_AsciiString GEOMImpl_IMeasureOperations::WhatIs (Handle(GEOM_Object)
   return Astr;
 }
 
-
-//=======================================================================
-//function : CheckSingularCase
-//purpose  : auxilary for GetMinDistance()
-//           workaround for bugs 19899, 19908 and 19910 from Mantis
-//=======================================================================
-static double CheckSingularCase(const TopoDS_Shape& aSh1,
-                                const TopoDS_Shape& aSh2,
-                                gp_Pnt& Ptmp1, gp_Pnt& Ptmp2)
-{
-  bool IsChange1 = false;
-  double AddDist1 = 0.0;
-  TopExp_Explorer anExp;
-  TopoDS_Shape tmpSh1, tmpSh2;
-  int nbf = 0;
-  for ( anExp.Init( aSh1, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-    nbf++;
-    tmpSh1 = anExp.Current();
-  }
-  if(nbf==1) {
-    TopoDS_Shape sh = aSh1;
-    while(sh.ShapeType()==TopAbs_COMPOUND) {
-      TopoDS_Iterator it(sh);
-      sh = it.Value();
-    }
-    Handle(Geom_Surface) S = BRep_Tool::Surface(TopoDS::Face(tmpSh1));
-    if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
-        S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
-      if( sh.ShapeType()==TopAbs_SHELL || sh.ShapeType()==TopAbs_FACE ) {
-        // non solid case
-        double U1,U2,V1,V2;
-        // changes for 0020677: EDF 1219 GEOM: MinDistance gives 0 instead of 20.88
-        //S->Bounds(U1,U2,V1,V2); changed by
-        ShapeAnalysis::GetFaceUVBounds(TopoDS::Face(tmpSh1),U1,U2,V1,V2);
-        // end of changes for 020677 (dmv)
-        Handle(Geom_RectangularTrimmedSurface) TrS1 =
-          new Geom_RectangularTrimmedSurface(S,U1,(U1+U2)/2.,V1,V2);
-        Handle(Geom_RectangularTrimmedSurface) TrS2 =
-          new Geom_RectangularTrimmedSurface(S,(U1+U2)/2.,U2,V1,V2);
-        BRep_Builder B;
-        TopoDS_Face F1,F2;
-        TopoDS_Compound Comp;
-        B.MakeCompound(Comp);
-        B.MakeFace(F1,TrS1,1.e-7);
-        B.Add(Comp,F1);
-        B.MakeFace(F2,TrS2,1.e-7);
-        B.Add(Comp,F2);
-        Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
-        sfs->Init(Comp);
-        sfs->SetPrecision(1.e-6);
-        sfs->SetMaxTolerance(1.0);
-        sfs->Perform();
-        tmpSh1 = sfs->Shape();
-        IsChange1 = true;
-      }
-      else {
-        if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
-          Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
-          gp_Pnt PC = SS->Location();
-          BRep_Builder B;
-          TopoDS_Vertex V;
-          B.MakeVertex(V,PC,1.e-7);
-          tmpSh1 = V;
-          AddDist1 = SS->Radius();
-          IsChange1 = true;
-        }
-        else {
-          Handle(Geom_ToroidalSurface) TS = Handle(Geom_ToroidalSurface)::DownCast(S);
-          gp_Ax3 ax3 = TS->Position();
-          Handle(Geom_Circle) C = new Geom_Circle(ax3.Ax2(),TS->MajorRadius());
-          BRep_Builder B;
-          TopoDS_Edge E;
-          B.MakeEdge(E,C,1.e-7);
-          tmpSh1 = E;
-          AddDist1 = TS->MinorRadius();
-          IsChange1 = true;
-        }
-      }
-    }
-    else
-      tmpSh1 = aSh1;
-  }
-  else
-    tmpSh1 = aSh1;
-  bool IsChange2 = false;
-  double AddDist2 = 0.0;
-  nbf = 0;
-  for ( anExp.Init( aSh2, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-    nbf++;
-    tmpSh2 = anExp.Current();
-  }
-  if(nbf==1) {
-    TopoDS_Shape sh = aSh2;
-    while(sh.ShapeType()==TopAbs_COMPOUND) {
-      TopoDS_Iterator it(sh);
-      sh = it.Value();
-    }
-    Handle(Geom_Surface) S = BRep_Tool::Surface(TopoDS::Face(tmpSh2));
-    if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
-        S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
-      if( sh.ShapeType()==TopAbs_SHELL || sh.ShapeType()==TopAbs_FACE ) {
-        // non solid case
-        double U1,U2,V1,V2;
-        //S->Bounds(U1,U2,V1,V2);
-        ShapeAnalysis::GetFaceUVBounds(TopoDS::Face(tmpSh2),U1,U2,V1,V2);
-        Handle(Geom_RectangularTrimmedSurface) TrS1 =
-          new Geom_RectangularTrimmedSurface(S,U1,(U1+U2)/2.,V1,V2);
-        Handle(Geom_RectangularTrimmedSurface) TrS2 =
-          new Geom_RectangularTrimmedSurface(S,(U1+U2)/2.,U2,V1,V2);
-        BRep_Builder B;
-        TopoDS_Face F1,F2;
-        TopoDS_Compound Comp;
-        B.MakeCompound(Comp);
-        B.MakeFace(F1,TrS1,1.e-7);
-        B.Add(Comp,F1);
-        B.MakeFace(F2,TrS2,1.e-7);
-        B.Add(Comp,F2);
-        Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
-        sfs->Init(Comp);
-        sfs->SetPrecision(1.e-6);
-        sfs->SetMaxTolerance(1.0);
-        sfs->Perform();
-        tmpSh2 = sfs->Shape();
-        IsChange2 = true;
-      }
-      else {
-        if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
-          Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
-          gp_Pnt PC = SS->Location();
-          BRep_Builder B;
-          TopoDS_Vertex V;
-          B.MakeVertex(V,PC,1.e-7);
-          tmpSh2 = V;
-          AddDist2 = SS->Radius();
-          IsChange2 = true;
-        }
-        else if( S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
-          Handle(Geom_ToroidalSurface) TS = Handle(Geom_ToroidalSurface)::DownCast(S);
-          gp_Ax3 ax3 = TS->Position();
-          Handle(Geom_Circle) C = new Geom_Circle(ax3.Ax2(),TS->MajorRadius());
-          BRep_Builder B;
-          TopoDS_Edge E;
-          B.MakeEdge(E,C,1.e-7);
-          tmpSh2 = E;
-          AddDist2 = TS->MinorRadius();
-          IsChange2 = true;
-        }
-      }
-    }
-    else
-      tmpSh2 = aSh2;
-  }
-  else
-    tmpSh2 = aSh2;
-
-  if( !IsChange1 && !IsChange2 )
-    return -2.0;
-
-  BRepExtrema_DistShapeShape dst(tmpSh1,tmpSh2);
-  if (dst.IsDone()) {
-    double MinDist = 1.e9;
-    gp_Pnt PMin1, PMin2, P1, P2;
-    for (int i = 1; i <= dst.NbSolution(); i++) {
-      P1 = dst.PointOnShape1(i);
-      P2 = dst.PointOnShape2(i);
-      Standard_Real Dist = P1.Distance(P2);
-      if (MinDist > Dist) {
-        MinDist = Dist;
-        PMin1 = P1;
-        PMin2 = P2;
-      }
-    }
-    if(MinDist<1.e-7) {
-      Ptmp1 = PMin1;
-      Ptmp2 = PMin2;
-    }
-    else {
-      gp_Dir aDir(gp_Vec(PMin1,PMin2));
-      if( MinDist > (AddDist1+AddDist2) ) {
-        Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
-                        PMin1.Y() + aDir.Y()*AddDist1,
-                        PMin1.Z() + aDir.Z()*AddDist1 );
-        Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
-                        PMin2.Y() - aDir.Y()*AddDist2,
-                        PMin2.Z() - aDir.Z()*AddDist2 );
-        return (MinDist - AddDist1 - AddDist2);
-      }
-      else {
-        if( AddDist1 > 0 ) {
-          Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
-                          PMin1.Y() + aDir.Y()*AddDist1,
-                          PMin1.Z() + aDir.Z()*AddDist1 );
-          Ptmp2 = Ptmp1;
-        }
-        else {
-          Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
-                          PMin2.Y() - aDir.Y()*AddDist2,
-                          PMin2.Z() - aDir.Z()*AddDist2 );
-          Ptmp1 = Ptmp2;
-        }
-      }
-    }
-    double res = MinDist - AddDist1 - AddDist2;
-    if(res<0.) res = 0.0;
-    return res;
-  }
-  return -2.0;
-}
-/* old variant
-static bool CheckSingularCase(const TopoDS_Shape& aSh1,
-                              const TopoDS_Shape& aSh2,
-                              gp_Pnt& Ptmp)
-{
-  TopExp_Explorer anExp;
-  TopoDS_Shape tmpSh1, tmpSh2;
-  int nbf = 0;
-  for ( anExp.Init( aSh1, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-    nbf++;
-    tmpSh1 = anExp.Current();
-  }
-  if(nbf==1) {
-    Handle(Geom_Surface) S1 = BRep_Tool::Surface(TopoDS::Face(tmpSh1));
-    if( S1->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
-        S1->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
-      nbf = 0;
-      for ( anExp.Init( aSh2, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-        nbf++;
-        tmpSh2 = anExp.Current();
-        Handle(Geom_Surface) S2 = BRep_Tool::Surface(TopoDS::Face(tmpSh2));
-        GeomAPI_IntSS ISS(S1,S2,1.e-7);
-        if(ISS.IsDone()) {
-          for(int i=1; i<=ISS.NbLines(); i++) {
-            Handle(Geom_Curve) C3d = ISS.Line(i);
-            BRep_Builder B;
-            TopoDS_Edge E;
-            B.MakeEdge(E,C3d,1.e-7);
-            BRepExtrema_DistShapeShape dst(tmpSh2,E);
-            if (dst.IsDone()) {
-              gp_Pnt PMin1, PMin2, P1, P2;
-              double MinDist = 1.e9;
-              for (int i = 1; i <= dst.NbSolution(); i++) {
-                P1 = dst.PointOnShape1(i);
-                P2 = dst.PointOnShape2(i);
-                Standard_Real Dist = P1.Distance(P2);
-                if (MinDist > Dist) {
-                  MinDist = Dist;
-                  Ptmp = P1;
-                }
-              }
-              if(MinDist<1.e-7)
-                return true;
-            }
-          }
-        }
-      }
-    }
-  }
-  nbf = 0;
-  for ( anExp.Init( aSh2, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-    nbf++;
-    tmpSh1 = anExp.Current();
-  }
-  if(nbf==1) {
-    Handle(Geom_Surface) S1 = BRep_Tool::Surface(TopoDS::Face(tmpSh1));
-    if( S1->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
-        S1->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
-      nbf = 0;
-      for ( anExp.Init( aSh1, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
-        nbf++;
-        tmpSh2 = anExp.Current();
-        Handle(Geom_Surface) S2 = BRep_Tool::Surface(TopoDS::Face(tmpSh2));
-        GeomAPI_IntSS ISS(S1,S2,1.e-7);
-        if(ISS.IsDone()) {
-          for(int i=1; i<=ISS.NbLines(); i++) {
-            Handle(Geom_Curve) C3d = ISS.Line(i);
-            BRep_Builder B;
-            TopoDS_Edge E;
-            B.MakeEdge(E,C3d,1.e-7);
-            BRepExtrema_DistShapeShape dst(tmpSh2,E);
-            if (dst.IsDone()) {
-              gp_Pnt P1,P2;
-              double MinDist = 1.e9;
-              for (int i = 1; i <= dst.NbSolution(); i++) {
-                P1 = dst.PointOnShape1(i);
-                P2 = dst.PointOnShape2(i);
-                Standard_Real Dist = P1.Distance(P2);
-                if (MinDist > Dist) {
-                  MinDist = Dist;
-                  Ptmp = P1;
-                }
-              }
-              if(MinDist<1.e-7)
-                return true;
-            }
-          }
-        }
-      }
-    }
-  }
-  return false;
-}
-*/
-
-
 //=============================================================================
 /*!
  *  AreCoordsInside
  */
 //=============================================================================
-std::vector<bool> GEOMImpl_IMeasureOperations::AreCoordsInside(Handle(GEOM_Object) theShape,
-                                                               const std::vector<double>& coords,
-                                                               double tolerance)
+std::vector<bool>
+GEOMImpl_IMeasureOperations::AreCoordsInside(Handle(GEOM_Object)        theShape,
+                                             const std::vector<double>& coords,
+                                             double                     tolerance)
 {
-  std::vector<bool> res;
+  std::vector<bool> isInsideRes;
   if (!theShape.IsNull()) {
     Handle(GEOM_Function) aRefShape = theShape->GetLastFunction();
     if (!aRefShape.IsNull()) {
       TopoDS_Shape aShape = aRefShape->GetValue();
-      if (!aShape.IsNull()) {
-        BRepClass3d_SolidClassifier SC(aShape);
-        unsigned int nb_points = coords.size()/3;
-        for (int i = 0; i < nb_points; i++) {
-          double x = coords[3*i];
-          double y = coords[3*i+1];
-          double z = coords[3*i+2];
-          gp_Pnt aPnt(x, y, z);
-          SC.Perform(aPnt, tolerance);
-          res.push_back( ( SC.State() == TopAbs_IN ) || ( SC.State() == TopAbs_ON ) );
+      if (!aShape.IsNull())
+      {
+        TopTools_IndexedMapOfShape mapShape;
+        {
+          TopExp_Explorer anExp;
+          for ( anExp.Init( aShape, TopAbs_SOLID ); anExp.More(); anExp.Next() )
+            mapShape.Add( anExp.Current() );
+          for ( anExp.Init( aShape, TopAbs_FACE, TopAbs_SOLID ); anExp.More(); anExp.Next() )
+            mapShape.Add( anExp.Current() );
+          for ( anExp.Init( aShape, TopAbs_EDGE, TopAbs_FACE ); anExp.More(); anExp.Next() )
+            mapShape.Add( anExp.Current() );
+          for ( anExp.Init( aShape, TopAbs_VERTEX, TopAbs_EDGE ); anExp.More(); anExp.Next() )
+            mapShape.Add( anExp.Current() ); //// ?????????
+        }
+        size_t nb_points = coords.size()/3, nb_points_inside = 0;
+        isInsideRes.resize( nb_points, false );
+
+        for ( int iS = 1; iS <= mapShape.Extent(); ++iS )
+        {
+          if ( nb_points_inside == nb_points )
+            break;
+          aShape = mapShape( iS );
+          switch ( aShape.ShapeType() ) {
+          case TopAbs_SOLID:
+          {
+            BRepClass3d_SolidClassifier SC( TopoDS::Solid( aShape ));
+            for ( size_t i = 0; i < nb_points; i++)
+            {
+              if ( isInsideRes[ i ]) continue;
+              gp_Pnt aPnt( coords[3*i], coords[3*i+1], coords[3*i+2] );
+              SC.Perform( aPnt, tolerance );
+              isInsideRes[ i ] = (( SC.State() == TopAbs_IN ) || ( SC.State() == TopAbs_ON ));
+              nb_points_inside += isInsideRes[ i ];
+            }
+            break;
+          }
+          case TopAbs_FACE:
+          {
+            Standard_Real u1,u2,v1,v2;
+            const TopoDS_Face&   face = TopoDS::Face( aShape );
+            Handle(Geom_Surface) surf = BRep_Tool::Surface( face );
+            surf->Bounds( u1,u2,v1,v2 );
+            GeomAPI_ProjectPointOnSurf project;
+            project.Init(surf, u1,u2, v1,v2, tolerance );
+            for ( size_t i = 0; i < nb_points; i++)
+            {
+              if ( isInsideRes[ i ]) continue;
+              gp_Pnt aPnt( coords[3*i], coords[3*i+1], coords[3*i+2] );
+              project.Perform( aPnt );
+              if ( project.IsDone() &&
+                   project.NbPoints() > 0 &&
+                   project.LowerDistance() <= tolerance )
+              {
+                Quantity_Parameter u, v;
+                project.LowerDistanceParameters(u, v);
+                gp_Pnt2d uv( u, v );
+                BRepClass_FaceClassifier FC ( face, uv, tolerance );
+                isInsideRes[ i ] = (( FC.State() == TopAbs_IN ) || ( FC.State() == TopAbs_ON ));
+                nb_points_inside += isInsideRes[ i ];
+              }
+            }
+            break;
+          }
+          case TopAbs_EDGE:
+          {
+            Standard_Real f,l;
+            const TopoDS_Edge&  edge = TopoDS::Edge( aShape );
+            Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f, l );
+            GeomAPI_ProjectPointOnCurve project;
+            project.Init( curve, f, l );
+            for ( size_t i = 0; i < nb_points; i++)
+            {
+              if ( isInsideRes[ i ]) continue;
+              gp_Pnt aPnt( coords[3*i], coords[3*i+1], coords[3*i+2] );
+              project.Perform( aPnt );
+              isInsideRes[ i ] = ( project.NbPoints() > 0 &&
+                                   project.LowerDistance() <= tolerance );
+              nb_points_inside += isInsideRes[ i ];
+            }
+            break;
+          }
+          case TopAbs_VERTEX:
+          {
+            gp_Pnt aVPnt = BRep_Tool::Pnt( TopoDS::Vertex( aShape ));
+            for ( size_t i = 0; i < nb_points; i++)
+            {
+              if ( isInsideRes[ i ]) continue;
+              gp_Pnt aPnt( coords[3*i], coords[3*i+1], coords[3*i+2] );
+              isInsideRes[ i ] = ( aPnt.SquareDistance( aVPnt ) <= tolerance * tolerance );
+              nb_points_inside += isInsideRes[ i ];
+            }
+            break;
+          }
+          default:;
+          } // switch ( aShape.ShapeType() )
         }
       }
     }
   }
-  return res;
+  return isInsideRes;
 }
 
 //=============================================================================
@@ -1881,10 +1670,15 @@ std::vector<bool> GEOMImpl_IMeasureOperations::AreCoordsInside(Handle(GEOM_Objec
  *  GetMinDistance
  */
 //=============================================================================
-Standard_Real GEOMImpl_IMeasureOperations::GetMinDistance
-  (Handle(GEOM_Object) theShape1, Handle(GEOM_Object) theShape2,
-   Standard_Real& X1, Standard_Real& Y1, Standard_Real& Z1,
-   Standard_Real& X2, Standard_Real& Y2, Standard_Real& Z2)
+Standard_Real
+GEOMImpl_IMeasureOperations::GetMinDistance (Handle(GEOM_Object) theShape1,
+                                             Handle(GEOM_Object) theShape2,
+                                             Standard_Real& X1,
+                                             Standard_Real& Y1,
+                                             Standard_Real& Z1,
+                                             Standard_Real& X2,
+                                             Standard_Real& Y2,
+                                             Standard_Real& Z2)
 {
   SetErrorCode(KO);
   Standard_Real MinDist = 1.e9;
@@ -1908,73 +1702,103 @@ Standard_Real GEOMImpl_IMeasureOperations::GetMinDistance
     OCC_CATCH_SIGNALS;
 #endif
 
-    // Issue 0020231: A min distance bug with torus and vertex.
-    // Make GetMinDistance() return zero if a sole VERTEX is inside any of SOLIDs
-
-    // which of shapes consists of only one vertex?
-    TopExp_Explorer exp1(aShape1,TopAbs_VERTEX), exp2(aShape2,TopAbs_VERTEX);
-    TopoDS_Shape V1 = exp1.More() ? exp1.Current() : TopoDS_Shape();
-    TopoDS_Shape V2 = exp2.More() ? exp2.Current() : TopoDS_Shape();
-    exp1.Next(); exp2.Next();
-    if ( exp1.More() ) V1.Nullify();
-    if ( exp2.More() ) V2.Nullify();
-    // vertex and container of solids
-    TopoDS_Shape V = V1.IsNull() ? V2 : V1;
-    TopoDS_Shape S = V1.IsNull() ? aShape1 : aShape2;
-    if ( !V.IsNull() ) {
-      // classify vertex against solids
-      gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( V ) );
-      for ( exp1.Init( S, TopAbs_SOLID ); exp1.More(); exp1.Next() ) {
-        BRepClass3d_SolidClassifier classifier( exp1.Current(), p, 1e-6);
-        if ( classifier.State() == TopAbs_IN ) {
-          p.Coord(X1, Y1, Z1);
-          p.Coord(X2, Y2, Z2);
-          SetErrorCode(OK);
-          return 0.0;
-        }
-      }
+    gp_Pnt aPnt1, aPnt2;
+
+    MinDist = GEOMUtils::GetMinDistance(aShape1, aShape2, aPnt1, aPnt2);
+
+    if (MinDist >= 0.0) {
+      aPnt1.Coord(X1, Y1, Z1);
+      aPnt2.Coord(X2, Y2, Z2);
+    } else {
+      return MinDist;
     }
-    // End Issue 0020231
+  }
+  catch (Standard_Failure) {
+    Handle(Standard_Failure) aFail = Standard_Failure::Caught();
+    SetErrorCode(aFail->GetMessageString());
+    return MinDist;
+  }
+
+  SetErrorCode(OK);
+  return MinDist;
+}
+
+//=======================================================================
+/*!
+ *  Get coordinates of closest points of two shapes
+ */
+//=======================================================================
+Standard_Integer GEOMImpl_IMeasureOperations::ClosestPoints (Handle(GEOM_Object) theShape1,
+                                                             Handle(GEOM_Object) theShape2,
+                                                             Handle(TColStd_HSequenceOfReal)& theDoubles)
+{
+  SetErrorCode(KO);
+  Standard_Integer nbSolutions = 0;
+
+  if (theShape1.IsNull() || theShape2.IsNull()) return nbSolutions;
+
+  Handle(GEOM_Function) aRefShape1 = theShape1->GetLastFunction();
+  Handle(GEOM_Function) aRefShape2 = theShape2->GetLastFunction();
+  if (aRefShape1.IsNull() || aRefShape2.IsNull()) return nbSolutions;
+
+  TopoDS_Shape aShape1 = aRefShape1->GetValue();
+  TopoDS_Shape aShape2 = aRefShape2->GetValue();
+  if (aShape1.IsNull() || aShape2.IsNull()) {
+    SetErrorCode("One of Objects has NULL Shape");
+    return nbSolutions;
+  }
+
+  // Compute the extremities
+  try {
+#if OCC_VERSION_LARGE > 0x06010000
+    OCC_CATCH_SIGNALS;
+#endif
 
     // skl 30.06.2008
     // additional workaround for bugs 19899, 19908 and 19910 from Mantis
-    gp_Pnt Ptmp1, Ptmp2;
-    double dist = CheckSingularCase(aShape1, aShape2, Ptmp1, Ptmp2);
-    if(dist>-1.0) {
-      Ptmp1.Coord(X1, Y1, Z1);
-      Ptmp2.Coord(X2, Y2, Z2);
+    gp_Pnt P1, P2;
+    double dist = GEOMUtils::GetMinDistanceSingular(aShape1, aShape2, P1, P2);
+    if (dist > -1.0) {
+      nbSolutions = 1;
+
+      theDoubles->Append(P1.X());
+      theDoubles->Append(P1.Y());
+      theDoubles->Append(P1.Z());
+      theDoubles->Append(P2.X());
+      theDoubles->Append(P2.Y());
+      theDoubles->Append(P2.Z());
+
       SetErrorCode(OK);
-      return dist;
+      return nbSolutions;
     }
 
     BRepExtrema_DistShapeShape dst (aShape1, aShape2);
     if (dst.IsDone()) {
-      gp_Pnt PMin1, PMin2, P1, P2;
+      nbSolutions = dst.NbSolution();
+      if (theDoubles.IsNull()) theDoubles = new TColStd_HSequenceOfReal;
 
-      for (int i = 1; i <= dst.NbSolution(); i++) {
+      gp_Pnt P1, P2;
+      for (int i = 1; i <= nbSolutions; i++) {
         P1 = dst.PointOnShape1(i);
         P2 = dst.PointOnShape2(i);
 
-        Standard_Real Dist = P1.Distance(P2);
-        if (MinDist > Dist) {
-          MinDist = Dist;
-          PMin1 = P1;
-          PMin2 = P2;
-        }
+        theDoubles->Append(P1.X());
+        theDoubles->Append(P1.Y());
+        theDoubles->Append(P1.Z());
+        theDoubles->Append(P2.X());
+        theDoubles->Append(P2.Y());
+        theDoubles->Append(P2.Z());
       }
-
-      PMin1.Coord(X1, Y1, Z1);
-      PMin2.Coord(X2, Y2, Z2);
     }
   }
   catch (Standard_Failure) {
     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
     SetErrorCode(aFail->GetMessageString());
-    return MinDist;
+    return nbSolutions;
   }
 
   SetErrorCode(OK);
-  return MinDist;
+  return nbSolutions;
 }
 
 //=======================================================================