DM-based Radius Compensation fixes (#4084)

* DM-based Radius Compensation fixes

* update

* remove debug output

* use fixed tolerance

Author: Grantim

source/MRMesh/MRRadiusCompensation.cpp

@@ -32,7 +32,7 @@ class RadiusCompensator
     // prepares raw distance map
     Expected<void> init();
 
-    // finds tool locations for each pixel and summary compensation volume for each
+    // finds tool locations for each pixel and summary compensation cost for each
     Expected<void> calcCompensations();
 
     // creates compensation distance map
@@ -61,8 +61,8 @@ class RadiusCompensator
     // if callback returns false iterations stops
     void iteratePixelsInRadius_( const Vector2i& pixelCoord, const std::function<bool( const Vector2i& )>& callback );
 
-    // calculates summary compensation volume for given tool location
-    float sumCompensationVolume_( const Vector3f& toolCenter );
+    // calculates summary compensation cost for given tool location
+    float sumCompensationCost_( const Vector3f& toolCenter );
 
     Mesh& mesh_;
     CompensateRadiusParams params_;
@@ -77,9 +77,13 @@ class RadiusCompensator
     Vector2i pixelsInRadius_;
 
     std::vector<Vector3f> toolCenters_; // per pixel
-    std::vector<std::pair<float,int>> volumes_; // <voulme, id> per pixel (this array will be sorted, thats why we need id as value and not only as key)
+    // cost is (approx compensation Volume)/(approx compensation Projected Area) - less is better
+    std::vector<std::pair<float,int>> costs_; // <cost, id> per pixel (this array will be sorted, thats why we need id as value and not only as key)
 
     DistanceMap compensatedDm_;
+
+    // tolerance of distance map height for "touching" compensations
+    static constexpr float cDMTolearance = 1e-5f;
 };
 
 Expected<void> RadiusCompensator::init()
@@ -113,20 +117,20 @@ Expected<void> RadiusCompensator::calcCompensations()
 {
     MR_TIMER;
     toolCenters_.resize( dm_.size(), Vector3f::diagonal( FLT_MAX ) );
-    volumes_.resize( dm_.size(), std::make_pair( -1.0f, -1 ) );
+    costs_.resize( dm_.size(), std::make_pair( -1.0f, -1 ) );
     bool keepGoing = ParallelFor( size_t( 0 ), dm_.size(), [&] ( size_t i )
     {
         auto pixelCoord = dm_.toPos( i );
         auto& toolCenter = toolCenters_[i];
         toolCenter = findToolCenterAtPixel_(pixelCoord);
         if ( toolCenter.x != FLT_MAX )
-            volumes_[i] = std::make_pair( sumCompensationVolume_( toolCenter ), int( i ) );
+            costs_[i] = std::make_pair( sumCompensationCost_( toolCenter ), int( i ) );
     }, subprogress( params_.callback, 0.05f, 0.15f ) );
 
     if ( !keepGoing )
         return unexpectedOperationCanceled();
 
-    tbb::parallel_sort( volumes_.begin(), volumes_.end(), [] ( const auto& l, const auto& r )
+    tbb::parallel_sort( costs_.begin(), costs_.end(), [] ( const auto& l, const auto& r )
     {
         return l.first < r.first;
     } );
@@ -143,12 +147,13 @@ Expected<void> RadiusCompensator::compensateDistanceMap()
     auto sb = subprogress( params_.callback, 0.2f, 0.5f );
     compensatedDm_ = DistanceMap( dm_.resX(), dm_.resY() );
     size_t i = 0;
-    for ( auto [volume, cId] : volumes_ )
+    const float cTolerance = params_.toolRadius * cDMTolearance;
+    for ( auto [cost, cId] : costs_ )
     {
-        if ( ( ++i % 1024 == 0 ) && !reportProgress( sb, float( i ) / volumes_.size() ) )
+        if ( ( ++i % 1024 == 0 ) && !reportProgress( sb, float( i ) / costs_.size() ) )
             return unexpectedOperationCanceled();
 
-        if ( cId < 0 || volume < 0.0f )
+        if ( cId < 0 || cost < 0.0f )
             continue;
 
         const auto& toolCenter = toolCenters_[cId];
@@ -169,10 +174,10 @@ Expected<void> RadiusCompensator::compensateDistanceMap()
         iteratePixelsInRadius_( pixelCoord, [&] ( const Vector2i& pixeli )->bool
         {
             auto realValue = dm_.getValue( pixeli.x, pixeli.y );
-            auto& value = compensatedDm_.getValue( pixeli.x, pixeli.y );
-            auto compValue = calcCompensatedHeightAtPixel_( pixeli, toolCenter );
-            if ( compValue > value && compValue >= realValue )
-                value = compValue;
+            auto& compValue = compensatedDm_.getValue( pixeli.x, pixeli.y );
+            auto newCompValue = calcCompensatedHeightAtPixel_( pixeli, toolCenter );
+            if ( newCompValue > compValue && newCompValue + cTolerance >= realValue )
+                compValue = std::max( newCompValue, realValue );
             return true;
         } );
     }
@@ -208,6 +213,7 @@ Expected<void> RadiusCompensator::applyCompensation()
     }
 
     // fix inverted faces (undercuts on original mesh)
+    for ( int iFlipped = 0; iFlipped < 10; ++iFlipped ) // repeat until no flipped faces left, 10 - max iters
     {
         // only fix flipped areas
         auto sumDirArea = Vector3f( mesh_.dirArea( faceRegion_ ).normalized() ); // we only care about direction
@@ -221,6 +227,9 @@ Expected<void> RadiusCompensator::applyCompensation()
         expand( mesh_.topology, flippedFaces, 4 );
         flippedFaces &= *faceRegion_;
 
+        if ( flippedFaces.none() )
+            break;
+
         auto equalizingVertRegion = getInnerVerts( mesh_.topology, flippedFaces );
         assert( ( equalizingVertRegion & bounds ).none() );
         MeshEqualizeTriAreasParams etParams;
@@ -254,7 +263,6 @@ Expected<void> RadiusCompensator::applyCompensation()
         vertRegion_.reset( v );
     }, subprogress( params_.callback, 0.65f, 0.8f ) );
 
