0025398: Modeling Algorithms - Provide shape proximity detector

Correction of test case for issue CR25398

src/BRepExtrema/BRepExtrema_ShapeProximity.cxx

@@ -0,0 +1,645 @@
+// Created on: 2014-10-20
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <BRepExtrema_ShapeProximity.hxx>
+
+#include <Precision.hxx>
+#include <TopExp_Explorer.hxx>
+
+//=======================================================================
+//function : BRepExtrema_ShapeProximity
+//purpose  : Creates empty proximity tool
+//=======================================================================
+BRepExtrema_ShapeProximity::BRepExtrema_ShapeProximity (const Standard_Real theTolerance)
+: myTolerance     (theTolerance),
+  myPrimitiveSet1 (new BRepExtrema_TriangleSet),
+  myPrimitiveSet2 (new BRepExtrema_TriangleSet)
+{
+  // Should be initialized later
+  myIsDone = myIsInitS1 = myIsInitS2 = Standard_False;
+}
+
+//=======================================================================
+//function : BRepExtrema_ShapeProximity
+//purpose  : Creates proximity tool for the given two shapes
+//=======================================================================
+BRepExtrema_ShapeProximity::BRepExtrema_ShapeProximity (const TopoDS_Shape& theShape1,
+                                                        const TopoDS_Shape& theShape2,
+                                                        const Standard_Real theTolerance)
+: myTolerance     (theTolerance),
+  myPrimitiveSet1 (new BRepExtrema_TriangleSet),
+  myPrimitiveSet2 (new BRepExtrema_TriangleSet)
+{
+  LoadShape1 (theShape1);
+  LoadShape2 (theShape2);
+}
+
+//=======================================================================
+//function : LoadShape1
+//purpose  : Loads 1st shape into proximity tool
+//=======================================================================
+Standard_Boolean BRepExtrema_ShapeProximity::LoadShape1 (const TopoDS_Shape& theShape1)
+{
+  myFaceList1.Clear();
+
+  for (TopExp_Explorer anIter (theShape1, TopAbs_FACE); anIter.More(); anIter.Next())
+  {
+    myFaceList1.Append (static_cast<const TopoDS_Face&> (anIter.Current()));
+  }
+
+  myIsDone = Standard_False;
+
+  return myIsInitS1 = myPrimitiveSet1->Init (myFaceList1);
+}
+
+//=======================================================================
+//function : LoadShape2
+//purpose  : Loads 2nd shape into proximity tool
+//=======================================================================
+Standard_Boolean BRepExtrema_ShapeProximity::LoadShape2 (const TopoDS_Shape& theShape2)
+{
+  myFaceList2.Clear();
+
+  for (TopExp_Explorer anIter (theShape2, TopAbs_FACE); anIter.More(); anIter.Next())
+  {
+    myFaceList2.Append (static_cast<const TopoDS_Face&> (anIter.Current()));
+  }
+
+  myIsDone = Standard_False;
+
+  return myIsInitS2 = myPrimitiveSet2->Init (myFaceList2);
+}
+
+namespace
+{
+  //! Tool class to describe stack item in traverse function.
+  struct BRepExtrema_StackItem
+  {
+    Standard_Integer Node1;
+    Standard_Integer Node2;
+
+    BRepExtrema_StackItem (const Standard_Integer theNode1 = 0,
+                           const Standard_Integer theNode2 = 0)
+    : Node1 (theNode1),
+      Node2 (theNode2)
+    {
+      //
+    }
+  };
+
+  //! Bounding triangular prism for specified triangle.
+  class BRepExtrema_BoundingPrism
+  {
+  public:
+
+    //! Vertices of the prism.
+    BVH_Vec3d Vertices[6];
+
+    //! Edges of the prism.
+    BVH_Vec3d Edges[3];
+
+    //! Normal to prism caps.
+    BVH_Vec3d Normal;
+
+    //! Normals to prism edges.
+    BVH_Vec3d EdgeNormals[3];
+
+    //! Is prism initialized?
+    Standard_Boolean IsInited;
+
+  public:
+
+    //! Creates uninitialized bounding prism.
+    BRepExtrema_BoundingPrism() : IsInited (Standard_False)
+    {
+      //
+    }
+
+    //! Creates new bounding prism for the given triangle.
+    BRepExtrema_BoundingPrism (const BVH_Vec3d&    theVertex0,
+                               const BVH_Vec3d&    theVertex1,
+                               const BVH_Vec3d&    theVertex2,
+                               const Standard_Real theDeflect)
+    {
+      Init (theVertex0,
+            theVertex1,
+            theVertex2,
+            theDeflect);
+    }
+
+    //! Calculates bounding prism for the given triangle.
+    void Init (const BVH_Vec3d&    theVertex0,
+               const BVH_Vec3d&    theVertex1,
+               const BVH_Vec3d&    theVertex2,
+               const Standard_Real theDeflect)
+    {
+      Edges[0] = theVertex1 - theVertex0;
+      Edges[1] = theVertex2 - theVertex0;
+      Edges[2] = theVertex2 - theVertex1;
+
+      Normal = BVH_Vec3d::Cross (Edges[0], Edges[1]);
+
+      EdgeNormals[0] = BVH_Vec3d::Cross (Edges[0], Normal);
+      EdgeNormals[1] = BVH_Vec3d::Cross (Edges[1], Normal);
+      EdgeNormals[2] = BVH_Vec3d::Cross (Edges[2], Normal);
+
+      EdgeNormals[0] *= 1.0 / Max (EdgeNormals[0].Modulus(), Precision::Confusion());
+      EdgeNormals[1] *= 1.0 / Max (EdgeNormals[1].Modulus(), Precision::Confusion());
+      EdgeNormals[2] *= 1.0 / Max (EdgeNormals[2].Modulus(), Precision::Confusion());
+
+      const BVH_Vec3d aDirect01 = EdgeNormals[0] - EdgeNormals[1];
+      const BVH_Vec3d aDirect02 = EdgeNormals[0] + EdgeNormals[2];
+      const BVH_Vec3d aDirect12 = EdgeNormals[2] - EdgeNormals[1];
+
+      Vertices[0] = Vertices[3] = theVertex0 + aDirect01 * (theDeflect / aDirect01.Dot (EdgeNormals[0]));
+      Vertices[1] = Vertices[4] = theVertex1 + aDirect02 * (theDeflect / aDirect02.Dot (EdgeNormals[2]));
+      Vertices[2] = Vertices[5] = theVertex2 + aDirect12 * (theDeflect / aDirect12.Dot (EdgeNormals[2]));
+
+      const BVH_Vec3d aNormOffset = Normal * (theDeflect / Max (Normal.Modulus(), Precision::Confusion()));
+
+      for (Standard_Integer aVertIdx = 0; aVertIdx < 3; ++aVertIdx)
+      {
+        Vertices[aVertIdx + 0] += aNormOffset;
+        Vertices[aVertIdx + 3] -= aNormOffset;
+      }
+
+      IsInited = Standard_True;
+    }
+
+    //! Checks if two prisms are separated along the given axis.
+    Standard_Boolean Separated (const BRepExtrema_BoundingPrism& thePrism, const BVH_Vec3d& theAxis) const
+    {
+      Standard_Real aMin1 =  DBL_MAX;
+      Standard_Real aMax1 = -DBL_MAX;
+
+      Standard_Real aMin2 =  DBL_MAX;
+      Standard_Real aMax2 = -DBL_MAX;
+
+      for (Standard_Integer aVertIdx = 0; aVertIdx < 6; ++aVertIdx)
+      {
+        const Standard_Real aProj1 = Vertices[aVertIdx].Dot (theAxis);
+
+        aMin1 = Min (aMin1, aProj1);
+        aMax1 = Max (aMax1, aProj1);
+
+        const Standard_Real aProj2 = thePrism.Vertices[aVertIdx].Dot (theAxis);
+
+        aMin2 = Min (aMin2, aProj2);
+        aMax2 = Max (aMax2, aProj2);
+
+        if (aMin1 <= aMax2 && aMax1 >= aMin2)
+        {
+          return Standard_False;
+        }
+      }
+
+      return aMin1 > aMax2 || aMax1 < aMin2;
+    }
+  };
+
+  // =======================================================================
+  // function : Separated
+  // purpose  : Checks if triangles can be separated along the given axis
+  //            (projects vertices on this axis and performs interval test)
+  // =======================================================================
+  inline Standard_Boolean SeparateTriangles (const BVH_Vec3d& theTrg1Vert0,
+                                             const BVH_Vec3d& theTrg1Vert1,
+                                             const BVH_Vec3d& theTrg1Vert2,
+                                             const BVH_Vec3d& theTrg2Vert0,
+                                             const BVH_Vec3d& theTrg2Vert1,
+                                             const BVH_Vec3d& theTrg2Vert2,
+                                             const BVH_Vec3d& theSplitAxis)
+  {
+    const Standard_Real aA1 = theTrg1Vert0.Dot (theSplitAxis);
+    const Standard_Real aB1 = theTrg1Vert1.Dot (theSplitAxis);
+    const Standard_Real aC1 = theTrg1Vert2.Dot (theSplitAxis);
+
+    const Standard_Real aA2 = theTrg2Vert0.Dot (theSplitAxis);
+    const Standard_Real aB2 = theTrg2Vert1.Dot (theSplitAxis);
+    const Standard_Real aC2 = theTrg2Vert2.Dot (theSplitAxis);
+
+    const Standard_Real aMin1 = Min (aA1, Min (aB1, aC1));
+    const Standard_Real aMax1 = Max (aA1, Max (aB1, aC1));
+
+    if (aMax1 < Min (aA2, Min (aB2, aC2)))
+    {
+      return Standard_True;
+    }
+
+    return aMin1 > Max (aA2, Max (aB2, aC2));
+  }
+
+  // =======================================================================
+  // function : TrianglesIntersected
+  // purpose  : Checks if two triangles are intersected
+  //            (test uses SAT - Separating Axis Theorem)
+  // =======================================================================
+  Standard_Boolean TrianglesIntersected (const BVH_Vec3d& theTrg1Vert0,
+                                         const BVH_Vec3d& theTrg1Vert1,
+                                         const BVH_Vec3d& theTrg1Vert2,
+                                         const BVH_Vec3d& theTrg2Vert0,
+                                         const BVH_Vec3d& theTrg2Vert1,
+                                         const BVH_Vec3d& theTrg2Vert2)
+  {
+    const BVH_Vec3d aEdges1[3] = { theTrg1Vert1 - theTrg1Vert0,
+                                   theTrg1Vert2 - theTrg1Vert0,
+                                   theTrg1Vert2 - theTrg1Vert1 };
+
+    const BVH_Vec3d aTrg1Normal = BVH_Vec3d::Cross (aEdges1[0], aEdges1[1]);
+
+    if (SeparateTriangles (theTrg1Vert0,
+                           theTrg1Vert1,
+                           theTrg1Vert2,
+                           theTrg2Vert0,
+                           theTrg2Vert1,
+                           theTrg2Vert2,
+                           aTrg1Normal))
+    {
+      return Standard_False;
+    }
+
+    const BVH_Vec3d aEdges2[3] = { theTrg2Vert1 - theTrg2Vert0,
+                                   theTrg2Vert2 - theTrg2Vert0,
+                                   theTrg2Vert2 - theTrg2Vert1 };
+
+    const BVH_Vec3d aTrg2Normal = BVH_Vec3d::Cross (aEdges2[0], aEdges2[1]);
+
+    if (SeparateTriangles (theTrg1Vert0,
+                           theTrg1Vert1,
+                           theTrg1Vert2,
+                           theTrg2Vert0,
+                           theTrg2Vert1,
+                           theTrg2Vert2,
+                           aTrg2Normal))
+    {
+      return Standard_False;
+    }
+
+    for (Standard_Integer anIdx1 = 0; anIdx1 < 3; ++anIdx1)
+    {
+      for (Standard_Integer anIdx2 = 0; anIdx2 < 3; ++anIdx2)
+      {
+        const BVH_Vec3d aSplitAxis = BVH_Vec3d::Cross (aEdges1[anIdx1], aEdges2[anIdx2]);
+
+        if (SeparateTriangles (theTrg1Vert0,
+                               theTrg1Vert1,
+                               theTrg1Vert2,
+                               theTrg2Vert0,
+                               theTrg2Vert1,
+                               theTrg2Vert2,
+                               aSplitAxis))
+        {
+          return Standard_False;
+        }
+      }
+    }
+
+    return Standard_True;
+  }
+
+  // =======================================================================
+  // function : PrismsIntersected
+  // purpose  : Checks if two triangular prisms are intersected
+  //            (test uses SAT - Separating Axis Theorem)
+  // =======================================================================
+  Standard_Boolean PrismsIntersected (const BRepExtrema_BoundingPrism& thePrism1,
+                                      const BRepExtrema_BoundingPrism& thePrism2)
+  {
+    if (thePrism1.Separated (thePrism2, thePrism1.Normal))
+    {
+      return Standard_False;
+    }
+
+    if (thePrism1.Separated (thePrism2, thePrism2.Normal))
+    {
+      return Standard_False;
+    }
+
+    for (Standard_Integer anIdx = 0; anIdx < 3; ++anIdx)
+    {
+      if (thePrism1.Separated (thePrism2, thePrism1.EdgeNormals[anIdx]))
+      {
+        return Standard_False;
+      }
+    }
+
+    for (Standard_Integer anIdx = 0; anIdx < 3; ++anIdx)
+    {
+      if (thePrism1.Separated (thePrism2, thePrism2.EdgeNormals[anIdx]))
+      {
+        return Standard_False;
+      }
+    }
+
+    for (Standard_Integer anIdx1 = 0; anIdx1 < 4; ++anIdx1)
+    {
+      const BVH_Vec3d& aEdge1 = (anIdx1 == 3) ? thePrism1.Normal : thePrism1.Edges[anIdx1];
+
+      for (Standard_Integer anIdx2 = 0; anIdx2 < 4; ++anIdx2)
+      {
+        const BVH_Vec3d& aEdge2 = (anIdx2 == 3) ? thePrism2.Normal : thePrism2.Edges[anIdx2];
+
+        if (thePrism1.Separated (thePrism2, BVH_Vec3d::Cross (aEdge1, aEdge2)))
+        {
+          return Standard_False;
+        }
+      }
+    }
+
+    return Standard_True;
+  }
+
+  // =======================================================================
+  // function : OverlapBoxes
+  // purpose  : Checks if two boxes (AABBs) are overlapped
+  // =======================================================================
+  inline Standard_Boolean OverlapBoxes (const BVH_Vec3d&    theBoxMin1,
+                                        const BVH_Vec3d&    theBoxMax1,
+                                        const BVH_Vec3d&    theBoxMin2,
+                                        const BVH_Vec3d&    theBoxMax2,
+                                        const Standard_Real theTolerance)
+  {
+    // Check for overlap
+    return !(theBoxMin1.x() > theBoxMax2.x() + theTolerance ||
+             theBoxMax1.x() < theBoxMin2.x() - theTolerance ||
+             theBoxMin1.y() > theBoxMax2.y() + theTolerance ||
+             theBoxMax1.y() < theBoxMin2.y() - theTolerance ||
+             theBoxMin1.z() > theBoxMax2.z() + theTolerance ||
+             theBoxMax1.z() < theBoxMin2.z() - theTolerance);
+  }
+
+  //=======================================================================
+  //function : getSetOfFaces
+  //purpose  :
+  //=======================================================================
+  TColStd_PackedMapOfInteger& getSetOfFaces (BRepExtrema_OverlappedSubShapes& theShapes,
+                                             const Standard_Integer           theFaceIdx)
+  {
+    if (!theShapes.IsBound (theFaceIdx))
+    {
+      theShapes.Bind (theFaceIdx, TColStd_PackedMapOfInteger());
+    }
+
+    return theShapes.ChangeFind (theFaceIdx);
+  }
+}
+
+//=======================================================================
+//function : IntersectLeavesExact
+//purpose  : Narrow-phase of overlap test (exact intersection)
+//=======================================================================
+void BRepExtrema_ShapeProximity::IntersectLeavesExact (const BVH_Vec4i& theLeaf1,
+                                                       const BVH_Vec4i& theLeaf2)
+{
+  for (Standard_Integer aTrgIdx1 = theLeaf1.y(); aTrgIdx1 <= theLeaf1.z(); ++aTrgIdx1)
+  {
+    const Standard_Integer aFaceIdx1 = myPrimitiveSet1->GetFaceID (aTrgIdx1);
+
+    BVH_Vec3d aTrg1Vert1;
+    BVH_Vec3d aTrg1Vert2;
+    BVH_Vec3d aTrg1Vert3;
+
+    myPrimitiveSet1->GetVertices (aTrgIdx1,
+                                  aTrg1Vert1,
+                                  aTrg1Vert2,
+                                  aTrg1Vert3);
+
+    const Standard_Boolean aIsInSet = myOverlapSubShapes1.IsBound (aFaceIdx1);
+
+    for (Standard_Integer aTrgIdx2 = theLeaf2.y(); aTrgIdx2 <= theLeaf2.z(); ++aTrgIdx2)
+    {
+      const Standard_Integer aFaceIdx2 = myPrimitiveSet2->GetFaceID (aTrgIdx2);
+
+      if (!aIsInSet || !myOverlapSubShapes1.Find (aFaceIdx1).Contains (aFaceIdx2))
+      {
+        BVH_Vec3d aTrg2Vert1;
+        BVH_Vec3d aTrg2Vert2;
+        BVH_Vec3d aTrg2Vert3;
+
+        myPrimitiveSet2->GetVertices (aTrgIdx2, aTrg2Vert1, aTrg2Vert2, aTrg2Vert3);
+
+        if (TrianglesIntersected (aTrg1Vert1,
+                                  aTrg1Vert2,
+                                  aTrg1Vert3,
+                                  aTrg2Vert1,
+                                  aTrg2Vert2,
+                                  aTrg2Vert3))
+        {
+          getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
+          getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
+        }
+      }
+    }
+  }
+}
+
+//=======================================================================
+//function : IntersectLeavesToler
+//purpose  : Narrow-phase of overlap test (with non-zero tolerance)
+//=======================================================================
+void BRepExtrema_ShapeProximity::IntersectLeavesToler (const BVH_Vec4i& theLeaf1,
+                                                       const BVH_Vec4i& theLeaf2)
+{
+  for (Standard_Integer aTrgIdx1 = theLeaf1.y(); aTrgIdx1 <= theLeaf1.z(); ++aTrgIdx1)
+  {
+    const Standard_Integer aFaceIdx1 = myPrimitiveSet1->GetFaceID (aTrgIdx1);
+
+    BVH_Vec3d aTrg1Vert1;
+    BVH_Vec3d aTrg1Vert2;
+    BVH_Vec3d aTrg1Vert3;
+
+    myPrimitiveSet1->GetVertices (aTrgIdx1,
+                                  aTrg1Vert1,
+                                  aTrg1Vert2,
+                                  aTrg1Vert3);
+
+    BRepExtrema_BoundingPrism aPrism1; // not initialized
+
+    const Standard_Boolean aIsInSet = myOverlapSubShapes1.IsBound (aFaceIdx1);
+
+    for (Standard_Integer aTrgIdx2 = theLeaf2.y(); aTrgIdx2 <= theLeaf2.z(); ++aTrgIdx2)
+    {
+      const Standard_Integer aFaceIdx2 = myPrimitiveSet2->GetFaceID (aTrgIdx2);
+
+      if (!aIsInSet || !myOverlapSubShapes1.Find (aFaceIdx1).Contains (aFaceIdx2))
+      {
+        if (!aPrism1.IsInited)
+        {
+          aPrism1.Init (aTrg1Vert1, aTrg1Vert2, aTrg1Vert3, myTolerance);
+        }
+
+        BVH_Vec3d aTrg2Vert1;
+        BVH_Vec3d aTrg2Vert2;
+        BVH_Vec3d aTrg2Vert3;
+
+        myPrimitiveSet2->GetVertices (aTrgIdx2,
+                                      aTrg2Vert1,
+                                      aTrg2Vert2,
+                                      aTrg2Vert3);
+
+        BRepExtrema_BoundingPrism aPrism2 (aTrg2Vert1,
+                                           aTrg2Vert2,
+                                           aTrg2Vert3,
+                                           myTolerance);
+
+        if (PrismsIntersected (aPrism1, aPrism2))
+        {
+          getSetOfFaces (myOverlapSubShapes1, aFaceIdx1).Add (aFaceIdx2);
+          getSetOfFaces (myOverlapSubShapes2, aFaceIdx2).Add (aFaceIdx1);
+        }
+      }
+    }
+  }
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : Performs search for overlapped faces
+//=======================================================================
+void BRepExtrema_ShapeProximity::Perform()
+{
+  if (myIsDone || !myIsInitS1 || !myIsInitS2)
+  {
+    return;
+  }
+
+  BRepExtrema_StackItem aStack[96];
+
+  const NCollection_Handle<BVH_Tree<Standard_Real, 3> >& aBVH1 = myPrimitiveSet1->BVH();
+  const NCollection_Handle<BVH_Tree<Standard_Real, 3> >& aBVH2 = myPrimitiveSet2->BVH();
+
+  if (aBVH1.IsNull() || aBVH2.IsNull())
+  {
+    return;
+  }
+
+  BRepExtrema_StackItem aNodes; // current pair of nodes
+
+  Standard_Integer aHead = -1; // stack head position
+
+  for (;;)
+  {
+    BVH_Vec4i aNodeData1 = aBVH1->NodeInfoBuffer()[aNodes.Node1];
+    BVH_Vec4i aNodeData2 = aBVH2->NodeInfoBuffer()[aNodes.Node2];
+
+    if (aNodeData1.x() != 0 && aNodeData2.x() != 0) // leaves
+    {
+      if (myTolerance == 0.0)
+      {
+        IntersectLeavesExact (aNodeData1, aNodeData2);
+      }
+      else
+      {
+        IntersectLeavesToler (aNodeData1, aNodeData2);
+      }
+
+      if (aHead < 0)
+        break;
+
+      aNodes = aStack[aHead--];
+    }
+    else
+    {
+      BRepExtrema_StackItem aPairsToProcess[4];
+
+      Standard_Integer aNbPairs = 0;
+
+      if (aNodeData1.x() == 0) // inner node
+      {
+        const BVH_Vec3d& aMinPntLft1 = aBVH1->MinPoint (aNodeData1.y());
+        const BVH_Vec3d& aMaxPntLft1 = aBVH1->MaxPoint (aNodeData1.y());
+        const BVH_Vec3d& aMinPntRgh1 = aBVH1->MinPoint (aNodeData1.z());
+        const BVH_Vec3d& aMaxPntRgh1 = aBVH1->MaxPoint (aNodeData1.z());
+
+        if (aNodeData2.x() == 0) // inner node
+        {
+          const BVH_Vec3d& aMinPntLft2 = aBVH2->MinPoint (aNodeData2.y());
+          const BVH_Vec3d& aMaxPntLft2 = aBVH2->MaxPoint (aNodeData2.y());
+          const BVH_Vec3d& aMinPntRgh2 = aBVH2->MinPoint (aNodeData2.z());
+          const BVH_Vec3d& aMaxPntRgh2 = aBVH2->MaxPoint (aNodeData2.z());
+
+          if (OverlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntLft2, aMaxPntLft2, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodeData2.y());
+          }
+          if (OverlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntRgh2, aMaxPntRgh2, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodeData2.z());
+          }
+          if (OverlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntLft2, aMaxPntLft2, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodeData2.y());
+          }
+          if (OverlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntRgh2, aMaxPntRgh2, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodeData2.z());
+          }
+        }
+        else
+        {
+          const BVH_Vec3d& aMinPntLeaf = aBVH2->MinPoint (aNodes.Node2);
+          const BVH_Vec3d& aMaxPntLeaf = aBVH2->MaxPoint (aNodes.Node2);
+
+          if (OverlapBoxes (aMinPntLft1, aMaxPntLft1, aMinPntLeaf, aMaxPntLeaf, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.y(), aNodes.Node2);
+          }
+          if (OverlapBoxes (aMinPntRgh1, aMaxPntRgh1, aMinPntLeaf, aMaxPntLeaf, myTolerance))
+          {
+            aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodeData1.z(), aNodes.Node2);
+          }
+        }
+      }
+      else
+      {
+        const BVH_Vec3d& aMinPntLeaf = aBVH1->MinPoint (aNodes.Node1);
+        const BVH_Vec3d& aMaxPntLeaf = aBVH1->MaxPoint (aNodes.Node1);
+
+        const BVH_Vec3d& aMinPntLft2 = aBVH2->MinPoint (aNodeData2.y());
+        const BVH_Vec3d& aMaxPntLft2 = aBVH2->MaxPoint (aNodeData2.y());
+        const BVH_Vec3d& aMinPntRgh2 = aBVH2->MinPoint (aNodeData2.z());
+        const BVH_Vec3d& aMaxPntRgh2 = aBVH2->MaxPoint (aNodeData2.z());
+
+        if (OverlapBoxes (aMinPntLft2, aMaxPntLft2, aMinPntLeaf, aMaxPntLeaf, myTolerance))
+        {
+          aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodes.Node1, aNodeData2.y());
+        }
+        if (OverlapBoxes (aMinPntRgh2, aMaxPntRgh2, aMinPntLeaf, aMaxPntLeaf, myTolerance))
+        {
+          aPairsToProcess[aNbPairs++] = BRepExtrema_StackItem (aNodes.Node1, aNodeData2.z());
+        }
+      }
+
+      if (aNbPairs > 0)
+      {
+        aNodes = aPairsToProcess[0];
+
+        for (Standard_Integer anIdx = 1; anIdx < aNbPairs; ++anIdx)
+        {
+          aStack[++aHead] = aPairsToProcess[anIdx];
+        }
+      }
+      else
+      {
+        if (aHead < 0)
+          break;
+
+        aNodes = aStack[aHead--];
+      }
+    }
+  }
+
+  myIsDone = Standard_True;
+}

src/BRepExtrema/BRepExtrema_ShapeProximity.hxx

@@ -0,0 +1,154 @@
+// Created on: 2014-10-20
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BRepExtrema_ShapeProximity_HeaderFile
+#define _BRepExtrema_ShapeProximity_HeaderFile
+
+#include <BVH_Geometry.hxx>
+#include <BRepExtrema_TriangleSet.hxx>
+#include <TColStd_PackedMapOfInteger.hxx>
+#include <NCollection_DataMap.hxx>
+
+//! Set of overlapped sub-shapes.
+typedef NCollection_DataMap<Standard_Integer, TColStd_PackedMapOfInteger > BRepExtrema_OverlappedSubShapes;
+
+//! Tool class for shape proximity detection.
+//! For two given shapes and given tolerance (offset from the mesh) the algorithm allows
+//! to determine whether or not they are overlapped. The algorithm input consists of any
+//! shapes which can be decomposed into individual faces (used as basic shape elements).
+//! High performance is achieved through the use of existing triangulation of faces. So
+//! poly triangulation (with the desired deflection) should already be built. Note that
+//! solution is approximate (and corresponds to the deflection used for triangulation).
+//!
+//! The algorithm can be run in two modes. If tolerance is set to zero, the algorithm
+//! will detect only intersecting faces (containing triangles with common points). If
+//! tolerance is set to positive value, the algorithm will also detect faces located
+//! on distance less than the given tolerance from each other.
+class BRepExtrema_ShapeProximity
+{
+ public:
+
+  //! Creates empty proximity tool.
+  Standard_EXPORT BRepExtrema_ShapeProximity (const Standard_Real theTolerance = 0.0);
+
+  //! Creates proximity tool for the given two shapes.
+  Standard_EXPORT BRepExtrema_ShapeProximity (const TopoDS_Shape& theShape1,
+                                              const TopoDS_Shape& theShape2,
+                                              const Standard_Real theTolerance = 0.0);
+
+public:
+
+  //! Returns tolerance value for overlap test (distance between shapes).
+  Standard_Real Tolerance() const
+  {
+    return myTolerance;
+  }
+
+  //! Sets tolerance value for overlap test (distance between shapes).
+  void SetTolerance (const Standard_Real theTolerance)
+  {
+    myTolerance = theTolerance;
+  }
+
+  //! Loads 1st shape into proximity tool.
+  Standard_EXPORT Standard_Boolean LoadShape1 (const TopoDS_Shape& theShape1);
+
+  //! Loads 2nd shape into proximity tool.
+  Standard_EXPORT Standard_Boolean LoadShape2 (const TopoDS_Shape& theShape2);
+
+  //! Performs search for overlapped faces.
+  Standard_EXPORT void Perform();
+
+  //! True if the search is completed.
+  Standard_Boolean IsDone() const
+  { 
+    return myIsDone;
+  }
+
+  //! Returns set of all the face triangles of the 1st shape.
+  const NCollection_Handle<BRepExtrema_TriangleSet>& PrimitiveSet1() const
+  {
+    return myPrimitiveSet1;
+  }
+
+  //! Returns set of all the face triangles of the 2nd shape.
+  const NCollection_Handle<BRepExtrema_TriangleSet>& PrimitiveSet2() const
+  {
+    return myPrimitiveSet2;
+  }
+
+  //! Returns set of IDs of overlapped faces of 1st shape.
+  const BRepExtrema_OverlappedSubShapes& OverlapSubShapes1() const
+  {
+    return myOverlapSubShapes1;
+  }
+
+  //! Returns set of IDs of overlapped faces of 2nd shape.
+  const BRepExtrema_OverlappedSubShapes& OverlapSubShapes2() const
+  {
+    return myOverlapSubShapes2;
+  }
+
+  //! Returns sub-shape from 1st shape with the given index.
+  const TopoDS_Face& GetSubShape1 (const Standard_Integer theID) const
+  {
+    return myFaceList1.Value (theID);
+  }
+
+  //! Returns sub-shape from 1st shape with the given index.
+  const TopoDS_Face& GetSubShape2 (const Standard_Integer theID) const
+  {
+    return myFaceList2.Value (theID);
+  }
+
+protected:
+
+  //! Performs narrow-phase of overlap test (exact intersection).
+  void IntersectLeavesExact (const BVH_Vec4i& theLeaf1, const BVH_Vec4i& theLeaf2);
+
+  //! Performs narrow-phase of overlap test (intersection with non-zero tolerance).
+  void IntersectLeavesToler (const BVH_Vec4i& theLeaf1, const BVH_Vec4i& theLeaf2);
+
+private:
+
+  //! Maximum overlapping distance.
+  Standard_Real myTolerance;
+
+  //! Is the 1st shape initialized?
+  Standard_Boolean myIsInitS1;
+  //! Is the 2nd shape initialized?
+  Standard_Boolean myIsInitS2;
+
+  //! List of faces of the 1st shape.
+  BRepExtrema_ShapeList myFaceList1;
+  //! List of faces of the 2nd shape.
+  BRepExtrema_ShapeList myFaceList2;
+
+  //! Set of all the face triangles of the 1st shape.
+  NCollection_Handle<BRepExtrema_TriangleSet> myPrimitiveSet1;
+  //! Set of all the face triangles of the 2nd shape.
+  NCollection_Handle<BRepExtrema_TriangleSet> myPrimitiveSet2;
+
+  //! Set of overlapped faces of 1st shape.
+  BRepExtrema_OverlappedSubShapes myOverlapSubShapes1;
+  //! Set of overlapped faces of 2nd shape.
+  BRepExtrema_OverlappedSubShapes myOverlapSubShapes2;
+
+  //! Is overlap test completed?
+  Standard_Boolean myIsDone;
+
+};
+
+#endif // _BRepExtrema_ShapeProximity_HeaderFile

