ReverseEngineering: approximation of a 2-degree polynomial surface and converting it into a Bezier surface

src/Mod/Mesh/App/Core/Approximation.cpp

@@ -412,6 +412,15 @@ std::vector<Base::Vector3f> PlaneFit::GetLocalPoints() const
     return localPoints;
 }
 
+Base::BoundBox3f PlaneFit::GetBoundings() const
+{
+    Base::BoundBox3f bbox;
+    std::vector<Base::Vector3f> 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<Base::Vector3f>& pts) const
+{
+    Base::Vector3d bs = Base::convertTo<Base::Vector3d>(this->_vBase);
+    Base::Vector3d ex = Base::convertTo<Base::Vector3d>(this->_vDirU);
+    Base::Vector3d ey = Base::convertTo<Base::Vector3d>(this->_vDirV);
+    Base::Vector3d ez = Base::convertTo<Base::Vector3d>(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<Base::Vector3d>& pts) const
+{
+    Base::Vector3d bs = Base::convertTo<Base::Vector3d>(this->_vBase);
+    Base::Vector3d ex = Base::convertTo<Base::Vector3d>(this->_vDirU);
+    Base::Vector3d ey = Base::convertTo<Base::Vector3d>(this->_vDirV);
+    Base::Vector3d ez = Base::convertTo<Base::Vector3d>(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<Base::Vector3d> SurfaceFit::toBezier(double umin, double umax, double vmin, double vmax) const
+{
+    std::vector<Base::Vector3d> 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 {

src/Mod/Mesh/App/Core/Approximation.h

@@ -35,6 +35,7 @@
 
 #include <Base/Vector3D.h>
 #include <Base/Matrix.h>
+#include <Base/BoundBox.h>
 
 namespace Wm4
 {
@@ -233,6 +234,12 @@ public:
      * array is returned.
      */
     std::vector<Base::Vector3f> 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<Base::Vector3f>&) const;
+    void Transform(std::vector<Base::Vector3d>&) 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<Base::Vector3d> toBezier(double umin=0.0, double umax=1.0, double vmin=0.0, double vmax=1.0) const;
 
 protected:
     double PolynomFit();

src/Mod/ReverseEngineering/Gui/Command.cpp

@@ -35,6 +35,8 @@
 #include <Mod/Part/App/TopoShape.h>
 #include <Mod/Part/App/PartFeature.h>
 #include <Mod/Part/App/FaceMakerCheese.h>
+#include <Mod/Part/App/Geometry.h>
+#include <Mod/Part/App/Tools.h>
 #include <Mod/Points/App/Structured.h>
 #include <Mod/Mesh/App/Mesh.h>
 #include <Mod/Mesh/App/MeshFeature.h>
@@ -50,6 +52,7 @@
 #include <Gui/MainWindow.h>
 #include <Gui/FileDialog.h>
 #include <Gui/Selection.h>
+#include <Base/Builder3D.h>
 #include <Base/CoordinateSystem.h>
 #include <Base/Converter.h>
 
@@ -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<Mesh::Feature*> sel = getSelection().getObjectsOfType<Mesh::Feature>();
+    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<Base::Vector3d> 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<gp_Pnt>(poles.at(0)));
+            grid.SetValue(2, 1, Base::convertTo<gp_Pnt>(poles.at(1)));
+            grid.SetValue(3, 1, Base::convertTo<gp_Pnt>(poles.at(2)));
+            grid.SetValue(1, 2, Base::convertTo<gp_Pnt>(poles.at(3)));
+            grid.SetValue(2, 2, Base::convertTo<gp_Pnt>(poles.at(4)));
+            grid.SetValue(3, 2, Base::convertTo<gp_Pnt>(poles.at(5)));
+            grid.SetValue(1, 3, Base::convertTo<gp_Pnt>(poles.at(6)));
+            grid.SetValue(2, 3, Base::convertTo<gp_Pnt>(poles.at(7)));
+            grid.SetValue(3, 3, Base::convertTo<gp_Pnt>(poles.at(8)));
+
+            Handle(Geom_BezierSurface) bezier(new Geom_BezierSurface(grid));
+            Part::Feature* part = static_cast<Part::Feature*>(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());

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;