src/GEOMImpl/GEOMImpl_WrappingDriver.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include <GEOMImpl_WrappingDriver.hxx>
  24
  25 #include <GEOMAlgo_Splitter.hxx>
  26 #include <GEOMUtils.hxx>
  27
  28 #include <GEOMImpl_ShapeDriver.hxx>
  29 #include <GEOMImpl_IWrap.hxx>
  30 #include <GEOMImpl_Types.hxx>
  31 #include <GEOMImpl_Block6Explorer.hxx>
  32 #include <GEOM_Object.hxx>
  33
  34 #include <BRep_Tool.hxx>
  35 #include <BRepBuilderAPI_MakeFace.hxx>
  36 #include <BRepPrimAPI_MakeSphere.hxx>
  37 #if OCC_VERSION_LARGE < 0x07070000
  38 #include <BRepAdaptor_HSurface.hxx>
  39 #else
  40 #include <BRepAdaptor_Surface.hxx>
  41 #endif
  42
  43 #include <Geom_Curve.hxx>
  44 #include <Geom_Surface.hxx>
  45 #include <Geom_Plane.hxx>
  46 #include <Geom_BSplineSurface.hxx>
  47 #include <GeomPlate_Surface.hxx>
  48 #include <GeomPlate_BuildPlateSurface.hxx>
  49 #include <GeomPlate_PointConstraint.hxx>
  50 #include <GeomPlate_MakeApprox.hxx>
  51 #include <GeomPlate_BuildAveragePlane.hxx>
  52
  53 #include <gp_Pnt.hxx>
  54 #include <gp_Sphere.hxx>
  55 #include <TopExp_Explorer.hxx>
  56 #include <TopoDS.hxx>
  57 #include <TopoDS_Shape.hxx>
  58 #include <TopoDS_Wire.hxx>
  59 #include <TopoDS_Vertex.hxx>
  60 #include <TColgp_HArray1OfPnt.hxx>
  61 #include <TColStd_HSequenceOfTransient.hxx>
  62 #include <ShapeFix_Face.hxx>
  63
  64 #include <Standard_NullObject.hxx>
  65
  66 #include <numeric>
  67 #include <vector>
  68
  69 #include <Eigen/Dense>
  70
  71 //=======================================================================
  72 // function : GetID
  73 // purpose  :
  74 //=======================================================================
  75 const Standard_GUID &GEOMImpl_WrappingDriver::GetID()
  76 {
  77   static Standard_GUID aWrapDriver("FF1BBB71-729D-4E83-8232-78E74FC5637C");
  78   return aWrapDriver;
  79 }
  80
  81 //=======================================================================
  82 // function : GEOMImpl_WrappingDriver
  83 // purpose  :
  84 //=======================================================================
  85 GEOMImpl_WrappingDriver::GEOMImpl_WrappingDriver()
  86 {
  87 }
  88
  89 //=======================================================================
  90 // function : createPointsOnEdges
  91 // purpose  : create points on edges
  92 //=======================================================================
  93 static void createPointsOnEdges(const Handle(TColStd_HSequenceOfTransient) & theEdgesFuncs,
  94                                 std::vector<gp_Pnt> &thePoints)
  95 {
  96   for (int i = 1; i <= theEdgesFuncs->Length(); i++)
  97   {
  98     Handle(GEOM_Function) aRefShape = Handle(GEOM_Function)::DownCast(theEdgesFuncs->Value(i));
  99     TopoDS_Shape aShape_i = aRefShape->GetValue();
 100     if (aShape_i.IsNull())
 101     {
 102       Standard_NullObject::Raise("Edge for face construction is null");
 103     }
 104     TopExp_Explorer exp(aShape_i, TopAbs_EDGE);
 105     for (; exp.More(); exp.Next())
 106     {
 107       Standard_Real aFP, aLP, aP;
 108       Handle(Geom_Curve) aCurve = BRep_Tool::Curve(TopoDS::Edge(exp.Current()), aFP, aLP);
 109       if (!aCurve.IsNull())
 110       {
 111         for (int aPar = 1; aPar <= 9; aPar++)
 112         {
 113           aP = aFP + (aLP - aFP) * (aPar * 0.1);
 114           gp_Pnt aPnt = aCurve->Value(aP);
 115           thePoints.push_back(aPnt);
 116         }
 117       }
 118     }
 119   }
 120 }
 121
 122 //================================================================================
 123 // function : maxDistanceToFace
 124 // purpose  : finds max distanse between points and a face
 125 //================================================================================
 126 static Standard_Real maxDistanceToFace(const std::vector<gp_Pnt>& thePoints,
 127                                        const TopoDS_Face& theFace,
 128                                        const Standard_Real theTolerance)
 129 {
 130   Standard_Real U, V;
 131   Standard_Real aMaxDist = 0.;
 132   for(auto& aPnt : thePoints)
 133   {
 134     gp_Pnt aProj = GEOMUtils::ProjectPointOnFace(aPnt, theFace,U,V,theTolerance);
 135     Standard_Real aDist = aProj.Distance(aPnt);
 136     if(aDist > aMaxDist)
 137     {
 138       aMaxDist = aDist;
 139     }
 140   }
 141   return aMaxDist;
 142 }
 143
 144 //================================================================================
 145 // function : divideSphericalShape
 146 // purpose  : divide spherical shape into two piece
 147 //            and choose that part that has the points
 148 //================================================================================
 149 static void divideSphericalShape(TopoDS_Shape &theShape,
 150                                  const TopoDS_Wire &theWire,
 151                                  const std::vector<gp_Pnt> &thePoints,
 152                                  const Standard_Real theTolerance)
 153 {
 154     TopExp_Explorer anExp(theShape, TopAbs_FACE);
 155     const TopoDS_Face& aFace = TopoDS::Face(anExp.Current());
 156     TopoDS_Face aToolFace;
 157     GEOMImpl_Block6Explorer::MakeFace(theWire, false, aToolFace);
 158     if(!aToolFace.IsNull())
 159     {
 160         //split sphere and choose right part
 161         GEOMAlgo_Splitter PS;
 162         PS.AddArgument(aFace);
 163         PS.AddTool(aToolFace);
 164         PS.SetLimit(TopAbs_FACE);
 165         PS.Perform();
 166         TopoDS_Shape aResultShape = PS.Shape();
 167         if(!aResultShape.IsNull())
 168         {
 169             anExp.Init(aResultShape, TopAbs_FACE);
 170             for (; anExp.More(); anExp.Next())
 171             {
 172                 Standard_Real aDist = maxDistanceToFace(thePoints, TopoDS::Face(anExp.Current()), theTolerance);
 173                 if(aDist < theTolerance)
 174                 {
 175                     theShape = TopoDS::Face(anExp.Current());
 176                     break;
 177                 }
 178             }
 179         }
 180     }
 181 }
 182
 183 //================================================================================
 184 // function : makeFaceFromPointsAndWire
 185 // purpose  : Create face from set of points with the same approach as for
 186 //            geompy.makeSmoothingSurface. Cut resulting surface with a wire
 187 //================================================================================
 188 static TopoDS_Shape makeFaceFromPointsAndWire(const TopoDS_Wire &theWire,
 189                                               const std::vector<gp_Pnt> &theAllPoints)
 190 {
 191   int nbPoints = theAllPoints.size();
 192   Handle(TColgp_HArray1OfPnt) anArrayofPnt = new TColgp_HArray1OfPnt(1, nbPoints);
 193   for (int i = 0; i < nbPoints; i++) {
 194     gp_Pnt aP = theAllPoints[i];
 195     anArrayofPnt->SetValue(i+1, aP);
 196   }
 197
 198   // Try to build smoothing surface
 199   GeomPlate_BuildAveragePlane gpbap(anArrayofPnt,anArrayofPnt->Length(),Precision::Confusion(),1,1);
 200   Handle(Geom_Plane) plane(gpbap.Plane());
 201
 202   Standard_Real Umin, Umax, Vmin, Vmax;
 203   gpbap.MinMaxBox(Umin,Umax,Vmin,Vmax);
 204
 205   TopoDS_Face aInitShape;
 206   BRepBuilderAPI_MakeFace mf(plane,Umin,Umax,Vmin,Vmax,Precision::Confusion());
 207   aInitShape =  mf.Face();
 208
 209   GeomPlate_BuildPlateSurface aBuilder(3,10);
 210   // ** Initialization of surface
 211 #if OCC_VERSION_LARGE < 0x07070000
 212   Handle(BRepAdaptor_HSurface) HSI = new BRepAdaptor_HSurface();
 213   HSI->ChangeSurface().Initialize(aInitShape);
 214   aBuilder.LoadInitSurface( BRep_Tool::Surface(HSI->ChangeSurface().Face()));
 215 #else
 216   Handle(BRepAdaptor_Surface) HSI = new BRepAdaptor_Surface();
 217   HSI->Initialize(aInitShape);
 218   aBuilder.LoadInitSurface( BRep_Tool::Surface(HSI->Face()) );
 219 #endif
 220
 221   Standard_Integer j, j1, j2;
 222   j1 = anArrayofPnt->Lower();
 223   j2 = anArrayofPnt->Upper();
 224   for (j=j1; j<=j2 ; j++)
 225   {
 226     gp_Pnt aPnt = anArrayofPnt->Value(j);
 227     Handle(GeomPlate_PointConstraint) PCont = new GeomPlate_PointConstraint(aPnt,0);
 228     aBuilder.Add(PCont);
 229   }
 230   aBuilder.Perform();
 231   Handle(GeomPlate_Surface) gpPlate = aBuilder.Surface();
 232   if(gpPlate.IsNull())
 233   {
 234     Standard_ConstructionError::Raise("Not possible to build a face with given input, GeomPlate_Surface failed");
 235   }
 236
 237   Standard_Real seuil = Max(0.0001,10*aBuilder.G0Error());
 238   GeomPlate_MakeApprox Mapp(gpPlate,0.0001,2,8,seuil);
 239   Handle(Geom_Surface) aSurf(Mapp.Surface());
 240
 241   //cut surface with a face
 242   TopoDS_Shape aShape;
 243   BRepBuilderAPI_MakeFace aMkFace(aSurf, theWire);
 244   if (aMkFace.IsDone())
 245   {
 246     aShape = aMkFace.Shape();
 247   }
 248
 249   TopoDS_Face aFace = TopoDS::Face(aShape);
 250   Handle(ShapeFix_Face) aFix = new ShapeFix_Face(aFace);
 251   aFix->Perform();
 252   aFix->FixOrientation();
 253   aFace = aFix->Face();
 254
 255   return aFace;
 256 }
 257
 258 //=======================================================================
 259 // function : Execute
 260 // purpose  :
 261 //=======================================================================
 262 Standard_Integer GEOMImpl_WrappingDriver::Execute(Handle(TFunction_Logbook) & log) const
 263 {
 264   if (Label().IsNull())
 265     return 0;
 266   Handle(GEOM_Function) aFunction = GEOM_Function::GetFunction(Label());
 267
 268   GEOMImpl_IWrap aCI(aFunction);
 269   Standard_Integer aType = aFunction->GetType();
 270
 271   TopoDS_Shape aShape;
 272
 273   if (aType == WRAPPING_FACE)
 274   {
 275     Handle(TColStd_HSequenceOfTransient) anEdges = aCI.GetEdges();
 276     Handle(TColStd_HSequenceOfTransient) aVertices = aCI.GetVertices();
 277
 278     int nbEdge = anEdges->Length();
 279     int nbVertices = aVertices->Length();
 280
 281     if (nbEdge < 1 || nbVertices < 1)
 282     {
 283       Standard_ConstructionError::Raise("No edges or vertices given");
 284     }
 285
 286     Standard_Real aTolerance = aCI.GetTolerance();
 287     if (aTolerance < Precision::Confusion())
 288       aTolerance = Precision::Confusion();
 289
 290     std::vector<gp_Pnt> anAllPoints, aPoints;
 291     createPointsOnEdges(anEdges, anAllPoints);
 292     for (int i = 1; i <= aVertices->Length(); i++)
 293     {
 294       Handle(GEOM_Function) aRefShape = Handle(GEOM_Function)::DownCast(aVertices->Value(i));
 295       TopoDS_Shape aShape_i = aRefShape->GetValue();
 296       if (aShape_i.IsNull())
 297       {
 298         Standard_NullObject::Raise("Vertex for face construction is null");
 299       }
 300       TopExp_Explorer exp (aShape_i, TopAbs_VERTEX);
 301       for (; exp.More(); exp.Next())
 302       {
 303         gp_Pnt aP = BRep_Tool::Pnt(TopoDS::Vertex(exp.Current()));
 304         anAllPoints.push_back(aP);
 305         aPoints.push_back(aP);
 306       }
 307     }
 308
 309     if(anAllPoints.size() == 0)
 310     {
 311       Standard_NullObject::Raise("Shape creation imposible, no valid objects were given");
 312     }
 313
 314     //create wire from edges
 315     TopoDS_Wire aWire = GEOMImpl_ShapeDriver::MakeWireFromEdges(anEdges, aTolerance);
 316     if(aWire.IsNull())
 317     {
 318       Standard_NullObject::Raise("Given edges does not create closed contour");
 319     }
 320
 321     aShape = createWrappedFace(aWire, anAllPoints, aPoints, aTolerance);
 322   }
 323
 324   if (aShape.IsNull())
 325     return 0;
 326
 327   aFunction->SetValue(aShape);
 328
 329   log->SetTouched(Label());
 330
 331   return 1;
 332 }
 333
 334 //================================================================================
 335 /*!
 336  * \brief Create wrapped face from vector of points and wire
 337  */
 338 //================================================================================
 339 TopoDS_Shape GEOMImpl_WrappingDriver::createWrappedFace(const TopoDS_Wire &theWire,
 340                                                         const std::vector<gp_Pnt> &theAllPoints,
 341                                                         const std::vector<gp_Pnt> &thePoints,
 342                                                         const Standard_Real theTolerance) const
 343 {
 344   TopoDS_Shape aShape;
 345
 346   gp_Pnt aCenter(0,0,0);
 347   Standard_Real aRadius = 100;
 348   if(isSphere(theAllPoints, aCenter, aRadius, theTolerance))
 349   {
 350     aShape = BRepPrimAPI_MakeSphere(aCenter, aRadius).Shape();
 351     divideSphericalShape(aShape, theWire, thePoints, theTolerance);
 352     return aShape;
 353   }
 354
 355   aShape = makeFaceFromPointsAndWire(theWire, theAllPoints);
 356   return aShape;
 357 }
 358
 359 //================================================================================
 360 /*!
 361  * \brief Returns a name of creation operation and names and values of creation parameters
 362  */
 363 //================================================================================
 364 Standard_Boolean GEOMImpl_WrappingDriver::isSphere(const std::vector<gp_Pnt>& thePoints,
 365                                  gp_Pnt& theCenter,
 366                                  Standard_Real& theRadius,
 367                                  const Standard_Real theTolerance)const
 368 {
 369   int aNumPoints = thePoints.size();
 370   if(aNumPoints == 0)
 371   {
 372     return false;
 373   }
 374
 375   // Create coefficient matrix A and right-hand side vector f
 376   Eigen::MatrixXd A(aNumPoints, 4);
 377   Eigen::VectorXd f(aNumPoints);
 378
 379   Standard_Real X(.0);
 380   Standard_Real Y(.0);
 381   Standard_Real Z(.0);
 382
 383   for (int i = 0; i < aNumPoints; ++i)
 384   {
 385     X = thePoints[i].X();
 386     Y = thePoints[i].Y();
 387     Z = thePoints[i].Z();
 388     A(i, 0) = X * 2;
 389     A(i, 1) = Y * 2;
 390     A(i, 2) = Z * 2;
 391     A(i, 3) = 1.0;
 392     f(i) = X * X + Y * Y + Z * Z;
 393   }
 394
 395   // Solve linear equations to get coefficients
 396   Eigen::VectorXd c = A.colPivHouseholderQr().solve(f);
 397
 398   Standard_Real t = c[0] * c[0] + c[1] * c[1] + c[2] * c[2] + c[3];
 399   theRadius = std::sqrt(t);
 400   theCenter.SetCoord(c[0], c[1], c[2]);
 401
 402   //check that all points belong to the sphere within theTolerance
 403   std::vector<Standard_Real> aDists;
 404   for(const auto& aPnt : thePoints)
 405   {
 406     aDists.push_back(aPnt.Distance(theCenter));
 407   }
 408   if (!aDists.empty())
 409   {
 410     Standard_Real sumDistances = std::accumulate(aDists.begin(), aDists.end(), 0.0);
 411     Standard_Real averageDistance = sumDistances / aDists.size();
 412
 413     if (std::fabs(averageDistance - theRadius) > theTolerance) {
 414         return false;
 415     }
 416   }
 417
 418   return true;
 419 }
 420
 421 //================================================================================
 422 /*!
 423  * \brief Returns a name of creation operation and names and values of creation parameters
 424  */
 425 //================================================================================
 426
 427 bool GEOMImpl_WrappingDriver::
 428     GetCreationInformation(std::string &theOperationName,
 429                            std::vector<GEOM_Param> &theParams)
 430 {
 431   if (Label().IsNull()) return 0;
 432   Handle(GEOM_Function) function = GEOM_Function::GetFunction(Label());
 433
 434   GEOMImpl_IWrap aCI( function );
 435   Standard_Integer aType = function->GetType();
 436
 437   theOperationName = "WRAPPEDFACE";
 438
 439   switch ( aType ) {
 440   case WRAPPING_FACE:
 441     AddParam(theParams, "Edges", aCI.GetEdges());
 442     AddParam(theParams, "Vertices", aCI.GetVertices());
 443     AddParam(theParams, "Tolerance", aCI.GetTolerance());
 444     break;
 445   default:
 446     return false;
 447   }
 448
 449   return true;
 450 }
 451
 452 IMPLEMENT_STANDARD_RTTIEXT(GEOMImpl_WrappingDriver, GEOM_BaseDriver)