From 4671993a5dcd751848e39b81fe970825738afcdf Mon Sep 17 00:00:00 2001 From: wmayer Date: Sat, 7 Mar 2020 14:21:57 +0100 Subject: [PATCH] ReverseEngineering: approximation of a 2-degree polynomial surface and converting it into a Bezier surface --- src/Mod/Mesh/App/Core/Approximation.cpp | 132 +++++++++++++++++++ src/Mod/Mesh/App/Core/Approximation.h | 22 ++++ src/Mod/ReverseEngineering/Gui/Command.cpp | 59 +++++++++ src/Mod/ReverseEngineering/Gui/Workbench.cpp | 1 + 4 files changed, 214 insertions(+) diff --git a/src/Mod/Mesh/App/Core/Approximation.cpp b/src/Mod/Mesh/App/Core/Approximation.cpp index 85efdcbef1..97658d35fa 100644 --- a/src/Mod/Mesh/App/Core/Approximation.cpp +++ b/src/Mod/Mesh/App/Core/Approximation.cpp @@ -412,6 +412,15 @@ std::vector PlaneFit::GetLocalPoints() const return localPoints; } +Base::BoundBox3f PlaneFit::GetBoundings() const +{ + Base::BoundBox3f bbox; + std::vector pts = GetLocalPoints(); + for (auto it : pts) + bbox.Add(it); + return bbox; +} + // ------------------------------------------------------------------------------- bool QuadraticFit::GetCurvatureInfo(double x, double y, double z, @@ -795,6 +804,129 @@ void SurfaceFit::GetCoefficients(double& a,double& b,double& c,double& d,double& f = _fCoeff[0]; } +void SurfaceFit::Transform(std::vector& pts) const +{ + Base::Vector3d bs = Base::convertTo(this->_vBase); + Base::Vector3d ex = Base::convertTo(this->_vDirU); + Base::Vector3d ey = Base::convertTo(this->_vDirV); + Base::Vector3d ez = Base::convertTo(this->_vDirW); + + Base::Matrix4D mat; + mat[0][0] = ex.x; + mat[0][1] = ey.x; + mat[0][2] = ez.x; + mat[0][3] = bs.x; + + mat[1][0] = ex.y; + mat[1][1] = ey.y; + mat[1][2] = ez.y; + mat[1][3] = bs.y; + + mat[2][0] = ex.z; + mat[2][1] = ey.z; + mat[2][2] = ez.z; + mat[2][3] = bs.z; + + std::transform(pts.begin(), pts.end(), pts.begin(), [&mat](const Base::Vector3f& pt) { + Base::Vector3f v(pt); + mat.multVec(v, v); + return v; + }); +} + +void SurfaceFit::Transform(std::vector& pts) const +{ + Base::Vector3d bs = Base::convertTo(this->_vBase); + Base::Vector3d ex = Base::convertTo(this->_vDirU); + Base::Vector3d ey = Base::convertTo(this->_vDirV); + Base::Vector3d ez = Base::convertTo(this->_vDirW); + + Base::Matrix4D mat; + mat[0][0] = ex.x; + mat[0][1] = ey.x; + mat[0][2] = ez.x; + mat[0][3] = bs.x; + + mat[1][0] = ex.y; + mat[1][1] = ey.y; + mat[1][2] = ez.y; + mat[1][3] = bs.y; + + mat[2][0] = ex.z; + mat[2][1] = ey.z; + mat[2][2] = ez.z; + mat[2][3] = bs.z; + + std::transform(pts.begin(), pts.end(), pts.begin(), [&mat](const Base::Vector3d& pt) { + Base::Vector3d v(pt); + mat.multVec(v, v); + return v; + }); +} + +/*! + * \brief SurfaceFit::toBezier + * This function computes the Bezier representation of the polynomial surface of the form + * f(x,y) = a*x*x + b*y*y + c*x*y + d*y + e*f + f + * by getting the 3x3 control points. + */ +std::vector SurfaceFit::toBezier(double umin, double umax, double vmin, double vmax) const +{ + std::vector pts; + pts.reserve(9); + + // the Bezier surface is defined by the 3x3 control points + // P11 P21 P31 + // P12 P22 P32 + // P13 P23 P33 + // + // The surface goes through the points P11, P31, P31 and P33 + // To get the four control points P21, P12, P32 and P23 we inverse + // the de-Casteljau algorithm used for Bezier curves of degree 2 + // as we already know the points for the parameters + // (0, 0.5), (0.5, 0), (0.5, 1.0) and (1.0, 0.5) + // To get the control point P22 we inverse the de-Casteljau algorithm + // for the surface point on (0.5, 0.5) + double umid = 0.5 * (umin + umax); + double vmid = 0.5 * (vmin + vmax); + + // first row + double z11 = Value(umin, vmin); + double v21 = Value(umid, vmin); + double z31 = Value(umax, vmin); + double z21 = 2.0 * v21 - 0.5 * (z11 + z31); + + // third row + double z13 = Value(umin, vmax); + double v23 = Value(umid, vmax); + double z33 = Value(umax, vmax); + double z23 = 2.0 * v23 - 0.5 * (z13 + z33); + + // second row + double v12 = Value(umin, vmid); + double z12 = 2.0 * v12 - 0.5 * (z11 + z13); + double v32 = Value(umax, vmid); + double z32 = 2.0 * v32 - 0.5 * (z31 + z33); + double v22 = Value(umid, vmid); + double z22 = 4.0 * v22 - 0.25 * (z11 + z31 + z13 + z33 + 2.0 * (z12 + z21 + z32 + z23)); + + // first row + pts.emplace_back(umin, vmin, z11); + pts.emplace_back(umid, vmin, z21); + pts.emplace_back(umax, vmin, z31); + + // second row + pts.emplace_back(umin, vmid, z12); + pts.emplace_back(umid, vmid, z22); + pts.emplace_back(umax, vmid, z32); + + // third row + pts.emplace_back(umin, vmax, z13); + pts.emplace_back(umid, vmax, z23); + pts.emplace_back(umax, vmax, z33); + return pts; +} + // ------------------------------------------------------------------------------- namespace MeshCore { diff --git a/src/Mod/Mesh/App/Core/Approximation.h b/src/Mod/Mesh/App/Core/Approximation.h index 168eb06e70..67a6760548 100644 --- a/src/Mod/Mesh/App/Core/Approximation.h +++ b/src/Mod/Mesh/App/Core/Approximation.h @@ -35,6 +35,7 @@ #include #include +#include namespace Wm4 { @@ -233,6 +234,12 @@ public: * array is returned. */ std::vector GetLocalPoints() const; + /** + * Returns the local bounding box of the transformed points releative to the + * coordinate system of the plane. If this method is called before the plane is + * computed an invalid bounding box is returned. + */ + Base::BoundBox3f GetBoundings() const; protected: Base::Vector3f _vBase; /**< Base vector of the plane. */ @@ -347,6 +354,21 @@ public: float Fit(); double Value(double x, double y) const; void GetCoefficients(double& a,double& b,double& c,double& d,double& e,double& f) const; + /** + * @brief Transform + * Transforms points from the local coordinate system to the world coordinate system + */ + void Transform(std::vector&) const; + void Transform(std::vector&) const; + /** + * @brief toBezier + * @param umin Parameter range + * @param umax Parameter range + * @param vmin Parameter range + * @param vmax Parameter range + * @return control points of the Bezier surface + */ + std::vector toBezier(double umin=0.0, double umax=1.0, double vmin=0.0, double vmax=1.0) const; protected: double PolynomFit(); diff --git a/src/Mod/ReverseEngineering/Gui/Command.cpp b/src/Mod/ReverseEngineering/Gui/Command.cpp index b6c8bc7295..1c35d21361 100644 --- a/src/Mod/ReverseEngineering/Gui/Command.cpp +++ b/src/Mod/ReverseEngineering/Gui/Command.cpp @@ -35,6 +35,8 @@ #include #include #include +#include +#include #include #include #include @@ -50,6 +52,7 @@ #include #include #include +#include #include #include @@ -296,6 +299,61 @@ bool CmdApproxSphere::isActive(void) return false; } +DEF_STD_CMD_A(CmdApproxPolynomial) + +CmdApproxPolynomial::CmdApproxPolynomial() + : Command("Reen_ApproxPolynomial") +{ + sAppModule = "Reen"; + sGroup = QT_TR_NOOP("Reverse Engineering"); + sMenuText = QT_TR_NOOP("Polynomial surface"); + sToolTipText = QT_TR_NOOP("Approximate a polynomial surface"); + sWhatsThis = "Reen_ApproxPolynomial"; + sStatusTip = sToolTipText; +} + +void CmdApproxPolynomial::activated(int) +{ + std::vector sel = getSelection().getObjectsOfType(); + App::Document* doc = App::GetApplication().getActiveDocument(); + openCommand("Fit polynomial surface"); + for (auto it : sel) { + const Mesh::MeshObject& mesh = it->Mesh.getValue(); + const MeshCore::MeshKernel& kernel = mesh.getKernel(); + MeshCore::SurfaceFit fit; + fit.AddPoints(kernel.GetPoints()); + if (fit.Fit() < FLOAT_MAX) { + Base::BoundBox3f bbox = fit.GetBoundings(); + std::vector poles = fit.toBezier(bbox.MinX, bbox.MaxX, bbox.MinY, bbox.MaxY); + fit.Transform(poles); + + TColgp_Array2OfPnt grid(1, 3, 1, 3); + grid.SetValue(1, 1, Base::convertTo(poles.at(0))); + grid.SetValue(2, 1, Base::convertTo(poles.at(1))); + grid.SetValue(3, 1, Base::convertTo(poles.at(2))); + grid.SetValue(1, 2, Base::convertTo(poles.at(3))); + grid.SetValue(2, 2, Base::convertTo(poles.at(4))); + grid.SetValue(3, 2, Base::convertTo(poles.at(5))); + grid.SetValue(1, 3, Base::convertTo(poles.at(6))); + grid.SetValue(2, 3, Base::convertTo(poles.at(7))); + grid.SetValue(3, 3, Base::convertTo(poles.at(8))); + + Handle(Geom_BezierSurface) bezier(new Geom_BezierSurface(grid)); + Part::Feature* part = static_cast(doc->addObject("Part::Spline", "Bezier")); + part->Shape.setValue(Part::GeomBezierSurface(bezier).toShape()); + } + } + commitCommand(); + updateActive(); +} + +bool CmdApproxPolynomial::isActive(void) +{ + if (getSelection().countObjectsOfType(Mesh::Feature::getClassTypeId()) > 0) + return true; + return false; +} + DEF_STD_CMD_A(CmdSegmentation) CmdSegmentation::CmdSegmentation() @@ -566,6 +624,7 @@ void CreateReverseEngineeringCommands(void) rcCmdMgr.addCommand(new CmdApproxPlane()); rcCmdMgr.addCommand(new CmdApproxCylinder()); rcCmdMgr.addCommand(new CmdApproxSphere()); + rcCmdMgr.addCommand(new CmdApproxPolynomial()); rcCmdMgr.addCommand(new CmdSegmentation()); rcCmdMgr.addCommand(new CmdSegmentationManual()); rcCmdMgr.addCommand(new CmdSegmentationFromComponents()); diff --git a/src/Mod/ReverseEngineering/Gui/Workbench.cpp b/src/Mod/ReverseEngineering/Gui/Workbench.cpp index 0d1b18682d..ab4d343ae9 100644 --- a/src/Mod/ReverseEngineering/Gui/Workbench.cpp +++ b/src/Mod/ReverseEngineering/Gui/Workbench.cpp @@ -79,6 +79,7 @@ Gui::MenuItem* Workbench::setupMenuBar() const *approx << "Reen_ApproxPlane" << "Reen_ApproxCylinder" << "Reen_ApproxSphere" + << "Reen_ApproxPolynomial" << "Separator" << "Reen_ApproxSurface"; *reen << approx;