gkjohnson / three-mesh-bvh

A BVH implementation to speed up raycasting and enable spatial queries against three.js meshes.
https://gkjohnson.github.io/three-mesh-bvh/example/bundle/raycast.html
MIT License
2.56k stars 268 forks source link

Potential bug raycasting on mesh vertices. #169

Closed obsius closed 3 years ago

obsius commented 3 years ago

Hi again, I'm getting an issue when a raycast exactly lines up with one of the input mesh vertices. The raycaster returns nothing when it should be returning an intersection.

In the below image, it just so happens that a height map grid I'm constructing aligns with a row of those cylinders (each cylinder has a value that I am using as a Z for that XY), and you can see the contour lines think that row has no height.

image

If on a raycast failure I try again and offset my ray origin point by .01, I get the following contour lines:

image2

The relevant code is:

let geometry = new THREE.BufferGeometry();
let positions = new Float32Array(refs.length * 3); //data is attached

geometry.setAttribute('position', new THREE.BufferAttribute(positions, 3));

geometry.boundsTree = new THREE.Bvh(geometry, {
    strategy: THREE.Bvh.SPLIT_TYPES.AVERAGE,
    maxDepth: 64,
    maxLeafTris: 16
});

this.mesh = new THREE.Mesh(geometry, {
    side: THREE.DoubleSide
});

this.raycaster = new THREE.Raycaster();
this.raycaster.firstHitOnly = true;

// this seems to be one of the bad raycasts
this.raycaster.ray.origin.x = -49;
this.raycaster.ray.origin.y = -4.139999999999418;
this.raycaster.ray.origin.z = -1000;

// raycast straight up
this.raycaster.ray.direction.x = 0;
this.raycaster.ray.direction.y = 0;
this.raycaster.ray.direction.z = 1;

this.raycaster.intersectObject(this.mesh);

Data (Float32Array serialized as an object, each key is the index):

data.json.txt

Thanks for making this library, I use it a lot for this project. I appreciate the time you've put into this.

gkjohnson commented 3 years ago

Great thanks for the report and the repro! I'm able to see the problem and it looks like the ray you provided is sitting exactly on the edge boundary of the AABB of one of the nodes. Artificially inflating the boxes by a small amount seems to fix the issue but I'll have to look into a bit more when I have time over the next couple weeks to see exactly what the issue is.

you can see the contour lines think that row has no height.

Just to clarify I'm seeing that raycaster.intersectObject( mesh ) is returning no results rather than a result with "no height" (which I interpret to be a hit where z === 0) is that correct?

Thank you again!

Notes - Some triangles are not completely encapsulated by the leaf bounding box - All child bounds _are_ encapsulated by parent bounding box - Splits lose precision when using float32 array for center and half extents. **Solutions** - Adjust bounds array to store just the min and max bounds: - `computeSAHPlanes` - `getCentroidBounds` (reads center point directly) - `getBounds` (reads center point directly) - `getAverage` uses the center points directly - Adjust the bounds array to store the min max and center of the bounds.
obsius commented 3 years ago

Correct, no results. As a fallback I set any grid cell in my height map that has no ray-casted intersection to zero in order to generate those contour lines.

Thanks for the quick response.

gkjohnson commented 3 years ago

Took a moment to look at where the triangle bounds are generated and it looks like we're losing precision because we store the the bounds center and half extents in a Float32Array which in some cases is changing the stored bounding center value by a magnitude up to 0.00001525.... Storing the min and max values rather than the half extents would prevent the precision loss because those values are directly from the position Float32 position buffer.

Using a Float64Array also fixes it but means more memory overheard on creation so I'm inclined to store the min / max values instead but I'll have to check what other functions that effects.

obsius commented 3 years ago

I'd go with the extents increase assuming that this is only being missed by the indexing and not the final calculation of the intersection.

gkjohnson commented 3 years ago

@obsius

I'd go with the extents increase assuming that this is only being missed by the indexing and not the final calculation of the intersection.

I think I'm now leaning towards this approach. It doesn't affect the memory and performance is minimally if at all effected.

And I've confirmed it's an error in the indexing -- I've made PR #170 that fixes the issue by adding an epsilon that's scaled by the half extents and center of the triangle bounds being stored. I've added a helper that can be used to make sure that all triangles are contained within the bounds after generation, as well. I'm not an expert in floating point number representation (I don't suppose you are? 😁) but the approach I've used works for the cases I was testing (5000 radius sphere sphere, your model). Take a look at the PR (specifically where the FLOAT32_EPSILON is applied) and let me know what you think.

I might do a little reading on floating point error estimation at different number scales at some point to understand this more.

obsius commented 3 years ago

I wish, haha, but I too have no idea how to handle the floating point arithmetic issues. I guess the drift depends on the size of the number. Here is a good website that shows the floating point error based on the original number: https://www.h-schmidt.net/FloatConverter/IEEE754.html

For example, 20.244 has an error of -5.18798828125E-7 but 2,000,000.244 has an error of 0.006

obsius commented 3 years ago

Oh, and the fix looks good, I deal with huge numbers since the mines that provide the data often have XYs at a million or more, but because of the issue I had last year I now offset those to a local origin so I never have coords greater than a thousand or so.

triangleBounds[ tri6 + el2 + 1 ] = halfExtents + ( Math.abs( min ) + halfExtents ) * FLOAT32_EPSILON;

Always nice to have a one liner.

gkjohnson commented 3 years ago

I think I'm happy with what's there for now. I don't love the change because we could have exact bounds if we put the triangle min / max directly in the buffer but that would come with a small performance or memory hit. I'll make another issue to track solutions for this in case another option comes up but I think this is good enough for now and is otherwise a pretty minimal change.

but because of the issue I had last year I now offset those to a local origin so I never have coords greater than a thousand or so

I was thinking about that -- it's possible that this would fix #112, as well.

Thanks again for reporting the issue! I'll have a new version published in the next couple weeks. I might try to get a few other changes in beforehand.

And a bit of an aside I don't know your exact use case but if you're trying to make contour or topographic lines I made a shader here awhile ago and recently updated it that renders crisp lines on the surface of a mesh. If you're dealing with meshes that have caves or other concave features like that a shader based approach may help handle those cases more easily.

obsius commented 3 years ago

Sounds good. Not relevant in this particular use case, but I do show topos and I may integrate that at some point. For now I use a shader to color topos by elevation:

test

In the case that caused this bug, the raycaster is actually built on a fake topography that is only used to populate a grid. This grid feeds a marching squares algorithm which produces those contour lines, but there isn't really an elevation when displayed since the Z value fed into the raycaster is actually each blast hole's detonation time. The coloring you saw in my first post is showing the timing delays between hole detonations.

gkjohnson commented 3 years ago

@obsius Just published the fix in v0.3.2:

https://github.com/gkjohnson/three-mesh-bvh/releases/tag/v0.3.2

obsius commented 3 years ago

Great, thanks for the heads up. I will update right away.