RebuildMesh: restore original performance (#3920)

source/MRCuda/MRCudaFastWindingNumber.cpp

@@ -145,7 +145,7 @@ Expected<void> FastWindingNumber::calcFromGrid( std::vector<float>& res, const V
     return {};
 }
 
-Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb )
+Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb )
 {
     MR_TIMER
     prepareData_( {} );
@@ -171,7 +171,7 @@ Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>&
         int3{ dims.x, dims.y, dims.z },
         cudaGridToMeshXf,
         data_->dipoles.data(), data_->cudaNodes.data(), data_->cudaMeshPoints.data(), data_->cudaFaces.data(),
-        cudaResult.data(), windingNumberThreshold, beta, maxDistSq, minDistSq );
+        cudaResult.data(), options );
 
     if ( auto code = CUDA_EXEC( cudaGetLastError() ) )
         return unexpected( Cuda::getError( code ) );

source/MRCuda/MRCudaFastWindingNumber.cu

@@ -1,6 +1,7 @@
 #include "MRCudaFastWindingNumber.cuh"
 #include "MRMesh/MRAABBTree.h"
 #include "MRMesh/MRConstants.h"
+#include "MRMesh/MRDistanceToMeshOptions.h"
 #include "device_launch_parameters.h"
 #include <limits>
 
@@ -200,7 +201,7 @@ static constexpr float cQuietNan = std::numeric_limits<float>::quiet_NaN();
 
 __global__ void signedDistanceKernel( int3 dims, Matrix4 gridToMeshXf,
     const Dipole* __restrict__ dipoles, const Node3* __restrict__ nodes, const float3* __restrict__ meshPoints, const FaceToThreeVerts* __restrict__ faces,
-    float* resVec, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, size_t size )
+    float* resVec, DistanceToMeshOptions options, size_t size ) // pass options by value to avoid reference on CPU memory
 {
     if ( size == 0 )
     {
@@ -218,17 +219,17 @@ __global__ void signedDistanceKernel( int3 dims, Matrix4 gridToMeshXf,
     const float3 point{ float( voxel.x ), float( voxel.y ), float( voxel.z ) };
     const float3 transformedPoint = gridToMeshXf.isIdentity ? point : gridToMeshXf.transform( point );
 
-    float resSq = calcDistanceSq( transformedPoint, nodes, meshPoints, faces, maxDistSq, minDistSq );
-    if ( resSq < minDistSq || resSq >= maxDistSq ) // note that resSq == minDistSq (e.g. == 0) is a valid situation
+    float resSq = calcDistanceSq( transformedPoint, nodes, meshPoints, faces, options.maxDistSq, options.minDistSq );
+    if ( options.nullOutsideMinMax && ( resSq < options.minDistSq || resSq >= options.maxDistSq ) ) // note that resSq == minDistSq (e.g. == 0) is a valid situation
     {
         resVec[index] = cQuietNan;
         return;
     }
 
     float fwn{ 0 };
-    processPoint( transformedPoint, fwn, dipoles, nodes, meshPoints, faces, beta, index );
+    processPoint( transformedPoint, fwn, dipoles, nodes, meshPoints, faces, options.windingNumberBeta, index );
     float res = sqrt( resSq );
-    if ( fwn > windingNumberThreshold )
+    if ( fwn > options.windingNumberThreshold )
         res = -res;
     resVec[index] = res;
 }
@@ -260,11 +261,11 @@ void fastWindingNumberFromGrid( int3 dims, Matrix4 gridToMeshXf,
 
 void signedDistance( int3 dims, Matrix4 gridToMeshXf,
                                         const Dipole* dipoles, const Node3* nodes, const float3* meshPoints, const FaceToThreeVerts* faces,
-                                        float* resVec, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq )
+                                        float* resVec, const DistanceToMeshOptions& options )
 {
     const size_t size = size_t( dims.x ) * dims.y * dims.z;
     int numBlocks = ( int( size ) + maxThreadsPerBlock - 1 ) / maxThreadsPerBlock;
-    signedDistanceKernel<<< numBlocks, maxThreadsPerBlock >>>( dims, gridToMeshXf, dipoles, nodes, meshPoints, faces, resVec, windingNumberThreshold, beta, maxDistSq, minDistSq, size );
+    signedDistanceKernel<<< numBlocks, maxThreadsPerBlock >>>( dims, gridToMeshXf, dipoles, nodes, meshPoints, faces, resVec, options, size );
 }
 
 } //namespece Cuda

source/MRCuda/MRCudaFastWindingNumber.cuh

@@ -1,11 +1,14 @@
 #pragma once
+
 #include "MRCudaBasic.cuh"
 #include "MRCudaMath.cuh"
 #include "MRCudaFloat.cuh"
 
 namespace MR
 {
 
+struct DistanceToMeshOptions;
+
 namespace Cuda
 {
 // GPU analog of CPU Dipole struct
@@ -48,7 +51,7 @@ void fastWindingNumberFromGrid( int3 gridSize, Matrix4 gridToMeshXf,
 /// calls fast winding number for each point in three-dimensional grid to get sign
 void signedDistance( int3 gridSize, Matrix4 gridToMeshXf,
                                       const Dipole* dipoles, const Node3* nodes, const float3* meshPoints, const FaceToThreeVerts* faces,
-                                      float* resVec, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq );
+                                      float* resVec, const DistanceToMeshOptions& options );
 
 
 } //namespece Cuda

source/MRCuda/MRCudaFastWindingNumber.h

@@ -2,6 +2,7 @@
 #include "exports.h"
 #include "MRMesh/MRFastWindingNumber.h"
 #include "MRMesh/MRMesh.h"
+#include "MRMesh/MRDistanceToMeshOptions.h" // only for bindings generation
 
 namespace MR
 {
@@ -26,7 +27,7 @@ class MRCUDA_CLASS FastWindingNumber : public IFastWindingNumber
     MRCUDA_API void calcFromVector( std::vector<float>& res, const std::vector<Vector3f>& points, float beta, FaceId skipFace = {} ) override;
     MRCUDA_API bool calcSelfIntersections( FaceBitSet& res, float beta, ProgressCallback cb ) override;
     MRCUDA_API Expected<void> calcFromGrid( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float beta, ProgressCallback cb ) override;
-    MRCUDA_API Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) override;
+    MRCUDA_API Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) override;
 
 private:
     bool prepareData_( ProgressCallback cb );

source/MRMesh/MRDistanceToMeshOptions.h

@@ -0,0 +1,42 @@
+#pragma once
+
+#include "MRSignDetectionMode.h"
+#include <cfloat>
+
+namespace MR
+{
+
+/// options determining computation of distance from a point to a mesh
+struct DistanceToMeshOptions
+{
+    /// minimum squared distance from a point to mesh to be computed precisely
+    float minDistSq{ 0 };
+
+    /// maximum squared distance from a point to mesh to be computed precisely
+    float maxDistSq{ FLT_MAX };
+
+    /// what to do if actual distance is outside [min, max) range:
+    /// true - return std::nullopt for std::optional<float> or NaN for float,
+    /// false - return approximate value of the distance (with correct sign in case of SignDetectionMode::HoleWindingRule);
+    /// please note that in HoleWindingRule the sign can change even for too small or too large distances,
+    /// so if you would like to get closed mesh from marching cubes, set false here
+    bool nullOutsideMinMax = true;
+
+    /// only for SignDetectionMode::HoleWindingRule:
+    /// positive distance if winding number below or equal this threshold;
+    /// ideal threshold: 0.5 for closed meshes; 0.0 for planar meshes
+    float windingNumberThreshold = 0.5f;
+
+    /// only for SignDetectionMode::HoleWindingRule:
+    /// determines the precision of fast approximation: the more the better, minimum value is 1
+    float windingNumberBeta = 2;
+};
+
+/// options determining computation of signed distance from a point to a mesh
+struct SignedDistanceToMeshOptions : DistanceToMeshOptions
+{
+    /// the method to compute distance sign
+    SignDetectionMode signMode{ SignDetectionMode::ProjectionNormal };
+};
+
+} //namespace MR

source/MRMesh/MRFastWindingNumber.cpp

@@ -8,6 +8,7 @@
 #include "MRAABBTree.h"
 #include "MRDipole.h"
 #include "MRIsNaN.h"
+#include "MRDistanceToMeshOptions.h"
 
 namespace MR
 {
@@ -59,16 +60,16 @@ Expected<void> FastWindingNumber::calcFromGrid( std::vector<float>& res, const V
     return {};
 }
 
-float FastWindingNumber::calcWithDistances( const Vector3f& p, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq )
+float FastWindingNumber::calcWithDistances( const Vector3f& p, const DistanceToMeshOptions& options )
 {
-    auto resSq = findProjection( p, mesh_, maxDistSq, nullptr, minDistSq ).distSq;
-    if ( resSq < minDistSq || resSq >= maxDistSq ) // note that resSq == minDistSq (e.g. == 0) is a valid situation
+    auto resSq = findProjection( p, mesh_, options.maxDistSq, nullptr, options.minDistSq ).distSq;
+    if ( options.nullOutsideMinMax && ( resSq < options.minDistSq || resSq >= options.maxDistSq ) ) // note that resSq == minDistSq (e.g. == 0) is a valid situation
         return cQuietNan;
-    const auto sign = calc_( p, beta ) > windingNumberThreshold ? -1.f : +1.f;
+    const auto sign = calc_( p, options.windingNumberBeta ) > options.windingNumberThreshold ? -1.f : +1.f;
     return sign * std::sqrt( resSq );
 }
 
-Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb )
+Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb )
 {
     MR_TIMER
 
@@ -80,7 +81,7 @@ Expected<void> FastWindingNumber::calcFromGridWithDistances( std::vector<float>&
     if ( !ParallelFor( 0_vox, indexer.endId(), [&]( VoxelId i )
         {
             const auto transformedPoint = gridToMeshXf( Vector3f( indexer.toPos( i ) ) );
-            res[i] = calcWithDistances( transformedPoint, windingNumberThreshold, beta, maxDistSq, minDistSq );
+            res[i] = calcWithDistances( transformedPoint, options );
         }, cb ) )
         return unexpectedOperationCanceled();
     return {};

source/MRMesh/MRFastWindingNumber.h

@@ -45,10 +45,7 @@ class IFastWindingNumber
     /// </summary>
     /// <param name="res">resulting signed distances, will be resized automatically</param>
     /// <param name="dims">dimensions of the grid</param>
-    /// <param name="gridToMeshXf">transform from integer grid locations to voxel's centers in mesh reference frame</param>
-    /// <param name="windingNumberThreshold">positive distance if winding number below or equal this threshold</param>
-    /// <param name="beta">determines the precision of the approximation: the more the better, recommended value 2 or more</param>
-    virtual Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) = 0;
+    virtual Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) = 0;
 };
 
 /// the class for fast approximate computation of winding number for a mesh (using its AABB tree)
