From 746997e48411f22a150b95b6fdbb9b02cbe741d5 Mon Sep 17 00:00:00 2001 From: wmayer Date: Fri, 12 Apr 2019 21:31:46 +0200 Subject: [PATCH] improve mesh repair functions --- src/Mod/Mesh/App/Core/Degeneration.cpp | 73 +++++++++++++++++++++++++- src/Mod/Mesh/App/Core/Degeneration.h | 42 +++++++++++++-- src/Mod/Mesh/App/Core/Elements.cpp | 67 ++++++++++++++++++----- src/Mod/Mesh/App/Core/Elements.h | 12 ++++- src/Mod/Mesh/App/FacetPy.xml | 20 ++++++- src/Mod/Mesh/App/FacetPyImp.cpp | 41 +++++++++++++++ src/Mod/Mesh/App/Mesh.cpp | 19 ++++++- src/Mod/Mesh/App/Mesh.h | 4 +- src/Mod/Mesh/App/MeshPy.xml | 21 ++++++-- src/Mod/Mesh/App/MeshPyImp.cpp | 43 ++++++++++++++- 10 files changed, 312 insertions(+), 30 deletions(-) diff --git a/src/Mod/Mesh/App/Core/Degeneration.cpp b/src/Mod/Mesh/App/Core/Degeneration.cpp index 378bc0b66a..699d186955 100644 --- a/src/Mod/Mesh/App/Core/Degeneration.cpp +++ b/src/Mod/Mesh/App/Core/Degeneration.cpp @@ -506,7 +506,7 @@ bool MeshFixDegeneratedFacets::Fixup() return true; } -bool MeshRemoveSmallEdges::Fixup() +bool MeshRemoveNeedles::Fixup() { typedef std::pair FaceEdge; // (face, edge) pair typedef std::pair FaceEdgePriority; @@ -685,6 +685,77 @@ bool MeshRemoveSmallEdges::Fixup() // ---------------------------------------------------------------------- +bool MeshFixCaps::Fixup() +{ + typedef std::pair FaceVertex; // (face, vertex) pair + typedef std::pair FaceVertexPriority; + + MeshTopoAlgorithm topAlg(_rclMesh); + const MeshFacetArray &rclFAry = _rclMesh.GetFacets(); + const MeshPointArray &rclPAry = _rclMesh.GetPoints(); + std::size_t facetCount = rclFAry.size(); + + float fCosMaxAngle = static_cast(cos(fMaxAngle)); + + std::priority_queue, + std::greater > todo; + for (std::size_t index = 0; index < facetCount; index++) { + for (int i=0; i<3; i++) { + const MeshFacet& facet = rclFAry[index]; + const Base::Vector3f& p1 = rclPAry[facet._aulPoints[i]]; + const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(i+1)%3]]; + const Base::Vector3f& p3 = rclPAry[facet._aulPoints[(i+2)%3]]; + Base::Vector3f dir1(p2-p1); dir1.Normalize(); + Base::Vector3f dir2(p3-p1); dir2.Normalize(); + + float fCosAngle = dir1.Dot(dir2); + if (fCosAngle < fCosMaxAngle) { + unsigned long facetIndex = static_cast(index); + todo.push(std::make_pair(fCosAngle, std::make_pair(facetIndex, i))); + } + } + } + + while (!todo.empty()) { + FaceVertex facevertex = todo.top().second; + todo.pop(); + + // the facet points may have changed, so check the current distance again + const MeshFacet& facet = rclFAry[facevertex.first]; + const Base::Vector3f& p1 = rclPAry[facet._aulPoints[facevertex.second]]; + const Base::Vector3f& p2 = rclPAry[facet._aulPoints[(facevertex.second+1)%3]]; + const Base::Vector3f& p3 = rclPAry[facet._aulPoints[(facevertex.second+2)%3]]; + Base::Vector3f dir1(p2-p1); dir1.Normalize(); + Base::Vector3f dir2(p3-p1); dir2.Normalize(); + + // check that the criterion is still OK in case + // an earlier edge-swap has an impact + float fCosAngle = dir1.Dot(dir2); + if (fCosAngle >= fCosMaxAngle) + continue; + + // the triangle shouldn't be a needle, therefore the projection of the point with + // the maximum angle must have a clear distance to the other corner points + // as factor we choose a default value of 25% of the corresponding edge length + Base::Vector3f p4 = p1.Perpendicular(p2, p3-p2); + float distP2P3 = Base::Distance(p2,p3); + float distP2P4 = Base::Distance(p2,p4); + float distP3P4 = Base::Distance(p3,p4); + if (distP2P4/distP2P3 < fSplitFactor || distP3P4/distP2P3 < fSplitFactor) + continue; + + unsigned long facetpos = facevertex.first; + unsigned long neighbour = rclFAry[facetpos]._aulNeighbours[(facevertex.second+1)%3]; + if (neighbour != ULONG_MAX) + topAlg.SwapEdge(facetpos, neighbour); + } + + return true; +} + +// ---------------------------------------------------------------------- + bool MeshEvalDeformedFacets::Evaluate() { float fCosMinAngle = cos(fMinAngle); diff --git a/src/Mod/Mesh/App/Core/Degeneration.h b/src/Mod/Mesh/App/Core/Degeneration.h index 2d7f1351e3..442ab16344 100644 --- a/src/Mod/Mesh/App/Core/Degeneration.h +++ b/src/Mod/Mesh/App/Core/Degeneration.h @@ -330,22 +330,25 @@ private: }; /** - * The MeshRemoveSmallEdges class tries to fix degenerations by removing small edges. + * The MeshRemoveNeedles class tries to fix degenerations by removing needles. + * Needles are triangles where its longest edge is much longer than its shortest edge. + * https://graphics.uni-bielefeld.de/publications/vmv01.pdf * @see MeshFixDegeneratedFacets + * @see MeshFixCaps * @author Werner Mayer */ -class MeshExport MeshRemoveSmallEdges : public MeshValidation +class MeshExport MeshRemoveNeedles : public MeshValidation { public: /** * Construction. */ - MeshRemoveSmallEdges (MeshKernel &rclM, float fMinEdgeLen = MeshDefinitions::_fMinPointDistance) + MeshRemoveNeedles (MeshKernel &rclM, float fMinEdgeLen = MeshDefinitions::_fMinPointDistance) : MeshValidation(rclM), fMinEdgeLength(fMinEdgeLen) { } /** * Destruction. */ - ~MeshRemoveSmallEdges () { } + ~MeshRemoveNeedles () { } /** * Removes all facets with an edge smaller than \a fMinEdgeLength without leaving holes or gaps * in the mesh. @@ -356,6 +359,37 @@ private: float fMinEdgeLength; }; +/** + * The MeshFixCaps class tries to fix degenerations by swapping the common edge of a cap + * and its neighbour. + * Caps are triangles with one angle close to 180 degree. The definitions of caps and needles + * are not mutually exclusive but here we only consider triangles that are caps but no needles. + * https://graphics.uni-bielefeld.de/publications/vmv01.pdf + * @see MeshFixDegeneratedFacets + * @see MeshRemoveNeedles + * @author Werner Mayer + */ +class MeshExport MeshFixCaps : public MeshValidation +{ +public: + /** + * Construction. The \arg fFactor must be in the range of 0.0 and 0.5. + */ + MeshFixCaps (MeshKernel &rclM, float fMaxAng = 2.61f, float fFactor = 0.25f) // ~150 degree + : MeshValidation(rclM), fMaxAngle(fMaxAng), fSplitFactor(fFactor) { } + /** + * Destruction. + */ + ~MeshFixCaps () { } + /** + */ + bool Fixup (); + +private: + float fMaxAngle; + float fSplitFactor; +}; + /** * The MeshEvalDeformedFacets class searches for deformed facets. A facet is regarded as deformed * if an angle is < 30 deg or > 120 deg. diff --git a/src/Mod/Mesh/App/Core/Elements.cpp b/src/Mod/Mesh/App/Core/Elements.cpp index db2ea6bcdc..57781133d6 100644 --- a/src/Mod/Mesh/App/Core/Elements.cpp +++ b/src/Mod/Mesh/App/Core/Elements.cpp @@ -1187,6 +1187,21 @@ float MeshGeomFacet::MaximumAngle () const return fMaxAngle; } +float MeshGeomFacet::MinimumAngle () const +{ + float fMinAngle = Mathf::PI; + + for ( int i=0; i<3; i++ ) { + Base::Vector3f dir1(_aclPoints[(i+1)%3]-_aclPoints[i]); + Base::Vector3f dir2(_aclPoints[(i+2)%3]-_aclPoints[i]); + float fAngle = dir1.GetAngle(dir2); + if (fAngle < fMinAngle) + fMinAngle = fAngle; + } + + return fMinAngle; +} + bool MeshGeomFacet::IsPointOfSphere(const Base::Vector3f& rP) const { float radius; @@ -1216,20 +1231,44 @@ bool MeshGeomFacet::IsPointOfSphere(const MeshGeomFacet& rFacet) const float MeshGeomFacet::AspectRatio() const { - Base::Vector3f center; - float radius = CenterOfCircumCircle(center); + Base::Vector3f d0 = _aclPoints[0] - _aclPoints[1]; + Base::Vector3f d1 = _aclPoints[1] - _aclPoints[2]; + Base::Vector3f d2 = _aclPoints[2] - _aclPoints[0]; - float minLength = FLOAT_MAX; - for (int i=0; i<3; i++) { - const Base::Vector3f& p0 = _aclPoints[i]; - const Base::Vector3f& p1 = _aclPoints[(i+1)%3]; - float distance = Base::Distance(p0, p1); - if (distance < minLength) { - minLength = distance; - } - } + float l2, maxl2 = d0.Sqr(); + if ((l2=d1.Sqr()) > maxl2) + maxl2 = l2; - if (minLength == 0) - return FLOAT_MAX; - return radius / minLength; + d1 = d2; + if ((l2=d1.Sqr()) > maxl2) + maxl2 = l2; + + // squared area of the parallelogram spanned by d0 and d1 + float a2 = (d0 % d1).Sqr(); + return float(sqrt( (maxl2 * maxl2) / a2 )); +} + +float MeshGeomFacet::AspectRatio2() const +{ + const Base::Vector3f& rcP1 = _aclPoints[0]; + const Base::Vector3f& rcP2 = _aclPoints[1]; + const Base::Vector3f& rcP3 = _aclPoints[2]; + + float a = Base::Distance(rcP1, rcP2); + float b = Base::Distance(rcP2, rcP3); + float c = Base::Distance(rcP3, rcP1); + + // https://stackoverflow.com/questions/10289752/aspect-ratio-of-a-triangle-of-a-meshed-surface + return a * b * c / ((b + c - a) * (c + a - b) * (a + b - c)); +} + +float MeshGeomFacet::Roundness() const +{ + const double FOUR_ROOT3 = 6.928203230275509; + double area = Area(); + Base::Vector3f d0 = _aclPoints[0] - _aclPoints[1]; + Base::Vector3f d1 = _aclPoints[1] - _aclPoints[2]; + Base::Vector3f d2 = _aclPoints[2] - _aclPoints[0]; + + return (float) (FOUR_ROOT3 * area / (d0.Sqr() + d1.Sqr() + d2.Sqr())); } diff --git a/src/Mod/Mesh/App/Core/Elements.h b/src/Mod/Mesh/App/Core/Elements.h index b5a138d68f..ead0169d97 100644 --- a/src/Mod/Mesh/App/Core/Elements.h +++ b/src/Mod/Mesh/App/Core/Elements.h @@ -369,7 +369,7 @@ public: */ bool IsDegenerated(float epsilon) const; /** - * Checks whether the triangle is deformed. A triangle is deformed if the an angle + * Checks whether the triangle is deformed. A triangle is deformed if an angle * exceeds a given maximum angle or falls below a given minimum angle. * For performance reasons the cosine of minimum and maximum angle is expected. */ @@ -415,6 +415,8 @@ public: inline float Area () const; /** Calculates the maximum angle of a facet. */ float MaximumAngle () const; + /** Calculates the minimum angle of a facet. */ + float MinimumAngle () const; /** Checks if the facet is inside the bounding box or intersects with it. */ inline bool ContainedByOrIntersectBoundingBox (const Base::BoundBox3f &rcBB) const; /** Checks if the facet intersects with the given bounding box. */ @@ -484,9 +486,15 @@ public: * If one of the facet's points is inside the sphere true is returned, otherwise false. */ bool IsPointOfSphere(const MeshGeomFacet& rFacet) const; - /** The aspect ratio is the ratio of the radius of the circum-circle to the shortest edge length. + /** The aspect ratio is the longest edge length divided by its height. */ float AspectRatio() const; + /** The alternative aspect ration is the ratio of the radius of the circum-circle and twice the radius of the in-circle. + */ + float AspectRatio2() const; + /** The roundness is in the range between 0.0 (colinear) and 1.0 (equilateral). + */ + float Roundness() const; protected: Base::Vector3f _clNormal; /**< Normal of the facet. */ diff --git a/src/Mod/Mesh/App/FacetPy.xml b/src/Mod/Mesh/App/FacetPy.xml index 523b617f71..c20f970447 100644 --- a/src/Mod/Mesh/App/FacetPy.xml +++ b/src/Mod/Mesh/App/FacetPy.xml @@ -100,15 +100,33 @@ The two angles are given in radian. - The aspect ratio of the facet + The aspect ratio of the facet computed by longest edge and its height + + + The aspect ratio of the facet computed by radius of circum-circle and in-circle + + + + + + The roundness of the facet + + + The center and radius of the circum-circle + + + The center and radius of the in-circle + + + diff --git a/src/Mod/Mesh/App/FacetPyImp.cpp b/src/Mod/Mesh/App/FacetPyImp.cpp index 53402518c4..dba93fb20e 100644 --- a/src/Mod/Mesh/App/FacetPyImp.cpp +++ b/src/Mod/Mesh/App/FacetPyImp.cpp @@ -233,6 +233,30 @@ Py::Float FacetPy::getAspectRatio(void) const return Py::Float(tria.AspectRatio()); } +Py::Float FacetPy::getAspectRatio2(void) const +{ + FacetPy::PointerType face = this->getFacetPtr(); + if (!face->isBound()) { + return Py::Float(-1.0); + } + + const MeshCore::MeshKernel& kernel = face->Mesh->getKernel(); + MeshCore::MeshGeomFacet tria = kernel.GetFacet(face->Index); + return Py::Float(tria.AspectRatio2()); +} + +Py::Float FacetPy::getRoundness(void) const +{ + FacetPy::PointerType face = this->getFacetPtr(); + if (!face->isBound()) { + return Py::Float(-1.0); + } + + const MeshCore::MeshKernel& kernel = face->Mesh->getKernel(); + MeshCore::MeshGeomFacet tria = kernel.GetFacet(face->Index); + return Py::Float(tria.Roundness()); +} + Py::Tuple FacetPy::getCircumCircle(void) const { FacetPy::PointerType face = this->getFacetPtr(); @@ -250,6 +274,23 @@ Py::Tuple FacetPy::getCircumCircle(void) const return tuple; } +Py::Tuple FacetPy::getInCircle(void) const +{ + FacetPy::PointerType face = this->getFacetPtr(); + if (!face->isBound()) { + return Py::None(); + } + + const MeshCore::MeshKernel& kernel = face->Mesh->getKernel(); + MeshCore::MeshGeomFacet tria = kernel.GetFacet(face->Index); + Base::Vector3f center; + float radius = tria.CenterOfInscribedCircle(center); + Py::Tuple tuple(2); + tuple.setItem(0, Py::Vector(center)); + tuple.setItem(1, Py::Float(radius)); + return tuple; +} + PyObject *FacetPy::getCustomAttributes(const char* /*attr*/) const { return 0; diff --git a/src/Mod/Mesh/App/Mesh.cpp b/src/Mod/Mesh/App/Mesh.cpp index 27fe0ca998..70d21a0ab0 100644 --- a/src/Mod/Mesh/App/Mesh.cpp +++ b/src/Mod/Mesh/App/Mesh.cpp @@ -1128,15 +1128,21 @@ void MeshObject::refine() this->_segments.clear(); } -void MeshObject::removeSmallEdges(float length) +void MeshObject::removeNeedles(float length) { unsigned long count = _kernel.CountFacets(); - MeshCore::MeshRemoveSmallEdges eval(_kernel, length); + MeshCore::MeshRemoveNeedles eval(_kernel, length); eval.Fixup(); if (_kernel.CountFacets() < count) this->_segments.clear(); } +void MeshObject::validateCaps(float fMaxAngle, float fSplitFactor) +{ + MeshCore::MeshFixCaps eval(_kernel, fMaxAngle, fSplitFactor); + eval.Fixup(); +} + void MeshObject::optimizeTopology(float fMaxAngle) { MeshCore::MeshTopoAlgorithm topalg(_kernel); @@ -1385,6 +1391,15 @@ void MeshObject::removeInvalidPoints() deletePoints(nan.GetIndices()); } +void MeshObject::mergeFacets() +{ + unsigned long count = _kernel.CountFacets(); + MeshCore::MeshFixMergeFacets merge(_kernel); + merge.Fixup(); + if (_kernel.CountFacets() < count) + this->_segments.clear(); +} + void MeshObject::validateIndices() { unsigned long count = _kernel.CountFacets(); diff --git a/src/Mod/Mesh/App/Mesh.h b/src/Mod/Mesh/App/Mesh.h index 2d00ab9b27..cd151d9ed0 100644 --- a/src/Mod/Mesh/App/Mesh.h +++ b/src/Mod/Mesh/App/Mesh.h @@ -255,7 +255,7 @@ public: /** @name Topological operations */ //@{ void refine(); - void removeSmallEdges(float); + void removeNeedles(float); void optimizeTopology(float); void optimizeEdges(); void splitEdges(); @@ -275,6 +275,7 @@ public: void flipNormals(); void harmonizeNormals(); void validateIndices(); + void validateCaps(float fMaxAngle, float fSplitFactor); void validateDeformations(float fMaxAngle, float fEps); void validateDegenerations(float fEps); void removeDuplicatedPoints(); @@ -289,6 +290,7 @@ public: void removeFullBoundaryFacets(); bool hasInvalidPoints() const; void removeInvalidPoints(); + void mergeFacets(); //@} /** @name Mesh segments */ diff --git a/src/Mod/Mesh/App/MeshPy.xml b/src/Mod/Mesh/App/MeshPy.xml index b67f11fd6d..80df47a629 100644 --- a/src/Mod/Mesh/App/MeshPy.xml +++ b/src/Mod/Mesh/App/MeshPy.xml @@ -136,11 +136,16 @@ mesh.write(Stream=file,Format='STL',[Name='Object name',Material=colors])Remove a list of facet indices from the mesh - + Remove all edges that are smaller than a given length + + + Remove facets whose all three points are on the boundary + + Builds a list of facet indices with triangles that are inside a volume mesh @@ -294,7 +299,12 @@ for c in mesh.getSeparatecomponents(): Repair any invalid indices - + + + Repair caps by swapping the edge + + + Repair deformed facets @@ -427,7 +437,12 @@ smooth([iteration=1,maxError=FLT_MAX]) - + + + Merge facets to optimize topology + + + Optimize the edges to get nicer facets diff --git a/src/Mod/Mesh/App/MeshPyImp.cpp b/src/Mod/Mesh/App/MeshPyImp.cpp index dc4538465b..b885fda3ff 100644 --- a/src/Mod/Mesh/App/MeshPyImp.cpp +++ b/src/Mod/Mesh/App/MeshPyImp.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include "Mesh.h" #include "MeshPy.h" @@ -1191,6 +1192,20 @@ PyObject* MeshPy::fixIndices(PyObject *args) Py_Return; } +PyObject* MeshPy::fixCaps(PyObject *args) +{ + float fMaxAngle = Base::toRadians(150.0f); + float fSplitFactor = 0.25f; + if (!PyArg_ParseTuple(args, "|ff", &fMaxAngle, &fSplitFactor)) + return NULL; + + PY_TRY { + getMeshObjectPtr()->validateCaps(fMaxAngle, fSplitFactor); + } PY_CATCH; + + Py_Return; +} + PyObject* MeshPy::fixDeformations(PyObject *args) { float fMaxAngle; @@ -1254,14 +1269,38 @@ PyObject* MeshPy::refine(PyObject *args) Py_Return; } -PyObject* MeshPy::removeSmallEdges(PyObject *args) +PyObject* MeshPy::removeNeedles(PyObject *args) { float length; if (!PyArg_ParseTuple(args, "f", &length)) return NULL; PY_TRY { - getMeshObjectPtr()->removeSmallEdges(length); + getMeshObjectPtr()->removeNeedles(length); + } PY_CATCH; + + Py_Return; +} + +PyObject* MeshPy::removeFullBoundaryFacets(PyObject *args) +{ + if (!PyArg_ParseTuple(args, "")) + return NULL; + + PY_TRY { + getMeshObjectPtr()->removeFullBoundaryFacets(); + } PY_CATCH; + + Py_Return; +} + +PyObject* MeshPy::mergeFacets(PyObject *args) +{ + if (!PyArg_ParseTuple(args, "")) + return NULL; + + PY_TRY { + getMeshObjectPtr()->mergeFacets(); } PY_CATCH; Py_Return;