Update edge indicator function #1215 (#1229)

* Update edge indicator function

* use k

* meshDenoiseViaNormals

source/MRMesh/MRNormalDenoising.cpp

@@ -2,6 +2,9 @@
 #include "MRMesh.h"
 #include "MRParallelFor.h"
 #include "MRRingIterator.h"
+#include "MRMeshNormals.h"
+#include "MRNormalsToPoints.h"
+#include "MRBitSetParallelFor.h"
 #include "MRTimer.h"
 
 #pragma warning(push)
@@ -97,4 +100,144 @@ void denoiseNormals( const Mesh & mesh, FaceNormals & normals, const Vector<floa
     } );
 }
 
+void updateIndicator( const Mesh & mesh, Vector<float, UndirectedEdgeId> & v, const FaceNormals & normals, float beta, float gamma )
+{
+    MR_TIMER
+
+    const auto sz = v.size();
+    assert( sz == mesh.topology.undirectedEdgeSize() );
+    assert( normals.size() == mesh.topology.faceSize() );
+    if ( sz <= 0 )
+        return;
+
+    std::vector< Eigen::Triplet<double> > mTriplets;
+    Eigen::VectorXd rhs;
+    rhs.resize( sz );
+    constexpr float eps = 0.001f;
+    const float rh = beta / ( 2 * eps );
+    const float k = 2 * beta * eps;
+    for ( auto ue = 0_ue; ue < sz; ++ue )
+    {
+        const EdgeId e = ue;
+        float centralWeight = rh;
+        const auto l = mesh.topology.left( e );
+        const auto r = mesh.topology.right( e );
+        if ( l && r )
+            centralWeight += 2 * gamma * ( normals[l] - normals[r] ).lengthSq();
+        const auto lenE = mesh.edgeLength( e );
+        if ( lenE > 0 )
+        {
+            if ( l )
+            {
+                const auto c = mesh.triCenter( l );
+                {
+                    const auto a = mesh.topology.next( e );
+                    const auto lenL = ( c - mesh.orgPnt( e ) ).length();
+                    const auto x = k * lenL / lenE;
+                    centralWeight += x;
+                    mTriplets.emplace_back( ue, a.undirected(), -x );
+                }
+                {
+                    const auto b = mesh.topology.prev( e.sym() );
+                    const auto lenL = ( c - mesh.destPnt( e ) ).length();
+                    const auto x = k * lenL / lenE;
+                    centralWeight += x;
+                    mTriplets.emplace_back( ue, b.undirected(), -x );
+                }
+            }
+            if ( r )
+            {
+                const auto c = mesh.triCenter( r );
+                {
+                    const auto a = mesh.topology.prev( e );
+                    const auto lenL = ( c - mesh.orgPnt( e ) ).length();
+                    const auto x = k * lenL / lenE;
+                    centralWeight += x;
+                    mTriplets.emplace_back( ue, a.undirected(), -x );
+                }
+                {
+                    const auto b = mesh.topology.next( e.sym() );
+                    const auto lenL = ( c - mesh.destPnt( e ) ).length();
+                    const auto x = k * lenL / lenE;
+                    centralWeight += x;
+                    mTriplets.emplace_back( ue, b.undirected(), -x );
+                }
+            }
+        }
+        mTriplets.emplace_back( ue, ue, centralWeight );
+        rhs[ue] = rh;
+    }
+
+    using SparseMatrix = Eigen::SparseMatrix<double,Eigen::RowMajor>;
+    SparseMatrix A;
+    A.resize( sz, sz );
+    A.setFromTriplets( mTriplets.begin(), mTriplets.end() );
+    Eigen::SimplicialLDLT<SparseMatrix> solver;
+    solver.compute( A );
+
+    Eigen::VectorXd sol = solver.solve( rhs );
+
+    // copy solution back into v
+    ParallelFor( v, [&]( UndirectedEdgeId ue )
+    {
+        v[ue] = (float) sol[ue];
+    } );
+}
+
+VoidOrErrStr meshDenoiseViaNormals( Mesh & mesh, const DenoiseViaNormalsSettings & settings )
+{
+    MR_TIMER
+    if ( settings.normalIters <= 0 || settings.pointIters <= 0 )
+    {
+        assert( false );
+        return tl::make_unexpected( "Bad parameters" );
+    }
+
+    if ( !reportProgress( settings.cb, 0.0f ) )
+        return tl::make_unexpected( "Operation was canceled" );
+
+    auto fnormals0 = computePerFaceNormals( mesh );
+    Vector<float, UndirectedEdgeId> v( mesh.topology.undirectedEdgeSize(), 1 );
+
+    if ( !reportProgress( settings.cb, 0.05f ) )
+        return tl::make_unexpected( "Operation was canceled" );
+
+    auto sp = subprogress( settings.cb, 0.05f, 0.95f );
+    FaceNormals fnormals;
+    for ( int i = 0; i < settings.normalIters; ++i )
+    {
+        fnormals = fnormals0;
+        denoiseNormals( mesh, fnormals, v, settings.gamma );
+        if ( !reportProgress( sp, float( 2 * i ) / ( 2 * settings.normalIters ) ) )
+            return tl::make_unexpected( "Operation was canceled" );
+
+        updateIndicator( mesh, v, fnormals, settings.beta, settings.gamma );
+        if ( !reportProgress( sp, float( 2 * i + 1 ) / ( 2 * settings.normalIters ) ) )
+            return tl::make_unexpected( "Operation was canceled" );
+    }
+
+    if ( settings.outCreases )
+    {
+        settings.outCreases->clear();
+        settings.outCreases->resize( mesh.topology.undirectedEdgeSize() );
+        BitSetParallelForAll( *settings.outCreases, [&]( UndirectedEdgeId ue )
+        {
+            if ( v[ue] < 0.5f )
+                settings.outCreases->set( ue );
+        } );
+    }
+
+    if ( !reportProgress( settings.cb, 0.95f ) )
+        return tl::make_unexpected( "Operation was canceled" );
+
+    const auto guide = mesh.points;
+    NormalsToPoints n2p;
+    n2p.prepare( mesh.topology );
+    for ( int i = 0; i < settings.pointIters; ++i )
+        n2p.run( guide, fnormals, mesh.points );
+
+    reportProgress( settings.cb, 1.0f );
+    return {};
+}
+
 } //namespace MR

source/MRMesh/MRNormalDenoising.h

@@ -1,16 +1,47 @@
 #pragma once
 
 #include "MRMeshFwd.h"
+#include "MRProgressCallback.h"
+#include "MRExpected.h"
 
 namespace MR
 {
 
 /// Smooth face normals, given
-/// \param mesh contains topology information and edge lengths for weights
+/// \param mesh contains topology information and coordinates for equation weights
 /// \param normals input noisy normals and output smooth normals
-/// \param v edge indicator: 1 - smooth edge, 0 - crease edge
-/// \param gamma the more the value, the larger smoothing is (approximately in the range [1,30])
-/// see the article "Mesh Denoising via a Novel Mumford-Shah Framework"
+/// \param v edge indicator function (1 - smooth edge, 0 - crease edge)
+/// \param gamma the amount of smoothing: 0 - no smoothing, 1 - average smoothing, ...
+/// see the article "Mesh Denoising via a Novel Mumford-Shah Framework", equation (19)
 MRMESH_API void denoiseNormals( const Mesh & mesh, FaceNormals & normals, const Vector<float, UndirectedEdgeId> & v, float gamma );
 
+/// Update edge indicator function (1 - smooth edge, 0 - crease edge), given
+/// \param mesh contains topology information and coordinates for equation weights
+/// \param v noisy input and smooth output
+/// \param normals per-face normals
+/// \param beta 0.001 - sharp edges, 0.01 - moderate edges, 0.1 - smooth edges
+/// \param gamma the amount of smoothing: 0 - no smoothing, 1 - average smoothing, ...
+/// see the article "Mesh Denoising via a Novel Mumford-Shah Framework", equation (20)
+MRMESH_API void updateIndicator( const Mesh & mesh, Vector<float, UndirectedEdgeId> & v, const FaceNormals & normals, float beta, float gamma );
+
+struct DenoiseViaNormalsSettings
+{
+    /// 0.001 - sharp edges, 0.01 - moderate edges, 0.1 - smooth edges
+    float beta = 0.001f;
+    /// the amount of smoothing: 0 - no smoothing, 1 - average smoothing, ...
+    float gamma = 5.f;
+    /// the number of iterations to smooth normals and find creases; the more the better quality, but longer computation
+    int normalIters = 10;
+    /// the number of iterations to update vertex coordinates from found normals; the more the better quality, but longer computation
+    int pointIters = 20;
+    /// optionally returns creases found during smoothing
+    UndirectedEdgeBitSet * outCreases = nullptr;
+    /// to get the progress and optionally cancel
+    ProgressCallback cb = {};
+};
+
+/// Reduces noise in given mesh,
+/// see the article "Mesh Denoising via a Novel Mumford-Shah Framework"
+MRMESH_API VoidOrErrStr meshDenoiseViaNormals( Mesh & mesh, const DenoiseViaNormalsSettings & settings = {} );
+
 } //namespace MR