-
     if ( vertRegion_.any() )
         positionVertsSmoothlySharpBd( mesh_, vertRegion_ );
 
@@ -270,7 +278,7 @@ Expected<void> RadiusCompensator::postprocessMesh()
     auto edgeBounds = findRegionBoundaryUndirectedEdgesInsideMesh( mesh_.topology, *faceRegion_ );
     RemeshSettings rParams;
     rParams.finalRelaxIters = 2;
-    rParams.targetEdgeLen = mesh_.averageEdgeLength();
+    rParams.targetEdgeLen = params_.remeshTargetEdgeLength <= 0.0f ? mesh_.averageEdgeLength() : params_.remeshTargetEdgeLength;
     rParams.region = params_.region;
     rParams.notFlippable = &edgeBounds;
     rParams.progressCallback = subprogress( params_.callback, 0.8f, 1.0f );
@@ -387,19 +395,24 @@ void RadiusCompensator::iteratePixelsInRadius_( const Vector2i& pixelCoord, cons
     }
 }
 
-float RadiusCompensator::sumCompensationVolume_( const Vector3f& toolCenter )
+float RadiusCompensator::sumCompensationCost_( const Vector3f& toolCenter )
 {
-    float sumVolume = 0.0f;
+    double sumVolume = 0.0;
+    double sumArea = 0.0;
     auto toolPixel = Vector2i( to2dim( toDmXf_( toolCenter ) ) );
+    const float cTolerance = cDMTolearance * params_.toolRadius;
     iteratePixelsInRadius_( toolPixel, [&] ( const Vector2i& pixelCoord )->bool
     {
         auto value = dm_.getValue( pixelCoord.x, pixelCoord.y );
         auto height = calcCompensatedHeightAtPixel_( pixelCoord, toolCenter );
-        if ( height > value )
-            sumVolume += ( height - value );
+        if ( height + cTolerance >= value )
+        {
+            sumVolume += std::max( height - value, 0.0f );
+            sumArea += 1.0;
+        }
         return true;
     } );
-    return sumVolume;
+    return sumArea == 0.0 ? -1.0f : float( sumVolume / sumArea );
 }
 
 Expected<void> compensateRadius( Mesh& mesh, const CompensateRadiusParams& params )

source/MRMesh/MRRadiusCompensation.h

@@ -17,14 +17,18 @@ struct CompensateRadiusParams
     float toolRadius{ 0.0f };
 
     /// resolution of distance map that is used for compensation
-    Vector2i distanceMapResolution = Vector2i( 100, 100 );
+    Vector2i distanceMapResolution = Vector2i( 150, 150 );
 
     /// region of the mesh that will be compensated
     /// it should not contain closed components
     /// it is updated during algorithm
     /// also please note that boundaries of the region are fixed
     FaceBitSet* region{ nullptr };
 
+    /// this value will be used for post-process re-meshing
+    /// value less or equal to zero will use average mesh edge length
+    float remeshTargetEdgeLength{ -1.0f };
+
     ProgressCallback callback;
 };