@@ -66,8 +63,8 @@ class MRMESH_CLASS FastWindingNumber : public IFastWindingNumber
     MRMESH_API void calcFromVector( std::vector<float>& res, const std::vector<Vector3f>& points, float beta, FaceId skipFace = {} ) override;
     MRMESH_API bool calcSelfIntersections( FaceBitSet& res, float beta, ProgressCallback cb ) override;
     MRMESH_API Expected<void> calcFromGrid( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float beta, ProgressCallback cb ) override;
-    MRMESH_API float calcWithDistances( const Vector3f& p, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq );
-    MRMESH_API Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) override;
+    MRMESH_API float calcWithDistances( const Vector3f& p, const DistanceToMeshOptions& options );
+    MRMESH_API Expected<void> calcFromGridWithDistances( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) override;
 
 private:
     [[nodiscard]] float calc_( const Vector3f & q, float beta, FaceId skipFace = {} ) const;

source/MRMesh/MRMesh.vcxproj

@@ -26,6 +26,7 @@
     <ClInclude Include="MRCylinderApproximator.h" />
     <ClInclude Include="MRCylinderObject.h" />
     <ClInclude Include="MRDipole.h" />
+    <ClInclude Include="MRDistanceToMeshOptions.h" />
     <ClInclude Include="MREnums.h" />
     <ClInclude Include="MRFaceDistance.h" />
     <ClInclude Include="MRFeatureObject.h" />