src/BRepExtrema/BRepExtrema_TriangleSet.cxx

@@ -0,0 +1,226 @@
+// Created on: 2014-10-20
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <BRepExtrema_TriangleSet.hxx>
+
+#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BVH_LinearBuilder.hxx>
+#include <Poly_Triangulation.hxx>
+#include <TColgp_Array1OfPnt2d.hxx>
+
+//=======================================================================
+//function : BRepExtrema_TriangleSet
+//purpose  : Creates empty triangle set
+//=======================================================================
+BRepExtrema_TriangleSet::BRepExtrema_TriangleSet()
+: BVH_PrimitiveSet<Standard_Real, 3>()
+{
+  // Set default builder - linear BVH (LBVH)
+  myBuilder = new BVH_LinearBuilder<Standard_Real, 3> (5, 32);
+}
+
+//=======================================================================
+//function : BRepExtrema_TriangleSet
+//purpose  : Creates triangle set from the given face
+//=======================================================================
+BRepExtrema_TriangleSet::BRepExtrema_TriangleSet (const BRepExtrema_ShapeList& theFaces)
+: BVH_PrimitiveSet<Standard_Real, 3>()
+{
+  // Set default builder - linear BVH (LBVH)
+  myBuilder = new BVH_LinearBuilder<Standard_Real, 3> (5, 32);
+
+  Init (theFaces);
+}
+
+//=======================================================================
+//function : ~BRepExtrema_TriangleSet
+//purpose  : Releases resources of triangle set
+//=======================================================================
+BRepExtrema_TriangleSet::~BRepExtrema_TriangleSet()
+{
+  //
+}
+
+//=======================================================================
+//function : Size
+//purpose  : Returns total number of triangles
+//=======================================================================
+Standard_Integer BRepExtrema_TriangleSet::Size() const
+{
+  return static_cast<Standard_Integer> (myTriangles.size());
+}
+
+//=======================================================================
+//function : Box
+//purpose  : Returns AABB of the given triangle
+//=======================================================================
+BVH_Box<Standard_Real, 3> BRepExtrema_TriangleSet::Box (const Standard_Integer theIndex) const
+{
+  const BVH_Vec4i& aTriangle = myTriangles[theIndex];
+
+  BVH_Vec3d aMinPnt = myVertexArray[aTriangle.x()].cwiseMin (
+    myVertexArray[aTriangle.y()].cwiseMin (myVertexArray[aTriangle.z()]));
+
+  BVH_Vec3d aMaxPnt = myVertexArray[aTriangle.x()].cwiseMax (
+    myVertexArray[aTriangle.y()].cwiseMax (myVertexArray[aTriangle.z()]));
+
+  return BVH_Box<Standard_Real, 3> (aMinPnt, aMaxPnt);
+}
+
+//=======================================================================
+//function : Center
+//purpose  : Returns centroid position along specified axis
+//=======================================================================
+Standard_Real BRepExtrema_TriangleSet::Center (const Standard_Integer theIndex, const Standard_Integer theAxis) const
+{
+  const BVH_Vec4i& aTriangle = myTriangles[theIndex];
+
+  if (theAxis == 0)
+  {
+    return (1.0 / 3.0) * (myVertexArray[aTriangle.x()].x() +
+                          myVertexArray[aTriangle.y()].x() +
+                          myVertexArray[aTriangle.z()].x());
+  }
+  else if (theAxis == 1)
+  {
+    return (1.0 / 3.0) * (myVertexArray[aTriangle.x()].y() +
+                          myVertexArray[aTriangle.y()].y() +
+                          myVertexArray[aTriangle.z()].y());
+  }
+  else
+  {
+    return (1.0 / 3.0) * (myVertexArray[aTriangle.x()].z() +
+                          myVertexArray[aTriangle.y()].z() +
+                          myVertexArray[aTriangle.z()].z());
+  }
+}
+
+//=======================================================================
+//function : Swap
+//purpose  : Swaps indices of two specified triangles
+//=======================================================================
+void BRepExtrema_TriangleSet::Swap (const Standard_Integer theIndex1, const Standard_Integer theIndex2)
+{
+  std::swap (myTriangles[theIndex1],
+             myTriangles[theIndex2]);
+}
+
+//=======================================================================
+//function : GetFaceID
+//purpose  : Returns face ID of the given triangle
+//=======================================================================
+Standard_Integer BRepExtrema_TriangleSet::GetFaceID (const Standard_Integer theIndex) const
+{
+  return myTriangles[theIndex].w();
+}
+
+//=======================================================================
+//function : GetVertices
+//purpose  : Returns vertices of the given triangle
+//=======================================================================
+void BRepExtrema_TriangleSet::GetVertices (const Standard_Integer theIndex,
+                                           BVH_Vec3d&             theVertex1,
+                                           BVH_Vec3d&             theVertex2,
+                                           BVH_Vec3d&             theVertex3) const
+{
+  BVH_Vec4i aTriangle = myTriangles[theIndex];
+
+  theVertex1 = myVertexArray[aTriangle.x()];
+  theVertex2 = myVertexArray[aTriangle.y()];
+  theVertex3 = myVertexArray[aTriangle.z()];
+}
+
+//=======================================================================
+//function : Clear
+//purpose  : Clears triangle set data
+//=======================================================================
+void BRepExtrema_TriangleSet::Clear()
+{
+  BVH_Array4i anEmptyTriangles;
+  myTriangles.swap (anEmptyTriangles);
+
+  BVH_Array2d anEmptyVertUVArray;
+  myVertUVArray.swap (anEmptyVertUVArray);
+
+  BVH_Array3d anEmptyVertexArray;
+  myVertexArray.swap (anEmptyVertexArray);
+}
+
+//=======================================================================
+//function : Init
+//purpose  : Initializes triangle set
+//=======================================================================
+Standard_Boolean BRepExtrema_TriangleSet::Init (const BRepExtrema_ShapeList& theFaces)
+{
+  Clear();
+
+  for (Standard_Integer aFaceIdx = 0; aFaceIdx < theFaces.Size(); ++aFaceIdx)
+  {
+    TopLoc_Location aLocation;
+
+    Handle(Poly_Triangulation) aTriangulation =
+      BRep_Tool::Triangulation (theFaces (aFaceIdx), aLocation);
+
+    if (aTriangulation.IsNull())
+    {
+      return Standard_False;
+    }
+
+    BRepAdaptor_Surface aFaceAdaptor (theFaces (aFaceIdx), Standard_False);
+
+    const Standard_Integer aVertOffset =
+      static_cast<Standard_Integer> (myVertexArray.size()) - 1;
+
+    for (Standard_Integer aVertIdx = 1; aVertIdx <= aTriangulation->NbNodes(); ++aVertIdx)
+    {
+      gp_Pnt aVertex = aTriangulation->Nodes().Value (aVertIdx);
+
+      aVertex.Transform (aLocation.Transformation());
+
+      myVertexArray.push_back (BVH_Vec3d (aVertex.X(),
+                                          aVertex.Y(),
+                                          aVertex.Z()));
+
+      const Standard_Real aU = aTriangulation->UVNodes().Value (aVertIdx).X();
+      const Standard_Real aV = aTriangulation->UVNodes().Value (aVertIdx).Y();
+
+      myVertUVArray.push_back (BVH_Vec2d (aU, aV));
+    }
+
+    for (Standard_Integer aTriIdx = 1; aTriIdx <= aTriangulation->NbTriangles(); ++aTriIdx)
+    {
+      Standard_Integer aVertex1;
+      Standard_Integer aVertex2;
+      Standard_Integer aVertex3;
+
+      aTriangulation->Triangles().Value (aTriIdx).Get (aVertex1,
+                                                       aVertex2,
+                                                       aVertex3);
+
+      myTriangles.push_back (BVH_Vec4i (aVertex1 + aVertOffset,
+                                        aVertex2 + aVertOffset,
+                                        aVertex3 + aVertOffset,
+                                        aFaceIdx));
+    }
+  }
+
+  MarkDirty(); // needs BVH rebuilding
+
+  Standard_ASSERT_RETURN (!BVH().IsNull(),
+    "Error: Failed to build BVH for primitive set", Standard_False);
+
+  return Standard_True;
+}

