Decimate: new option angleWeightedDistToPlane (#3929)

source/MRMesh/MRMesh.cpp

@@ -751,7 +751,7 @@ float Mesh::leftCotan( EdgeId e ) const
     return nom / den;
 }
 
-QuadraticForm3f Mesh::quadraticForm( VertId v, const FaceBitSet * region, const UndirectedEdgeBitSet * creases ) const
+QuadraticForm3f Mesh::quadraticForm( VertId v, bool angleWeigted, const FaceBitSet * region, const UndirectedEdgeBitSet * creases ) const
 {
     QuadraticForm3f qf;
     for ( EdgeId e : orgRing( topology, v ) )
@@ -764,9 +764,20 @@ QuadraticForm3f Mesh::quadraticForm( VertId v, const FaceBitSet * region, const
         }
         if ( topology.left( e ) ) // intentionally do not check that left face is in region to respect its plane as well
         {
-            // zero-area triangle is treated as no triangle with no penalty at all,
-            // otherwise it penalizes the shift proportionally to the distance from the plane containing the triangle
-            qf.addDistToPlane( leftNormal( e ) );
+            if ( angleWeigted )
+            {
+                auto d0 = edgeVector( e );
+                auto d1 = edgeVector( topology.next( e ) );
+                auto angle = MR::angle( d0, d1 );
+                static constexpr float INV_PIF = 1 / PI_F;
+                qf.addDistToPlane( leftNormal( e ), angle * INV_PIF );
+            }
+            else
+            {
+                // zero-area triangle is treated as no triangle with no penalty at all,
+                // otherwise it penalizes the shift proportionally to the distance from the plane containing the triangle
+                qf.addDistToPlane( leftNormal( e ) );
+            }
         }
     }
     return qf;

source/MRMesh/MRMesh.h

@@ -301,9 +301,10 @@ struct [[nodiscard]] Mesh
     [[nodiscard]] float cotan( UndirectedEdgeId ue ) const { EdgeId e{ ue }; return leftCotan( e ) + leftCotan( e.sym() ); }
 
     /// computes quadratic form in the vertex as the sum of squared distances from
-    /// 1) planes of adjacent triangles
+    /// 1) planes of adjacent triangles, with the weight equal to the angle of adjacent triangle at this vertex divided on PI in case of angleWeigted=true;
     /// 2) lines of adjacent boundary and crease edges
-    [[nodiscard]] MRMESH_API QuadraticForm3f quadraticForm( VertId v, const FaceBitSet * region = nullptr, const UndirectedEdgeBitSet * creases = nullptr ) const;
+    [[nodiscard]] MRMESH_API QuadraticForm3f quadraticForm( VertId v, bool angleWeigted,
+        const FaceBitSet * region = nullptr, const UndirectedEdgeBitSet * creases = nullptr ) const;
 
     /// passes through all valid vertices and finds the minimal bounding box containing all of them;
     /// if toWorld transformation is given then returns minimal bounding box in world space

source/MRMesh/MRMeshDecimate.cpp

@@ -203,14 +203,14 @@ class MeshDecimator::EdgeMetricCalc
     std::vector<QueueElement> elems_;
 };
 
-QuadraticForm3f computeFormAtVertex( const MR::MeshPart & mp, MR::VertId v, float stabilizer, const UndirectedEdgeBitSet * creases )
+QuadraticForm3f computeFormAtVertex( const MR::MeshPart & mp, MR::VertId v, float stabilizer, bool angleWeigted, const UndirectedEdgeBitSet * creases )
 {
-    QuadraticForm3f qf = mp.mesh.quadraticForm( v, mp.region, creases );
+    QuadraticForm3f qf = mp.mesh.quadraticForm( v, angleWeigted, mp.region, creases );
     qf.addDistToOrigin( stabilizer );
     return qf;
 }
 
