resolveMeshDegenerations redesigned to reduce triangle aspect ratio

source/MRMesh/MRMeshDecimate.cpp

@@ -193,17 +193,26 @@ MR::QuadraticForm3f computeFormAtVertex( const MR::MeshPart & mp, MR::VertId v,
     return qf;
 }
 
-bool resolveMeshDegenerations( MR::Mesh& mesh, int maxIters, float maxDeviation )
+bool resolveMeshDegenerations( MR::Mesh& mesh, int maxIters, float maxDeviation, float maxAngleChange, float criticalAspectRatio )
 {
     MR_TIMER;
     bool meshChanged = false;
     for( int i = 0; i < maxIters; ++i )
     {
-        bool changedThisIter = makeDeloneEdgeFlips( mesh, { .maxDeviationAfterFlip = maxDeviation }, 5 ) > 0;
+        DeloneSettings delone
+        {
+            .maxDeviationAfterFlip = maxDeviation,
+            .maxAngleChange = maxAngleChange,
+            .criticalTriAspectRatio = criticalAspectRatio
+        };
+        bool changedThisIter = makeDeloneEdgeFlips( mesh, delone, 5 ) > 0;
 
-        DecimateSettings settings;
-        settings.maxError = maxDeviation;
-        changedThisIter = decimateMesh( mesh, settings ).vertsDeleted > 0 || changedThisIter;
+        DecimateSettings decimate
+        {
+            .maxError = maxDeviation,
+            .criticalTriAspectRatio = criticalAspectRatio
+        };
+        changedThisIter = decimateMesh( mesh, decimate ).vertsDeleted > 0 || changedThisIter;
         meshChanged = meshChanged || changedThisIter;
         if ( !changedThisIter )
             break;
@@ -393,12 +402,15 @@ VertId MeshDecimator::collapse_( EdgeId edgeToCollapse, const Vector3f & collaps
         return {}; // new triangle aspect ratio would be larger than all of old triangle aspect ratios and larger than allowed in settings
 
     // checks that all new normals are consistent (do not check for degenerate edges)
-    if ( ( po != pd ) || ( po != collapsePos ) )
+    if ( maxOldAspectRatio < settings_.criticalTriAspectRatio )
     {
-        auto n = Vector3f{ sumDblArea_.normalized() };
-        for ( const auto da : triDblAreas_ )
-            if ( dot( da, n ) < 0 )
-                return {};
+        if ( ( po != pd ) || ( po != collapsePos ) )
+        {
+            auto n = Vector3f{ sumDblArea_.normalized() };
+            for ( const auto da : triDblAreas_ )
+                if ( dot( da, n ) < 0 )
+                    return {};
+        }
     }
 
     if ( settings_.preCollapse && !settings_.preCollapse( edgeToCollapse, collapsePos ) )

source/MRMesh/MRMeshDecimate.h

@@ -3,6 +3,7 @@
 #include "MRMeshFwd.h"
 #include "MRProgressCallback.h"
 #include "MRConstants.h"
+#include <cfloat>
 #include <climits>
 #include <functional>
 
@@ -39,6 +40,8 @@ struct DecimateSettings
     float maxEdgeLen = 1;
     /// Maximal possible aspect ratio of a triangle introduced during decimation
     float maxTriangleAspectRatio = 20;
+    /// the algorithm will try to eliminate triangles with equal or larger aspect ratio, ignoring normal orientation checks
+    float criticalTriAspectRatio = FLT_MAX;
     /// Small stabilizer is important to achieve good results on completely planar mesh parts,
     /// if your mesh is not-planer everywhere, then you can set it to zero
     float stabilizer = 0.001f;
@@ -130,7 +133,7 @@ MRMESH_API QuadraticForm3f computeFormAtVertex( const MeshPart & mp, VertId v, f
  * 
  * \sa \ref decimateMesh
  */
-MRMESH_API bool resolveMeshDegenerations( Mesh& mesh, int maxIters = 1, float maxDeviation = 0 );
+MRMESH_API bool resolveMeshDegenerations( Mesh& mesh, int maxIters = 1, float maxDeviation = 0, float maxAngleChange = PI_F / 3, float criticalAspectRatio = 1e7 );
 
 struct RemeshSettings
 {

source/MRMesh/MRMeshDecimateParallel.cpp

@@ -19,6 +19,7 @@ DecimateResult decimateParallelMesh( MR::Mesh & mesh, const DecimateParallelSett
     seqSettings.maxError = settings.maxError;
     seqSettings.maxEdgeLen = settings.maxEdgeLen;
     seqSettings.maxTriangleAspectRatio = settings.maxTriangleAspectRatio;
+    seqSettings.criticalTriAspectRatio = settings.criticalTriAspectRatio;
     seqSettings.stabilizer = settings.stabilizer;
     seqSettings.optimizeVertexPos = settings.optimizeVertexPos;
     seqSettings.region = settings.region;

source/MRMesh/MRMeshDecimateParallel.h

@@ -24,6 +24,8 @@ struct DecimateParallelSettings
     float maxEdgeLen = 1;
     /// Maximal possible aspect ratio of a triangle introduced during decimation
     float maxTriangleAspectRatio = 20;
+    /// the algorithm will try to eliminate triangles with equal or larger aspect ratio, ignoring normal orientation checks
+    float criticalTriAspectRatio = FLT_MAX;
     /// Small stabilizer is important to achieve good results on completely planar mesh parts,
     /// if your mesh is not-planer everywhere, then you can set it to zero
     float stabilizer = 0.001f;

source/MRMesh/MRMeshDelone.cpp

@@ -12,6 +12,10 @@ inline auto dir( const auto& p, const auto& q, const auto& r )
 {
     return cross( q - p, r - p );
 }
+inline auto area( const auto& p, const auto& q, const auto& r )
+{
+    return dir( p, q, r ).length();
+}
 
 bool checkDeloneQuadrangle( const Vector3d& a, const Vector3d& b, const Vector3d& c, const Vector3d& d, double maxAngleChange )
 {
@@ -47,6 +51,8 @@ bool checkAspectRatiosInQuadrangleT( const Vector3<T>& a, const Vector3<T>& b, c
 {
     auto metricAC = std::max( triangleAspectRatio( a, c, d ), triangleAspectRatio( c, a, b ) );
     auto metricBD = std::max( triangleAspectRatio( b, d, a ), triangleAspectRatio( d, b, c ) );
+    if ( metricAC <= metricBD )
+        return true;
     if ( metricAC < criticalTriAspectRatio && maxAngleChange < NoAngleChangeLimit )
     {
         const auto dirABD = dir( a, b, d );
@@ -61,7 +67,19 @@ bool checkAspectRatiosInQuadrangleT( const Vector3<T>& a, const Vector3<T>& b, c
         if ( angleChange > maxAngleChange )
             return true;
     }
-    return metricAC <= metricBD;
+    else if ( metricAC >= criticalTriAspectRatio )
+    {
+        const auto sABC = area( a, b, c );
+        const auto sACD = area( a, c, d );
+
+        const auto sABD = area( a, b, d );
+        const auto sDBC = area( d, b, c );
+
+        // in case of degenerate triangles, select the subdivision with smaller total area
+        if ( sABC + sACD < sABD + sDBC )
+            return true;
+    }
+    return false;
 }
 
 bool checkAspectRatiosInQuadrangle( const Vector3d& a, const Vector3d& b, const Vector3d& c, const Vector3d& d, double maxAngleChange, double criticalTriAspectRatio )
@@ -127,6 +145,7 @@ bool checkDeloneQuadrangleInMesh( const Mesh & mesh, EdgeId edge, const DeloneSe
         auto diag1 = dp - bp;
         // distance between them
         double dist = fabs( dot( cross( diag0, diag1 ).normalized(), bp - ap ) );
+        // TODO: this actually does not work if one of the diagonals is collapsed in point, where the formula gives 0, which is wrong
         if ( dist > settings.maxDeviationAfterFlip )
             return true; // flipping of given edge will change the surface shape too much
     }

source/MRMesh/MRMeshDelone.h

@@ -39,7 +39,7 @@ MRMESH_API bool checkDeloneQuadrangle( const Vector3f& a, const Vector3f& b, con
 
 /// given quadrangle ABCD, checks whether its subdivision on 
 /// (triangles ABD and DBC) is better than on (triangles ABC and ACD),
-/// where "better" means smaller maximum triangle aspect ratio and the flip of edge from BD to AC will
+/// where "better" means smaller maximum triangle aspect ratio and the flip of edge from AC to BD will
 /// not change dihedral angle more than on given threshold
 MRMESH_API bool checkAspectRatiosInQuadrangle( const Vector3d& a, const Vector3d& b, const Vector3d& c, const Vector3d& d, double maxAngleChange, double criticalTriAspectRatio );
 MRMESH_API bool checkAspectRatiosInQuadrangle( const Vector3f& a, const Vector3f& b, const Vector3f& c, const Vector3f& d, float maxAngleChange, float criticalTriAspectRatio );