Add selectable inside/outside method for SDF sign determination (#1393)

Summary:

Add a `SignMethod` enum to `MeshToSdfConfig` with three options:

- `RayCasting` (default): existing 6-direction majority vote. Fast but fails on non-watertight meshes where rays pass through gaps.
- `WindingNumber`: signed winding number (wn > 0.5), assumes consistent outward-facing normals. Will reject meshes with accidentally flipped normals rather than silently accepting them.
- `WindingNumberPermissive`: absolute winding number (abs(wn) > 0.5), tolerant of either winding convention. Handles mesh soups, mixed orientations, and geometry from multiple sources.

Both winding number variants use the Van Oosterom-Strackee solid angle formula, computed per voxel over all triangles. No new dependencies — the computation is self-contained.

The sign determination step was already isolated in `applySignsToDistanceField`, so the change is minimal: the function now accepts a `SignMethod` parameter and dispatches to the appropriate classifier.

Reviewed By: cstollmeta

Differential Revision: D104091366

Author: cdtwigg

axel/axel/MeshToSdf.cpp

@@ -49,11 +49,8 @@ struct VoxelCandidate {
 };
 
 // ================================================================================================
-// STEP 2: FAST MARCHING PROPAGATION
-// UTILITY FUNCTIONS
-// UTILITY FUNCTIONS
-// UTILITY FUNCTIONS
-// =========================================================================================
+// UTILITY TYPES
+// ================================================================================================
 
 /**
  * Voxel states for fast marching algorithm.
@@ -137,7 +134,7 @@ void meshToSdfImpl(
   if (config.verbose) {
     std::cout << "Step 3: Applying signs..." << std::endl;
   }
-  applySignsToDistanceField(sdf, vertices, triangles);
+  applySignsToDistanceField(sdf, vertices, triangles, config.signMethod);
   if (config.verbose) {
     std::cout << "Signs applied" << std::endl;
   }
@@ -674,34 +671,100 @@ bool isPointInsideByRayCasting(
   return insideCount > (directions.size() / 2);
 }
 
+/**
+ * Compute the solid angle subtended by a triangle as seen from a point.
+ * Uses the Van Oosterom-Strackee formula.
+ */
+template <typename ScalarType>
+ScalarType triangleSolidAngle(
+    const Eigen::Vector3<ScalarType>& point,
+    const Eigen::Vector3<ScalarType>& v0,
+    const Eigen::Vector3<ScalarType>& v1,
+    const Eigen::Vector3<ScalarType>& v2) {
+  const Eigen::Vector3<ScalarType> a = v0 - point;
+  const Eigen::Vector3<ScalarType> b = v1 - point;
+  const Eigen::Vector3<ScalarType> c = v2 - point;
+
+  const ScalarType la = a.norm();
+  const ScalarType lb = b.norm();
+  const ScalarType lc = c.norm();
+
+  // Degenerate case: point is on a vertex
+  constexpr ScalarType kEps = std::numeric_limits<ScalarType>::epsilon() * ScalarType{100};
+  if (la < kEps || lb < kEps || lc < kEps) {
+    return ScalarType{0};
+  }
+
+  const ScalarType numerator = a.dot(b.cross(c));
+  const ScalarType denominator = la * lb * lc + a.dot(b) * lc + b.dot(c) * la + c.dot(a) * lb;
+
+  return ScalarType{2} * std::atan2(numerator, denominator);
+}
+
+/**
+ * Compute the generalized winding number at a point.
+ * Sums the solid angle of each triangle as seen from the query point,
+ * normalized by 4pi. Returns +1 for inside (outward normals), -1 for
+ * inside (inward normals), ~0 for outside.
+ */
+template <typename ScalarType>
+ScalarType computeWindingNumber(
+    const Eigen::Vector3<ScalarType>& point,
+    std::span<const Eigen::Vector3<ScalarType>> vertices,
+    std::span<const Eigen::Vector3i> triangles) {
+  ScalarType totalSolidAngle{0};
+  constexpr ScalarType kFourPi = ScalarType{4} * ScalarType{M_PI};
+
+  for (const auto& tri : triangles) {
+    totalSolidAngle +=
+        triangleSolidAngle(point, vertices[tri.x()], vertices[tri.y()], vertices[tri.z()]);
+  }
+
+  return totalSolidAngle / kFourPi;
+}
+
 template <typename ScalarType>
 void applySignsToDistanceField(
     SignedDistanceField<ScalarType>& sdf,
     std::span<const Eigen::Vector3<ScalarType>> vertices,
-    std::span<const Eigen::Vector3i> triangles) {
+    std::span<const Eigen::Vector3i> triangles,
+    SignMethod signMethod) {
   const auto& resolution = sdf.resolution();
 
-  // BUILD BVH ONCE for efficient ray casting
-  // Convert spans to matrices for TriBvh constructor
-  Eigen::MatrixX3<ScalarType> vertexMatrix(vertices.size(), 3);
-  for (size_t i = 0; i < vertices.size(); ++i) {
-    vertexMatrix.row(i) = vertices[i];
-  }
+  // Build BVH for ray casting (only needed for RayCasting method, but cheap to build)
+  std::unique_ptr<TriBvh<ScalarType>> bvh;
+  if (signMethod == SignMethod::RayCasting) {
+    Eigen::MatrixX3<ScalarType> vertexMatrix(vertices.size(), 3);
+    for (size_t i = 0; i < vertices.size(); ++i) {
+      vertexMatrix.row(i) = vertices[i];
+    }
 
-  Eigen::MatrixX3i triangleMatrix(triangles.size(), 3);
-  for (size_t i = 0; i < triangles.size(); ++i) {
-    triangleMatrix.row(i) = triangles[i];
-  }
+    Eigen::MatrixX3i triangleMatrix(triangles.size(), 3);
+    for (size_t i = 0; i < triangles.size(); ++i) {
+      triangleMatrix.row(i) = triangles[i];
+    }
 
-  const TriBvh<ScalarType> bvh(std::move(vertexMatrix), std::move(triangleMatrix));
+    bvh = std::make_unique<TriBvh<ScalarType>>(std::move(vertexMatrix), std::move(triangleMatrix));
+  }
 
   // Process each voxel
   const auto processVoxel = [&](Index i, Index j, Index k) {
     const Eigen::Vector3<ScalarType> gridPos(
         static_cast<ScalarType>(i), static_cast<ScalarType>(j), static_cast<ScalarType>(k));
     const auto worldPos = sdf.gridToWorld(gridPos);
 
-    bool isInside = isPointInsideByRayCasting(worldPos, bvh);
+    bool isInside = false;
+    switch (signMethod) {
+      case SignMethod::RayCasting:
+        isInside = isPointInsideByRayCasting(worldPos, *bvh);
+        break;
+      case SignMethod::WindingNumber:
+        isInside = computeWindingNumber(worldPos, vertices, triangles) > ScalarType{0.5};
+        break;
+      case SignMethod::WindingNumberPermissive:
+        isInside = std::abs(computeWindingNumber(worldPos, vertices, triangles)) > ScalarType{0.5};
+        break;
+    }
 
     // Apply sign: negative inside, positive outside
     if (isInside) {
@@ -836,12 +899,14 @@ template void fastMarchingPropagate<double>(SignedDistanceField<double>&);
 template void applySignsToDistanceField<float>(
     SignedDistanceField<float>&,
     std::span<const Eigen::Vector3<float>>,
-    std::span<const Eigen::Vector3i>);
+    std::span<const Eigen::Vector3i>,
+    SignMethod);
 
 template void applySignsToDistanceField<double>(
     SignedDistanceField<double>&,
     std::span<const Eigen::Vector3<double>>,
-    std::span<const Eigen::Vector3i>);
+    std::span<const Eigen::Vector3i>,
+    SignMethod);
 
 template BoundingBox<float> computeMeshBounds<float>(std::span<const Eigen::Vector3<float>>);
 

axel/axel/MeshToSdf.h

@@ -20,6 +20,24 @@
 
 namespace axel {
 
+/**
+ * Method for determining inside/outside classification in SDF sign determination.
+ */
+enum class SignMethod {
+  /// Ray casting with 6-direction majority vote. Fast but fails on non-watertight meshes.
+  RayCasting,
+
+  /// Winding number assuming consistent outward-facing normals (wn > 0.5).
+  /// Will reject meshes with accidentally flipped normals rather than silently
+  /// accepting them. Use when triangle orientation is known to be correct.
+  WindingNumber,
+
+  /// Winding number tolerant of either winding convention (abs(wn) > 0.5).
+  /// Handles mesh soups, mixed orientations, and geometry from multiple sources
+  /// where consistent winding cannot be guaranteed.
+  WindingNumberPermissive,
+};
+
 /**
  * Configuration parameters for mesh-to-SDF conversion.
  */
@@ -39,6 +57,9 @@ struct MeshToSdfConfig {
 
   /// Print progress messages to stdout
   bool verbose = false;
+
+  /// Method for inside/outside classification during sign determination.
+  SignMethod signMethod = SignMethod::RayCasting;
 };
 
 /**
@@ -141,7 +162,8 @@ template <typename ScalarType>
 void applySignsToDistanceField(
     SignedDistanceField<ScalarType>& sdf,
     std::span<const Eigen::Vector3<ScalarType>> vertices,
-    std::span<const Eigen::Vector3i> triangles);
+    std::span<const Eigen::Vector3i> triangles,
+    SignMethod signMethod = SignMethod::RayCasting);
 
 /**
  * Compute mesh bounding box from vertex spans.

axel/axel/test/MeshToSdfTest.cpp

@@ -448,37 +448,97 @@ TEST_F(MeshToSdfTest, Step3_SignDetermination_InsideOutsideAccuracy) {
   EXPECT_GT(clearlyOutsideVoxels, 0) << "Should have some clearly outside voxels";
 }
 
-TEST_F(MeshToSdfTest, Step3_SignDetermination_WindingNumbers) {
-  // Test sign determination using winding numbers
+TEST_F(MeshToSdfTest, Step3_WindingNumber_OutwardNormals) {
+  // Test WindingNumber (signed, wn > 0.5) with a cube that has outward-facing normals.
+  // Flip the test cube's winding order so normals point outward.
+  std::vector<Eigen::Vector3i> outwardFaces;
+  outwardFaces.reserve(cubeFaces.size());
+  for (const auto& f : cubeFaces) {
+    outwardFaces.emplace_back(f.x(), f.z(), f.y());
+  }
+
   const BoundingBoxf bounds(
       Eigen::Vector3f(-1.0f, -1.0f, -1.0f), Eigen::Vector3f(1.0f, 1.0f, 1.0f));
   const Eigen::Vector3<Index> resolution(10, 10, 10);
 
   MeshToSdfConfigf config;
   config.narrowBandWidth = 2.0f;
+  config.signMethod = SignMethod::WindingNumber;
 
-  // Generate complete SDF
+  const auto sdf = meshToSdf<float>(
+      std::span<const Eigen::Vector3f>(cubeVertices),
+      std::span<const Eigen::Vector3i>(outwardFaces),
+      bounds,
+      resolution,
+      config);
+
+  const float centerValue = sdf.sample(Eigen::Vector3f(0.0f, 0.0f, 0.0f));
+  EXPECT_LT(centerValue, 0.0f) << "Center should be inside with outward normals";
+
+  const float outsideValue = sdf.sample(Eigen::Vector3f(0.9f, 0.9f, 0.9f));
+  EXPECT_GT(outsideValue, 0.0f) << "Outside point should be positive";
+}
+
+TEST_F(MeshToSdfTest, Step3_WindingNumber_RejectsFlippedNormals) {
+  // WindingNumber (signed) should classify inside as outside when normals point inward,
+  // since winding number will be -1 (not > 0.5). This is the desired behavior —
+  // it catches orientation errors rather than silently accepting them.
+  const BoundingBoxf bounds(
+      Eigen::Vector3f(-1.0f, -1.0f, -1.0f), Eigen::Vector3f(1.0f, 1.0f, 1.0f));
+  const Eigen::Vector3<Index> resolution(10, 10, 10);
+
+  MeshToSdfConfigf config;
+  config.narrowBandWidth = 2.0f;
+  config.signMethod = SignMethod::WindingNumber;
+
+  // The test cube has inward-pointing normals
   const auto sdf = meshToSdf<float>(
       std::span<const Eigen::Vector3f>(cubeVertices),
       std::span<const Eigen::Vector3i>(cubeFaces),
       bounds,
       resolution,
       config);
 
+  // Center should be classified as OUTSIDE (positive) because normals are flipped
+  const float centerValue = sdf.sample(Eigen::Vector3f(0.0f, 0.0f, 0.0f));
+  EXPECT_GT(centerValue, 0.0f) << "Signed winding number should reject inward normals";
+}
+
+TEST_F(MeshToSdfTest, Step3_WindingNumberPermissive_AcceptsBothOrientations) {
+  // WindingNumberPermissive should work regardless of normal orientation.
+  const BoundingBoxf bounds(
+      Eigen::Vector3f(-1.0f, -1.0f, -1.0f), Eigen::Vector3f(1.0f, 1.0f, 1.0f));
+  const Eigen::Vector3<Index> resolution(10, 10, 10);
+
+  MeshToSdfConfigf config;
+  config.narrowBandWidth = 2.0f;
+  config.signMethod = SignMethod::WindingNumberPermissive;
+
+  const float boundaryThreshold = 0.05f;
   int correctSigns = 0;
   int totalVoxels = 0;
 
+  // Test with the inward-normal cube — Permissive should handle it
+  const auto sdf = meshToSdf<float>(
+      std::span<const Eigen::Vector3f>(cubeVertices),
+      std::span<const Eigen::Vector3i>(cubeFaces),
+      bounds,
+      resolution,
+      config);
+
   for (Index i = 0; i < resolution.x(); ++i) {
     for (Index j = 0; j < resolution.y(); ++j) {
       for (Index k = 0; k < resolution.z(); ++k) {
-        const float sdfValue = sdf.at(i, j, k);
-
         const Eigen::Vector3f gridPos(
             static_cast<float>(i), static_cast<float>(j), static_cast<float>(k));
         const Eigen::Vector3f worldPos = sdf.gridToWorld(gridPos);
 
+        if (exactDistanceToUnitCube(worldPos) < boundaryThreshold) {
+          continue;
+        }
+
         const bool shouldBeInside = isInsideUnitCube(worldPos);
-        const bool sdfSaysInside = sdfValue < 0.0f;
+        const bool sdfSaysInside = sdf.at(i, j, k) < 0.0f;
 
         if (shouldBeInside == sdfSaysInside) {
           correctSigns++;
@@ -489,10 +549,91 @@ TEST_F(MeshToSdfTest, Step3_SignDetermination_WindingNumbers) {
   }
 
   const float accuracy = static_cast<float>(correctSigns) / totalVoxels;
-  std::cout << "Step 3 (Winding): Sign accuracy: " << accuracy * 100.0f << "% (" << correctSigns
-            << "/" << totalVoxels << ")" << std::endl;
+  EXPECT_GT(accuracy, 0.95f) << "Permissive winding number should handle either orientation";
+}
+
+TEST_F(MeshToSdfTest, Step3_WindingNumberPermissive_NonWatertightMesh) {
+  // Test that permissive winding numbers handle a non-watertight mesh.
+  std::vector<Eigen::Vector3i> openCubeFaces(cubeFaces.begin(), cubeFaces.end());
+  openCubeFaces.erase(openCubeFaces.begin() + 2, openCubeFaces.begin() + 4);
+
+  const BoundingBoxf bounds(
+      Eigen::Vector3f(-1.0f, -1.0f, -1.0f), Eigen::Vector3f(1.0f, 1.0f, 1.0f));
+  const Eigen::Vector3<Index> resolution(10, 10, 10);
+
+  MeshToSdfConfigf config;
+  config.narrowBandWidth = 2.0f;
+  config.signMethod = SignMethod::WindingNumberPermissive;
+
+  const auto sdf = meshToSdf<float>(
+      std::span<const Eigen::Vector3f>(cubeVertices),
+      std::span<const Eigen::Vector3i>(openCubeFaces),
+      bounds,
+      resolution,
+      config);
+
+  const float centerValue = sdf.sample(Eigen::Vector3f(0.0f, 0.0f, 0.0f));
+  EXPECT_LT(centerValue, 0.0f) << "Center should be inside even with missing face";
+
+  const float outsideValue = sdf.sample(Eigen::Vector3f(0.9f, 0.9f, 0.9f));
+  EXPECT_GT(outsideValue, 0.0f) << "Clearly outside point should be positive";
+}
+
+TEST_F(MeshToSdfTest, Step3_WindingNumberPermissive_MatchesRayCastingOnWatertight) {
+  // On a watertight mesh, both methods should produce the same sign classification.
+  const BoundingBoxf bounds(
+      Eigen::Vector3f(-1.0f, -1.0f, -1.0f), Eigen::Vector3f(1.0f, 1.0f, 1.0f));
+  const Eigen::Vector3<Index> resolution(10, 10, 10);
+
+  MeshToSdfConfigf configRay;
+  configRay.narrowBandWidth = 2.0f;
+  configRay.signMethod = SignMethod::RayCasting;
+
+  MeshToSdfConfigf configWind;
+  configWind.narrowBandWidth = 2.0f;
+  configWind.signMethod = SignMethod::WindingNumberPermissive;
+
+  const auto sdfRay = meshToSdf<float>(
+      std::span<const Eigen::Vector3f>(cubeVertices),
+      std::span<const Eigen::Vector3i>(cubeFaces),
+      bounds,
+      resolution,
+      configRay);
+
+  const auto sdfWind = meshToSdf<float>(
+      std::span<const Eigen::Vector3f>(cubeVertices),
+      std::span<const Eigen::Vector3i>(cubeFaces),
+      bounds,
+      resolution,
+      configWind);
+
+  const float boundaryThreshold = 0.05f;
+  int agreements = 0;
+  int total = 0;
+
+  for (Index i = 0; i < resolution.x(); ++i) {
+    for (Index j = 0; j < resolution.y(); ++j) {
+      for (Index k = 0; k < resolution.z(); ++k) {
+        const Eigen::Vector3f gridPos(
+            static_cast<float>(i), static_cast<float>(j), static_cast<float>(k));
+        const Eigen::Vector3f worldPos = sdfRay.gridToWorld(gridPos);
+
+        if (exactDistanceToUnitCube(worldPos) < boundaryThreshold) {
+          continue;
+        }
+
+        const bool rayInside = sdfRay.at(i, j, k) < 0.0f;
+        const bool windInside = sdfWind.at(i, j, k) < 0.0f;
+        if (rayInside == windInside) {
+          agreements++;
+        }
+        total++;
+      }
+    }
+  }
 
-  EXPECT_GT(accuracy, 0.95f) << "Winding number sign determination should be very accurate";
+  const float agreement = static_cast<float>(agreements) / total;
+  EXPECT_GT(agreement, 0.99f) << "Both methods should agree on watertight mesh";
 }
 
 TEST_F(MeshToSdfTest, IntegratedTest_CubeSDFProperties) {