SDF interpolation (#866)

* start of test

* started FindSurface

* beautiful smoothing

* implemented ITP root finding

* improved ITP

* added precision to SDF API

* test comparison to linear interpolation

* fixed quad 4-manifold edge

* fixed bug

* optimized ITP

* fix NaN tangents

.vscode/settings.json

@@ -114,7 +114,11 @@
     "numbers": "cpp",
     "ranges": "cpp",
     "stop_token": "cpp",
-    "execution": "cpp"
+    "execution": "cpp",
+    "csetjmp": "cpp",
+    "cuchar": "cpp",
+    "propagate_const": "cpp",
+    "scoped_allocator": "cpp"
   },
   "C_Cpp.clang_format_fallbackStyle": "google",
   "editor.formatOnSave": true,

src/manifold/src/properties.cpp

@@ -146,19 +146,6 @@ struct CheckHalfedges {
   }
 };
 
-struct NoDuplicates {
-  VecView<const Halfedge> halfedges;
-
-  bool operator()(size_t edge) const {
-    const Halfedge halfedge = halfedges[edge];
-    if (halfedge.startVert == -1 && halfedge.endVert == -1 &&
-        halfedge.pairedHalfedge == -1)
-      return true;
-    return halfedge.startVert != halfedges[edge + 1].startVert ||
-           halfedge.endVert != halfedges[edge + 1].endVert;
-  }
-};
-
 struct CheckCCW {
   VecView<const Halfedge> halfedges;
   VecView<const glm::vec3> vertPos;
@@ -226,8 +213,14 @@ bool Manifold::Impl::Is2Manifold() const {
   Vec<Halfedge> halfedge(halfedge_);
   stable_sort(halfedge.begin(), halfedge.end());
 
-  return all_of(countAt(0_z), countAt(2 * NumEdge() - 1),
-                NoDuplicates({halfedge}));
+  return all_of(
+      countAt(0_z), countAt(2 * NumEdge() - 1), [halfedge](size_t edge) {
+        const Halfedge h = halfedge[edge];
+        if (h.startVert == -1 && h.endVert == -1 && h.pairedHalfedge == -1)
+          return true;
+        return h.startVert != halfedge[edge + 1].startVert ||
+               h.endVert != halfedge[edge + 1].endVert;
+      });
 }
 
 /**

src/manifold/src/smoothing.cpp

@@ -736,7 +736,7 @@ void Manifold::Impl::DistributeTangents(const Vec<bool>& fixedHalfedges) {
           lastTangent = thisTangent;
         } while (!fixedHalfedges[current]);
 
-        if (currentAngle.size() == 1) return;
+        if (currentAngle.size() == 1 || glm::dot(normal, normal) == 0) return;
 
         const float scale = currentAngle.back() / desiredAngle.back();
         float offset = 0;

src/manifold/src/subdivision.cpp

@@ -89,19 +89,19 @@ class Partition {
     return partition;
   }
 
-  Vec<glm::ivec3> Reindex(glm::ivec4 tri, glm::ivec4 edgeOffsets,
+  Vec<glm::ivec3> Reindex(glm::ivec4 triVerts, glm::ivec4 edgeOffsets,
                           glm::bvec4 edgeFwd, int interiorOffset) const {
     Vec<int> newVerts;
     newVerts.reserve(vertBary.size());
     glm::ivec4 triIdx = idx;
     glm::ivec4 outTri = {0, 1, 2, 3};
-    if (tri[3] < 0 && idx[1] != Next3(idx[0])) {
+    if (triVerts[3] < 0 && idx[1] != Next3(idx[0])) {
       triIdx = {idx[2], idx[0], idx[1], idx[3]};
       edgeFwd = glm::not_(edgeFwd);
       std::swap(outTri[0], outTri[1]);
     }
     for (const int i : {0, 1, 2, 3}) {
-      if (tri[triIdx[i]] >= 0) newVerts.push_back(tri[triIdx[i]]);
+      if (triVerts[triIdx[i]] >= 0) newVerts.push_back(triVerts[triIdx[i]]);
     }
     for (const int i : {0, 1, 2, 3}) {
       const int n = sortedDivisions[i] - 1;
@@ -413,20 +413,12 @@ glm::ivec4 Manifold::Impl::GetHalfedges(int tri) const {
     if (pair / 3 < tri) {
       return glm::ivec4(-1);  // only process lower tri index
     }
-    glm::ivec2 otherHalf;
-    otherHalf[0] = NextHalfedge(pair);
-    otherHalf[1] = NextHalfedge(otherHalf[0]);
-    halfedges[neighbor] = otherHalf[0];
-    if (neighbor == 2) {
-      halfedges[3] = otherHalf[1];
-    } else if (neighbor == 1) {
-      halfedges[3] = halfedges[2];
-      halfedges[2] = otherHalf[1];
-    } else {
-      halfedges[3] = halfedges[2];
-      halfedges[2] = halfedges[1];
-      halfedges[1] = otherHalf[1];
-    }
+    // The order here matters to keep small quads split the way they started, or
+    // else it can create a 4-manifold edge.
+    halfedges[2] = NextHalfedge(halfedges[neighbor]);
+    halfedges[3] = NextHalfedge(halfedges[2]);
+    halfedges[0] = NextHalfedge(pair);
+    halfedges[1] = NextHalfedge(halfedges[0]);
   }
   return halfedges;
 }
@@ -451,10 +443,9 @@ Manifold::Impl::BaryIndices Manifold::Impl::GetIndices(int halfedge) const {
     const int pair = halfedge_[3 * tri + neighbor].pairedHalfedge;
     if (pair / 3 < tri) {
       tri = pair / 3;
-      const int j = pair % 3;
-      idx = Next3(neighbor) == idx ? j : (j + 1) % 4;
-    } else if (idx > neighbor) {
-      ++idx;
+      idx = Next3(neighbor) == idx ? 0 : 1;
+    } else {
+      idx = Next3(neighbor) == idx ? 2 : 3;
     }
     return {tri, idx, (idx + 1) % 4};
   }

src/sdf/include/sdf.h

@@ -20,5 +20,6 @@
 
 namespace manifold {
 Mesh LevelSet(std::function<float(glm::vec3)> sdf, Box bounds, float edgeLength,
-              float level = 0, bool canParallel = true);
+              float level = 0, bool canParallel = true,
+              float precision = std::numeric_limits<float>::infinity());
 }

src/sdf/src/sdf.cpp

@@ -131,6 +131,7 @@ struct ComputeVerts {
   const glm::ivec3 gridSize;
   const glm::vec3 spacing;
   const float level;
+  const float tol;
 
   inline glm::vec3 Position(glm::ivec4 gridIndex) const {
     return origin +
@@ -150,6 +151,49 @@ struct ComputeVerts {
     return d;
   }
 
+  // Simplified ITP root finding algorithm - same worst-case performance as
+  // bisection, better average performance.
+  inline glm::vec3 FindSurface(glm::vec3 pos0, float d0, glm::vec3 pos1,
+                               float d1) const {
+    if (d0 == 0) {
+      return pos0;
+    } else if (d1 == 0) {
+      return pos1;
+    }
+
+    // Sole tuning parameter, k: (0, 1) - smaller value gets better median
+    // performance, but also hits the worst case more often.
+    const float k = 0.1;
+    const float tol2 = 4 * tol * tol;
+    const float dist2 = glm::dot(pos0 - pos1, pos0 - pos1);
+    float frac = 1;
+    float biFrac = 1;
+    do {
+      const float a = d0 / (d0 - d1);
+      if (dist2 * frac * frac < tol2) {
+        return glm::mix(pos0, pos1, a);
+      }
+
+      const float t = glm::mix(a, 0.5f, k);
+      const float r = biFrac / frac - 0.5;
+      const float x = glm::abs(t - 0.5) < r ? t : 0.5 - r * (t < 0.5 ? 1 : -1);
+
+      const glm::vec3 mid = glm::mix(pos0, pos1, x);
+      const float d = sdf(mid) - level;
+
+      if ((d > 0) == (d0 > 0)) {
+        d0 = d;
+        pos0 = mid;
+        frac *= 1 - x;
+      } else {
+        d1 = d;
+        pos1 = mid;
+        frac *= x;
+      }
+      biFrac /= 2;
+    } while (1);
+  }
+
   inline void operator()(Uint64 mortonCode) {
     ZoneScoped;
     if (gridVerts.Full()) return;
@@ -177,9 +221,8 @@ struct ComputeVerts {
       keep = true;
 
       const int idx = AtomicAdd(vertIndex[0], 1);
-      vertPos[idx] =
-          (val * position - gridVert.distance * Position(neighborIndex)) /
-          (val - gridVert.distance);
+      vertPos[idx] = FindSurface(position, gridVert.distance,
+                                 Position(neighborIndex), val);
       gridVert.edgeVerts[i] = idx;
     }
 
@@ -298,12 +341,16 @@ namespace manifold {
  * with runtime locks that expect to not be called back by unregistered threads.
  * This allows bindings use LevelSet despite being compiled with MANIFOLD_PAR
  * active.
+ * @param precision Ensure each vertex is within this distance of the true
+ * surface. Defaults to infinity, which will return the interpolated
+ * crossing-point based on the two nearest grid points. Smaller values will
+ * require more sdf evaluations per output vertex.
  * @return Mesh This class does not depend on Manifold, so it just returns a
  * Mesh, but it is guaranteed to be manifold and so can always be used as
  * input to the Manifold constructor for further operations.
  */
 Mesh LevelSet(std::function<float(glm::vec3)> sdf, Box bounds, float edgeLength,
-              float level, bool canParallel) {
+              float level, bool canParallel, float precision) {
   Mesh out;
 
   const glm::vec3 dim = bounds.Size();
@@ -327,7 +374,7 @@ Mesh LevelSet(std::function<float(glm::vec3)> sdf, Box bounds, float edgeLength,
     Vec<int> index(1, 0);
     for_each_n(pol, countAt(0_z), maxMorton + 1,
                ComputeVerts({vertPos, index, gridVerts.D(), sdf, bounds.min,
-                             gridSize + 1, spacing, level}));
+                             gridSize + 1, spacing, level, precision}));
 
     if (gridVerts.Full()) {  // Resize HashTable
       const glm::vec3 lastVert = vertPos[index[0] - 1];

test/smooth_test.cpp

@@ -373,3 +373,60 @@ TEST(Smooth, SineSurface) {
   }
 #endif
 }
