IFastWindingNumber::calcFromGridWithDistances returns NaN outside min…

…-max range (#3917)
MeshInspector · Dec 24, 2024 · 6e4b7b2 · 6e4b7b2
1 parent 9339fed
commit 6e4b7b2
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 6 deletions.
diff --git a/source/MRCuda/MRCudaFastWindingNumber.cu b/source/MRCuda/MRCudaFastWindingNumber.cu
@@ -2,6 +2,7 @@
 #include "MRMesh/MRAABBTree.h"
 #include "MRMesh/MRConstants.h"
 #include "device_launch_parameters.h"
+#include <limits>
 
 namespace MR
 {
@@ -59,7 +60,7 @@ __device__ void processPoint( const float3& q, float& res, const Dipole* dipoles
     res *= INV_4PI;
 }
 
-__device__ float calcDistance( const float3& pt,
+__device__ float calcDistanceSq( const float3& pt,
     const Node3* __restrict__ nodes, const float3* __restrict__ meshPoints, const FaceToThreeVerts* __restrict__ faces,
     float maxDistSq, float minDistSq )
 {
@@ -130,7 +131,7 @@ __device__ float calcDistance( const float3& pt,
         addSubTask( s1 ); // larger distance to look later
         addSubTask( s2 ); // smaller distance to look first
     }
-    return sqrt( resSq );
+    return resSq;
 }
 
 __global__ void fastWindingNumberFromVectorKernel( const float3* points,
@@ -195,6 +196,8 @@ __global__ void fastWindingNumberFromGridKernel( int3 dims, Matrix4 gridToMeshXf
     processPoint( transformedPoint, resVec[index], dipoles, nodes, meshPoints, faces, beta, index );
 }
 
+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 )
@@ -215,13 +218,19 @@ __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& res = resVec[index];
-    res = calcDistance( transformedPoint, nodes, meshPoints, faces, maxDistSq, minDistSq );
+    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
+    {
+        resVec[index] = cQuietNan;
+        return;
+    }
 
     float fwn{ 0 };
     processPoint( transformedPoint, fwn, dipoles, nodes, meshPoints, faces, beta, index );
+    float res = sqrt( resSq );
     if ( fwn > windingNumberThreshold )
         res = -res;
+    resVec[index] = res;
 }
 
 void fastWindingNumberFromVector( const float3* points, const Dipole* dipoles,

diff --git a/source/MRMesh/MRFastWindingNumber.cpp b/source/MRMesh/MRFastWindingNumber.cpp
@@ -7,6 +7,7 @@
 #include "MRMeshProject.h"
 #include "MRAABBTree.h"
 #include "MRDipole.h"
+#include "MRIsNaN.h"
 
 namespace MR
 {
@@ -60,8 +61,11 @@ Expected<void> FastWindingNumber::calcFromGrid( std::vector<float>& res, const V
 
 float FastWindingNumber::calcWithDistances( const Vector3f& p, float windingNumberThreshold, float beta, float maxDistSq, float minDistSq )
 {
+    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
+        return cQuietNan;
     const auto sign = calc_( p, beta ) > windingNumberThreshold ? -1.f : +1.f;
-    return sign * std::sqrt( findProjection( p, mesh_, maxDistSq, nullptr, minDistSq ).distSq );
+    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 )

diff --git a/source/MRMesh/MRFastWindingNumber.h b/source/MRMesh/MRFastWindingNumber.h
@@ -40,7 +40,8 @@ class IFastWindingNumber
     virtual Expected<void> calcFromGrid( std::vector<float>& res, const Vector3i& dims, const AffineXf3f& gridToMeshXf, float beta, ProgressCallback cb = {} ) = 0;
 
     /// <summary>
-    /// calculates distances with the sign obtained from winding number in each point from a three-dimensional grid
+    /// calculates distances with the sign obtained from generalized winding number in each point from a three-dimensional grid;
+    /// if sqr(res) < minDistSq or sqr(res) >= maxDistSq, then NaN is returned for such point
     /// </summary>
     /// <param name="res">resulting signed distances, will be resized automatically</param>
     /// <param name="dims">dimensions of the grid</param>