29470 - Point cloud on a face

idl/GEOM_Gen.idl

@@ -816,7 +816,8 @@ module GEOM
      *  \param theFace The referenced face.
      *  \return New GEOM_Object, containing the created point.
      */
-    GEOM_Object MakePointOnFace (in GEOM_Object theFace);
+    GEOM_Object MakePointOnFace (in GEOM_Object theFace,
+				 in long        theNumberOfPnts);
 
     /*!
      *  \brief Create a point, on two lines intersection.

src/GEOMAlgo/GEOMAlgo_AlgoTools.cxx

@@ -25,11 +25,17 @@
 
 #include <GEOMAlgo_AlgoTools.hxx>
 
+#include <Basics_OCCTVersion.hxx>
+
 #include <gp_Pnt.hxx>
 #include <gp_Pnt2d.hxx>
 #include <gp_Dir2d.hxx>
+#include <gp_Ax2.hxx>
 #include <Bnd_Box.hxx>
 
+#include <BRepAdaptor_Curve2d.hxx>
+#include <BRepTopAdaptor_FClass2d.hxx>
+
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
 #include <Geom2d_Line.hxx>
@@ -41,6 +47,7 @@
 #include <Geom2dAdaptor_Curve.hxx>
 #include <HatchGen_Domain.hxx>
 
+#include <GeomLib.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Surface.hxx>
 
@@ -49,6 +56,8 @@
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
 
+#include <GProp_GProps.hxx>
+
 #include <Poly_Triangulation.hxx>
 
 #include <TopAbs_Orientation.hxx>
@@ -59,20 +68,28 @@
 #include <TopoDS_Iterator.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Edge.hxx>
+#include <TopoDS_Compound.hxx>
 
+#include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 
 #include <BRep_Tool.hxx>
 #include <BRep_Builder.hxx>
+#include <BRepLib_MakeVertex.hxx>
 
 #include <BRepTools.hxx>
+#include <BRepTools_WireExplorer.hxx>
 #include <BRepBndLib.hxx>
 #include <BRepMesh_IncrementalMesh.hxx>
+#include <BRepGProp.hxx>
 
 #include <IntTools_Tools.hxx>
 
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
 #include <TopTools_ListOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopTools_DataMapOfShapeReal.hxx>
+#include <TColgp_SequenceOfPnt2d.hxx>
 
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
@@ -83,8 +100,12 @@
 #include <BOPTools_AlgoTools3D.hxx>
 #include <BOPTools_AlgoTools2D.hxx>
 
+#include <ShapeUpgrade_ShapeDivideArea.hxx>
+#include <ShapeUpgrade_UnifySameDomain.hxx>
+
 #include <GEOMAlgo_PassKeyShape.hxx>
 
+#include <algorithm>
 
 static
   void GetCount(const TopoDS_Shape& aS,
@@ -94,6 +115,30 @@ static
     TopTools_IndexedDataMapOfShapeShape& aMapSS,
     TopoDS_Shape& aSC);
 
+static Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
+                             const std::pair<TopoDS_Shape, Standard_Real>& theB);
+
+static Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
+                                const TopoDS_Face& theFace);
+
+static void CorrectShell (const TopoDS_Shape& theShell,
+                          const TopoDS_Face&  theFace);
+
+static gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
+                          const Standard_Boolean theIsNaturalRestrictions);
+
+static void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
+                                       const Standard_Real    theAverageArea,
+                                       const Standard_Boolean theIsToAddFaces,
+                                       Standard_Integer&      theNbExtremalFaces,
+                                       TopTools_MapOfShape&   theExtremalFaces,
+                                       const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
+                                       const TopTools_DataMapOfShapeReal& theFaceAreaMap,
+                                       const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
+                                       TopoDS_Shape&          theRes,
+                                       TopoDS_Shape&          theGlobalRes,
+                                       TopTools_MapOfShape&   theRemovedFaces);
+
 //=======================================================================
 //function : CopyShape
 //purpose  :
@@ -883,6 +928,7 @@ void GetCount(const TopoDS_Shape& aS,
     GetCount(aSx, iCnt);
   }
 }
+
 //=======================================================================
 //function : PntInFace
 //purpose  :
@@ -1004,3 +1050,592 @@ Standard_Integer GEOMAlgo_AlgoTools::PntInFace(const TopoDS_Face& aF,
   //
   return iErr;
 }
+
+//=======================================================================
+//function : PointCloudInFace
+//purpose  :
+//=======================================================================
+Standard_Integer GEOMAlgo_AlgoTools::PointCloudInFace(const TopoDS_Face& theFace,
+                                                      const int          theNbPnts,
+                                                      TopoDS_Compound&   theCompound)
+{
+#if OCC_VERSION_LARGE < 0x07050304
+  return -1;
+#else
+  ShapeUpgrade_ShapeDivideArea tool (theFace);
+  tool.SetSplittingByNumber (Standard_True);
+  tool.NbParts() = theNbPnts;
+  tool.Perform();
+  TopoDS_Shape res = tool.Result();
+
+  BRep_Builder aBB;
+  TopoDS_Compound aGlobalRes;
+  aBB.MakeCompound (aGlobalRes);
+
+  TopTools_IndexedMapOfShape aFaceMap;
+  TopExp::MapShapes (res, TopAbs_FACE, aFaceMap);
+  Standard_Integer aNbFaces = aFaceMap.Extent();
+
+  TopTools_IndexedDataMapOfShapeListOfShape aEFmap;
+  TopExp::MapShapesAndAncestors (res, TopAbs_EDGE, TopAbs_FACE, aEFmap);
+  
+  TopTools_MapOfShape aBiggestFaces, aSmallestFaces;
+  Standard_Boolean aUseTriangulation = Standard_True;
+  Standard_Boolean aSkipShared = Standard_False;
+  if (aNbFaces != theNbPnts)
+  {
+    Standard_Real aTotalArea = 0.;
+    std::vector<std::pair<TopoDS_Shape, Standard_Real> > aFacesAndAreas (aNbFaces);
+    for (Standard_Integer ii = 1; ii <= aNbFaces; ii++)
+    {
+      GProp_GProps aProps;
+      BRepGProp::SurfaceProperties (aFaceMap(ii), aProps, aSkipShared, aUseTriangulation);
+      Standard_Real anArea = aProps.Mass();
+      aTotalArea += anArea;
+      std::pair<TopoDS_Shape, Standard_Real> aFaceWithArea (aFaceMap(ii), anArea);
+      aFacesAndAreas[ii-1] = aFaceWithArea;
+    }
+    std::sort (aFacesAndAreas.begin(), aFacesAndAreas.end(), comp);
+
+    Standard_Real anAverageArea = aTotalArea / theNbPnts;
+
+    TopTools_DataMapOfShapeReal aFaceAreaMap;
+    for (Standard_Integer ii = 0; ii < aNbFaces; ii++)
+      aFaceAreaMap.Bind (aFacesAndAreas[ii].first, aFacesAndAreas[ii].second);
+    
+    TopTools_MapOfShape aRemovedFaces;
+    
+    if (aNbFaces < theNbPnts)
+    {
+      Standard_Integer aNbMissingFaces = theNbPnts - aNbFaces;
+      for (Standard_Integer ii = aNbFaces-1; ii > aNbFaces - aNbMissingFaces - 1; ii--)
+        aBiggestFaces.Add (aFacesAndAreas[ii].first);
+
+      ModifyFacesForGlobalResult (theFace, anAverageArea,
+                                  Standard_True, //to add faces
+                                  aNbMissingFaces, aBiggestFaces,
+                                  aFacesAndAreas, aFaceAreaMap, aEFmap,
+                                  res, aGlobalRes,
+                                  aRemovedFaces);
+    }
+    else //aNbFaces > theNbPnts
+    {
+      Standard_Integer aNbExcessFaces = aNbFaces - theNbPnts;
+      for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
+        aSmallestFaces.Add (aFacesAndAreas[ii].first);
+
+      TopTools_IndexedDataMapOfShapeListOfShape aVFmap;
+      TopExp::MapShapesAndAncestors (res, TopAbs_VERTEX, TopAbs_FACE, aVFmap);
+
+      //Remove smallest faces with free boundaries
+      for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
+      {
+        const TopoDS_Face& aFace = TopoDS::Face (aFacesAndAreas[ii].first);
+        Standard_Boolean anIsFreeBoundFound = Standard_False;
+        TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
+        for (; anExplo.More(); anExplo.Next())
+        {
+          const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
+          if (!BRep_Tool::Degenerated (anEdge) &&
+              aEFmap.FindFromKey(anEdge).Extent() < 2)
+          {
+            anIsFreeBoundFound = Standard_True;
+            break;
+          }
+        }
+        if (anIsFreeBoundFound)
+        {
+          Standard_Real aMaxArea = 0.;
+          for (anExplo.Init(aFace, TopAbs_VERTEX); anExplo.More(); anExplo.Next())
+          {
+            const TopoDS_Shape& aVertex = anExplo.Current();
+            const TopTools_ListOfShape& aFaceList = aVFmap.FindFromKey (aVertex);
+            TopTools_ListIteratorOfListOfShape anItl (aFaceList);
+            for (; anItl.More(); anItl.Next())
+            {
+              Standard_Real anArea = aFaceAreaMap (anItl.Value());
+              if (anArea > aMaxArea)
+                aMaxArea = anArea;
+            }
+          }
+          Standard_Real anArreaOfSmallestFace = aFaceAreaMap (aFace);
+          if (anArreaOfSmallestFace < aMaxArea / 16)
+          {
+            aBB.Remove (res, aFace);
+            aRemovedFaces.Add (aFace);
+          }
+        }
+      }
+
+      ModifyFacesForGlobalResult (theFace, anAverageArea,
+                                  Standard_False, //to decrease number of faces
+                                  aNbExcessFaces, aSmallestFaces,
+                                  aFacesAndAreas, aFaceAreaMap, aEFmap,
+                                  res, aGlobalRes,
+                                  aRemovedFaces);
+    }
+  }
+
+  aBB.Add (aGlobalRes, res);
+
+  aBB.MakeCompound (theCompound);
+  for (TopExp_Explorer aGlobalExplo(aGlobalRes, TopAbs_FACE); aGlobalExplo.More(); aGlobalExplo.Next())
+  {
+    const TopoDS_Face& aFace = TopoDS::Face (aGlobalExplo.Current());
+    Standard_Boolean anIsNaturalRestrictions = Standard_True;
+    TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
+    for (; anExplo.More(); anExplo.Next())
+    {
+      const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
+      if (BRep_Tool::Degenerated (anEdge))
+        continue;
+      if (!aEFmap.Contains(anEdge) ||
+          aEFmap.FindFromKey(anEdge).Extent() < 2)
+      {
+        anIsNaturalRestrictions = Standard_False;
+        break;
+      }
+    }
+
+    gp_Pnt aPnt = GetMidPnt2d (aFace, anIsNaturalRestrictions);
+    TopoDS_Vertex aVertex = BRepLib_MakeVertex (aPnt);
+    aBB.Add (theCompound, aVertex);
+  }
+
+  return 0;
+#endif
+}
+
+Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
+                      const std::pair<TopoDS_Shape, Standard_Real>& theB)
+{
+  return (theA.second < theB.second);
+}
+
+Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
+                         const TopoDS_Face& theFace)
+{
+  BRepAdaptor_Curve2d aBAcurve2d (theEdge, theFace);
+  gp_Pnt2d aP2d;
+  gp_Vec2d aVec;
+  aBAcurve2d.D1 (aBAcurve2d.FirstParameter(), aP2d, aVec);
+  return (Abs(aVec.Y()) > Abs(aVec.X()));
+}
+
+void CorrectShell (const TopoDS_Shape& theShell,
+                   const TopoDS_Face&  theFace)
+{
+  BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
+  GeomAbs_SurfaceType aType = aBAsurf.GetType();
+  if (aType <= GeomAbs_Torus) //elementary surfaces
+    return;
+
+  TopLoc_Location anInputLoc;
+  const Handle(Geom_Surface)& anInputSurf = BRep_Tool::Surface (theFace, anInputLoc);
+
+  BRep_Builder aBB;
+  
+  TopoDS_Iterator anIter (theShell);
+  for (; anIter.More(); anIter.Next())
+  {
+    const TopoDS_Face& aFace = TopoDS::Face (anIter.Value());
+    TopLoc_Location aLoc;
+    const Handle(Geom_Surface)& aSurf = BRep_Tool::Surface (aFace, aLoc);
+    if (aSurf == anInputSurf)
+      continue;
+
+    TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
+    for (; anExplo.More(); anExplo.Next())
+    {
+      const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
+      Standard_Real aFirst, aLast;
+      Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
+      aBB.UpdateEdge (anEdge, aPCurve, anInputSurf, anInputLoc, 0.);
+    }
+    Standard_Real aTol = BRep_Tool::Tolerance (aFace);
+    aBB.UpdateFace (aFace, anInputSurf, anInputLoc, aTol);
+  }
+}
+
+gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
+                   const Standard_Boolean theIsNaturalRestrictions)
+{
+  gp_Pnt aResPnt;
+  
+  if (theIsNaturalRestrictions)
+  {
+    BRepAdaptor_Surface aBAsurf (theFace);
+    Standard_Real aUmin, aUmax, aVmin, aVmax;
+    aUmin = aBAsurf.FirstUParameter();
+    aUmax = aBAsurf.LastUParameter();
+    aVmin = aBAsurf.FirstVParameter();
+    aVmax = aBAsurf.LastVParameter();
+    aResPnt = aBAsurf.Value ((aUmin + aUmax)/2, (aVmin + aVmax)/2);
+  }
+  else
+  {
+    const Standard_Integer aNbSamples = 4;
+    TopoDS_Wire aWire = BRepTools::OuterWire (theFace);
+    TopTools_IndexedMapOfShape aEmap;
+    TopExp::MapShapes (aWire, TopAbs_EDGE, aEmap);
+    Standard_Integer aNbPointsOnContour = aNbSamples * aEmap.Extent();
+    TColgp_Array1OfPnt anArray (1, aNbPointsOnContour);
+    
+    BRepTools_WireExplorer aWexp (aWire, theFace);
+    Standard_Integer anInd = 0;
+    for (; aWexp.More(); aWexp.Next())
+    {
+      const TopoDS_Edge& anEdge = aWexp.Current();
+      BRepAdaptor_Curve2d aBAcurve2d (anEdge, theFace);
+      Standard_Real aDelta = (aBAcurve2d.LastParameter() - aBAcurve2d.FirstParameter())/aNbSamples;
+      for (Standard_Integer ii = 0; ii < aNbSamples; ii++)
+      {
+        Standard_Real aParam = aBAcurve2d.FirstParameter() + ii * aDelta;
+        gp_Pnt2d aP2d = aBAcurve2d.Value (aParam);
+        gp_Pnt aPnt (aP2d.X(), aP2d.Y(), 0.);
+        anArray (++anInd) = aPnt;
+      }
+    }
+    
+    gp_Ax2 anAxis;
+    Standard_Boolean anIsSingular;
+    GeomLib::AxeOfInertia (anArray, anAxis, anIsSingular);
+    gp_Pnt aBaryCentre = anAxis.Location();
+    gp_Pnt2d aCentre2d (aBaryCentre.X(), aBaryCentre.Y());
+    BRepTopAdaptor_FClass2d aClassifier (theFace, Precision::Confusion());
+    BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
+    
+    TopAbs_State aStatus = aClassifier.Perform (aCentre2d);
+    gp_Pnt2d aP2d = aCentre2d;
+    Standard_Integer anIndVertex = 0;
+    const Standard_Integer aNbIter = 10;
+    while (aStatus != TopAbs_IN && anIndVertex < aNbPointsOnContour)
+    {
+      gp_Pnt aVertexPnt = anArray (anIndVertex+1);
+      gp_Pnt2d aVertexP2d (aVertexPnt.X(), aVertexPnt.Y());
+      TColgp_SequenceOfPnt2d aPseq;
+      aPseq.Append (aCentre2d);
+      aPseq.Append (aVertexP2d);
+      for (Standard_Integer ii = 1; ii <= aNbIter; ii++)
+      {
+        for (Standard_Integer jj = 1; jj < aPseq.Length(); jj++)
+        {
+          aP2d.SetXY ((aPseq(jj).XY() + aPseq(jj+1).XY())/2);
+          aStatus = aClassifier.Perform (aP2d);
+          if (aStatus == TopAbs_IN)
+            break;
+          else
+          {
+            aPseq.InsertAfter (jj, aP2d);
+            jj++;
+          }
+        }
+        if (aStatus == TopAbs_IN)
+          break;
+      }
+      anIndVertex += aNbSamples;
+    }
+    aResPnt = aBAsurf.Value (aP2d.X(), aP2d.Y());
+  } //case of complex boundaries
+
+  return aResPnt;
+}
+
+void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
+                                const Standard_Real    theAverageArea,
+                                const Standard_Boolean theIsToAddFaces,
+                                Standard_Integer&      theNbExtremalFaces,
+                                TopTools_MapOfShape&   theExtremalFaces,
+                                const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
+                                const TopTools_DataMapOfShapeReal& theFaceAreaMap,
+                                const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
+                                TopoDS_Shape&          theRes,
+                                TopoDS_Shape&          theGlobalRes,
+                                TopTools_MapOfShape&   theRemovedFaces)
+{
+  BRep_Builder aBB;
+  const Standard_Integer aNbFaces = (Standard_Integer) theFacesAndAreas.size();
+
+  const Standard_Integer aDiff = theNbExtremalFaces - theRemovedFaces.Extent();
+
+  Standard_Integer aSum = 0;
+  while (aSum < aDiff) //global loop
+  {
+    Standard_Integer aNbFacesDone = 0, aNbFacesInTape = 0;
+    TopoDS_Face aStartFace;
+    
+    Standard_Integer aStartIndex = (theIsToAddFaces)? aNbFaces-1 : 0;
+    Standard_Integer anEndIndex  = (theIsToAddFaces)? 0 : aNbFaces-1;
+    Standard_Integer aStep = (theIsToAddFaces)? -1 : 1;
+    
+    for (Standard_Integer ii = aStartIndex; ii != anEndIndex; ii += aStep)
+    {
+      const TopoDS_Face& aFace = TopoDS::Face (theFacesAndAreas[ii].first);
+      if (!theRemovedFaces.Contains(aFace))
+      {
+        aStartFace = aFace;
+        break;
+      }
+    }
+    if (aStartFace.IsNull())
+      break;
+
+    theRemovedFaces.Add (aStartFace);
+    
+    TopoDS_Edge aCommonEdge;
+    TopoDS_Face aNextFace;
+    Standard_Real anExtremalArea = (theIsToAddFaces)? 0. : Precision::Infinite();
+    for (TopExp_Explorer anExplo(aStartFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
+    {
+      const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
+      const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
+      TopTools_ListIteratorOfListOfShape anItl (aFaceList);
+      for (; anItl.More(); anItl.Next())
+      {
+        const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
+        if (aFace.IsSame (aStartFace) ||
+            theRemovedFaces.Contains(aFace))
+          continue;
+        Standard_Real anArea = theFaceAreaMap(aFace);
+        Standard_Boolean anIsToExchange = (theIsToAddFaces)? (anArea > anExtremalArea) : (anArea < anExtremalArea);
+        if (anIsToExchange)
+        {
+          anExtremalArea = anArea;
+          aCommonEdge = anEdge;
+          aNextFace = aFace;
+        }
+      }
+    }
+    if (aCommonEdge.IsNull()) //all adjacent faces are already removed
+    {
+      theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
+      theNbExtremalFaces++;
+      continue;
+    }
+
+    //Start filling the shell
+    aBB.Remove (theRes, aStartFace);
+    aNbFacesDone++;
+    TopoDS_Shell aShell;
+    aBB.MakeShell (aShell);
+    Standard_Real anAreaOfTape = 0.;
+    aBB.Add (aShell, aStartFace);
+    aNbFacesInTape++;
+    anAreaOfTape += theFaceAreaMap (aStartFace);
+    
+    Standard_Boolean anIsUiso = IsUiso (aCommonEdge, aStartFace);
+    //Find another faces on this level
+    TopoDS_Face aCurrentFace = aNextFace;
+    TopoDS_Edge aCurrentEdge = aCommonEdge;
+    Standard_Boolean anIsFirstDirection = Standard_True;
+    aBB.Remove (theRes, aCurrentFace);
+    theRemovedFaces.Add (aCurrentFace);
+    if (theExtremalFaces.Contains (aCurrentFace))
+    {
+      aNbFacesDone++;
+    }
+    aBB.Add (aShell, aCurrentFace);
+    aNbFacesInTape++;
+    anAreaOfTape += theFaceAreaMap (aCurrentFace);
+    Standard_Boolean anIsRound = Standard_False;
+    for (;;) //local loop
+    {
+      TopoDS_Edge aNextEdge;
+      for (TopExp_Explorer anExplo(aCurrentFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
+      {
+        const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
+        if (anEdge.IsSame (aCurrentEdge))
+          continue;
+        const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
+        TopTools_ListIteratorOfListOfShape anItl (aFaceList);
+        for (; anItl.More(); anItl.Next())
+        {
+          const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
+          if (aFace.IsSame (aCurrentFace))
+            continue;
+          if (aFace.IsSame (aStartFace))
+          {
+            anIsRound = Standard_True;
+            break;
+          }
+          if (theRemovedFaces.Contains(aFace))
+            continue;
+          if (anIsUiso == IsUiso (anEdge, aFace))
+          {
+            aNextEdge = anEdge;
+            aNextFace = aFace;
+            break;
+          }
+        }
+        if (anIsRound || !aNextEdge.IsNull())
+          break;
+      }
+      if (anIsRound) //round tape: returned to start face
+        break;
+      if (aNextEdge.IsNull())
+      {
+        if (anIsFirstDirection)
+        {
+          aCurrentFace = aStartFace;
+          aCurrentEdge = aCommonEdge;
+          anIsFirstDirection = Standard_False;
+          continue;
+        }
+        else
+          break;
+      }
+      
+      aBB.Add (aShell, aNextFace);
+      aNbFacesInTape++;
+      anAreaOfTape += theFaceAreaMap (aNextFace);
+      aBB.Remove (theRes, aNextFace);
+      theRemovedFaces.Add (aNextFace);
+      if (theExtremalFaces.Contains (aNextFace))
+      {
+        aNbFacesDone++;
+      }
+      aCurrentEdge = aNextEdge;
+      aNextEdge.Nullify();
+      aCurrentFace = aNextFace;
+    } //end of local loop
+    
+    //Tape is formed
+    Standard_Integer aNumberToSplit = (theIsToAddFaces)? aNbFacesInTape + aNbFacesDone : aNbFacesInTape - aNbFacesDone;
+    if (!theIsToAddFaces && aNbFacesDone > 1)
+    {
+      Standard_Integer aRealNumberToSplit = (aNumberToSplit > 0)? aNumberToSplit : 1;
+      Standard_Real anAverageAreaInTape = anAreaOfTape / aRealNumberToSplit;
+      if (anAverageAreaInTape > theAverageArea)
+      {
+        Standard_Integer aNewNumberToSplit = RealToInt(round(anAreaOfTape / theAverageArea));
+        if (aNewNumberToSplit < aNbFacesInTape)
+        {
+          Standard_Integer aNumberToIncrease = aNewNumberToSplit - aNumberToSplit;
+          for (Standard_Integer jj = theNbExtremalFaces; jj < theNbExtremalFaces + aNumberToIncrease; jj++)
+            theExtremalFaces.Add (theFacesAndAreas[jj].first);
+          theNbExtremalFaces += aNumberToIncrease;
+          aNumberToSplit = aNewNumberToSplit;
+        }
+      }
+    }
+    if (anIsRound && aNumberToSplit <= 1)
+    {
+      Standard_Integer aNumberToIncrease = 3 - aNumberToSplit;
+      for (Standard_Integer jj = theNbExtremalFaces; jj < theNbExtremalFaces + aNumberToIncrease; jj++)
+        theExtremalFaces.Add (theFacesAndAreas[jj].first);
+      theNbExtremalFaces += aNumberToIncrease;
+      aNumberToSplit = 3;
+    }
+    CorrectShell (aShell, theInputFace);
+    ShapeUpgrade_UnifySameDomain aUnifier;
+    aUnifier.Initialize (aShell, Standard_True, Standard_True);
+    aUnifier.Build();
+    TopoDS_Shape aUnifiedShape = aUnifier.Shape();
+    //Splitting
+    TopoDS_Shape aLocalResult = aUnifiedShape;
+    Standard_Integer aNbFacesInLocalResult;
+    if (aNumberToSplit > 1)
+    {
+#if OCC_VERSION_LARGE < 0x07050304
+      aNbFacesInLocalResult = 0;
+#else
+      ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
+      aLocalTool.SetSplittingByNumber (Standard_True);
+      aLocalTool.MaxArea() = -1;
+      if (anIsUiso)
+        aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
+      else
+        aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
+      aLocalTool.Perform();
+      aLocalResult = aLocalTool.Result();
+      aNbFacesInLocalResult = aNumberToSplit;
+#endif
+    }
+    else
+    {
+      aNbFacesInLocalResult = 1;
+      if (aNumberToSplit == 0)
+      {
+        theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
+        theNbExtremalFaces++;
+      }
+    }
+    aBB.Add (theGlobalRes, aLocalResult);
+
+    aSum += Abs(aNbFacesInTape - aNbFacesInLocalResult);
+  } //end of global loop
+
+  //Second global loop
+  TopoDS_Compound aSecondComp;
+  aBB.MakeCompound (aSecondComp);
+  while (aSum < aDiff)
+  {
+    TopoDS_Shape aMaxShell;
+    Standard_Integer aMaxNbFaces = 0;
+    TopoDS_Iterator anIter (theGlobalRes);
+    for (; anIter.More(); anIter.Next())
+    {
+      const TopoDS_Shape& aShell = anIter.Value();
+      TopTools_IndexedMapOfShape aFaceMap;
+      TopExp::MapShapes (aShell, TopAbs_FACE, aFaceMap);
+      if (aFaceMap.Extent() > aMaxNbFaces)
+      {
+        aMaxNbFaces = aFaceMap.Extent();
+        aMaxShell = aShell;
+      }
+    }
+
+    if (aMaxNbFaces == 1)
+      break;
+    
+    aBB.Remove (theGlobalRes, aMaxShell);
+    //Find iso
+    Standard_Boolean anIsUiso = Standard_True;
+    TopTools_IndexedDataMapOfShapeListOfShape aLocalEFmap;
+    TopExp::MapShapesAndAncestors (aMaxShell, TopAbs_EDGE, TopAbs_FACE, aLocalEFmap);
+    for (Standard_Integer jj = 1; jj <= aLocalEFmap.Extent(); jj++)
+    {
+      const TopoDS_Edge& anEdge = TopoDS::Edge (aLocalEFmap.FindKey(jj));
+      const TopTools_ListOfShape& aFaceList = aLocalEFmap(jj);
+      if (aFaceList.Extent() == 2)
+      {
+        const TopoDS_Face& aFace = TopoDS::Face (aFaceList.First());
+        anIsUiso = IsUiso (anEdge, aFace);
+        break;
+      }
+    }
+    CorrectShell (aMaxShell, theInputFace);
+    ShapeUpgrade_UnifySameDomain aUnifier;
+    aUnifier.Initialize (aMaxShell, Standard_True, Standard_True);
+    aUnifier.Build();
+    TopoDS_Shape aUnifiedShape = aUnifier.Shape();
+    TopoDS_Shape aLocalResult = aUnifiedShape;
+    
+    Standard_Integer aNumberToSplit = (theIsToAddFaces)? aMaxNbFaces + (aDiff-aSum) : aMaxNbFaces - (aDiff-aSum);
+    if (aNumberToSplit > 1)
+    {
+#if OCC_VERSION_LARGE < 0x07050304
+      aNumberToSplit = 1;
+#else
+      ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
+      aLocalTool.SetSplittingByNumber (Standard_True);
+      aLocalTool.MaxArea() = -1;
+      if (anIsUiso)
+        aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
+      else
+        aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
+      aLocalTool.Perform();
+      aLocalResult = aLocalTool.Result();
+#endif
+    }
+    else
+      aNumberToSplit = 1;
+
+    aBB.Add (aSecondComp, aLocalResult);
+    
+    if (theIsToAddFaces)
+      break;
+    aSum += aMaxNbFaces - aNumberToSplit;
+  }
+  aBB.Add (theGlobalRes, aSecondComp);
+}

src/GEOMAlgo/GEOMAlgo_AlgoTools.hxx

@@ -42,6 +42,7 @@
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Shape.hxx>
 #include <TopoDS_Edge.hxx>
+#include <TopoDS_Compound.hxx>
 
 #include <TopTools_ListOfShape.hxx>
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
@@ -75,6 +76,13 @@ class GEOMAlgo_AlgoTools  {
 				       gp_Pnt& theP,
 				       gp_Pnt2d& theP2D) ;
 
+  //! Computes a set of points inside the face <theF>. <br>
+  //!          Returns 0 in case of success. <br>
+  Standard_EXPORT
+    static  Standard_Integer PointCloudInFace(const TopoDS_Face& theF,
+                                              const int          theNbPnts,
+                                              TopoDS_Compound&   theCompound) ;
+
   Standard_EXPORT
     static  Standard_Boolean IsCompositeShape(const TopoDS_Shape& aS) ;
 

src/GEOMImpl/GEOMImpl_IBasicOperations.cxx

@@ -184,6 +184,7 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::makePointOnGeom
                      double              theParam2,
                      double              theParam3,
                      const PointLocation theLocation,
+                     int                 theNumberOfPnts,
                      const bool          takeOrientationIntoAccount,
                      Handle(GEOM_Object) theRefPoint)
 {
@@ -250,6 +251,7 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::makePointOnGeom
       break;
     case PointOn_Face:
       aPI.SetSurface(aRefFunction);
+      aPI.SetNumberOfPoints(theNumberOfPnts);
       break;
     default: break;
     }
@@ -296,7 +298,7 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::makePointOnGeom
       break;
     case PointOn_Face:
       GEOM::TPythonDump(aFunction) << aPoint << " = geompy.MakeVertexInsideFace("
-                                   << theGeomObj << ")";
+                                   << theGeomObj << ", " << theNumberOfPnts << ")";
       break;
     default: break;
     }
@@ -316,7 +318,7 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::MakePointOnCurve
                              bool                takeOrientationIntoAccount)
 {
   return makePointOnGeom(theCurve, theParameter, 0.0, 0.0, PointOn_CurveByParam,
-                         takeOrientationIntoAccount);
+                         1, takeOrientationIntoAccount);
 }
 
 //=============================================================================
@@ -344,7 +346,7 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::MakePointOnCurveByLength
                      Handle(GEOM_Object) theStartPoint)
 {
   return makePointOnGeom(theCurve, theLength, 0.0, 0.0, PointOn_CurveByLength,
-                         false, theStartPoint);
+                         1, false, theStartPoint);
 }
 
 //=============================================================================
@@ -379,9 +381,10 @@ Handle(GEOM_Object) GEOMImpl_IBasicOperations::MakePointOnSurfaceByCoord
  *  MakePointOnFace
  */
 //=============================================================================