+
+TEST(Smooth, SDF) {
+  const float r = 10;
+  const float extra = 2;
+
+  auto sphericalGyroid = [r](glm::vec3 p) {
+    const float gyroid =
+        cos(p.x) * sin(p.y) + cos(p.y) * sin(p.z) + cos(p.z) * sin(p.x);
+    const float d = glm::min(0.0f, r - glm::length(p));
+    return gyroid - d * d;
+  };
+
+  auto gradient = [r](glm::vec3 pos) {
+    const float rad = glm::length(pos);
+    const float d = glm::min(0.0f, r - rad) / (rad > 0 ? rad : 1);
+    const glm::vec3 sphereGrad = 2 * d * pos;
+    const glm::vec3 gyroidGrad(
+        cos(pos.z) * cos(pos.x) - sin(pos.x) * sin(pos.y),
+        cos(pos.x) * cos(pos.y) - sin(pos.y) * sin(pos.z),
+        cos(pos.y) * cos(pos.z) - sin(pos.z) * sin(pos.x));
+    return gyroidGrad + sphereGrad;
+  };
+
+  auto error = [sphericalGyroid](float* newProp, glm::vec3 pos,
+                                 const float* oldProp) {
+    newProp[0] = glm::abs(sphericalGyroid(pos));
+  };
+
+  Manifold gyroid = Manifold(
+      LevelSet(sphericalGyroid, {glm::vec3(-r - extra), glm::vec3(r + extra)},
+               0.5, 0, true, 0.00001));
+
+  Manifold interpolated = gyroid.Refine(3).SetProperties(1, error);
+
+  Manifold smoothed =
+      gyroid
+          .SetProperties(
+              3,
+              [gradient](float* newProp, glm::vec3 pos, const float* oldProp) {
+                const glm::vec3 normal = -glm::normalize(gradient(pos));
+                for (const int i : {0, 1, 2}) newProp[i] = normal[i];
+              })
+          .SmoothByNormals(0)
+          .RefineToLength(0.1)
+          .SetProperties(1, error);
+
+  MeshGL out = smoothed.GetMeshGL();
+  EXPECT_NEAR(GetMaxProperty(out, 3), 0, 0.016);
+  EXPECT_NEAR(GetMaxProperty(interpolated.GetMeshGL(), 3), 0, 0.068);
+
+#ifdef MANIFOLD_EXPORT
+  if (options.exportModels) {
+    ExportOptions options2;
+    ExportMesh("smoothGyroid.glb", out, options2);
+  }
+#endif
+}