source/MRMesh/MRMesh.vcxproj.filters

@@ -1248,6 +1248,9 @@
     <ClInclude Include="MRMeshTopologyDiff.h">
       <Filter>Source Files\Mesh</Filter>
     </ClInclude>
+    <ClInclude Include="MRDistanceToMeshOptions.h">
+      <Filter>Source Files\AABBTree</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="MRParallelProgressReporter.cpp">

source/MRMesh/MRMeshDistance.cpp

@@ -9,14 +9,12 @@
 namespace MR
 {
 
-std::optional<float> signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const DistanceToMeshOptions& op )
+std::optional<float> signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const SignedDistanceToMeshOptions& op )
 {
     assert( op.signMode != SignDetectionMode::OpenVDB );
     const auto proj = findProjection( p, mp, op.maxDistSq, nullptr, op.minDistSq );
 
-    // please note that in HoleWindingRule the sign can change even for too small or too large distances,
-    // so if you would like to get closed mesh from marching cubes, do not change op.minDistSq and op.maxDistSq
-    if ( proj.distSq < op.minDistSq || proj.distSq >= op.maxDistSq ) // note that proj.distSq == op.minDistSq (e.g. == 0) is a valid situation
+    if ( op.nullOutsideMinMax && ( proj.distSq < op.minDistSq || proj.distSq >= op.maxDistSq ) ) // note that proj.distSq == op.minDistSq (e.g. == 0) is a valid situation
         return {}; // distance is too small or too large, discard them
 
     float dist = std::sqrt( proj.distSq );

source/MRMesh/MRMeshDistance.h

@@ -3,8 +3,7 @@
 // distance queries to one mesh only, please see MRMeshMeshDistance.h for queries involving two meshes
 
 #include "MRMeshPart.h"
-#include "MRSignDetectionMode.h"
-#include <cfloat>
+#include "MRDistanceToMeshOptions.h"
 #include <functional>
 #include <optional>
 
@@ -32,30 +31,9 @@ using TriangleCallback = std::function<ProcessOneResult( const Vector3f & p, Fac
 /// given squared distance from t-triangle
 MRMESH_API void processCloseTriangles( const MeshPart& mp, const Triangle3f & t, float rangeSq, const TriangleCallback & call );
 
-struct DistanceToMeshOptions
-{
-    /// minimum squared distance from a point to mesh
-    float minDistSq{ 0 };
-
-    /// maximum squared distance from a point to mesh
-    float maxDistSq{ FLT_MAX };
-
-    /// the method to compute distance sign
-    SignDetectionMode signMode{ SignDetectionMode::ProjectionNormal };
-
-    /// only for SignDetectionMode::HoleWindingRule:
-    /// positive distance if winding number below or equal this threshold;
-    /// ideal threshold: 0.5 for closed meshes; 0.0 for planar meshes
-    float windingNumberThreshold = 0.5f;
-
-    /// only for SignDetectionMode::HoleWindingRule:
-    /// determines the precision of fast approximation: the more the better, minimum value is 1
-    float windingNumberBeta = 2;
-};
-
 /// computes signed distance from point (p) to mesh part (mp) following options (op);
 /// returns std::nullopt if distance is smaller than op.minDist or larger than op.maxDist (except for op.signMode == HoleWindingRule)
-[[nodiscard]] MRMESH_API std::optional<float> signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const DistanceToMeshOptions& op );
+[[nodiscard]] MRMESH_API std::optional<float> signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const SignedDistanceToMeshOptions& op );
 
 /// \}
 