-Handle(GEOM_Object) GEOMImpl_IBasicOperations::MakePointOnFace (Handle(GEOM_Object) theFace)
+Handle(GEOM_Object) GEOMImpl_IBasicOperations::MakePointOnFace (Handle(GEOM_Object) theFace,
+                                                                int                 theNumberOfPnts)
 {
-  return makePointOnGeom(theFace, 0., 0., 0., PointOn_Face);
+  return makePointOnGeom(theFace, 0., 0., 0., PointOn_Face, theNumberOfPnts);
 }
 
 //=============================================================================

src/GEOMImpl/GEOMImpl_IBasicOperations.hxx

@@ -66,7 +66,8 @@ class GEOMImpl_IBasicOperations : public GEOM_IOperations {
                                                                  double theYParam,
                                                                  double theZParam);
 
-  Standard_EXPORT Handle(GEOM_Object) MakePointOnFace (Handle(GEOM_Object) theFace);
+  Standard_EXPORT Handle(GEOM_Object) MakePointOnFace (Handle(GEOM_Object) theFace,
+                                                       int                 theNumberOfPnts);
 
   // Vector
   Standard_EXPORT Handle(GEOM_Object) MakeVectorDXDYDZ (double theDX, double theDY, double theDZ);
@@ -142,6 +143,7 @@ class GEOMImpl_IBasicOperations : public GEOM_IOperations {
                       double              theParam2,
                       double              theParam3,
                       const PointLocation theLocation,
+                      int                 theNumberOfPnts = 1,
                       const bool          takeOrientationIntoAccount = false,
                       Handle(GEOM_Object) theRefPoint = 0);
 };

