diff --git a/source/MRMesh/MRMeshSubdivide.cpp b/source/MRMesh/MRMeshSubdivide.cpp index e1d2cb21a311..3ad5faad24fb 100644 --- a/source/MRMesh/MRMeshSubdivide.cpp +++ b/source/MRMesh/MRMeshSubdivide.cpp @@ -12,6 +12,7 @@ #include "MRRegionBoundary.h" #include "MRBitSetParallelFor.h" #include "MRParallelFor.h" +#include "MRBuffer.h" #include namespace MR @@ -211,7 +212,44 @@ int subdivideMesh( Mesh & mesh, const SubdivideSettings & settings ) return splitsDone; } -TEST(MRMesh, SubdivideMesh) +Expected 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 }, diff --git a/source/MRMesh/MRMeshSubdivide.h b/source/MRMesh/MRMeshSubdivide.h index f318ce4669b1..6853c8b4bb7a 100644 --- a/source/MRMesh/MRMeshSubdivide.h +++ b/source/MRMesh/MRMeshSubdivide.h @@ -2,6 +2,7 @@ #include "MRMeshFwd.h" #include "MRProgressCallback.h" +#include "MRExpected.h" #include "MRConstants.h" #include #include @@ -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 beforeEdgeSplit; + /// this function is called each time a new vertex has been created, but before the ring is made Delone std::function onVertCreated; + /// this function is called each time edge (e) is split into (e1->e), but before the ring is made Delone std::function 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 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 copySubdividePackMesh( const MeshPart & mp, float voxelSize, const ProgressCallback & cb = {} ); + /// \} -} +} //namespace MR diff --git a/source/MRVoxels/MRRebuildMesh.cpp b/source/MRVoxels/MRRebuildMesh.cpp index ef1554f02c78..fc94b7c8d2c2 100644 --- a/source/MRVoxels/MRRebuildMesh.cpp +++ b/source/MRVoxels/MRRebuildMesh.cpp @@ -6,6 +6,7 @@ #include "MRMesh/MRMeshDecimate.h" #include "MRMesh/MRMeshCollide.h" #include "MRMesh/MRTimer.h" +#include "MRMesh/MRMeshSubdivide.h" namespace MR { @@ -13,8 +14,10 @@ namespace MR Expected 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 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 rebuildMesh( const MeshPart& mp, const RebuildMeshSettings& setti if ( settings.onSignDetectionModeSelected ) settings.onSignDetectionModeSelected( genOffsetParams.signDetectionMode ); + std::optional 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 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 ) diff --git a/source/MRVoxels/MRRebuildMesh.h b/source/MRVoxels/MRRebuildMesh.h index cde3fd5a7b24..4e2a9e5986ef 100644 --- a/source/MRVoxels/MRRebuildMesh.h +++ b/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;