kindred/create
1
1
Fork 0
Code Issues 30 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

improve mesh repair functions

Browse Source
This commit is contained in:
wmayer
2019-04-12 21:31:46 +02:00
parent 57287e8ca6
commit 746997e484
10 changed files with 312 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
src/Mod/Mesh/App/Core/Degeneration.cpp
View File

@@ -506,7 +506,7 @@ bool MeshFixDegeneratedFacets::Fixup()
return true;
}
bool MeshRemoveSmallEdges::Fixup()
bool MeshRemoveNeedles::Fixup()
{
typedef std::pair<unsigned long, int> FaceEdge; // (face, edge) pair
typedef std::pair<float, FaceEdge> FaceEdgePriority;
@@ -685,6 +685,77 @@ bool MeshRemoveSmallEdges::Fixup()
// ----------------------------------------------------------------------
bool MeshFixCaps::Fixup()
{
typedef std::pair<unsigned long, int> FaceVertex; // (face, vertex) pair
typedef std::pair<float, FaceVertex> FaceVertexPriority;
MeshTopoAlgorithm topAlg(_rclMesh);
const MeshFacetArray &rclFAry = _rclMesh.GetFacets();
const MeshPointArray &rclPAry = _rclMesh.GetPoints();
std::size_t facetCount = rclFAry.size();
float fCosMaxAngle = static_cast<float>(cos(fMaxAngle));
std::priority_queue<FaceVertexPriority,
std::vector<FaceVertexPriority>,
std::greater<FaceVertexPriority> > 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<unsigned long>(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);

42
src/Mod/Mesh/App/Core/Degeneration.h
View File

@@ -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.

67
src/Mod/Mesh/App/Core/Elements.cpp
View File

@@ -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()));
}

12
src/Mod/Mesh/App/Core/Elements.h
View File

@@ -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. */

20
src/Mod/Mesh/App/FacetPy.xml
View File

@@ -100,15 +100,33 @@ The two angles are given in radian.
</Attribute>
<Attribute Name="AspectRatio" ReadOnly="true">
<Documentation>
<UserDocu>The aspect ratio of the facet</UserDocu>
<UserDocu>The aspect ratio of the facet computed by longest edge and its height</UserDocu>
</Documentation>
<Parameter Name="AspectRatio" Type="Float"/>
</Attribute>
<Attribute Name="AspectRatio2" ReadOnly="true">
<Documentation>
<UserDocu>The aspect ratio of the facet computed by radius of circum-circle and in-circle</UserDocu>
</Documentation>
<Parameter Name="AspectRatio2" Type="Float"/>
</Attribute>
<Attribute Name="Roundness" ReadOnly="true">
<Documentation>
<UserDocu>The roundness of the facet</UserDocu>
</Documentation>
<Parameter Name="Roundness" Type="Float"/>
</Attribute>
<Attribute Name="CircumCircle" ReadOnly="true">
<Documentation>
<UserDocu>The center and radius of the circum-circle</UserDocu>
</Documentation>
<Parameter Name="CircumCircle" Type="Tuple"/>
</Attribute>
<Attribute Name="InCircle" ReadOnly="true">
<Documentation>
<UserDocu>The center and radius of the in-circle</UserDocu>
</Documentation>
<Parameter Name="InCircle" Type="Tuple"/>
</Attribute>
</PythonExport>
</GenerateModel>

41
src/Mod/Mesh/App/FacetPyImp.cpp
View File

@@ -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;

19
src/Mod/Mesh/App/Mesh.cpp
View File

@@ -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();

4
src/Mod/Mesh/App/Mesh.h
View File

@@ -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 */

21
src/Mod/Mesh/App/MeshPy.xml
View File

@@ -136,11 +136,16 @@ mesh.write(Stream=file,Format='STL',[Name='Object name',Material=colors])</UserD
<UserDocu>Remove a list of facet indices from the mesh</UserDocu>
</Documentation>
</Methode>
<Methode Name="removeSmallEdges">
<Methode Name="removeNeedles">
<Documentation>
<UserDocu>Remove all edges that are smaller than a given length</UserDocu>
</Documentation>
</Methode>
<Methode Name="removeFullBoundaryFacets">
<Documentation>
<UserDocu>Remove facets whose all three points are on the boundary</UserDocu>
</Documentation>
</Methode>
<Methode Name="getInternalFacets" Const="true">
<Documentation>
<UserDocu>Builds a list of facet indices with triangles that are inside a volume mesh</UserDocu>
@@ -294,7 +299,12 @@ for c in mesh.getSeparatecomponents():
<UserDocu>Repair any invalid indices</UserDocu>
</Documentation>
</Methode>
<Methode Name="fixDeformations">
<Methode Name="fixCaps">
<Documentation>
<UserDocu>Repair caps by swapping the edge</UserDocu>
</Documentation>
</Methode>
<Methode Name="fixDeformations">
<Documentation>
<UserDocu>Repair deformed facets</UserDocu>
</Documentation>
@@ -427,7 +437,12 @@ smooth([iteration=1,maxError=FLT_MAX])</UserDocu>
</UserDocu>
</Documentation>
</Methode>
<Methode Name="optimizeTopology" Const="true">
<Methode Name="mergeFacets">
<Documentation>
<UserDocu>Merge facets to optimize topology</UserDocu>
</Documentation>
</Methode>
<Methode Name="optimizeTopology" Const="true">
<Documentation>
<UserDocu>Optimize the edges to get nicer facets</UserDocu>
</Documentation>

43
src/Mod/Mesh/App/MeshPyImp.cpp
View File

@@ -28,6 +28,7 @@
#include <Base/Builder3D.h>
#include <Base/GeometryPyCXX.h>
#include <Base/MatrixPy.h>
#include <Base/Tools.h>
#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<float>(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;