src/GeomAlgoAPI/GeomAlgoAPI_PointCloudOnFace.cpp

   1 // Copyright (C) 2014-2023  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "GeomAlgoAPI_PointCloudOnFace.h"
  21
  22 #include <TopoDS_Shape.hxx>
  23 #include <TopoDS_Face.hxx>
  24 #include <TopoDS.hxx>
  25 #include <BRep_Tool.hxx>
  26 #include <BRep_Builder.hxx>
  27 #include <BRepTools.hxx>
  28 #include <BRepTools_WireExplorer.hxx>
  29 #include <BRepTopAdaptor_FClass2d.hxx>
  30 #include <BRepLib_MakeVertex.hxx>
  31 #include <GeomLib.hxx>
  32 #include <gp_Ax3.hxx>
  33 #include <Geom_Plane.hxx>
  34 #include <gp_Pln.hxx>
  35 #include <TopoDS_Iterator.hxx>
  36 #include <TopExp_Explorer.hxx>
  37 #include <TopExp.hxx>
  38 #include <GProp_GProps.hxx>
  39 #include <BRepGProp.hxx>
  40 #include <ShapeAnalysis.hxx>
  41 #include <ShapeAnalysis_Surface.hxx>
  42 #include <ShapeUpgrade_UnifySameDomain.hxx>
  43 #include <ShapeUpgrade_ShapeDivideArea.hxx>
  44 #include <BRepAdaptor_Surface.hxx>
  45 #include <BRepAdaptor_Curve2d.hxx>
  46 #include <BRepBndLib.hxx>
  47 #include <GeomLProp_SLProps.hxx>
  48 #include <BRepBuilderAPI_MakeEdge.hxx>
  49 #include <TopTools_DataMapOfShapeReal.hxx>
  50
  51 #include <Standard_Version.hxx>
  52 // code from KERNEL_SRC/src/Basics/Basics_OCCTVersion.hxx
  53 #ifdef OCC_VERSION_SERVICEPACK
  54 #  define OCC_VERSION_LARGE (OCC_VERSION_MAJOR << 24 | OCC_VERSION_MINOR << 16 | OCC_VERSION_MAINTENANCE << 8 | OCC_VERSION_SERVICEPACK)
  55 #else
  56 #  ifdef OCC_VERSION_DEVELOPMENT
  57 #    define OCC_VERSION_LARGE ((OCC_VERSION_MAJOR << 24 | OCC_VERSION_MINOR << 16 | OCC_VERSION_MAINTENANCE << 8)-1)
  58 #  else
  59 #    define OCC_VERSION_LARGE (OCC_VERSION_MAJOR << 24 | OCC_VERSION_MINOR << 16 | OCC_VERSION_MAINTENANCE << 8)
  60 #  endif
  61 #endif
  62
  63 #include <algorithm>
  64
  65 Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
  66                       const std::pair<TopoDS_Shape, Standard_Real>& theB)
  67 {
  68   return (theA.second < theB.second);
  69 }
  70
  71 Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
  72                          const TopoDS_Face& theFace)
  73 {
  74   BRepAdaptor_Curve2d aBAcurve2d (theEdge, theFace);
  75   gp_Pnt2d aP2d;
  76   gp_Vec2d aVec;
  77   aBAcurve2d.D1 (aBAcurve2d.FirstParameter(), aP2d, aVec);
  78   return (Abs(aVec.Y()) > Abs(aVec.X()));
  79 }
  80
  81 void CorrectShell (const TopoDS_Shape& theShell,
  82                    const TopoDS_Face&  theFace)
  83 {
  84   BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
  85   GeomAbs_SurfaceType aType = aBAsurf.GetType();
  86   if (aType <= GeomAbs_Torus) //elementary surfaces
  87     return;
  88
  89   TopLoc_Location anInputLoc;
  90   const Handle(Geom_Surface)& anInputSurf = BRep_Tool::Surface (theFace, anInputLoc);
  91
  92   BRep_Builder aBB;
  93
  94   TopoDS_Iterator anIter (theShell);
  95   for (; anIter.More(); anIter.Next())
  96   {
  97     const TopoDS_Face& aFace = TopoDS::Face (anIter.Value());
  98     TopLoc_Location aLoc;
  99     const Handle(Geom_Surface)& aSurf = BRep_Tool::Surface (aFace, aLoc);
 100     if (aSurf == anInputSurf)
 101       continue;
 102
 103     TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
 104     for (; anExplo.More(); anExplo.Next())
 105     {
 106       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
 107       Standard_Real aFirst, aLast;
 108       Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
 109       aBB.UpdateEdge (anEdge, aPCurve, anInputSurf, anInputLoc, 0.);
 110     }
 111     Standard_Real aTol = BRep_Tool::Tolerance (aFace);
 112     aBB.UpdateFace (aFace, anInputSurf, anInputLoc, aTol);
 113   }
 114 }
 115
 116 gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
 117                    const Standard_Boolean theIsNaturalRestrictions)
 118 {
 119   gp_Pnt aResPnt;
 120
 121   if (theIsNaturalRestrictions)
 122   {
 123     BRepAdaptor_Surface aBAsurf (theFace);
 124     Standard_Real aUmin, aUmax, aVmin, aVmax;
 125     aUmin = aBAsurf.FirstUParameter();
 126     aUmax = aBAsurf.LastUParameter();
 127     aVmin = aBAsurf.FirstVParameter();
 128     aVmax = aBAsurf.LastVParameter();
 129     aResPnt = aBAsurf.Value ((aUmin + aUmax)/2, (aVmin + aVmax)/2);
 130   }
 131   else
 132   {
 133     const Standard_Integer aNbSamples = 4;
 134     TopoDS_Wire aWire = BRepTools::OuterWire (theFace);
 135     TopTools_IndexedMapOfShape aEmap;
 136     TopExp::MapShapes (aWire, TopAbs_EDGE, aEmap);
 137     Standard_Integer aNbPointsOnContour = aNbSamples * aEmap.Extent();
 138     TColgp_Array1OfPnt anArray (1, aNbPointsOnContour);
 139
 140     BRepTools_WireExplorer aWexp (aWire, theFace);
 141     Standard_Integer anInd = 0;
 142     TopTools_MapOfShape aUsedEmap;
 143     for (; aWexp.More(); aWexp.Next())
 144     {
 145       const TopoDS_Edge& anEdge = aWexp.Current();
 146       if (!aUsedEmap.Add(anEdge)) continue;
 147       BRepAdaptor_Curve2d aBAcurve2d (anEdge, theFace);
 148       Standard_Real aDelta = (aBAcurve2d.LastParameter() - aBAcurve2d.FirstParameter())/aNbSamples;
 149       for (Standard_Integer ii = 0; ii < aNbSamples; ii++)
 150       {
 151         Standard_Real aParam = aBAcurve2d.FirstParameter() + ii * aDelta;
 152         gp_Pnt2d aP2d = aBAcurve2d.Value (aParam);
 153         gp_Pnt aPnt (aP2d.X(), aP2d.Y(), 0.);
 154         anArray (++anInd) = aPnt;
 155       }
 156     }
 157
 158     gp_Ax2 anAxis;
 159     Standard_Boolean anIsSingular;
 160     GeomLib::AxeOfInertia (anArray, anAxis, anIsSingular);
 161     gp_Pnt aBaryCentre = anAxis.Location();
 162     gp_Pnt2d aCentre2d (aBaryCentre.X(), aBaryCentre.Y());
 163     BRepTopAdaptor_FClass2d aClassifier (theFace, Precision::Confusion());
 164     BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
 165
 166     TopAbs_State aStatus = aClassifier.Perform (aCentre2d);
 167     gp_Pnt2d aP2d = aCentre2d;
 168     Standard_Integer anIndVertex = 0;
 169     const Standard_Integer aNbIter = 10;
 170     while (aStatus != TopAbs_IN && anIndVertex < aNbPointsOnContour)
 171     {
 172       gp_Pnt aVertexPnt = anArray (anIndVertex+1);
 173       gp_Pnt2d aVertexP2d (aVertexPnt.X(), aVertexPnt.Y());
 174       TColgp_SequenceOfPnt2d aPseq;
 175       aPseq.Append (aCentre2d);
 176       aPseq.Append (aVertexP2d);
 177       for (Standard_Integer ii = 1; ii <= aNbIter; ii++)
 178       {
 179         for (Standard_Integer jj = 1; jj < aPseq.Length(); jj++)
 180         {
 181           aP2d.SetXY ((aPseq(jj).XY() + aPseq(jj+1).XY())/2);
 182           aStatus = aClassifier.Perform (aP2d);
 183           if (aStatus == TopAbs_IN)
 184             break;
 185           else
 186           {
 187             aPseq.InsertAfter (jj, aP2d);
 188             jj++;
 189           }
 190         }
 191         if (aStatus == TopAbs_IN)
 192           break;
 193       }
 194       anIndVertex += aNbSamples;
 195     }
 196     aResPnt = aBAsurf.Value (aP2d.X(), aP2d.Y());
 197   } //case of complex boundaries
 198
 199   return aResPnt;
 200 }
 201
 202 void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
 203                                 const Standard_Real    theAverageArea,
 204                                 const Standard_Boolean theIsToAddFaces,
 205                                 Standard_Integer&      theNbExtremalFaces,
 206                                 TopTools_MapOfShape&   theExtremalFaces,
 207                                 const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
 208                                 const TopTools_DataMapOfShapeReal& theFaceAreaMap,
 209                                 const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
 210                                 TopoDS_Shape&          theRes,
 211                                 TopoDS_Shape&          theGlobalRes,
 212                                 TopTools_MapOfShape&   theRemovedFaces)
 213 {
 214   BRepAdaptor_Surface aBAsurf (theInputFace, Standard_False);
 215   GeomAbs_SurfaceType aType = aBAsurf.GetType();
 216
 217   BRep_Builder aBB;
 218   const Standard_Integer aNbFaces = (Standard_Integer) theFacesAndAreas.size();
 219
 220   const Standard_Integer aDiff = theNbExtremalFaces - theRemovedFaces.Extent();
 221
 222   Standard_Integer aSum = 0;
 223   while (aSum < aDiff) //global loop
 224   {
 225     Standard_Integer aNbFacesDone = 0, aNbFacesInTape = 0;
 226     TopoDS_Face aStartFace;
 227
 228     Standard_Integer aStartIndex = (theIsToAddFaces)? aNbFaces-1 : 0;
 229     Standard_Integer anEndIndex  = (theIsToAddFaces)? 0 : aNbFaces-1;
 230     Standard_Integer aStep = (theIsToAddFaces)? -1 : 1;
 231
 232     for (Standard_Integer ii = aStartIndex; ii != anEndIndex; ii += aStep)
 233     {
 234       const TopoDS_Face& aFace = TopoDS::Face (theFacesAndAreas[ii].first);
 235       if (!theRemovedFaces.Contains(aFace))
 236       {
 237         aStartFace = aFace;
 238         break;
 239       }
 240     }
 241     if (aStartFace.IsNull())
 242       break;
 243
 244     theRemovedFaces.Add (aStartFace);
 245
 246     TopoDS_Edge aCommonEdge;
 247     TopoDS_Face aNextFace;
 248     Standard_Real anExtremalArea = (theIsToAddFaces)? 0. : Precision::Infinite();
 249     for (TopExp_Explorer anExplo(aStartFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
 250     {
 251       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
 252       const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
 253       TopTools_ListIteratorOfListOfShape anItl (aFaceList);
 254       for (; anItl.More(); anItl.Next())
 255       {
 256         const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
 257         if (aFace.IsSame (aStartFace) ||
 258             theRemovedFaces.Contains(aFace))
 259           continue;
 260         Standard_Real anArea = theFaceAreaMap(aFace);
 261         Standard_Boolean anIsToExchange = (theIsToAddFaces)? (anArea > anExtremalArea) : (anArea < anExtremalArea);
 262         if (anIsToExchange)
 263         {
 264           anExtremalArea = anArea;
 265           aCommonEdge = anEdge;
 266           aNextFace = aFace;
 267         }
 268       }
 269     }
 270     if (aCommonEdge.IsNull()) //all adjacent faces are already removed
 271     {
 272       theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
 273       theNbExtremalFaces++;
 274       continue;
 275     }
 276
 277     //Start filling the shell
 278     aBB.Remove (theRes, aStartFace);
 279     aNbFacesDone++;
 280     TopoDS_Shell aShell;
 281     aBB.MakeShell (aShell);
 282     Standard_Real anAreaOfTape = 0.;
 283     aBB.Add (aShell, aStartFace);
 284     aNbFacesInTape++;
 285     anAreaOfTape += theFaceAreaMap (aStartFace);
 286
 287     Standard_Boolean anIsUiso = IsUiso (aCommonEdge, aStartFace);
 288     //Find another faces on this level
 289     TopoDS_Face aCurrentFace = aNextFace;
 290     TopoDS_Edge aCurrentEdge = aCommonEdge;
 291     Standard_Boolean anIsFirstDirection = Standard_True;
 292     aBB.Remove (theRes, aCurrentFace);
 293     theRemovedFaces.Add (aCurrentFace);
 294     if (theExtremalFaces.Contains (aCurrentFace))
 295     {
 296       aNbFacesDone++;
 297     }
 298     aBB.Add (aShell, aCurrentFace);
 299     aNbFacesInTape++;
 300     anAreaOfTape += theFaceAreaMap (aCurrentFace);
 301     Standard_Boolean anIsRound = Standard_False;
 302     for (;;) //local loop
 303     {
 304       TopoDS_Edge aNextEdge;
 305       for (TopExp_Explorer anExplo(aCurrentFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
 306       {
 307         const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
 308         if (anEdge.IsSame (aCurrentEdge))
 309           continue;
 310         const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
 311         TopTools_ListIteratorOfListOfShape anItl (aFaceList);
 312         for (; anItl.More(); anItl.Next())
 313         {
 314           const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
 315           if (aFace.IsSame (aCurrentFace))
 316             continue;
 317           if (aFace.IsSame (aStartFace))
 318           {
 319             anIsRound = Standard_True;
 320             if (aType > GeomAbs_Torus) { // non-elementary surfaces
 321               // remove last face to prevent close tape creation
 322               // it is a workaround for Tulip bos #26791
 323               // as there is a problem with closed tape on some surface types
 324               aBB.Remove (aShell, aCurrentFace);
 325               aNbFacesInTape--;
 326               anAreaOfTape -= theFaceAreaMap(aCurrentFace);
 327               aBB.Add(theRes, aCurrentFace); // aaajfa ???
 328               theRemovedFaces.Remove(aCurrentFace);
 329               if (theExtremalFaces.Contains(aCurrentFace)) {
 330                 aNbFacesDone--;
 331               }
 332             }
 333             break;
 334           }
 335           if (theRemovedFaces.Contains(aFace))
 336             continue;
 337           if (anIsUiso == IsUiso (anEdge, aFace))
 338           {
 339             aNextEdge = anEdge;
 340             aNextFace = aFace;
 341             break;
 342           }
 343         }
 344         if (anIsRound || !aNextEdge.IsNull())
 345           break;
 346       }
 347       if (anIsRound) //round tape: returned to start face
 348         break;
 349       if (aNextEdge.IsNull())
 350       {
 351         if (anIsFirstDirection)
 352         {
 353           aCurrentFace = aStartFace;
 354           aCurrentEdge = aCommonEdge;
 355           anIsFirstDirection = Standard_False;
 356           continue;
 357         }
 358         else
 359           break;
 360       }
 361
 362       aBB.Add (aShell, aNextFace);
 363       aNbFacesInTape++;
 364       anAreaOfTape += theFaceAreaMap (aNextFace);
 365       aBB.Remove (theRes, aNextFace);
 366       theRemovedFaces.Add (aNextFace);
 367       if (theExtremalFaces.Contains (aNextFace))
 368       {
 369         aNbFacesDone++;
 370       }
 371       aCurrentEdge = aNextEdge;
 372       aNextEdge.Nullify();
 373       aCurrentFace = aNextFace;
 374     } //end of local loop
 375
 376     //Tape is formed
 377     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aNbFacesInTape + aNbFacesDone : aNbFacesInTape - aNbFacesDone;
 378     if (!theIsToAddFaces && aNbFacesDone > 1)
 379     {
 380       Standard_Integer aRealNumberToSplit = (aNumberToSplit > 0)? aNumberToSplit : 1;
 381       Standard_Real anAverageAreaInTape = anAreaOfTape / aRealNumberToSplit;
 382       if (anAverageAreaInTape > theAverageArea)
 383       {
 384         Standard_Integer aNewNumberToSplit = RealToInt(round(anAreaOfTape / theAverageArea));
 385         if (aNewNumberToSplit < aNbFacesInTape)
 386         {
 387           Standard_Integer aNumberToIncrease = aNewNumberToSplit - aNumberToSplit;
 388           for (Standard_Integer jj = theNbExtremalFaces;
 389                jj < theNbExtremalFaces + aNumberToIncrease && jj < aNbFaces;
 390                jj++)
 391             theExtremalFaces.Add (theFacesAndAreas[jj].first);
 392           theNbExtremalFaces += aNumberToIncrease;
 393           aNumberToSplit = aNewNumberToSplit;
 394         }
 395       }
 396     }
 397     if (anIsRound && aNumberToSplit <= 1)
 398     {
 399       Standard_Integer aNumberToIncrease = 3 - aNumberToSplit;
 400       for (Standard_Integer jj = theNbExtremalFaces;
 401            jj < theNbExtremalFaces + aNumberToIncrease && jj < aNbFaces;
 402            jj++)
 403         theExtremalFaces.Add (theFacesAndAreas[jj].first);
 404       theNbExtremalFaces += aNumberToIncrease;
 405       aNumberToSplit = 3;
 406     }
 407     CorrectShell (aShell, theInputFace);
 408     ShapeUpgrade_UnifySameDomain aUnifier;
 409     aUnifier.Initialize (aShell, Standard_True, Standard_True);
 410     aUnifier.Build();
 411     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
 412     //Splitting
 413     TopoDS_Shape aLocalResult = aUnifiedShape;
 414     Standard_Integer aNbFacesInLocalResult = 1;
 415     if (aNumberToSplit > 1)
 416     {
 417 #if OCC_VERSION_LARGE < 0x07050304
 418       aNbFacesInLocalResult = 1;
 419 #else
 420       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
 421       aLocalTool.SetSplittingByNumber (Standard_True);
 422       aLocalTool.MaxArea() = -1;
 423       if (anIsUiso)
 424         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
 425       else
 426         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
 427       aLocalTool.Perform();
 428       aLocalResult = aLocalTool.Result();
 429       aNbFacesInLocalResult = aNumberToSplit;
 430 #endif
 431     }
 432     else
 433     {
 434       if (aNumberToSplit == 0)
 435       {
 436         if (theNbExtremalFaces < aNbFaces) {
 437           theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
 438           theNbExtremalFaces++;
 439         }
 440       }
 441     }
 442     aBB.Add (theGlobalRes, aLocalResult);
 443
 444     aSum += Abs(aNbFacesInTape - aNbFacesInLocalResult);
 445   } //end of global loop
 446
 447   //Second global loop
 448   TopoDS_Compound aSecondComp;
 449   aBB.MakeCompound (aSecondComp);
 450   while (aSum < aDiff)
 451   {
 452     TopoDS_Shape aMaxShell;
 453     Standard_Integer aMaxNbFaces = 0;
 454     TopoDS_Iterator anIter (theGlobalRes);
 455     for (; anIter.More(); anIter.Next())
 456     {
 457       const TopoDS_Shape& aShell = anIter.Value();
 458       TopTools_IndexedMapOfShape aFaceMap;
 459       TopExp::MapShapes (aShell, TopAbs_FACE, aFaceMap);
 460       if (aFaceMap.Extent() > aMaxNbFaces)
 461       {
 462         aMaxNbFaces = aFaceMap.Extent();
 463         aMaxShell = aShell;
 464       }
 465     }
 466
 467     if (aMaxNbFaces == 1)
 468       break;
 469
 470     aBB.Remove (theGlobalRes, aMaxShell);
 471     //Find iso
 472     Standard_Boolean anIsUiso = Standard_True;
 473     TopTools_IndexedDataMapOfShapeListOfShape aLocalEFmap;
 474     TopExp::MapShapesAndAncestors (aMaxShell, TopAbs_EDGE, TopAbs_FACE, aLocalEFmap);
 475     for (Standard_Integer jj = 1; jj <= aLocalEFmap.Extent(); jj++)
 476     {
 477       const TopoDS_Edge& anEdge = TopoDS::Edge (aLocalEFmap.FindKey(jj));
 478       const TopTools_ListOfShape& aFaceList = aLocalEFmap(jj);
 479       if (aFaceList.Extent() == 2)
 480       {
 481         const TopoDS_Face& aFace = TopoDS::Face (aFaceList.First());
 482         anIsUiso = IsUiso (anEdge, aFace);
 483         break;
 484       }
 485     }
 486     CorrectShell (aMaxShell, theInputFace);
 487     ShapeUpgrade_UnifySameDomain aUnifier;
 488     aUnifier.Initialize (aMaxShell, Standard_True, Standard_True);
 489     aUnifier.Build();
 490     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
 491     TopoDS_Shape aLocalResult = aUnifiedShape;
 492
 493     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aMaxNbFaces + (aDiff-aSum) : aMaxNbFaces - (aDiff-aSum);
 494     if (aNumberToSplit > 1)
 495     {
 496 #if OCC_VERSION_LARGE < 0x07050304
 497       aNumberToSplit = 1;
 498 #else
 499       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
 500       aLocalTool.SetSplittingByNumber (Standard_True);
 501       aLocalTool.MaxArea() = -1;
 502       if (anIsUiso)
 503         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
 504       else
 505         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
 506       aLocalTool.Perform();
 507       aLocalResult = aLocalTool.Result();
 508 #endif
 509     }
 510     else
 511       aNumberToSplit = 1;
 512
 513     aBB.Add (aSecondComp, aLocalResult);
 514
 515     if (theIsToAddFaces)
 516       break;
 517     aSum += aMaxNbFaces - aNumberToSplit;
 518   }
 519   aBB.Add (theGlobalRes, aSecondComp);
 520 }
 521
 522 //=================================================================================================
 523 bool GeomAlgoAPI_PointCloudOnFace::PointCloud(GeomShapePtr theFace,
 524                                               const int    theNumberOfPoints,
 525                                               GeomShapePtr thePoints,
 526                                               std::string& theError)
 527 {
 528 #ifdef _DEBUG
 529   std::cout << "GeomAlgoAPI_PointCloudOnFace::PointCloud" << std::endl;
 530 #endif
 531
 532 #if OCC_VERSION_LARGE < 0x07050304
 533   theError = "Improper OCCT version: please, use OCCT 7.5.3p4 or newer.";
 534   return false;
 535 #else
 536
 537   if (!theFace.get()) {
 538     theError = "Face for point cloud calculation is null";
 539     return false;
 540   }
 541
 542   TopoDS_Shape anInputShape = theFace->impl<TopoDS_Shape>();
 543
 544   if (anInputShape.ShapeType() != TopAbs_FACE) {
 545     theError = "Shape for normale calculation is not a face";
 546     return false;
 547   }
 548
 549   TopoDS_Face anInputFace = TopoDS::Face (anInputShape);
 550
 551   ShapeUpgrade_ShapeDivideArea tool (anInputFace);
 552   tool.SetSplittingByNumber (Standard_True);
 553   tool.NbParts() = theNumberOfPoints;
 554   tool.Perform();
 555   TopoDS_Shape res = tool.Result();
 556
 557   BRep_Builder aBB;
 558   TopoDS_Compound aGlobalRes;
 559   aBB.MakeCompound (aGlobalRes);
 560
 561   TopTools_IndexedMapOfShape aFaceMap;
 562   TopExp::MapShapes (res, TopAbs_FACE, aFaceMap);
 563   Standard_Integer aNbFaces = aFaceMap.Extent();
 564
 565   TopTools_IndexedDataMapOfShapeListOfShape aEFmap;
 566   TopExp::MapShapesAndAncestors (res, TopAbs_EDGE, TopAbs_FACE, aEFmap);
 567
 568   TopTools_MapOfShape aBiggestFaces, aSmallestFaces;
 569   Standard_Boolean aUseTriangulation = Standard_True;
 570   Standard_Boolean aSkipShared = Standard_False;
 571   if (aNbFaces != theNumberOfPoints)
 572   {
 573     Standard_Real aTotalArea = 0.;
 574     std::vector<std::pair<TopoDS_Shape, Standard_Real>> aFacesAndAreas (aNbFaces);
 575     for (Standard_Integer ii = 1; ii <= aNbFaces; ii++)
 576     {
 577       GProp_GProps aProps;
 578       BRepGProp::SurfaceProperties (aFaceMap(ii), aProps, aSkipShared, aUseTriangulation);
 579       Standard_Real anArea = aProps.Mass();
 580       aTotalArea += anArea;
 581       std::pair<TopoDS_Shape, Standard_Real> aFaceWithArea (aFaceMap(ii), anArea);
 582       aFacesAndAreas[ii-1] = aFaceWithArea;
 583     }
 584     std::sort (aFacesAndAreas.begin(), aFacesAndAreas.end(), comp);
 585
 586     Standard_Real anAverageArea = aTotalArea / theNumberOfPoints;
 587
 588     TopTools_DataMapOfShapeReal aFaceAreaMap;
 589     for (Standard_Integer ii = 0; ii < aNbFaces; ii++)
 590       aFaceAreaMap.Bind (aFacesAndAreas[ii].first, aFacesAndAreas[ii].second);
 591
 592     TopTools_MapOfShape aRemovedFaces;
 593
 594     if (aNbFaces < theNumberOfPoints)
 595     {
 596       Standard_Integer aNbMissingFaces = theNumberOfPoints - aNbFaces;
 597       for (Standard_Integer ii = aNbFaces-1; ii > aNbFaces - aNbMissingFaces - 1; ii--)
 598         aBiggestFaces.Add (aFacesAndAreas[ii].first);
 599
 600       ModifyFacesForGlobalResult (anInputFace, anAverageArea,
 601                                   Standard_True, //to add faces
 602                                   aNbMissingFaces, aBiggestFaces,
 603                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
 604                                   res, aGlobalRes,
 605                                   aRemovedFaces);
 606     }
 607     else //aNbFaces > theNumberOfPoints
 608     {
 609       Standard_Integer aNbExcessFaces = aNbFaces - theNumberOfPoints;
 610       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
 611         aSmallestFaces.Add (aFacesAndAreas[ii].first);
 612
 613       TopTools_IndexedDataMapOfShapeListOfShape aVFmap;
 614       TopExp::MapShapesAndAncestors (res, TopAbs_VERTEX, TopAbs_FACE, aVFmap);
 615
 616       //Remove smallest faces with free boundaries
 617       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
 618       {
 619         const TopoDS_Face& aFace = TopoDS::Face (aFacesAndAreas[ii].first);
 620         Standard_Boolean anIsFreeBoundFound = Standard_False;
 621         TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
 622         for (; anExplo.More(); anExplo.Next())
 623         {
 624           const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
 625           if (!BRep_Tool::Degenerated (anEdge) &&
 626               aEFmap.FindFromKey(anEdge).Extent() < 2)
 627           {
 628             anIsFreeBoundFound = Standard_True;
 629             break;
 630           }
 631         }
 632         if (anIsFreeBoundFound)
 633         {
 634           Standard_Real aMaxArea = 0.;
 635           for (anExplo.Init(aFace, TopAbs_VERTEX); anExplo.More(); anExplo.Next())
 636           {
 637             const TopoDS_Shape& aVertex = anExplo.Current();
 638             const TopTools_ListOfShape& aFaceList = aVFmap.FindFromKey (aVertex);
 639             TopTools_ListIteratorOfListOfShape anItl (aFaceList);
 640             for (; anItl.More(); anItl.Next())
 641             {
 642               Standard_Real anArea = aFaceAreaMap (anItl.Value());
 643               if (anArea > aMaxArea)
 644                 aMaxArea = anArea;
 645             }
 646           }
 647           Standard_Real anArreaOfSmallestFace = aFaceAreaMap (aFace);
 648           if (anArreaOfSmallestFace < aMaxArea / 16)
 649           {
 650             aBB.Remove (res, aFace);
 651             aRemovedFaces.Add (aFace);
 652           }
 653         }
 654       }
 655
 656       ModifyFacesForGlobalResult (anInputFace, anAverageArea,
 657                                   Standard_False, //to decrease number of faces
 658                                   aNbExcessFaces, aSmallestFaces,
 659                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
 660                                   res, aGlobalRes,
 661                                   aRemovedFaces);
 662     }
 663   }
 664
 665   aBB.Add (aGlobalRes, res);
 666
 667   int iface = 1;
 668   TopoDS_Compound aCompound;
 669   aBB.MakeCompound (aCompound);
 670   for (TopExp_Explorer aGlobalExplo(aGlobalRes, TopAbs_FACE);
 671        aGlobalExplo.More(), iface <= theNumberOfPoints;
 672        aGlobalExplo.Next(), iface++)
 673   {
 674     const TopoDS_Face& aFace = TopoDS::Face (aGlobalExplo.Current());
 675     Standard_Boolean anIsNaturalRestrictions = Standard_True;
 676     TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
 677     for (; anExplo.More(); anExplo.Next())
 678     {
 679       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
 680       if (BRep_Tool::Degenerated (anEdge))
 681         continue;
 682       if (!aEFmap.Contains(anEdge) ||
 683           aEFmap.FindFromKey(anEdge).Extent() < 2)
 684       {
 685         anIsNaturalRestrictions = Standard_False;
 686         break;
 687       }
 688     }
 689
 690     gp_Pnt aPnt = GetMidPnt2d (aFace, anIsNaturalRestrictions);
 691     TopoDS_Vertex aVertex = BRepLib_MakeVertex (aPnt);
 692     aBB.Add (aCompound, aVertex);
 693   }
 694
 695   thePoints->setImpl(new TopoDS_Shape(aCompound));
 696
 697   return true;
 698 #endif
 699 }