SDF interpolation

[elalish]

This tests out a new use-case for our smooth mesh interpolation: evaluating SDFs on a coarse grid, then smoothly interpolating them to get a cleaner shape with less computation. The SDF grid size becomes about feature extraction (what is the thinnest feature desired, etc), but the smoothness of the surface can be refined arbitrarily after the fact.

In order to not introduce wiggles to the interpolated mesh, the vertices must be very precise, so I've added a precision parameter to LevelSet and an improved version of bisection root-finding. This API is not final, as I'm making a follow-on PR to pull LevelSet into Manifold as a static constructor. It is also important to use precise normals, so I evaluated them analytically, rather than using our pseudonormal calculation.

This isn't quite ready to merge, as I think I found a bug when I was playing around with the new test.

Original SDF result:

With smoothing:

[fire]

Can I use this to create manifold meshes out of non manifold input?

This was done in coacd via sony's vdb and also the same idea was implemented in Unreal Engine's mesh decimation tool.

• https://github.com/SarahWeiii/CoACD/blob/2389bd4d5882b60c6d4a3236e0811056701d9208/src/preprocess.cpp#L7

[elalish]

Yes, an SDF is a great way to guarantee manifoldness. Of course, making "good" output from non manifold input will always be a hard problem since it's so arbitrary. And the SDF will not retain sharp corners. So, as usual, it's all tradeoffs.

[elalish]

@pca006132 this is odd - we're getting Vec out of Range errors on some of the CI bots just for Smooth.SDF, but my Mac is fine - even the address sanitizer shows no problems. Any chance you can repro the problem on your machine?

[pca006132]

Is there any chance this is something related to concurrency? I am a bit busy this week (for some reason I need to do a 40 minutes presentation to a bunch of random peolpe now), so it will take some time.

[pca006132]

Ah, seems something that can be reproduced on any machine. Will check.

[pca006132]

Apparently it is somewhere in GetFaceBoxMorton:

manifold::VecView<glm::vec<3, float, (glm::qualifier)0> >::operator[] (this=0x251481b0a30, i=649655) at /home/pca006132/code/manifold/src/utilities/include/vec_view.h:50
50	   return ptr_[i];
(gdb) bt
#0  manifold::VecView<glm::vec<3, float, (glm::qualifier)0> >::operator[] (this=0x251481b0a30, i=649655) at /home/pca006132/code/manifold/src/utilities/include/vec_view.h:50
#1  manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0::operator()(int) const (this=0x7ffffffefdd8, 
    face=face@entry=757329) at /home/pca006132/code/manifold/src/manifold/src/sort.cpp:291
#2  0x00007ffff7d424cb in manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1}::operator()(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&) const (this=0x7ffff690b550, range=...) at /home/pca006132/code/manifold/src/utilities/include/par.h:377
#3  std::__invoke_impl<void, manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&>(std::__invoke_other, manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&) (__f=..., __args=...) at /nix/store/llmjvk4i2yncv8xqdvs4382wr3kgdmvp-gcc-13.2.0/include/c++/13.2.0/bits/invoke.h:61
#4  std::__invoke<manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&>(manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&) (__fn=..., __args=...) at /nix/store/llmjvk4i2yncv8xqdvs4382wr3kgdmvp-gcc-13.2.0/include/c++/13.2.0/bits/invoke.h:96
#5  std::invoke<manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&>(manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0)::{lambda(tbb::detail::d1::blocked_range<manifold::CountingIterator<int> > const&)#1} const&, tbb::detail::d1::blocked_range<manifold::CountingIterator<int> >&) (__fn=..., __args=...) at /nix/store/llmjvk4i2yncv8xqdvs4382wr3kgdmvp-gcc-13.2.0/include/c++/13.2.0/functional:113

i = 649655, size = 649647

Compiling without parallelization also has this bug:

manifold::VecView<glm::vec<3, float, (glm::qualifier)0> >::operator[] (this=0x598081a0a30, i=649639, i@entry=140737488289008)
    at /home/pca006132/code/manifold/src/utilities/include/vec_view.h:50
50	   return ptr_[i];
(gdb) bt
#0  manifold::VecView<glm::vec<3, float, (glm::qualifier)0> >::operator[] (this=0x598081a0a30, i=649639, i@entry=140737488289008)
    at /home/pca006132/code/manifold/src/utilities/include/vec_view.h:50
#1  manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0::operator()(int) const (face=233400, this=<optimized out>)
    at /home/pca006132/code/manifold/src/manifold/src/sort.cpp:291
#2  std::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0) (
    __first=..., __last=..., __f=...) at /nix/store/llmjvk4i2yncv8xqdvs4382wr3kgdmvp-gcc-13.2.0/include/c++/13.2.0/bits/stl_algo.h:3833
#3  manifold::for_each<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0) (policy=manifold::ExecutionPolicy::Par, first=..., last=..., f=...) at /home/pca006132/code/manifold/src/utilities/include/par.h:382
#4  manifold::for_each_n<manifold::CountingIterator<int>, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0>(manifold::ExecutionPolicy, manifold::CountingIterator<int>, unsigned long, manifold::Manifold::Impl::GetFaceBoxMorton(manifold::Vec<manifold::Box>&, manifold::Vec<unsigned int>&) const::$_0) (policy=manifold::ExecutionPolicy::Par, first=..., n=1299556, f=...) at /home/pca006132/code/manifold/src/utilities/include/par.h:392
#5  manifold::Manifold::Impl::GetFaceBoxMorton (this=this@entry=0x598081a0a10, faceBox=..., faceMorton=...) at /home/pca006132/code/manifold/src/manifold/src/sort.cpp:277
#6  0x00007ffff7d8c1d5 in manifold::Manifold::Impl::Finish (this=0x598081a0a10) at /home/pca006132/code/manifold/src/manifold/src/sort.cpp:137
#7  0x00007ffff7d8846b in manifold::Manifold::Impl::Refine(std::function<int (glm::vec<3, float, (glm::qualifier)0>)>) (this=0x598081a0a10, edgeDivisions=...)
    at /home/pca006132/code/manifold/src/manifold/src/smoothing.cpp:1019