src/GEOMImpl/GEOMImpl_IPoint.hxx

@@ -42,6 +42,8 @@
 
 #define ARG_USE_ORIENTATION 12
 
+#define ARG_NBPNTS 13
+
 class GEOMImpl_IPoint
 {
  public:
@@ -72,12 +74,14 @@ class GEOMImpl_IPoint
 
   void SetParameter(double theParam) { _func->SetReal(ARG_PARAM, theParam); }
   void SetParameter2(double theParam) { _func->SetReal(ARG_PARAM2, theParam); }
+  void SetNumberOfPoints(int theNumberOfPnts) { _func->SetInteger(ARG_NBPNTS, theNumberOfPnts); }
   void SetLength(double theLength) { _func->SetReal(ARG_LENGTH, theLength); }
   void SetTakeOrientationIntoAccount(bool takeOrientationIntoAccount)
         { _func->SetInteger(ARG_USE_ORIENTATION, takeOrientationIntoAccount); }
 
   double GetParameter() { return _func->GetReal(ARG_PARAM); }
   double GetParameter2() { return _func->GetReal(ARG_PARAM2); }
+  int    GetNumberOfPoints() { return _func->GetInteger(ARG_NBPNTS); }
   double GetLength() { return _func->GetReal(ARG_LENGTH); }
   bool   GetTakeOrientationIntoAccount() { return _func->GetInteger(ARG_USE_ORIENTATION); }
 

src/GEOMImpl/GEOMImpl_PointDriver.cxx

@@ -22,6 +22,8 @@
 
 #include <Standard_Stream.hxx>
 
+#include <Basics_OCCTVersion.hxx>
+
 #include <GEOMImpl_PointDriver.hxx>
 #include <GEOMImpl_IPoint.hxx>
 #include <GEOMImpl_Types.hxx>
@@ -57,6 +59,7 @@
 #include <Precision.hxx>
 
 #include <Standard_NullObject.hxx>
+#include <Standard_NotImplemented.hxx>
 
 //=======================================================================
 //function : GetID
@@ -301,13 +304,31 @@ Standard_Integer GEOMImpl_PointDriver::Execute(Handle(TFunction_Logbook)& log) c
   else if (aType == POINT_FACE_ANY) {
     Handle(GEOM_Function) aRefFunc = aPI.GetSurface();
     TopoDS_Shape aRefShape = aRefFunc->GetValue();
+    int aNbPnts = aPI.GetNumberOfPoints();
+    if (aNbPnts < 1) {
+      Standard_TypeMismatch::Raise
+        ("Point On Surface creation aborted : number of points is zero or negative");
+    }
     if (aRefShape.ShapeType() != TopAbs_FACE) {
       Standard_TypeMismatch::Raise
         ("Point On Surface creation aborted : surface shape is not a face");
     }
     TopoDS_Face F = TopoDS::Face(aRefShape);
-    gp_Pnt2d aP2d;
+    if (aNbPnts == 1)
-    GEOMAlgo_AlgoTools::PntInFace(F, aPnt, aP2d);
+    {
+      gp_Pnt2d aP2d;
+      GEOMAlgo_AlgoTools::PntInFace(F, aPnt, aP2d);
+    }
+    else
+    {
+#if OCC_VERSION_LARGE < 0x07050304
+      Standard_NotImplemented::Raise("Point cloud creation aborted. Improper OCCT version: please, use OCCT 7.5.3p4 or newer.");
+#else
+      if (GEOMAlgo_AlgoTools::PointCloudInFace(F, aNbPnts, aCompound) < 0)
+        Standard_ConstructionError::Raise("Point cloud creation aborted : algorithm failed");
+      retCompound = true;
+#endif
+    }
   }
   else if (aType == POINT_LINES_INTERSECTION) {
     Handle(GEOM_Function) aRef1 = aPI.GetLine1();

src/GEOM_I/GEOM_IBasicOperations_i.cc

@@ -284,7 +284,8 @@ GEOM::GEOM_Object_ptr GEOM_IBasicOperations_i::MakePointOnSurfaceByCoord
  *  MakePointOnFace
  */
 //=============================================================================
-GEOM::GEOM_Object_ptr GEOM_IBasicOperations_i::MakePointOnFace (GEOM::GEOM_Object_ptr theFace)
+GEOM::GEOM_Object_ptr GEOM_IBasicOperations_i::MakePointOnFace (GEOM::GEOM_Object_ptr theFace,
+                                                                CORBA::Long           theNumberOfPnts)
 {
   GEOM::GEOM_Object_var aGEOMObject;
 
@@ -296,7 +297,7 @@ GEOM::GEOM_Object_ptr GEOM_IBasicOperations_i::MakePointOnFace (GEOM::GEOM_Objec
   if (aReference.IsNull()) return aGEOMObject._retn();
 
   //Create the point
-  Handle(::GEOM_Object) anObject = GetOperations()->MakePointOnFace(aReference);
+  Handle(::GEOM_Object) anObject = GetOperations()->MakePointOnFace(aReference, theNumberOfPnts);
   if (!GetOperations()->IsDone() || anObject.IsNull())
     return aGEOMObject._retn();
 

src/GEOM_I/GEOM_IBasicOperations_i.hh

@@ -75,7 +75,8 @@ class GEOM_I_EXPORT GEOM_IBasicOperations_i :
                                                     CORBA::Double theYParameter,
                                                     CORBA::Double theZParameter);
 
-   GEOM::GEOM_Object_ptr MakePointOnFace (GEOM::GEOM_Object_ptr theFace);
+   GEOM::GEOM_Object_ptr MakePointOnFace (GEOM::GEOM_Object_ptr theFace,
+                                          CORBA::Long           theNumberOfPnts);
 
    GEOM::GEOM_Object_ptr MakePointOnLinesIntersection (GEOM::GEOM_Object_ptr theLine1,
 						       GEOM::GEOM_Object_ptr theLine2);

src/GEOM_SWIG/geomBuilder.py

@@ -1441,7 +1441,7 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
         #
         #  @ref swig_MakeVertexInsideFace "Example"
         @ManageTransactions("BasicOp")
-        def MakeVertexInsideFace (self, theFace, theName=None):
+        def MakeVertexInsideFace (self, theFace, theNumberOfPnts=1, theName=None):
             """
             Create a point, which lays on the given face.
             The point will lay in arbitrary place of the face.
@@ -1451,6 +1451,7 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
 
             Parameters:
                 theFace The referenced face.
+                theNumberOfPnts The number of points we want to get, 1 by default.
                 theName Object name; when specified, this parameter is used
                         for result publication in the study. Otherwise, if automatic
                         publication is switched on, default value is used for result name.
@@ -1462,7 +1463,7 @@ class geomBuilder(GEOM._objref_GEOM_Gen):
                 p_on_face = geompy.MakeVertexInsideFace(Face)
             """
             # Example: see GEOM_TestAll.py
-            anObj = self.BasicOp.MakePointOnFace(theFace)
+            anObj = self.BasicOp.MakePointOnFace(theFace, theNumberOfPnts)
             RaiseIfFailed("MakeVertexInsideFace", self.BasicOp)
             self._autoPublish(anObj, theName, "vertex")
             return anObj

test/test_point_cloud_on_face.py

@@ -0,0 +1,44 @@
+# Point Cloud on Face
+
+import math
+import salome
+salome.salome_init_without_session()
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+
+# create spherical face
+Sphere = geompy.MakeSphereR(10, "Sphere")
+[Face] = geompy.ExtractShapes(Sphere, geompy.ShapeType["FACE"], True, "Sphere_face")
+
+# make a cloud of 100 points on the spherical face
+try:
+    CompoundOfVertices = geompy.MakeVertexInsideFace(Face, 100, "CompoundOfVertices")
+except Exception as err:
+    print(err)
+    # this test should not fail in case of "Improper OCCT version"
+    assert("Improper OCCT version" in str(err))
+else:
+    # check result
+    assert(geompy.NumberOfSubShapes(CompoundOfVertices, geompy.ShapeType["VERTEX"]) == 100)
+
+# test point cloud on a "Horse saddle"
+OX = geompy.MakeVectorDXDYDZ(1, 0, 0, 'OX')
+OY = geompy.MakeVectorDXDYDZ(0, 1, 0, 'OY')
+[Edge_1,Edge_2,Edge_3] = geompy.ExtractShapes(Sphere, geompy.ShapeType["EDGE"], True, "Edge")
+Rotation_1 = geompy.MakeRotation(Edge_3, OX, 90*math.pi/180.0, 'Rotation_1')
+Rotation_2 = geompy.MakeRotation(Rotation_1, OY, 180*math.pi/180.0, 'Rotation_2')
+Translation_1 = geompy.MakeTranslation(Rotation_2, 200, 0, 0, 'Translation_1')
+Translation_2 = geompy.MakeTranslation(Edge_3, 100, 100, 0, 'Translation_2')
+Translation_3 = geompy.MakeTranslation(Translation_2, 0, -200, 0, 'Translation_3')
+Filling_1 = geompy.MakeFilling([Translation_2, Edge_3, Translation_3])
+geompy.addToStudy(Filling_1, 'Filling_1')
+
+try:
+    PointCloudOnFilling = geompy.MakeVertexInsideFace(Filling_1, 30, "PointCloudOnFilling")
+except Exception as err:
+    print(err)
+    # this test should not fail in case of "Improper OCCT version"
+    assert("Improper OCCT version" in str(err))
+else:
+    assert(geompy.NumberOfSubShapes(PointCloudOnFilling, geompy.ShapeType["VERTEX"]) == 30)

test/tests.set

@@ -19,4 +19,5 @@
 
 SET(ALL_TESTS
   test_perf_01.py
+  test_point_cloud_on_face.py
   )