implement trimming of a mesh with a plane

src/Mod/Mesh/App/CMakeLists.txt

@@ -88,6 +88,8 @@ SET(Core_SRCS
     Core/Triangulation.h
     Core/Trim.cpp
     Core/Trim.h
+    Core/TrimByPlane.cpp
+    Core/TrimByPlane.h
     Core/tritritest.h
     Core/Visitor.cpp
     Core/Visitor.h

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

@@ -268,6 +268,25 @@ bool MeshGeomEdge::IntersectWithLine (const Base::Vector3f &rclPt,
     return dist2 + dist3 <= dist1 + eps;
 }
 
+bool MeshGeomEdge::IntersectWithPlane (const Base::Vector3f &rclPt,
+                                       const Base::Vector3f &rclDir,
+                                       Base::Vector3f &rclRes) const
+{
+    float dist1 = _aclPoints[0].DistanceToPlane(rclPt, rclDir);
+    float dist2 = _aclPoints[1].DistanceToPlane(rclPt, rclDir);
+
+    // either both points are below or above the plane
+    if (dist1 * dist2 >= 0.0f)
+        return false;
+
+    Base::Vector3f u = _aclPoints[1] - _aclPoints[0];
+    Base::Vector3f b = rclPt - _aclPoints[0];
+    float t = b.Dot(rclDir) / u.Dot(rclDir);
+    rclRes = _aclPoints[0] + t * u;
+
+    return true;
+}
+
 void MeshGeomEdge::ProjectPointToLine (const Base::Vector3f &rclPoint,
                                        Base::Vector3f &rclProj) const
 {

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

@@ -168,6 +168,10 @@ class MeshExport MeshGeomEdge
    * with the edge. The intersection must be inside the edge. If there is no intersection false is returned.
    */
   bool IntersectWithLine (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, Base::Vector3f &rclRes) const;
+  /** Calculates the intersection point of the plane defined by the base \a rclPt and the direction \a rclDir
+   * with the edge. The intersection must be inside the edge. If there is no intersection false is returned.
+   */
+  bool IntersectWithPlane (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, Base::Vector3f &rclRes) const;
   /**
    * Calculates the projection of a point onto the line defined by the edge. The caller must check if
    * the projection point is inside the edge.

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

@@ -0,0 +1,168 @@
+/***************************************************************************
+ *   Copyright (c) 2019 Werner Mayer <wmayer[at]users.sourceforge.net>     *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 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 Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "PreCompiled.h"
+#include <algorithm>
+
+#include "TrimByPlane.h"
+#include "Grid.h"
+#include "Iterator.h"
+
+using namespace MeshCore;
+
+MeshTrimByPlane::MeshTrimByPlane(MeshKernel &rclM)
+  : myMesh(rclM)
+{
+}
+
+MeshTrimByPlane::~MeshTrimByPlane()
+{
+}
+
+void MeshTrimByPlane::CheckFacets(const MeshFacetGrid& rclGrid, const Base::Vector3f& base, const Base::Vector3f& normal,
+                                  std::vector<unsigned long> &trimFacets, std::vector<unsigned long>& removeFacets) const
+{
+    // Go through the grid and check for each cell if its bounding box intersects the plane.
+    // If the box is completely below the plane all facets will be kept, if it's above the
+    // plane all triangles will be removed.
+    std::vector<unsigned long> checkElements;
+    MeshGridIterator clGridIter(rclGrid);
+    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
+        Base::BoundBox3f clBBox3d = clGridIter.GetBoundBox();
+        if (clBBox3d.IsCutPlane(base, normal)) {
+            // save all elements in checkElements
+            clGridIter.GetElements(checkElements);
+        }
+        else if (clBBox3d.CalcPoint(0).DistanceToPlane(base, normal) > 0.0f) {
+            // save all elements in removeFacets
+            clGridIter.GetElements(removeFacets);
+        }
+    }
+
+    // remove double elements
+    std::sort(checkElements.begin(), checkElements.end());
+    checkElements.erase(std::unique(checkElements.begin(), checkElements.end()), checkElements.end());
+
+    trimFacets.reserve(checkElements.size()/2); // reserve some memory
+    for (auto it = checkElements.begin(); it != checkElements.end(); ++it) {
+        MeshGeomFacet clFacet = myMesh.GetFacet(*it);
+        if (clFacet.IntersectWithPlane(base, normal)) {
+            trimFacets.push_back(*it);
+            removeFacets.push_back(*it);
+        }
+        else if (clFacet._aclPoints[0].DistanceToPlane(base, normal) > 0.0f) {
+            removeFacets.push_back(*it);
+        }
+    }
+
+    // remove double elements
+    std::sort(removeFacets.begin(), removeFacets.end());
+    removeFacets.erase(std::unique(removeFacets.begin(), removeFacets.end()), removeFacets.end());
+}
+
+void MeshTrimByPlane::CreateOneFacet(const Base::Vector3f& base, const Base::Vector3f& normal, unsigned short shift,
+                                     const MeshGeomFacet& facet, std::vector<MeshGeomFacet>& trimmedFacets) const
+{
+    unsigned short nul = shift % 3;
+    unsigned short one = (shift + 1) % 3;
+    unsigned short two = (shift + 2) % 3;
+
+    Base::Vector3f p1, p2;
+    MeshGeomEdge edge;
+
+    edge._aclPoints[0] = facet._aclPoints[nul];
+    edge._aclPoints[1] = facet._aclPoints[one];
+    edge.IntersectWithPlane(base, normal, p1);
+
+    edge._aclPoints[0] = facet._aclPoints[nul];
+    edge._aclPoints[1] = facet._aclPoints[two];
+    edge.IntersectWithPlane(base, normal, p2);
+
+    MeshGeomFacet create;
+    create._aclPoints[0] = facet._aclPoints[nul];
+    create._aclPoints[1] = p1;
+    create._aclPoints[2] = p2;
+    trimmedFacets.push_back(create);
+}
+
+void MeshTrimByPlane::CreateTwoFacet(const Base::Vector3f& base, const Base::Vector3f& normal, unsigned short shift,
+                                     const MeshGeomFacet& facet, std::vector<MeshGeomFacet>& trimmedFacets) const
+{
+    unsigned short nul = shift % 3;
+    unsigned short one = (shift + 1) % 3;
+    unsigned short two = (shift + 2) % 3;
+
+    Base::Vector3f p1, p2;
+    MeshGeomEdge edge;
+
+    edge._aclPoints[0] = facet._aclPoints[nul];
+    edge._aclPoints[1] = facet._aclPoints[two];
+    edge.IntersectWithPlane(base, normal, p1);
+
+    edge._aclPoints[0] = facet._aclPoints[one];
+    edge._aclPoints[1] = facet._aclPoints[two];
+    edge.IntersectWithPlane(base, normal, p2);
+
+    MeshGeomFacet create;
+    create._aclPoints[0] = facet._aclPoints[nul];
+    create._aclPoints[1] = facet._aclPoints[one];
+    create._aclPoints[2] = p1;
+    trimmedFacets.push_back(create);
+
+    create._aclPoints[0] = facet._aclPoints[one];
+    create._aclPoints[1] = p2;
+    create._aclPoints[2] = p1;
+    trimmedFacets.push_back(create);
+}
+
+void MeshTrimByPlane::TrimFacets(const std::vector<unsigned long>& trimFacets, const Base::Vector3f& base,
+                                 const Base::Vector3f& normal, std::vector<MeshGeomFacet>& trimmedFacets)
+{
+    trimmedFacets.reserve(2 * trimFacets.size());
+    for (auto it = trimFacets.begin(); it != trimFacets.end(); ++it) {
+        MeshGeomFacet facet = myMesh.GetFacet(*it);
+        float dist1 = facet._aclPoints[0].DistanceToPlane(base, normal);
+        float dist2 = facet._aclPoints[1].DistanceToPlane(base, normal);
+        float dist3 = facet._aclPoints[2].DistanceToPlane(base, normal);
+
+        // only one point below
+        if (dist1 < 0.0f && dist2 > 0.0f && dist3 > 0.0f) {
+            CreateOneFacet(base, normal, 0, facet, trimmedFacets);
+        }
+        else if (dist1 > 0.0f && dist2 < 0.0f && dist3 > 0.0f) {
+            CreateOneFacet(base, normal, 1, facet, trimmedFacets);
+        }
+        else if (dist1 > 0.0f && dist2 > 0.0f && dist3 < 0.0f) {
+            CreateOneFacet(base, normal, 2, facet, trimmedFacets);
+        }
+        // two points below
+        else if (dist1 < 0.0f && dist2 < 0.0f && dist3 > 0.0f) {
+            CreateTwoFacet(base, normal, 0, facet, trimmedFacets);
+        }
+        else if (dist1 > 0.0f && dist2 < 0.0f && dist3 < 0.0f) {
+            CreateTwoFacet(base, normal, 1, facet, trimmedFacets);
+        }
+        else if (dist1 < 0.0f && dist2 > 0.0f && dist3 < 0.0f) {
+            CreateTwoFacet(base, normal, 2, facet, trimmedFacets);
+        }
+    }
+}

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

@@ -0,0 +1,66 @@
+/***************************************************************************
+ *   Copyright (c) 2019 Werner Mayer <wmayer[at]users.sourceforge.net>     *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 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 Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef MESHTRIM_BY_PLANE_H
+#define MESHTRIM_BY_PLANE_H
+
+#include <Mod/Mesh/App/Core/Elements.h>
+#include <Mod/Mesh/App/Core/MeshKernel.h>
+
+namespace MeshCore {
+
+/**
+ * Trim the the facets in 3D with a plane
+ * \author Werner Mayer
+ */
+class MeshExport MeshTrimByPlane
+{
+public:
+    MeshTrimByPlane(MeshKernel& mesh);
+    ~MeshTrimByPlane();
+
+public:
+    /**
+     * Checks all facets for intersection with the plane and writes all touched facets into the vector
+     */
+    void CheckFacets(const MeshFacetGrid& rclGrid, const Base::Vector3f& base, const Base::Vector3f& normal,
+                     std::vector<unsigned long>& trimFacets, std::vector<unsigned long>& removeFacets) const;
+
+    /**
+     * The facets from \a trimFacets will be trimmed or deleted and \a trimmedFacets holds the newly generated facets
+     */
+    void TrimFacets(const std::vector<unsigned long>& trimFacets, const Base::Vector3f& base,
+                    const Base::Vector3f& normal, std::vector<MeshGeomFacet>& trimmedFacets);
+
+private:
+    void CreateOneFacet(const Base::Vector3f& base, const Base::Vector3f& normal, unsigned short shift,
+                        const MeshGeomFacet& facet, std::vector<MeshGeomFacet>& trimmedFacets) const;
+    void CreateTwoFacet(const Base::Vector3f& base, const Base::Vector3f& normal, unsigned short shift,
+                        const MeshGeomFacet& facet, std::vector<MeshGeomFacet>& trimmedFacets) const;
+
+private:
+    MeshKernel& myMesh;
+};
+
+} //namespace MeshCore
+
+#endif //MESHTRIM_BY_PLANE_H

src/Mod/Mesh/App/Mesh.cpp

@@ -50,6 +50,7 @@
 #include "Core/SetOperations.h"
 #include "Core/Triangulation.h"
 #include "Core/Trim.h"
+#include "Core/TrimByPlane.h"
 #include "Core/Visitor.h"
 #include "Core/Decimation.h"
 
@@ -60,7 +61,7 @@ using namespace Mesh;
 
 float MeshObject::Epsilon = 1.0e-5f;
 
-TYPESYSTEM_SOURCE(Mesh::MeshObject, Data::ComplexGeoData);
+TYPESYSTEM_SOURCE(Mesh::MeshObject, Data::ComplexGeoData)
 
 MeshObject::MeshObject()
 {
@@ -1053,6 +1054,21 @@ void MeshObject::trim(const Base::Polygon2d& polygon2d,
         this->_kernel.AddFacets(triangle);
 }
 
+void MeshObject::trim(const Base::Vector3f& base, const Base::Vector3f& normal)
+{
+    MeshCore::MeshTrimByPlane trim(this->_kernel);
+    std::vector<unsigned long> trimFacets, removeFacets;
+    std::vector<MeshCore::MeshGeomFacet> triangle;
+
+    MeshCore::MeshFacetGrid meshGrid(this->_kernel);
+    trim.CheckFacets(meshGrid, base, normal, trimFacets, removeFacets);
+    trim.TrimFacets(trimFacets, base, normal, triangle);
+    if (!removeFacets.empty())
+        this->deleteFacets(removeFacets);
+    if (!triangle.empty())
+        this->_kernel.AddFacets(triangle);
+}
+
 MeshObject* MeshObject::unite(const MeshObject& mesh) const
 {
     MeshCore::MeshKernel result;

src/Mod/Mesh/App/Mesh.h

@@ -223,6 +223,7 @@ class MeshExport MeshObject : public Data::ComplexGeoData
                        float fMinEps = 1.0e-2f, bool bConnectPolygons = false) const;
     void cut(const Base::Polygon2d& polygon, const Base::ViewProjMethod& proj, CutType);
     void trim(const Base::Polygon2d& polygon, const Base::ViewProjMethod& proj, CutType);
+    void trim(const Base::Vector3f& base, const Base::Vector3f& normal);
     //@}
 
     /** @name Selection */