src/BRepExtrema/BRepExtrema_TriangleSet.hxx

@@ -0,0 +1,84 @@
+// Created on: 2014-10-20
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BRepExtrema_TriangleSet_HeaderFile
+#define _BRepExtrema_TriangleSet_HeaderFile
+
+#include <TopoDS_Face.hxx>
+
+#include <BVH_PrimitiveSet.hxx>
+
+//! List of shapes and their IDs for collision detection.
+typedef NCollection_Vector<TopoDS_Face> BRepExtrema_ShapeList;
+
+//! Triangle set corresponding to specific face.
+class BRepExtrema_TriangleSet : public BVH_PrimitiveSet<Standard_Real, 3>
+{
+public:
+
+  //! Creates empty triangle set.
+  Standard_EXPORT BRepExtrema_TriangleSet();
+
+  //! Creates triangle set from the given face.
+  Standard_EXPORT BRepExtrema_TriangleSet (const BRepExtrema_ShapeList& theFaces);
+
+  //! Releases resources of triangle set.
+  Standard_EXPORT ~BRepExtrema_TriangleSet();
+
+public: //! @name methods implementing BVH set interface
+
+  //! Returns total number of triangles.
+  Standard_Integer Size() const;
+
+  //! Returns AABB of the given triangle.
+  BVH_Box<Standard_Real, 3> Box (const Standard_Integer theIndex) const;
+
+  //! Returns centroid position along specified axis.
+  Standard_Real Center (const Standard_Integer theIndex, const Standard_Integer theAxis) const;
+
+  //! Swaps indices of two specified triangles.
+  void Swap (const Standard_Integer theIndex1, const Standard_Integer theIndex2);
+
+public:
+
+  //! Clears triangle set data.
+  Standard_EXPORT void Clear();
+
+  //! Initializes triangle set.
+  Standard_EXPORT Standard_Boolean Init (const BRepExtrema_ShapeList& theFaces);
+
+  //! Returns vertices of the given triangle.
+  Standard_EXPORT void GetVertices (const Standard_Integer theIndex,
+                                    BVH_Vec3d&             theVertex1,
+                                    BVH_Vec3d&             theVertex2,
+                                    BVH_Vec3d&             theVertex3) const;
+
+  //! Returns face ID of the given triangle.
+  Standard_EXPORT Standard_Integer GetFaceID (const Standard_Integer theIndex) const;
+
+protected:
+
+  //! Array of vertex indices.
+  BVH_Array4i myTriangles;
+
+  //! Array of vertex UV params.
+  BVH_Array2d myVertUVArray;
+
+  //! Array of vertex coordinates.
+  BVH_Array3d myVertexArray;
+
+};
+
+#endif // _BRepExtrema_TriangleSet_HeaderFile

