rebuildMesh: preSubdivide option (#3936)

source/MRMesh/MRMeshSubdivide.cpp

@@ -12,6 +12,7 @@
 #include "MRRegionBoundary.h"
 #include "MRBitSetParallelFor.h"
 #include "MRParallelFor.h"
+#include "MRBuffer.h"
 #include <queue>
 
 namespace MR
@@ -211,7 +212,44 @@ int subdivideMesh( Mesh & mesh, const SubdivideSettings & settings )
     return splitsDone;
 }
 
-TEST(MRMesh, SubdivideMesh) 
+Expected<Mesh> copySubdividePackMesh( const MeshPart & mp, float voxelSize, const ProgressCallback & cb )
+{
+    MR_TIMER
+    Mesh subMesh;
+
+    // copy mesh
+    if ( mp.region )
+        subMesh.addMeshPart( mp );
+    else
+        subMesh = mp.mesh;
+    if ( !reportProgress( cb, 0.25f ) )
+        return unexpectedOperationCanceled();
+
+    // smaller edge lengths require too much space and time for subdivision,
+    // larger edge lengths do not give significant benefits in AABB-tree optimized searches
+    const float approxMaxEdgeLen = 5 * voxelSize;
+
+    // subdivide copied mesh
+    SubdivideSettings subSettings
+    {
+        .maxEdgeLen = approxMaxEdgeLen,
+        .maxEdgeSplits = int( subMesh.area() / sqr( approxMaxEdgeLen ) ),
+        .maxDeviationAfterFlip = 0.1f * voxelSize,
+        .progressCallback = subprogress( cb, 0.25f, 0.75f )
+    };
+    subdivideMesh( subMesh, subSettings );
+    if ( !reportProgress( cb, 0.75f ) )
+        return unexpectedOperationCanceled();
+
+    // pack subdivided mesh
+    subMesh.packOptimally();
+    if ( !reportProgress( cb, 1.0f ) )
+        return unexpectedOperationCanceled();
+
+    return subMesh;
+}
+
+TEST(MRMesh, SubdivideMesh)
 {
     Triangulation t{
         { 0_v, 1_v, 2_v },

source/MRMesh/MRMeshSubdivide.h

@@ -2,6 +2,7 @@
 
 #include "MRMeshFwd.h"
 #include "MRProgressCallback.h"
+#include "MRExpected.h"
 #include "MRConstants.h"
 #include <cfloat>
 #include <functional>
@@ -17,53 +18,74 @@ struct SubdivideSettings
 {
     /// Subdivision is stopped when all edges inside or on the boundary of the region are not longer than this value
     float maxEdgeLen = 0;
+
     /// Maximum number of edge splits allowed
     int maxEdgeSplits = 1000;
+
     /// Improves local mesh triangulation by doing edge flips if it does not make too big surface deviation
     float maxDeviationAfterFlip = 1;
+
     /// Improves local mesh triangulation by doing edge flips if it does not change dihedral angle more than on this value (in radians)
     float maxAngleChangeAfterFlip = FLT_MAX;
+
     /// If this value is less than FLT_MAX then edge flips will
     /// ignore dihedral angle check if one of triangles has aspect ratio more than this value
     /// Unit: rad
     float criticalAspectRatioFlip = 1000.0f;
+
     /// Region on mesh to be subdivided, it is updated during the operation
     FaceBitSet * region = nullptr;
+
     /// Edges specified by this bit-set will never be flipped, but they can be split so it is updated during the operation
     UndirectedEdgeBitSet* notFlippable = nullptr;
+
     /// New vertices appeared during subdivision will be added here
     VertBitSet * newVerts = nullptr;
+
     /// If false do not touch border edges (cannot subdivide lone faces)\n
     /// use \ref MR::findRegionOuterFaces to find boundary faces
     bool subdivideBorder = true;
+
     /// The subdivision stops as soon as all triangles (in the region) have aspect ratio below or equal to this value
     float maxTriAspectRatio = 0;
+
     /// An edge is subdivided only if both its left and right triangles have aspect ratio below or equal to this value.
     /// So this is a maximum aspect ratio of a triangle that can be split on two before Delone optimization.
     /// Please set it to a smaller value only if subdivideBorder==false, otherwise many narrow triangles can appear near border
     float maxSplittableTriAspectRatio = FLT_MAX;
+
     /// Puts new vertices so that they form a smooth surface together with existing vertices.
     /// This option works best for natural surfaces without sharp edges in between triangles
     bool smoothMode = false;
+
     /// In case of activated smoothMode, the smoothness is locally deactivated at the edges having
     /// dihedral angle at least this value
     float minSharpDihedralAngle = PI_F / 6; // 30 degrees
+
     /// if true, then every new vertex will be projected on the original mesh (before smoothing)
     bool projectOnOriginalMesh = false;
+
+    /// this function is called each time edge (e) is going to split, if it returns false then this split will be skipped
+    std::function<bool(EdgeId e)> beforeEdgeSplit;
+
     /// this function is called each time a new vertex has been created, but before the ring is made Delone
     std::function<void(VertId)> onVertCreated;
+
     /// this function is called each time edge (e) is split into (e1->e), but before the ring is made Delone
     std::function<void(EdgeId e1, EdgeId e)> onEdgeSplit;
-    /// this function is called each time edge (e) is going to split, if it returns false then this split will be skipped
-    std::function<bool(EdgeId e)> beforeEdgeSplit;
+
     /// callback to report algorithm progress and cancel it by user request
-    ProgressCallback progressCallback = {};
+    ProgressCallback progressCallback;
 };
 
-/// Split edges in mesh region according to the settings;\n
+/// splits edges in mesh region according to the settings;\n
 /// \return The total number of edge splits performed
 MRMESH_API int subdivideMesh( Mesh & mesh, const SubdivideSettings & settings = {} );
 
+/// creates a copy of given mesh part, subdivides it to get rid of too long edges compared with voxelSize, then packs resulting mesh,
+/// this is called typically in preparation for 3D space sampling with voxelSize step, and subdivision is important for making leaves of AABB tree not too big compared with voxelSize
+[[nodiscard]] MRMESH_API Expected<Mesh> copySubdividePackMesh( const MeshPart & mp, float voxelSize, const ProgressCallback & cb = {} );
+
 /// \}
 
-}
+} //namespace MR