#8  0x00007ffff7d77ed5 in manifold::Manifold::RefineToLength (this=<optimized out>, length=15.5284863) at /home/pca006132/code/manifold/src/manifold/src/manifold.cpp:757
#9  0x00005555555f7384 in Smooth_SDF_Test::TestBody (this=<optimized out>) at /home/pca006132/code/manifold/test/smooth_test.cpp:419
#10 0x00007ffff7f11f77 in testing::internal::HandleSehExceptionsInMethodIfSupported<testing::Test, void> (method=<optimized out>, location=0x7ffff7f208a7 "the test body", 
    object=<optimized out>) at _deps/googletest-src/googletest/src/gtest.cc:2612
#11 testing::internal::HandleExceptionsInMethodIfSupported<testing::Test, void> (object=object@entry=0x598081401a0, method=<optimized out>, location=0x7ffff7f208a7 "the test body")
    at _deps/googletest-src/googletest/src/gtest.cc:2648
#12 0x00007ffff7ef6788 in testing::Test::Run (this=this@entry=0x598081401a0) at _deps/googletest-src/googletest/src/gtest.cc:2687
#13 0x00007ffff7ef7800 in testing::TestInfo::Run (this=0x5980813f600) at _deps/googletest-src/googletest/src/gtest.cc:2836
#14 0x00007ffff7ef8405 in testing::TestSuite::Run (this=0x5980813e840) at _deps/googletest-src/googletest/src/gtest.cc:3015
#15 0x00007ffff7f0865d in testing::internal::UnitTestImpl::RunAllTests (this=this@entry=0x598080e0600) at _deps/googletest-src/googletest/src/gtest.cc:5920
#16 0x00007ffff7f12e87 in testing::internal::HandleSehExceptionsInMethodIfSupported<testing::internal::UnitTestImpl, bool> (method=<optimized out>, 
    location=0x7ffff7f210fe "auxiliary test code (environments or event listeners)", object=<optimized out>) at _deps/googletest-src/googletest/src/gtest.cc:2612
#17 testing::internal::HandleExceptionsInMethodIfSupported<testing::internal::UnitTestImpl, bool> (object=<optimized out>, method=<optimized out>, 
    location=0x7ffff7f210fe "auxiliary test code (environments or event listeners)") at _deps/googletest-src/googletest/src/gtest.cc:2648
#18 0x00007ffff7f081aa in testing::UnitTest::Run (this=0x7ffff7f351a8 <testing::UnitTest::GetInstance()::instance>) at _deps/googletest-src/googletest/src/gtest.cc:5484
#19 0x00005555555a1709 in RUN_ALL_TESTS () at _deps/googletest-src/googletest/include/gtest/gtest.h:2317
#20 main (argc=1, argv=0x7fffffff0a78) at /home/pca006132/code/manifold/test/test_main.cpp:83

(gdb) p i
$1 = 649639
(gdb) p size_
$2 = 649639

[elalish]

@pca006132 Thanks, it looks like I'm getting a few NaN verts out of smoothing, but not on my machine. I must have missed a divide-by-zero somewhere. Since you have a repro, is there a chance you can expand on the last commit I made to hunt down where the first NaN creeps in? It's a bit tedious doing a debug cycle through our CI, but that's what I'll do if you don't have a chance to look at this first.

[pca006132]

I think the simplest way may be to bring up a linux vm? I am outside almost everyday now so I have very limited debugging capacity

[pca006132]

Got some time and did some debugging.

diff --git a/src/manifold/src/smoothing.cpp b/src/manifold/src/smoothing.cpp
index 2fa994cf..280de296 100644
--- a/src/manifold/src/smoothing.cpp
+++ b/src/manifold/src/smoothing.cpp
@@ -769,6 +769,10 @@ void Manifold::Impl::DistributeTangents(const Vec<bool>& fixedHalfedges) {
           glm::vec3 tangent(halfedgeTangent_[current]);
           halfedgeTangent_[current] = glm::vec4(
               glm::rotate(tangent, angle, normal), halfedgeTangent_[current].w);
+          if (!isfinite(halfedgeTangent_[current].x)) {
+            std::cout << "angle = " << angle << ", normal = " << normal << std::endl;
+            __builtin_debugtrap();
+          }
           ++i;
         } while (!fixedHalfedges[current]);
       });

The normal is somehow zero here, which makes the rotation ill-defined and causes the tangent to become NaN.

[pca006132]

Thinking if automatic differentiation can be useful for computing a precise normal automatically, when users are given a set of function combinators...

[elalish]

@pca006132 Thank you for narrowing that down for me! This was an existing bug, just brought to the surface by this new test. After some checking, the one-line fix should be fine. The reason the normal was zero was that the only valid tangents were exactly opposed, in which case there is nothing DistributeTangents needs to do anyway, so the early-out is the simple answer.