src/BRepExtrema/FILES

@@ -17,3 +17,7 @@ BRepExtrema_Poly.cxx
 BRepExtrema_SeqOfSolution.hxx
 BRepExtrema_SolutionElem.hxx
 BRepExtrema_SupportType.hxx
+BRepExtrema_TriangleSet.hxx
+BRepExtrema_TriangleSet.cxx
+BRepExtrema_ShapeProximity.hxx
+BRepExtrema_ShapeProximity.cxx

src/BRepTest/BRepTest_ExtremaCommands.cxx

@@ -20,9 +20,16 @@
 #include <BRepTest.hxx>
 #include <BRepExtrema_Poly.hxx>
 #include <BRepExtrema_DistShapeShape.hxx>
+#include <BRepExtrema_ShapeProximity.hxx>
 #include <BRepLib_MakeEdge.hxx>
 #include <BRepLib_MakeVertex.hxx>
+#include <TopoDS_Builder.hxx>
+#include <TopoDS_Compound.hxx>
 #include <Draw.hxx>
+#include <OSD_Timer.hxx>
+#include <TCollection_AsciiString.hxx>
+#include <TColStd_MapIteratorOfPackedMapOfInteger.hxx>
+
 
 //#ifdef WNT
 #include <stdio.h>
@@ -131,6 +138,143 @@ static Standard_Integer distmini(Draw_Interpretor& di, Standard_Integer n, const
   return 0;
 }
 
