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

axel/axel/MeshToSdf.cpp

@@ -13,6 +13,7 @@
 #include <iostream>
 #include <limits>
 #include <memory>
+#include <numbers>
 #include <random>
 #include <unordered_map>
 
@@ -49,11 +50,8 @@ struct VoxelCandidate {
 };
 
 // ================================================================================================
-// STEP 2: FAST MARCHING PROPAGATION
-// UTILITY FUNCTIONS
-// UTILITY FUNCTIONS
-// UTILITY FUNCTIONS
-// =========================================================================================
+// UTILITY TYPES
+// ================================================================================================
 
 /**
  * Voxel states for fast marching algorithm.
@@ -137,7 +135,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 +672,102 @@ 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.
+ *
+ * O(T) per query — brute-force over all triangles.
+ */
+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};
+  const ScalarType kFourPi = ScalarType{4} * std::numbers::pi_v<ScalarType>;
+
+  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 (only needed for RayCasting)
+  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 +902,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;
 };
 
 /**
@@ -60,7 +81,7 @@ struct MeshToSdfPadding {
  * Convert a triangle mesh to a signed distance field using modern 3-step approach:
  * 1. Narrow band initialization with exact triangle distances
  * 2. Fast marching propagation using Eikonal equation
- * 3. Sign determination using ray casting
+ * 3. Sign determination (configurable via `SignMethod`)
  *
  * @param vertices Vertex positions as span (works with std::vector, arrays, subranges)
  * @param triangles Triangle indices as span (indices must be valid within vertices)
@@ -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.