source/MRVoxels/MRRebuildMesh.cpp

@@ -6,15 +6,18 @@
 #include "MRMesh/MRMeshDecimate.h"
 #include "MRMesh/MRMeshCollide.h"
 #include "MRMesh/MRTimer.h"
+#include "MRMesh/MRMeshSubdivide.h"
 
 namespace MR
 {
 
 Expected<Mesh> rebuildMesh( const MeshPart& mp, const RebuildMeshSettings& settings )
 {
     MR_TIMER
-    GeneralOffsetParameters genOffsetParams;
 
+    auto progress = settings.progress;
+
+    GeneralOffsetParameters genOffsetParams;
     switch ( settings.signMode )
     {
     default:
@@ -23,7 +26,8 @@ Expected<Mesh> rebuildMesh( const MeshPart& mp, const RebuildMeshSettings& setti
     case SignDetectionModeShort::Auto:
         if ( mp.mesh.topology.isClosed( mp.region ) )
         {
-            auto expSelfy = findSelfCollidingTriangles( mp, nullptr, subprogress( settings.progress, 0.0f, 0.1f ) );
+            auto expSelfy = findSelfCollidingTriangles( mp, nullptr, subprogress( progress, 0.0f, 0.1f ) );
+            progress = subprogress( progress, 0.1f, 1.0f );
             if ( !expSelfy )
                 return unexpected( std::move( expSelfy.error() ) );
             if ( *expSelfy )
@@ -48,25 +52,35 @@ Expected<Mesh> rebuildMesh( const MeshPart& mp, const RebuildMeshSettings& setti
     if ( settings.onSignDetectionModeSelected )
         settings.onSignDetectionModeSelected( genOffsetParams.signDetectionMode );
 
+    std::optional<Mesh> subMesh;
+    if ( settings.preSubdivide && genOffsetParams.signDetectionMode != SignDetectionMode::OpenVDB ) // OpenVDB slows down with more input triangles
+    {
+        if ( auto maybeMesh = copySubdividePackMesh( mp, settings.voxelSize, subprogress( progress, 0.0f, 0.1f ) ) )
+            subMesh = std::move( *maybeMesh );
+        else
+            return unexpected( std::move( maybeMesh.error() ) );
+        progress = subprogress( progress, 0.1f, 1.0f );
+    }
+
     genOffsetParams.closeHolesInHoleWindingNumber = settings.closeHolesInHoleWindingNumber;
     genOffsetParams.voxelSize = settings.voxelSize;
     genOffsetParams.mode = settings.offsetMode;
     genOffsetParams.windingNumberThreshold = settings.windingNumberThreshold;
     genOffsetParams.windingNumberBeta = settings.windingNumberBeta;
     genOffsetParams.fwn = settings.fwn;
-    genOffsetParams.callBack = subprogress( settings.progress, 0.1f, ( settings.decimate ? 0.7f : 1.0f ) );
+    genOffsetParams.callBack = subprogress( progress, 0.0f, ( settings.decimate ? 0.7f : 1.0f ) );
 
     UndirectedEdgeBitSet sharpEdges;
     genOffsetParams.outSharpEdges = &sharpEdges;
 
-    auto resMesh = generalOffsetMesh( mp, 0.0f, genOffsetParams );
+    auto resMesh = generalOffsetMesh( subMesh ? *subMesh : mp, 0.0f, genOffsetParams );
     if ( !resMesh.has_value() )
         return resMesh;
 
     if ( settings.decimate && resMesh->topology.numValidFaces() > 0 )
     {
         const auto map = resMesh->packOptimally( false );
-        if ( !reportProgress( settings.progress, 0.75f ) )
+        if ( !reportProgress( progress, 0.75f ) )
             return unexpectedOperationCanceled();
 
         sharpEdges = mapEdges( map.e, sharpEdges );
@@ -78,7 +92,7 @@ Expected<Mesh> rebuildMesh( const MeshPart& mp, const RebuildMeshSettings& setti
             .stabilizer = 1e-5f, // 1e-6 here resulted in a bit worse mesh
             .notFlippable = sharpEdges.any() ? &sharpEdges : nullptr,
             .packMesh = true,
-            .progressCallback = subprogress( settings.progress, 0.75f, 1.0f ),
+            .progressCallback = subprogress( progress, 0.75f, 1.0f ),
             .subdivideParts = 64
         };
         if ( decimateMesh( *resMesh, decimSettings ).cancelled )

source/MRVoxels/MRRebuildMesh.h

@@ -11,6 +11,13 @@ namespace MR
 
 struct RebuildMeshSettings
 {
+    /// whether to make subdivision of initial mesh before conversion to voxels,
+    /// despite time and memory required for the subdivision, it typically makes the whole rebuilding faster (or even much faster in case of large initial triangles),
+    /// because AABB tree contains small triangles in leaves, which is good for both
+    /// 1) search for closest triangle because the closest box more frequently contains the closest triangle,
+    /// 2) and winding number approximation because of more frequent usage of approximation for distant dipoles
+    bool preSubdivide = true;
+
     /// Size of voxel in grid conversions;
     /// The user is responsible for setting some positive value here
     float voxelSize = 0;