-Vector<QuadraticForm3f, VertId> computeFormsAtVertices( const MeshPart & mp, float stabilizer, const UndirectedEdgeBitSet * creases )
+Vector<QuadraticForm3f, VertId> computeFormsAtVertices( const MeshPart & mp, float stabilizer, bool angleWeigted, const UndirectedEdgeBitSet * creases )
 {
     MR_TIMER;
 
@@ -220,7 +220,7 @@ Vector<QuadraticForm3f, VertId> computeFormsAtVertices( const MeshPart & mp, flo
     Vector<QuadraticForm3f, VertId> res( regionVertices.find_last() + 1 );
     BitSetParallelFor( regionVertices, [&]( VertId v )
     {
-        res[v] = computeFormAtVertex( mp, v, stabilizer, creases );
+        res[v] = computeFormAtVertex( mp, v, stabilizer, angleWeigted, creases );
     } );
 
     return res;
@@ -264,7 +264,7 @@ bool MeshDecimator::initializeQueue_()
         pVertForms_ = &myVertForms_;
 
     if ( pVertForms_->empty() )
-        *pVertForms_ = computeFormsAtVertices( MeshPart{ mesh_, settings_.region }, settings_.stabilizer, settings_.notFlippable );
+        *pVertForms_ = computeFormsAtVertices( MeshPart{ mesh_, settings_.region }, settings_.stabilizer, settings_.angleWeightedDistToPlane, settings_.notFlippable );
 
     if ( settings_.progressCallback && !settings_.progressCallback( 0.1f ) )
         return false;
@@ -940,7 +940,7 @@ static DecimateResult decimateMeshParallelInplace( MR::Mesh & mesh, const Decima
     if ( settings.vertForms )
         mVertForms = std::move( *settings.vertForms );
     if ( mVertForms.empty() )
-        mVertForms = computeFormsAtVertices( MeshPart{ mesh, settings.region }, settings.stabilizer, settings.notFlippable );
+        mVertForms = computeFormsAtVertices( MeshPart{ mesh, settings.region }, settings.stabilizer, settings.angleWeightedDistToPlane, settings.notFlippable );
     if ( settings.progressCallback && !settings.progressCallback( 0.2f ) )
         return res;
 

source/MRMesh/MRMeshDecimate.h

@@ -58,6 +58,10 @@ struct DecimateSettings
     /// if your mesh is not-planer everywhere, then you can set it to zero
     float stabilizer = 0.001f;
 
+    /// if false, then quadratic error metric is equal to the sum of distances to the planes of original mesh triangles;
+    /// if true, then the sum is weighted, and the weight is equal to the angle of adjacent triangle at the vertex divided on PI (to get one after summing all 3 vertices of the triangle)
+    bool angleWeightedDistToPlane = false;
+
     /// if true then after each edge collapse the position of remaining vertex is optimized to
     /// minimize local shape change, if false then the edge is collapsed in one of its vertices, which keeps its position
     bool optimizeVertexPos = true;
@@ -193,13 +197,13 @@ MRMESH_API DecimateResult decimateMesh( Mesh & mesh, const DecimateSettings & se
  * \brief Computes quadratic form at given vertex of the initial surface before decimation
  * \ingroup DecimateGroup
  */
-[[nodiscard]] MRMESH_API QuadraticForm3f computeFormAtVertex( const MeshPart & mp, VertId v, float stabilizer, const UndirectedEdgeBitSet * creases = nullptr );
+[[nodiscard]] MRMESH_API QuadraticForm3f computeFormAtVertex( const MeshPart & mp, VertId v, float stabilizer, bool angleWeigted, const UndirectedEdgeBitSet * creases = nullptr );
 
 /**
  * \brief Computes quadratic forms at every vertex of mesh part before decimation
  * \ingroup DecimateGroup
  */
-[[nodiscard]] MRMESH_API Vector<QuadraticForm3f, VertId> computeFormsAtVertices( const MeshPart & mp, float stabilizer, const UndirectedEdgeBitSet * creases = nullptr );
+[[nodiscard]] MRMESH_API Vector<QuadraticForm3f, VertId> computeFormsAtVertices( const MeshPart & mp, float stabilizer, bool angleWeigted, const UndirectedEdgeBitSet * creases = nullptr );
 
 /**
  * \brief returns given subdivision part of all valid faces;