Wrong distance sign computation in close-by geometric configurations
I found what I think is a bug in openVDB 7.2.
I created a minimum working example that you will find at the end of this post.
Bug description
The MWE converts into a vdb grid the meshes of two spheres of radius 1 centered in (-1-epsilon/2,0,0) and (1+epsilon/2,0,0), so that they are separated by exactly epsilon=0.1 along x. Now, when reading the interpolated value of the grid at the origin, the expected result is 0.05, as is lies exactly midway between both spheres.
But if I trace this value against voxel size, here is what i get :

It starts ok, but then around 0.062, it suddenly flips sign before stepping down (at about 0.072) and proceeding to continuously decrease. The expected sampling issue where samples merely miss the gap should not arise before the voxel size gets greater than the gap itself, so 0.1. Below that, I would also expect some variation of the value (based on the relative location of the sample in the gap), but not a wrong sign entirely.
Here is a meshed iso-surface at 0 with a voxel size of 0.07:

As expected from the above, the negative value at the origin translates into this point belonging to the inside of the resulting object, and hence merging both spheres whereas they should be separated.
Minimum Working example
Here is the MWE (code + mesh file): project.zip
Is sphere.obj attached anywhere here? Not being an actual sphere, but a polygonal representation, I suspect tessellation will affect reproducibility.
I think I have reproduced it in Houdini by constructing two spheres and using a vdb from polygons. A Volume Slice set to show a sudden change on sign change will show the origin flipping sign as you adjust the voxelsize. The voxelsize that triggers the flip is 0.13 in the default tesselated sphere, using a denser tesselation seems to cause the sign flip to occur with smaller voxel sizes.
Thank you for your answer !
Is sphere.obj attached anywhere here?
I forgot to attach it, so I just edited my initial comment to include both the code of the MWE and the mesh file. Note that I tried with a much finer mesh (27 MB vs 1.3 MB) and obtained the exact same result, so I think that this issue is independent of resolution.
We discussed this in the last TSC meeting, notes can be found at: https://github.com/AcademySoftwareFoundation/openvdb/pull/925/files
The problem we suspect is with intentional smoothing we apply to the levelset fields. We do not treat the mesh as ground-truth because it is often a discretized representation of something else. And when it comes from an artist, it often isn't manifold, watertight, non-self-intersecting, or any other adjective you'd like to apply. So, in particular, between two spheres such as this the gradient of the level set hits zero, so there are some issues trying to build a valid levelset in any case at this location.
I think there is an RFE to add some flags to suppress the clean-up code for cases when one knows the incoming polygonal mesh is "correct".
I was unaware of this smoothing step, and it does seem to explain the strange behavior I noticed.
Being able to disable this step would then be very convenient for me.