source/MRMesh/MRMeshFwd.h

@@ -540,6 +540,8 @@ MR_CANONICAL_TYPEDEFS( (template<typename V> struct [[nodiscard]]), PolylineProj
 )
 
 class DistanceMap;
+struct DistanceToMeshOptions;
+struct SignedDistanceToMeshOptions;
 
 using GcodeSource = std::vector<std::string>;
 

source/MRVoxels/MRMeshToDistanceVolume.cpp

@@ -52,9 +52,7 @@ Expected<SimpleVolumeMinMax> meshToDistanceVolume( const MeshPart& mp, const Mes
             params.fwn = std::make_shared<FastWindingNumber>( mp.mesh );
         assert( !mp.region ); // only whole mesh is supported for now
         auto basis = AffineXf3f( Matrix3f::scale( params.vol.voxelSize ), params.vol.origin + 0.5f * params.vol.voxelSize );
-        if ( auto d = params.fwn->calcFromGridWithDistances( res.data, res.dims, basis,
-            params.dist.windingNumberThreshold, params.dist.windingNumberBeta,
-            params.dist.maxDistSq, params.dist.minDistSq, params.vol.cb ); !d )
+        if ( auto d = params.fwn->calcFromGridWithDistances( res.data, res.dims, basis, params.dist, params.vol.cb ); !d )
         {
             return unexpected( std::move( d.error() ) );
         }

source/MRVoxels/MRMeshToDistanceVolume.h

@@ -13,7 +13,7 @@ struct MeshToDistanceVolumeParams
 {
     DistanceVolumeParams vol;
 
-    DistanceToMeshOptions dist;
+    SignedDistanceToMeshOptions dist;
 
     std::shared_ptr<IFastWindingNumber> fwn;
 };

source/MRVoxels/MROffset.cpp

@@ -155,11 +155,9 @@ Expected<Mesh> mcOffsetMesh( const MeshPart& mp, float offset,
         msParams.vol.voxelSize = Vector3f::diagonal( params.voxelSize );
         msParams.vol.dimensions = dimensions;
         msParams.dist.signMode = params.signDetectionMode;
-        if ( params.signDetectionMode != SignDetectionMode::HoleWindingRule || !params.closeHolesInHoleWindingNumber )
-        {
-            msParams.dist.maxDistSq = sqr( absOffset + 1.001f * params.voxelSize ); // we multiply by 1.001f to be sure not to have rounding errors (which may lead to unexpected NaN values )
-            msParams.dist.minDistSq = sqr( std::max( absOffset - 1.001f * params.voxelSize, 0.0f ) ); // we multiply by 1.001f to be sure not to have rounding errors (which may lead to unexpected NaN values )
-        }
+        msParams.dist.maxDistSq = sqr( absOffset + 1.001f * params.voxelSize ); // we multiply by 1.001f to be sure not to have rounding errors (which may lead to unexpected NaN values )
+        msParams.dist.minDistSq = sqr( std::max( absOffset - 1.001f * params.voxelSize, 0.0f ) ); // we multiply by 1.001f to be sure not to have rounding errors (which may lead to unexpected NaN values )
+        msParams.dist.nullOutsideMinMax = params.signDetectionMode != SignDetectionMode::HoleWindingRule || !params.closeHolesInHoleWindingNumber;
         msParams.dist.windingNumberThreshold = params.windingNumberThreshold;
         msParams.dist.windingNumberBeta = params.windingNumberBeta;
         msParams.fwn = params.fwn;