+//==============================================================================
+//function : ShapeProximity
+//purpose  :
+//==============================================================================
+static int ShapeProximity (Draw_Interpretor& theDI, Standard_Integer theNbArgs, const char** theArgs)
+{
+  if (theNbArgs < 3 || theNbArgs > 6)
+  {
+    std::cout << "Usage: " << theArgs[0] <<
+      " Shape1 Shape2 [-tol <value>] [-profile]" << std::endl;
+
+    return 1;
+  }
+
+  TopoDS_Shape aShape1 = DBRep::Get (theArgs[1]);
+  TopoDS_Shape aShape2 = DBRep::Get (theArgs[2]);
+
+  if (aShape1.IsNull() || aShape2.IsNull())
+  {
+    std::cout << "Error: Failed to find specified shapes" << std::endl;
+    return 1;
+  }
+
+  BRepExtrema_ShapeProximity aTool;
+
+  Standard_Boolean aProfile = Standard_False;
+
+  for (Standard_Integer anArgIdx = 3; anArgIdx < theNbArgs; ++anArgIdx)
+  {
+    TCollection_AsciiString aFlag (theArgs[anArgIdx]);
+    aFlag.LowerCase();
+
+    if (aFlag == "-tol")
+    {
+      if (++anArgIdx >= theNbArgs)
+      {
+        std::cout << "Error: wrong syntax at argument '" << aFlag << std::endl;
+        return 1;
+      }
+
+      const Standard_Real aTolerance = Draw::Atof (theArgs[anArgIdx]);
+      if (aTolerance < 0.0)
+      {
+        std::cout << "Error: Tolerance value should be non-negative" << std::endl;
+        return 1;
+      }
+      else
+      {
+        aTool.SetTolerance (aTolerance);
+      }
+    }
+
+    if (aFlag == "-profile")
+    {
+      aProfile = Standard_True;
+    }
+  }
+
+  Standard_Real aInitTime = 0.0;
+  Standard_Real aWorkTime = 0.0;
+
+  OSD_Timer aTimer;
+
+  if (aProfile)
+  {
+    aTimer.Start();
+  }
+
+  aTool.LoadShape1 (aShape1);
+  aTool.LoadShape2 (aShape2);
+
+  if (aProfile)
+  {
+    aInitTime = aTimer.ElapsedTime();
+    aTimer.Reset();
+    aTimer.Start();
+  }
+
+  // Perform shape proximity test
+  aTool.Perform();
+
+  if (aProfile)
+  {
+    aWorkTime = aTimer.ElapsedTime();
+    aTimer.Stop();
+  }
+
+  if (!aTool.IsDone())
+  {
+    std::cout << "Error: Failed to perform proximity test" << std::endl;
+    return 1;
+  }
+
+  if (aProfile)
+  {
+    theDI << "Number of primitives in shape 1: " << aTool.PrimitiveSet1()->Size() << "\n";
+    theDI << "Number of primitives in shape 2: " << aTool.PrimitiveSet2()->Size() << "\n";
+    theDI << "Building data structures: " << aInitTime << "\n";
+    theDI << "Executing proximity test: " << aWorkTime << "\n";
+  }
+
+  TopoDS_Builder  aCompBuilder;
+
+  TopoDS_Compound aFaceCompound1;
+  aCompBuilder.MakeCompound (aFaceCompound1);
+
+  for (BRepExtrema_OverlappedSubShapes::Iterator anIt1 (aTool.OverlapSubShapes1()); anIt1.More(); anIt1.Next())
+  {
+    TCollection_AsciiString aStr = TCollection_AsciiString (theArgs[1]) + "_" + (anIt1.Key() + 1);
+
+    const TopoDS_Face& aFace = aTool.GetSubShape1 (anIt1.Key());
+    aCompBuilder.Add (aFaceCompound1, aFace);
+    DBRep::Set (aStr.ToCString(), aFace);
+
+    theDI << aStr << " \n";
+  }
+
+  TopoDS_Compound aFaceCompound2;
+  aCompBuilder.MakeCompound (aFaceCompound2);
+
+  for (BRepExtrema_OverlappedSubShapes::Iterator anIt2 (aTool.OverlapSubShapes2()); anIt2.More(); anIt2.Next())
+  {
+    TCollection_AsciiString aStr = TCollection_AsciiString (theArgs[2]) + "_" + (anIt2.Key() + 1);
+
+    const TopoDS_Face& aFace = aTool.GetSubShape2 (anIt2.Key());
+    aCompBuilder.Add (aFaceCompound2, aFace);
+    DBRep::Set (aStr.ToCString(), aFace);
+
+    theDI << aStr << " \n";
+  }
+
+  DBRep::Set ((TCollection_AsciiString (theArgs[1]) + "_" + "overlapped").ToCString(), aFaceCompound1);
+  DBRep::Set ((TCollection_AsciiString (theArgs[2]) + "_" + "overlapped").ToCString(), aFaceCompound2);
+
+  return 0;
+}
+
 //=======================================================================
 //function : ExtremaCommands
 //purpose  : 
