geom/src/GEOMImpl/GEOMImpl_WrappingDriver.cxx

457 lines
15 KiB
C++

// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
// These GeomPlate_* headers must be included before Salome ones
// because of conflicting OK variable/constant between
// OCCT Plate_Plate.hxx and GEOM GEOM_IOperations.hxx
#include <GeomPlate_Surface.hxx>
#include <GeomPlate_BuildPlateSurface.hxx>
#include <GeomPlate_PointConstraint.hxx>
#include <GeomPlate_MakeApprox.hxx>
#include <GeomPlate_BuildAveragePlane.hxx>
#include <GEOMImpl_WrappingDriver.hxx>
#include <GEOMAlgo_Splitter.hxx>
#include <GEOMUtils.hxx>
#include <GEOMImpl_ShapeDriver.hxx>
#include <GEOMImpl_IWrap.hxx>
#include <GEOMImpl_Types.hxx>
#include <GEOMImpl_Block6Explorer.hxx>
#include <GEOM_Object.hxx>
#include <BRep_Tool.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
#if OCC_VERSION_LARGE < 0x07070000
#include <BRepAdaptor_HSurface.hxx>
#else
#include <BRepAdaptor_Surface.hxx>
#endif
#include <Geom_Curve.hxx>
#include <Geom_Surface.hxx>
#include <Geom_Plane.hxx>
#include <Geom_BSplineSurface.hxx>
#include <gp_Pnt.hxx>
#include <gp_Sphere.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Wire.hxx>
#include <TopoDS_Vertex.hxx>
#include <TColgp_HArray1OfPnt.hxx>
#include <TColStd_HSequenceOfTransient.hxx>
#include <ShapeFix_Face.hxx>
#include <Standard_NullObject.hxx>
#include <numeric>
#include <vector>
#include <Eigen/Dense>
//=======================================================================
// function : GetID
// purpose :
//=======================================================================
const Standard_GUID &GEOMImpl_WrappingDriver::GetID()
{
static Standard_GUID aWrapDriver("FF1BBB71-729D-4E83-8232-78E74FC5637C");
return aWrapDriver;
}
//=======================================================================
// function : GEOMImpl_WrappingDriver
// purpose :
//=======================================================================
GEOMImpl_WrappingDriver::GEOMImpl_WrappingDriver()
{
}
//=======================================================================
// function : createPointsOnEdges
// purpose : create points on edges
//=======================================================================
static void createPointsOnEdges(const Handle(TColStd_HSequenceOfTransient) & theEdgesFuncs,
std::vector<gp_Pnt> &thePoints)
{
for (int i = 1; i <= theEdgesFuncs->Length(); i++)
{
Handle(GEOM_Function) aRefShape = Handle(GEOM_Function)::DownCast(theEdgesFuncs->Value(i));
TopoDS_Shape aShape_i = aRefShape->GetValue();
if (aShape_i.IsNull())
{
Standard_NullObject::Raise("Edge for face construction is null");
}
TopExp_Explorer exp(aShape_i, TopAbs_EDGE);
for (; exp.More(); exp.Next())
{
Standard_Real aFP, aLP, aP;
Handle(Geom_Curve) aCurve = BRep_Tool::Curve(TopoDS::Edge(exp.Current()), aFP, aLP);
if (!aCurve.IsNull())
{
for (int aPar = 1; aPar <= 9; aPar++)
{
aP = aFP + (aLP - aFP) * (aPar * 0.1);
gp_Pnt aPnt = aCurve->Value(aP);
thePoints.push_back(aPnt);
}
}
}
}
}
//================================================================================
// function : maxDistanceToFace
// purpose : finds max distanse between points and a face
//================================================================================
static Standard_Real maxDistanceToFace(const std::vector<gp_Pnt>& thePoints,
const TopoDS_Face& theFace,
const Standard_Real theTolerance)
{
Standard_Real U, V;
Standard_Real aMaxDist = 0.;
for(auto& aPnt : thePoints)
{
gp_Pnt aProj = GEOMUtils::ProjectPointOnFace(aPnt, theFace,U,V,theTolerance);
Standard_Real aDist = aProj.Distance(aPnt);
if(aDist > aMaxDist)
{
aMaxDist = aDist;
}
}
return aMaxDist;
}
//================================================================================
// function : divideSphericalShape
// purpose : divide spherical shape into two piece
// and choose that part that has the points
//================================================================================
static void divideSphericalShape(TopoDS_Shape &theShape,
const TopoDS_Wire &theWire,
const std::vector<gp_Pnt> &thePoints,
const Standard_Real theTolerance)
{
TopExp_Explorer anExp(theShape, TopAbs_FACE);
const TopoDS_Face& aFace = TopoDS::Face(anExp.Current());
TopoDS_Face aToolFace;
GEOMImpl_Block6Explorer::MakeFace(theWire, false, aToolFace);
if(!aToolFace.IsNull())
{
//split sphere and choose right part
GEOMAlgo_Splitter PS;
PS.AddArgument(aFace);
PS.AddTool(aToolFace);
PS.SetLimit(TopAbs_FACE);
PS.Perform();
TopoDS_Shape aResultShape = PS.Shape();
if(!aResultShape.IsNull())
{
anExp.Init(aResultShape, TopAbs_FACE);
for (; anExp.More(); anExp.Next())
{
Standard_Real aDist = maxDistanceToFace(thePoints, TopoDS::Face(anExp.Current()), theTolerance);
if(aDist < theTolerance)
{
theShape = TopoDS::Face(anExp.Current());
break;
}
}
}
}
}
//================================================================================
// function : makeFaceFromPointsAndWire
// purpose : Create face from set of points with the same approach as for
// geompy.makeSmoothingSurface. Cut resulting surface with a wire
//================================================================================
static TopoDS_Shape makeFaceFromPointsAndWire(const TopoDS_Wire &theWire,
const std::vector<gp_Pnt> &theAllPoints)
{
int nbPoints = theAllPoints.size();
Handle(TColgp_HArray1OfPnt) anArrayofPnt = new TColgp_HArray1OfPnt(1, nbPoints);
for (int i = 0; i < nbPoints; i++) {
gp_Pnt aP = theAllPoints[i];
anArrayofPnt->SetValue(i+1, aP);
}
// Try to build smoothing surface
GeomPlate_BuildAveragePlane gpbap(anArrayofPnt,anArrayofPnt->Length(),Precision::Confusion(),1,1);
Handle(Geom_Plane) plane(gpbap.Plane());
Standard_Real Umin, Umax, Vmin, Vmax;
gpbap.MinMaxBox(Umin,Umax,Vmin,Vmax);
TopoDS_Face aInitShape;
BRepBuilderAPI_MakeFace mf(plane,Umin,Umax,Vmin,Vmax,Precision::Confusion());
aInitShape = mf.Face();
GeomPlate_BuildPlateSurface aBuilder(3,10);
// ** Initialization of surface
#if OCC_VERSION_LARGE < 0x07070000
Handle(BRepAdaptor_HSurface) HSI = new BRepAdaptor_HSurface();
HSI->ChangeSurface().Initialize(aInitShape);
aBuilder.LoadInitSurface( BRep_Tool::Surface(HSI->ChangeSurface().Face()));
#else
Handle(BRepAdaptor_Surface) HSI = new BRepAdaptor_Surface();
HSI->Initialize(aInitShape);
aBuilder.LoadInitSurface( BRep_Tool::Surface(HSI->Face()) );
#endif
Standard_Integer j, j1, j2;
j1 = anArrayofPnt->Lower();
j2 = anArrayofPnt->Upper();
for (j=j1; j<=j2 ; j++)
{
gp_Pnt aPnt = anArrayofPnt->Value(j);
Handle(GeomPlate_PointConstraint) PCont = new GeomPlate_PointConstraint(aPnt,0);
aBuilder.Add(PCont);
}
aBuilder.Perform();
Handle(GeomPlate_Surface) gpPlate = aBuilder.Surface();
if(gpPlate.IsNull())
{
Standard_ConstructionError::Raise("Not possible to build a face with given input, GeomPlate_Surface failed");
}
Standard_Real seuil = Max(0.0001,10*aBuilder.G0Error());
GeomPlate_MakeApprox Mapp(gpPlate,0.0001,2,8,seuil);
Handle(Geom_Surface) aSurf(Mapp.Surface());
//cut surface with a face
TopoDS_Shape aShape;
BRepBuilderAPI_MakeFace aMkFace(aSurf, theWire);
if (aMkFace.IsDone())
{
aShape = aMkFace.Shape();
}
TopoDS_Face aFace = TopoDS::Face(aShape);
Handle(ShapeFix_Face) aFix = new ShapeFix_Face(aFace);
aFix->Perform();
aFix->FixOrientation();
aFace = aFix->Face();
return aFace;
}
//=======================================================================
// function : Execute
// purpose :
//=======================================================================
Standard_Integer GEOMImpl_WrappingDriver::Execute(Handle(TFunction_Logbook) & log) const
{
if (Label().IsNull())
return 0;
Handle(GEOM_Function) aFunction = GEOM_Function::GetFunction(Label());
GEOMImpl_IWrap aCI(aFunction);
Standard_Integer aType = aFunction->GetType();
TopoDS_Shape aShape;
if (aType == WRAPPING_FACE)
{
Handle(TColStd_HSequenceOfTransient) anEdges = aCI.GetEdges();
Handle(TColStd_HSequenceOfTransient) aVertices = aCI.GetVertices();
int nbEdge = anEdges->Length();
int nbVertices = aVertices->Length();
if (nbEdge < 1 || nbVertices < 1)
{
Standard_ConstructionError::Raise("No edges or vertices given");
}
Standard_Real aTolerance = aCI.GetTolerance();
if (aTolerance < Precision::Confusion())
aTolerance = Precision::Confusion();
std::vector<gp_Pnt> anAllPoints, aPoints;
createPointsOnEdges(anEdges, anAllPoints);
for (int i = 1; i <= aVertices->Length(); i++)
{
Handle(GEOM_Function) aRefShape = Handle(GEOM_Function)::DownCast(aVertices->Value(i));
TopoDS_Shape aShape_i = aRefShape->GetValue();
if (aShape_i.IsNull())
{
Standard_NullObject::Raise("Vertex for face construction is null");
}
TopExp_Explorer exp (aShape_i, TopAbs_VERTEX);
for (; exp.More(); exp.Next())
{
gp_Pnt aP = BRep_Tool::Pnt(TopoDS::Vertex(exp.Current()));
anAllPoints.push_back(aP);
aPoints.push_back(aP);
}
}
if(anAllPoints.size() == 0)
{
Standard_NullObject::Raise("Shape creation imposible, no valid objects were given");
}
//create wire from edges
TopoDS_Wire aWire = GEOMImpl_ShapeDriver::MakeWireFromEdges(anEdges, aTolerance);
if(aWire.IsNull())
{
Standard_NullObject::Raise("Given edges does not create closed contour");
}
aShape = createWrappedFace(aWire, anAllPoints, aPoints, aTolerance);
}
if (aShape.IsNull())
return 0;
aFunction->SetValue(aShape);
log->SetTouched(Label());
return 1;
}
//================================================================================
/*!
* \brief Create wrapped face from vector of points and wire
*/
//================================================================================
TopoDS_Shape GEOMImpl_WrappingDriver::createWrappedFace(const TopoDS_Wire &theWire,
const std::vector<gp_Pnt> &theAllPoints,
const std::vector<gp_Pnt> &thePoints,
const Standard_Real theTolerance) const
{
TopoDS_Shape aShape;
gp_Pnt aCenter(0,0,0);
Standard_Real aRadius = 100;
if(isSphere(theAllPoints, aCenter, aRadius, theTolerance))
{
aShape = BRepPrimAPI_MakeSphere(aCenter, aRadius).Shape();
divideSphericalShape(aShape, theWire, thePoints, theTolerance);
return aShape;
}
aShape = makeFaceFromPointsAndWire(theWire, theAllPoints);
return aShape;
}
//================================================================================
/*!
* \brief Returns a name of creation operation and names and values of creation parameters
*/
//================================================================================
Standard_Boolean GEOMImpl_WrappingDriver::isSphere(const std::vector<gp_Pnt>& thePoints,
gp_Pnt& theCenter,
Standard_Real& theRadius,
const Standard_Real theTolerance)const
{
int aNumPoints = thePoints.size();
if(aNumPoints == 0)
{
return false;
}
// Create coefficient matrix A and right-hand side vector f
Eigen::MatrixXd A(aNumPoints, 4);
Eigen::VectorXd f(aNumPoints);
Standard_Real X(.0);
Standard_Real Y(.0);
Standard_Real Z(.0);
for (int i = 0; i < aNumPoints; ++i)
{
X = thePoints[i].X();
Y = thePoints[i].Y();
Z = thePoints[i].Z();
A(i, 0) = X * 2;
A(i, 1) = Y * 2;
A(i, 2) = Z * 2;
A(i, 3) = 1.0;
f(i) = X * X + Y * Y + Z * Z;
}
// Solve linear equations to get coefficients
Eigen::VectorXd c = A.colPivHouseholderQr().solve(f);
Standard_Real t = c[0] * c[0] + c[1] * c[1] + c[2] * c[2] + c[3];
theRadius = std::sqrt(t);
theCenter.SetCoord(c[0], c[1], c[2]);
//check that all points belong to the sphere within theTolerance
std::vector<Standard_Real> aDists;
for(const auto& aPnt : thePoints)
{
aDists.push_back(aPnt.Distance(theCenter));
}
if (!aDists.empty())
{
Standard_Real sumDistances = std::accumulate(aDists.begin(), aDists.end(), 0.0);
Standard_Real averageDistance = sumDistances / aDists.size();
if (std::fabs(averageDistance - theRadius) > theTolerance) {
return false;
}
}
return true;
}
//================================================================================
/*!
* \brief Returns a name of creation operation and names and values of creation parameters
*/
//================================================================================
bool GEOMImpl_WrappingDriver::
GetCreationInformation(std::string &theOperationName,
std::vector<GEOM_Param> &theParams)
{
if (Label().IsNull()) return 0;
Handle(GEOM_Function) function = GEOM_Function::GetFunction(Label());
GEOMImpl_IWrap aCI( function );
Standard_Integer aType = function->GetType();
theOperationName = "WRAPPEDFACE";
switch ( aType ) {
case WRAPPING_FACE:
AddParam(theParams, "Edges", aCI.GetEdges());
AddParam(theParams, "Vertices", aCI.GetVertices());
AddParam(theParams, "Tolerance", aCI.GetTolerance());
break;
default:
return false;
}
return true;
}
IMPLEMENT_STANDARD_RTTIEXT(GEOMImpl_WrappingDriver, GEOM_BaseDriver)