From 4739ccd3557fda52edad4c8a43fa2e057703630d Mon Sep 17 00:00:00 2001 From: Fedor Chelnokov Date: Thu, 26 Dec 2024 14:10:54 +0300 Subject: [PATCH] RebuildMesh: restore original performance (#3920) --- source/MRCuda/MRCudaFastWindingNumber.cpp | 4 +-- source/MRCuda/MRCudaFastWindingNumber.cu | 15 ++++---- source/MRCuda/MRCudaFastWindingNumber.cuh | 5 ++- source/MRCuda/MRCudaFastWindingNumber.h | 3 +- source/MRMesh/MRDistanceToMeshOptions.h | 42 ++++++++++++++++++++++ source/MRMesh/MRFastWindingNumber.cpp | 13 +++---- source/MRMesh/MRFastWindingNumber.h | 9 ++--- source/MRMesh/MRMesh.vcxproj | 1 + source/MRMesh/MRMesh.vcxproj.filters | 3 ++ source/MRMesh/MRMeshDistance.cpp | 6 ++-- source/MRMesh/MRMeshDistance.h | 26 ++------------ source/MRMesh/MRMeshFwd.h | 2 ++ source/MRVoxels/MRMeshToDistanceVolume.cpp | 4 +-- source/MRVoxels/MRMeshToDistanceVolume.h | 2 +- source/MRVoxels/MROffset.cpp | 8 ++--- 15 files changed, 83 insertions(+), 60 deletions(-) create mode 100644 source/MRMesh/MRDistanceToMeshOptions.h diff --git a/source/MRCuda/MRCudaFastWindingNumber.cpp b/source/MRCuda/MRCudaFastWindingNumber.cpp index 303c8c1b6b41..df5aef52b959 100644 --- a/source/MRCuda/MRCudaFastWindingNumber.cpp +++ b/source/MRCuda/MRCudaFastWindingNumber.cpp @@ -145,7 +145,7 @@ Expected FastWindingNumber::calcFromGrid( std::vector& res, const V return {}; } -Expected FastWindingNumber::calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) +Expected FastWindingNumber::calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) { MR_TIMER prepareData_( {} ); @@ -171,7 +171,7 @@ Expected FastWindingNumber::calcFromGridWithDistances( std::vector& 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 ) ); diff --git a/source/MRCuda/MRCudaFastWindingNumber.cu b/source/MRCuda/MRCudaFastWindingNumber.cu index dfd48cad6d80..d9d7f5c7bec0 100644 --- a/source/MRCuda/MRCudaFastWindingNumber.cu +++ b/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 @@ -200,7 +201,7 @@ static constexpr float cQuietNan = std::numeric_limits::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 diff --git a/source/MRCuda/MRCudaFastWindingNumber.cuh b/source/MRCuda/MRCudaFastWindingNumber.cuh index d23db49a0d2f..6cd8a63b7094 100644 --- a/source/MRCuda/MRCudaFastWindingNumber.cuh +++ b/source/MRCuda/MRCudaFastWindingNumber.cuh @@ -1,4 +1,5 @@ #pragma once + #include "MRCudaBasic.cuh" #include "MRCudaMath.cuh" #include "MRCudaFloat.cuh" @@ -6,6 +7,8 @@ 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 diff --git a/source/MRCuda/MRCudaFastWindingNumber.h b/source/MRCuda/MRCudaFastWindingNumber.h index 9f2876ab431e..8ec477791321 100644 --- a/source/MRCuda/MRCudaFastWindingNumber.h +++ b/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& res, const std::vector& points, float beta, FaceId skipFace = {} ) override; MRCUDA_API bool calcSelfIntersections( FaceBitSet& res, float beta, ProgressCallback cb ) override; MRCUDA_API Expected calcFromGrid( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float beta, ProgressCallback cb ) override; - MRCUDA_API Expected calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) override; + MRCUDA_API Expected calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) override; private: bool prepareData_( ProgressCallback cb ); diff --git a/source/MRMesh/MRDistanceToMeshOptions.h b/source/MRMesh/MRDistanceToMeshOptions.h new file mode 100644 index 000000000000..c8e8bfb1210c --- /dev/null +++ b/source/MRMesh/MRDistanceToMeshOptions.h @@ -0,0 +1,42 @@ +#pragma once + +#include "MRSignDetectionMode.h" +#include + +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 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 diff --git a/source/MRMesh/MRFastWindingNumber.cpp b/source/MRMesh/MRFastWindingNumber.cpp index 1c90c2020d7d..6b2602e5642a 100644 --- a/source/MRMesh/MRFastWindingNumber.cpp +++ b/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 FastWindingNumber::calcFromGrid( std::vector& 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 FastWindingNumber::calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) +Expected FastWindingNumber::calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, const DistanceToMeshOptions& options, const ProgressCallback& cb ) { MR_TIMER @@ -80,7 +81,7 @@ Expected FastWindingNumber::calcFromGridWithDistances( std::vector& 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 {}; diff --git a/source/MRMesh/MRFastWindingNumber.h b/source/MRMesh/MRFastWindingNumber.h index 62bbbc151bcf..468140799d3b 100644 --- a/source/MRMesh/MRFastWindingNumber.h +++ b/source/MRMesh/MRFastWindingNumber.h @@ -45,10 +45,7 @@ class IFastWindingNumber /// /// resulting signed distances, will be resized automatically /// dimensions of the grid - /// transform from integer grid locations to voxel's centers in mesh reference frame - /// positive distance if winding number below or equal this threshold - /// determines the precision of the approximation: the more the better, recommended value 2 or more - virtual Expected calcFromGridWithDistances( std::vector& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq, ProgressCallback cb ) = 0; + virtual Expected calcFromGridWithDistances( std::vector& 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& res, const std::vector& points, float beta, FaceId skipFace = {} ) override; MRMESH_API bool calcSelfIntersections( FaceBitSet& res, float beta, ProgressCallback cb ) override; MRMESH_API Expected calcFromGrid( std::vector& 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 calcFromGridWithDistances( std::vector& 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 calcFromGridWithDistances( std::vector& 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; diff --git a/source/MRMesh/MRMesh.vcxproj b/source/MRMesh/MRMesh.vcxproj index ef9c060c8653..2d984d3a74ee 100644 --- a/source/MRMesh/MRMesh.vcxproj +++ b/source/MRMesh/MRMesh.vcxproj @@ -26,6 +26,7 @@ + diff --git a/source/MRMesh/MRMesh.vcxproj.filters b/source/MRMesh/MRMesh.vcxproj.filters index 7e21f19dd6d7..fc78e09d16f6 100644 --- a/source/MRMesh/MRMesh.vcxproj.filters +++ b/source/MRMesh/MRMesh.vcxproj.filters @@ -1248,6 +1248,9 @@ Source Files\Mesh + + Source Files\AABBTree + diff --git a/source/MRMesh/MRMeshDistance.cpp b/source/MRMesh/MRMeshDistance.cpp index be161273fdda..0b8964dd5a51 100644 --- a/source/MRMesh/MRMeshDistance.cpp +++ b/source/MRMesh/MRMeshDistance.cpp @@ -9,14 +9,12 @@ namespace MR { -std::optional signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const DistanceToMeshOptions& op ) +std::optional 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 ); diff --git a/source/MRMesh/MRMeshDistance.h b/source/MRMesh/MRMeshDistance.h index 91d8745763cb..fcf9809fed1e 100644 --- a/source/MRMesh/MRMeshDistance.h +++ b/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 +#include "MRDistanceToMeshOptions.h" #include #include @@ -32,30 +31,9 @@ using TriangleCallback = std::function signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const DistanceToMeshOptions& op ); +[[nodiscard]] MRMESH_API std::optional signedDistanceToMesh( const MeshPart& mp, const Vector3f& p, const SignedDistanceToMeshOptions& op ); /// \} diff --git a/source/MRMesh/MRMeshFwd.h b/source/MRMesh/MRMeshFwd.h index 9c0d2ea3fef2..d7e7eb7d9bc1 100644 --- a/source/MRMesh/MRMeshFwd.h +++ b/source/MRMesh/MRMeshFwd.h @@ -540,6 +540,8 @@ MR_CANONICAL_TYPEDEFS( (template struct [[nodiscard]]), PolylineProj ) class DistanceMap; +struct DistanceToMeshOptions; +struct SignedDistanceToMeshOptions; using GcodeSource = std::vector; diff --git a/source/MRVoxels/MRMeshToDistanceVolume.cpp b/source/MRVoxels/MRMeshToDistanceVolume.cpp index 92b7bf395afa..bc43e94cdd16 100644 --- a/source/MRVoxels/MRMeshToDistanceVolume.cpp +++ b/source/MRVoxels/MRMeshToDistanceVolume.cpp @@ -52,9 +52,7 @@ Expected meshToDistanceVolume( const MeshPart& mp, const Mes params.fwn = std::make_shared( 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() ) ); } diff --git a/source/MRVoxels/MRMeshToDistanceVolume.h b/source/MRVoxels/MRMeshToDistanceVolume.h index 95b44df873ad..8061ec8e25a0 100644 --- a/source/MRVoxels/MRMeshToDistanceVolume.h +++ b/source/MRVoxels/MRMeshToDistanceVolume.h @@ -13,7 +13,7 @@ struct MeshToDistanceVolumeParams { DistanceVolumeParams vol; - DistanceToMeshOptions dist; + SignedDistanceToMeshOptions dist; std::shared_ptr fwn; }; diff --git a/source/MRVoxels/MROffset.cpp b/source/MRVoxels/MROffset.cpp index bb151b1741dc..346c147214ab 100644 --- a/source/MRVoxels/MROffset.cpp +++ b/source/MRVoxels/MROffset.cpp @@ -155,11 +155,9 @@ Expected 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;