@@ -157,4 +301,16 @@ void BRepTest::ExtremaCommands (Draw_Interpretor& theCommands)
                    __FILE__,
                    distmini,
                    aGroup);
+
+  theCommands.Add ("proximity",
+                   "proximity Shape1 Shape2 [-tol <value>] [-profile]"
+                   "\n\t\t: Searches for pairs of overlapping faces of the given shapes."
+                   "\n\t\t: The options are:"
+                   "\n\t\t:   -tol     : non-negative tolerance value used for overlapping"
+                   "\n\t\t:              test (for zero tolerance, the strict intersection"
+                   "\n\t\t:              test will be performed)"
+                   "\n\t\t:   -profile : outputs execution time for main algorithm stages",
+                   __FILE__,
+                   ShapeProximity,
+                   aGroup);
 }

tests/bugs/begin

@@ -277,3 +277,29 @@ proc checkList {List Tolerance D_good Limit_Tol} {
       }
    }
 }
+
+# Procedure to check result of nbshapes command
+proc checknbshapes { res nbshapes_expected_s count_locations message} {
+
+    upvar $res shape
+    if { ${count_locations} == 0 } {
+      set nb_info [nbshapes shape]
+    } else {
+      set nb_info [nbshapes shape -t]
+    }
+
+    set EntityList {VERTEX EDGE WIRE FACE SHELL SOLID COMPSOLID COMPOUND SHAPE}
+
+    foreach Entity ${EntityList} {
+       set expr_string "${Entity} +: +(\[-0-9.+eE\]+)"
+       if { [regexp "${expr_string}" ${nbshapes_expected_s} full nb_entity1] > 0 } {
+       	  if { [regexp "${expr_string}" ${nb_info} full nb_entity2] > 0 } {
+    	     if { ${nb_entity2} != ${nb_entity1} } {
+             	puts "Error : ${message} is WRONG because number of ${Entity} entities is ${nb_entity2} while ${nb_entity1} is expected"
+    	     } else {
+             	puts "OK : ${message} is GOOD because number of ${Entity} entities is equal to number of expected ${Entity} entities"
+    	     }
+       	  }
+       }
+    }
+}

tests/bugs/modalg_5/bug25398

@@ -0,0 +1,38 @@
+puts "========"
+puts "OCC25398"
+puts "========"
+puts ""
+######################################################
+# Provide shape proximity detector
+######################################################
+
+psphere s 10
+tclean s
+incmesh s 0.001
+
+pcone c 5 20 30
+tclean c
+incmesh c 0.001
+
+proximity s c -tol 0.01 -profile
+
+set nbshapes_expected "
+Number of shapes in shape
+ VERTEX    : 2
+ EDGE      : 3
+ WIRE      : 1
+ FACE      : 1
+ SHELL     : 0
+ SOLID     : 0
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 8
+"
+checknbshapes s_overlapped ${nbshapes_expected} 1 "Overlapped part of shere"
+checknbshapes c_overlapped ${nbshapes_expected} 1 "Overlapped part of shere"
+
+vinit
+vdisplay s_overlapped
+vdisplay c_overlapped
+vfit
+set only_screen 1