getTriangleHitPointInfo: Compute Interpolated Normal #705
Conversation
This improves the insidedness heuristic for the ClosestPointToPoint function...
|
Using interpolated normals is much more of a situational solution and feels bit like a hack to to fix this case. If vertex normals aren't provided then you're back with the same issue. Likewise vertex normals shouldn't be necessary for determining containment and there's no guarantee that the vertex normals are coherent with the facing direction of the triangle which could cause more issues. For context these incorrect distances are caused by this extraneous geometry with averaged normals by the rabbits foot: 
This does seem to be an issue caused by what you've laid out in #704. I would think there would be some way to choose a "best" triangle from a set if there multiple triangles that are at the "closest point". We should create a simple repro cases (not just the rabbit model) that can be used to reproduce this issue. A simple pyramid like the one in the rabbit model may work just fine. |
|
A pyramid should work for this case; I can make that. I don’t think there’s a “best triangle” heuristic that works unless you already know whether you’re inside or outside the mesh 😅 I think you need the local information from all the triangles on the vertex/edge to properly make an inside/outside judgement. This non-local(?) information is why interpolating the normals out to the edges helps so much. But, even then, as you say, it’s easy to fool area weighted normals. If I’m imagining correctly, a trapezoid with a small top will have funky normals at the top. I can add a test case for that too. I wrote a more principled solution with Generalized Winding Numbers on the locally closest triangles master...zalo:three-mesh-bvh:master but it still didn’t seem to help Suzanne… which, you’ve now told me isn’t even manifold anyway 🫠 (so maybe it does work…) |
Looking at the image you posted in the other issue I can intuit that the left edge is the one that we would want to be considered "closer" even though the actual closest point is shared by both edges. The question is how codify that. Using something like the face normals and the delta edge from point to closest point to compare the triangles might be the right track - ie taking the triangle that has a normal closest to parallel. I don't think this should have to be a heuristic. 
Signed winding numbers to compute volume would work but you risk numeric instability when the point being checked is really close to a geometry edge. From #704 I mention that the ray-crossing algorithm or raycast to check if a back face is hit would work, as well, which may be more robust though still subject to some numeric issues near edges. |
|
Hrm, taking the normal closest to parallel/antiparallel with the (normalized) delta might work; will have to test it. |
|
Hrm, it seems like neither the most orthogonal triangle or my "Generalized Winding Number" solution fix the Rabbit SDF's invalid zone: https://github.com/zalo/three-mesh-bvh/commit/65a35c47ed984b855263f976b6fbfb5bea8f3da9 Either I have bugs, or I feel like there's something funny I'm misunderstanding about the way that shapecast is interacting with the variables outside of its scope... but the way If the closest point is on an edge or a vertex, is shapecast guaranteed to include all of the adjacent triangles? Or would I need to build a half-edge structure to traverse all of the adjacent triangles? Another heuristic I'm curious about is clamping the baryCoords to 0.99 when finding the closest point on the triangle (which could act as a tie-breaker between adjacent triangles...) EDIT: Ohohoho, could this option in the bunny's generator have something to do with it? EDIT2: Came up with a torture test mesh
import trimesh
import manifold3d
square : manifold3d.CrossSection = manifold3d.CrossSection.square([0.5, 0.5], True)
square -= manifold3d.CrossSection.circle(0.5, 100).translate([0.24, 0.24])
extrusion = square.extrude(0.7, scale_top=(0.0, 0.0))
mesh = extrusion.rotate([-90, 0.0, 0.0]).to_mesh()
meshOut = trimesh.Trimesh(vertices=mesh.vert_properties, faces=mesh.tri_verts)
trimesh.exchange.export.export_mesh(meshOut, 'torture_tester.obj', 'obj')v -0.25000000 0.00000000 0.25000000
v -0.14525677 0.00000000 0.07871194
v -0.12448426 0.00000000 0.10227354
v -0.10227363 0.00000000 0.12448426
v -0.07871203 0.00000000 0.14525659
v -0.21241347 0.00000000 -0.02711038
v -0.18216400 0.00000000 0.02791335
v -0.16450851 0.00000000 0.05389263
v -0.19815339 0.00000000 0.00087675
v -0.23552828 0.00000000 -0.08549149
v -0.22488834 0.00000000 -0.05593780
v -0.24429159 0.00000000 -0.11565510
v -0.25000000 0.00000000 -0.14119309
v -0.05389263 0.00000000 0.16450851
v -0.02791359 0.00000000 0.18216391
v 0.02711025 0.00000000 0.21241347
v -0.00087675 0.00000000 0.19815333
v 0.05593768 0.00000000 0.22488822
v 0.08549149 0.00000000 0.23552828
v 0.11565484 0.00000000 0.24429156
v 0.14119309 0.00000000 0.25000000
v 0.00000000 0.69999999 0.00000000
f 1 3 4
f 11 6 1
f 1 9 7
f 1 6 9
f 1 7 8
f 1 8 2
f 1 2 3
f 1 4 5
f 1 5 14
f 1 14 15
f 1 15 17
f 1 16 18
f 1 17 16
f 22 3 2
f 22 5 4
f 22 4 3
f 11 1 10
f 10 1 12
f 12 1 13
f 22 13 1
f 22 11 10
f 22 10 12
f 22 12 13
f 22 8 7
f 22 7 9
f 22 9 6
f 22 6 11
f 22 2 8
f 1 18 19
f 1 19 20
f 1 20 21
f 22 1 21
f 22 15 14
f 22 17 15
f 22 14 5
f 22 16 17
f 22 18 16
f 22 19 18
f 22 20 19
f 22 21 20The default behavior and my solid angle "technique" (super bad!) With most-orthogonal triangle With interpolated normals |
|
Ah, when you mentioned the extraneous geometry on the rabbit's foot, I didn't realize they were disconnected, floating, non-manifold triangles... 😅
I have a heuristic that works for the torture test (pick the triangle with a normal furthest along the delta, no absolute value), but it fails badly for the extra features I added to it later… still thinking I need a heuristic that uses the info from all the adjacent triangles somehow…I’ll keep trying things as I have time… |
Shapecast will only iterate over triangles that are within bounds that determined to be intersecting a shape as returned by "intersectsBounds". The current implementation of "closestPointToPoint" doesn't use an "equal to" boolean check with the current best point so it's possible a bounds with an equivalent closest point are skipped. Floating point error may also play a role. Checking the approach by checking every triangle by brute force, first, to ensure it works before using the optimizations from shapecast would be best. You can do this by always returning "true" from "intersectsBounds" which will cause every triangle to be checked.
No this just impacts how deep the BVH is generated to.
I would probably use a more basic pyramid with 4 faces to see if it can be reproduced that way. It would make things easier to debug.
I'd be interested to see the code for this and be able to look at the two or three triangles that meet at a point causing this the issue.
I'm not sure it is non-manifold. The image I shared above shows that there are multiple triangles at that spot. I haven't looked at the topology in detail but it's possible it's a small pyramid shape. |
Here's the code I was using; I might be able to test it more tonight, but does this look reasonable? |
|
Thanks. The logic looks right but it's not clear what other degenerate cases there might be here. I've repro'd it using a much simpler cone geometry with 3 sides with cpu generation:
Looking at the generated SDF layers you can see where we're getting incorrect results: 
By performing a ray cast after collecting the closest distance you can determine whether it's inside the model and the results look good:
It just requires adding these lines after performing the distance check: const ray = new THREE.Ray();
ray.origin.copy( point );
ray.direction.set( 0, 0, 1 );
const hit = bvh.raycastFirst( ray, DoubleSide );
sdfTex.image.data[ index ] = hit && hit.face.normal.dot( ray.direction ) > 0.0 ? - dist : dist;Of course this incurs more overhead, though. I would be interested to see what is happening at one of those incorrect points to try to see if it's another detectable case. This raycasting method could be used as a ground truth to find a point that's producing an incorrect result if you'd like try to identify a failure case and visualize the transformations of the triangles and point. |
That's a good idea; the fact that the output field is so noisy suggests there is a conceptual failure I'm not understanding yet... I'd bet that the adjacent triangles aren't always getting sampled for some reason... If that turns out to be the case, I don't suppose you'd let me make a half-edge data structure in the BVH construction step? I can also add the adjacent triangles to the ExtendedTriangle representation to facilitate this enumeration🤔 |
In the above images I have ensured that every triangle is always sampled by returning |
|
Bah I've tackled this a couple times since, and haven't been able to make any progress. The one "good" looking estimator I thought I had ended up always returning positive distances. Perhaps it's worth it to bite the bullet and just use |
I've toyed around with this a bit more, as well, and while I can get this working well in cases where an edge is shared there are still issues in some cases where there's a common vertex shared by 3 or more faces. I suspect this is a more complicated case that it not handled well by testing the normal. I'm seeing cases where floating point precision is causing issues close to triangle edges and surfaces, as well. Overall I think the raycast method is going to be significantly more robust. If you'd like to make a PR into the sdfGeneration example that updates it to use the raycast method for detecting whether a point is inside or outside the mesh in both the CPU and GPU generation case that would be great.
The GPU variant will be using the same approach as the CPU version currently which means it will also be incorrectly using the normal for containment detection.
If you're interested in seeing if the GPU variation can be pushed even further in terms of speed by converting the binary BVH to a compressed one, which can apparently help improve performance on the GPU (#519). Likewise running this logic in a compute shader vs a rendering shader would likely provide improvements, as well, though this requires WebGPU. This would both be great improvements to the project if you'd like to help with either. |
|
Alrighty; bullet time. |
|
I came across a paper discussing the 2d-variation of the issue we ran into when measuring distance to vertices. They also use orthogonality but it may be discussed more rigorously. I'll note it here for future reference. It's mentioned in section 2.4 of the paper: https://github.com/Chlumsky/msdfgen/files/3050967/thesis.pdf |
|
This version of the heuristic might be viable for it too… from the GeometryProcessingWorld discord server |
|
I would think using winding numbers for a complex geometry will always be slower than just performing a raycast to determine inside-ness. That might be a good discord to ask about 3d equivalents of the 2d approach mentioned in #705 (comment) and #705 (comment) so we can reuse the calculations used in computing closest point values. |
@gkjohnson This adds normal computation to
getTriangleHitPointInfoin the spirit of my sister PR to three.js: mrdoob/three.js#25566Using an interpolated normal significantly improves the insidedness-via-normal heuristic for the ClosestPointToPoint function, especially when the closest point lies along an edge or vertex (and is generally just more correct).
This is starting to address #704 with relatively little extra cost.
Before:
After:
Also a comparison shot from my pet project.
Before:
After:
