Caching, performance improvements, and three-mesh-bvh?

[gkjohnson]

Hey! Love seeing work happen in this space and this looks great -- I think a good performant, memory efficient CSG library would be a huge win for the three.js community.

As you've seen in https://github.com/gkjohnson/three-mesh-bvh/issues/281 - I've thought a bit about how I'd implement CSG using three-mesh-bvh but there are too many projects I want to work on and it's looking more and more like I won't ever get to it 😂 Depending on how much time someone is interested in putting time into this and building out a more robust CSG solution I'm happy to share some of my thoughts, ideas, and research on how I was imagining it working. Some of the newer CSG solutions out there now actually don't use BSPs. And I'd be happy to prioritize some changes to the three-mesh-bvh project to enable this work more cleanly if that's going to be worked on.

In terms of additional features I was thinking about - maintaining a CSG hierarchy of operations, caching intermediate results and data structures to help improve memory allocation and run time performance when regenerating a mesh, and longer term possibly enabling it to be run in a Web Worker (or multiple) would all be greatt. I haven't taken a deep look at the code here, yet, though so I'm not sure if any of these are already here!

And I'm generally interested in seeing three-mesh-bvh used under the hood because then with one spatial data structure representation you could run CSG, path traced rendering, physics, performant raycasts, geometric intersections, etc all with a memory-compact hierarchy representation backing it!

cc @manthrax

[giladdarshan]

Hey @gkjohnson, agreed!

I thought about using BVH as the data structure for CSG but when running it in my head I could see it running into similar issues that we see in BSP (such as recursion overflow), it will be interesting to compare BSP vs BVH vs Octree for CSG.

I didn't fully test it but I think there is also a drawback when using only raytracing for triangle intersection for consecutive CSG operations, for example, raytracing can see a triangle as not intersecting when checking intersection with a hollow cube (subtract a 8mm cube from a 10mm cube), if you check with triangle center and the center is outside the object or inside the hollow area it will not show the triangle as intersecting, same goes if you check with the triangle points.

The main logic behind this library is:

1. Mark polygons (triangles) that are intersecting with the other Octree
2. Loop over the intersecting polygons and handle the intersection and for each polygon: a. Store polygons from the other Octree that are passing the triangle-triangle intersection test with the current polygon b. If none found, use raytracing to check intersection type (inside, outside, spanning), if spanning it will check intersection type against all planes in the other Octree since the regular triangle-triangle intersection test failed c. Loop over all the stored polygons from 2.a and use their planes for the split by plane function (taken from @manthrax's THREE-CSGMesh which is based on csg.js)
3. Delete polygons based on the CSG operation and states (types)
4. Combine the two Octrees into a new Octree based on the CSG operation

Additional features - all on the (still unwritten) todo list :) Web workers in general is something I need to be more familiar with. Regarding the caching intermediate results, how would you approach it? If I'm not mistaken, raytracing will require you to rebuild the geometry between CSG operations so it will work properly, no?

Adding CSG to three-mesh-bvh to have it as a all-in-one package I think will have performance impact on the other awesome features, please correct me if I'm wrong but currently you are not storing the triangles inside a BVH tree but only the information about them (offset, count, split axis, etc) as a buffer and re-arranging the index, adding CSG functionality would require you to start storing the triangles which can affect memory consumption and performance.

I will add a test option to the demo @manthrax prepared to see what kind of impact adding BVH to the mix will have (performance impact of adding BVH tree vs. raytracing improvement thanks to three-mesh-bvh).

There is still a lot of performance improvement to do in this library, my main goal was to get a POC working first and then improve the performance (for example, right now all the CSG operations are iterating over the polygons multiple times during the deletion for each state, I plan to reduce it to a single iteration).

Fun times ahead :)

[gkjohnson]

To be clear I don't expect raycasting to play a big role in the CSG generation. A BVH is just a spatial acceleration structure like an octree that's organized a bit differently. And three-mesh-bvh provides a number of spatial queries in addition to raycasting including geometry / geometry, box / geometry, sphere / geometry, etc intersections.

To ground the discussion a bit on of the use cases that originally prompted some interest was IFC.js which currently uses three-mesh-bvh to partition their geometry to accelerate mouse ray interaction, voxelize the building models, clip the geometry, and generate floor plans. Interactive CSG was one use case they mentioned for modifying or extending their models. That might be a ways out but it's a good use case, I think, and targeting interactive frame rates makes everything more fun.

I'm actually seeing now that your approach does not use BSPs so maybe it would be useful for you to walk through the approach used in this project so I can understand at a higher level. I'll walk through what I was imagining, though, and there may be a lot of overlap. Of course this assumes that all geometry in play is two-manifold otherwise there will likely be issues:

1. Collect all triangles that intersect between geometry A and B.
2. Split all intersecting triangles on the intersection planes with the other triangles. The geometries are now split into disconnected pieces along the intersecting edges.
3. For each piece in the new split sets of triangles use something like the ray-crossing algorithm determine whether or not a vertex is inside or outside the other mesh and include or discard it depending on whether union, difference, or intersection are being used.
4. Construct a new geometry from the triangles.
5. Repeat!

This should support the nested box example you mention. Originally inspired by this use case I put together this demo showing the edges resulting from the intersection of two torus geometries which hopefully helps illustrate the separate "chunks" of geometry formed by the intersections.

To find these "chunks" of geometry you can use a half edge data structure to cache sibling triangle connections and find the coherent set of triangles to group. I'm imagining that this can be stored in a vertex buffer attribute to save memory, as well.

And lastly, regarding caching and the operation hierarchy, CSG is often performed in a hierarchical fashion so multiple groups of operations can be reused as a larger group of operations.

Here you could optimize things by caching precomputed geometry, BVHs (or octrees), and half edge structures at every step to help accelerate interactive adjustment of subtrees and only recompute the objects necessary. You can also save a lot of memory allocation by caching data and reusing what you have. And three-mesh-bvh is very compact being stored in tight array buffers rather than JS objects.

Of course the downsides to three-mesh-bvh is that it's a fairly rigid structure - it can't be modified too extensively in terms of adding removing new triangles from it so you'd need to regenerate a new BVH for every new step in the hierarchy which may be an issue at some point? But again this can be done in web workers off the main thread for an interactive delay. There are definitely some trades to be made. Memory management for stuff like this is difficult in JS since the only real way to be compact is to use fixed length array buffers so there may be no getting around that rigidity if memory footprint is considered.

[giladdarshan]

To be clear I don't expect raycasting to play a big role in the CSG generation. A BVH is just a spatial acceleration structure like an octree that's organized a bit differently. And three-mesh-bvh provides a number of spatial queries in addition to raycasting including geometry / geometry, box / geometry, sphere / geometry, etc intersections.

Understood, I took another look at the three-mesh-bvh file and the demo you referenced and it seems I missed SeparatingAxisTriangle.prototype.intersectsTriangle, I added SeparatingAxisTriangle to the list of exports and tested it quickly with a couple of (coplanar) triangles that caused me issues and the intersection algorithm fails the intersection test -

let bvhTriangle1 = new BVH.SeparatingAxisTriangle(new THREE.Vector3(1, 5, -1), new THREE.Vector3(5, 5, -5), new THREE.Vector3(1, 5, -5));
let bvhTriangle2 = new BVH.SeparatingAxisTriangle(new THREE.Vector3(1, 5, -5), new THREE.Vector3(1, 5, 5), new THREE.Vector3(11, 5, -5));
console.log("BVH Intersect", bvhTriangle1.intersectsTriangle(bvhTriangle2)); 
// Expected: true
// Actual: false

I based the triangle intersection algorithm I'm using on this library which also has issues with coplanar triangles so I added another check of line intersections. I'm sure it's still not bulletproof so it is on the todo list to find a triangle intersection algorithm that performs better.

To ground the discussion a bit on of the use cases that originally prompted some interest was IFC.js which currently uses three-mesh-bvh to partition their geometry to accelerate mouse ray interaction, voxelize the building models, clip the geometry, and generate floor plans. Interactive CSG was one use case they mentioned for modifying or extending their models. That might be a ways out but it's a good use case, I think, and targeting interactive frame rates makes everything more fun.

What do you mean when you say "Interactive CSG"? CSG operations in real time and not while building the scene?

I'm actually seeing now that your approach does not use BSPs so maybe it would be useful for you to walk through the approach used in this project so I can understand at a higher level. I'll walk through what I was imagining, though, and there may be a lot of overlap.

I tried going for Octree-embedded BSPs at first but didn't like the results so I started working on an Octree only library. The approach in high level is pretty much the same, I'm working on a flowchart that can help.

1. Create Octrees from the two meshes
2. Mark intersecting triangles
3. Handle intersecting triangles a. Collect all intersecting triangles and split by plane if needed b. If there are no intersecting triangles, use raycasting to determine if the triangle is inside or outside the other octree / geometry
4. Discard triangles based on criteria for the CSG operation
5. Add all relevant triangles to a new Octree based on the CSG operation
6. Create a new mesh / geometry
   

Of course this assumes that all geometry in play is two-manifold otherwise there will likely be issues

What are your thoughts about T-junctions? Currently my approach is creating T-junctions in order to keep the number of triangles low.

And lastly, regarding caching and the operation hierarchy, CSG is often performed in a hierarchical fashion so multiple groups of operations can be reused as a larger group of operations. Here you could optimize things by caching precomputed geometry, BVHs (or octrees), and half edge structures at every step to help accelerate interactive adjustment of subtrees and only recompute the objects necessary.

It sure is something to look into, the logic itself shouldn't be too complicated, currently my only limitation is that I will need to build the mesh after each step for the raycasting to work. I plan to work on the rayIntersect function which is currently commented out and not used, if I can get it to work properly I believe it will remove the requirement to build a new mesh after each step. This is what I have in mind:

// 1. build the base 5 meshes (cube, sphere and 3 cylinders)
// 2. position them in their location
// 3. do the CSG operations
let result = OctreeCSG.operation({
    op: "subtract",
    objA: {
        op: "intersect",
        objA: cube,
        objB: sphere
    },
    objB: {
        op: "union",
        objA: {
            op: "union",
            objA: cylinder1,
            objB: cylinder2,
        },
        objB: cylinder3
    }
});

You can also save a lot of memory allocation by caching data and reusing what you have. And three-mesh-bvh is very compact being stored in tight array buffers rather than JS objects.

I've been meaning to ask you why you chose to use buffers and not JS objects, is it only memory allocation or also performance?

[gkjohnson]

quickly with a couple of (coplanar) triangles that caused me issues and the intersection algorithm fails the intersection test

I'm not seeing the issue you're mentioning. I'm seeing that the function returns "true" and logs an warning, though. The function does not return a proper result for coplanar triangles so that is something that could be added. That aside, though, the bvh "shapecast" function used for that demo is flexible enough to let you use any triangle / triangle intersection function you want during traversal.

What do you mean when you say "Interactive CSG"? CSG operations in real time and not while building the scene?

I just mean some approach to CSG that enables more interactive framerates so a user can move something around and edit the CSG. At the least this means keeping performance good for CSG operations but it could also mean caching data where possible so it doesn't need to be generated, using a generator / web worker to run the CSG operations asynchronously, and / or debouncing CSG updates until the user stops moving. Some of these things might involve consideration to core code while others (like debouncing) are better considered at the application level, I'd think.

One of the most compelling use cases for CSG, I'd imagine, are things like editors for blocking out level designs for games, etc so I think the interactive nature of the libray is good to consider.

What are your thoughts about T-junctions? Currently my approach is creating T-junctions in order to keep the number of triangles low.

I think T-junctions are fine -- otherwise topology could get crazy. Lower triangle count for memory and improved performance / lower complexity on subsequent CSG operations seems more important. A model can be post processed to have more complete connectivity if needed. This is something I'd wait for someone to complain about before changing, too 😁

if I can get it to work properly I believe it will remove the requirement to build a new mesh after each step.

It's definitely possible to do this all in one step as you've suggested but that doesn't lend itself to interactive interfaces. I'm imagining each level of the hierarchy in the above pictures being transformable and updatable. The fastest way to regenerate the mesh when only a small piece has moved is to have as many other static operation results cached so they don't have to be rerun.

You can see in this video on the Unity plugin "RealTimeCSG" that the user can drag around a pre-built window component and have it update: https://youtu.be/bChuOnELDvQ?t=136. Not sure if enabling this kind of thing is outside the scope of what you were imagining but it's where mind went. And I'm sure there are more tricks going on in enabling that to work with more complex cases.

I've been meaning to ask you why you chose to use buffers and not JS objects, is it only memory allocation or also performance?

The BVH project started off with just a hierarchy of JS objects and memory usage was astronomical. Sometimes nearing (or over) a GB for a single reasonably complex model. JS numbers are all 64 bit and JS objects are effectively hashmaps and require string information to store keys, pointers for arrays / objects, etc. It's all a lot of hidden cost that adds up very quickly when building out a recursive structure. Typed arrays are really the only way to have complete control over memory allocation in JS outside of something like WASM. Once the BVH to switched to use that data started being stored in bitmasked fields when possible the memory allocation was significantly improved and now it's a fairly tight and optimal representation. Array buffers also are significantly faster to transfer between web workers and can be easily serialized and reloaded into the app so that was a big win, as well. I'm not sure if there are any intrinsic performance benefits to them over using basic arrays and objects but there could be. I just haven't tested that.

[giladdarshan]

I'm not seeing the issue you're mentioning. I'm seeing that the function returns "true" and logs an warning, though. The function does not return a proper result for coplanar triangles so that is something that could be added.

My bad, I tested with an older version of the library, will do some tests with the new version and see how it goes.

That aside, though, the bvh "shapecast" function used for that demo is flexible enough to let you use any triangle / triangle intersection function you want during traversal.

I'm still not 100% following the logic behind bvhcast and shapecast, I will get there :)

One of the most compelling use cases for CSG, I'd imagine, are things like editors for blocking out level designs for games, etc so I think the interactive nature of the libray is good to consider.

Sure does sound interesting, will also add this one to the todo list.

I think T-junctions are fine -- otherwise topology could get crazy. Lower triangle count for memory and improved performance / lower complexity on subsequent CSG operations seems more important. A model can be post processed to have more complete connectivity if needed. This is something I'd wait for someone to complain about before changing, too 😁

Agreed! There is some interesting logic in the SimplifyModifier class that I plan to look into in the future, should help minimize the number of triangles even more.

It's definitely possible to do this all in one step as you've suggested but that doesn't lend itself to interactive interfaces. I'm imagining each level of the hierarchy in the above pictures being transformable and updatable. The fastest way to regenerate the mesh when only a small piece has moved is to have as many other static operation results cached so they don't have to be rerun.

This basically means keeping "history" of the objects (octrees) so when one piece changes there is no need to recalculate the rest of the pieces, correct?

The BVH project started off with just a hierarchy of JS objects and memory usage was astronomical. Sometimes nearing (or over) a GB for a single reasonably complex model. JS numbers are all 64 bit and JS objects are effectively hashmaps and require string information to store keys, pointers for arrays / objects, etc. It's all a lot of hidden cost that adds up very quickly when building out a recursive structure. Typed arrays are really the only way to have complete control over memory allocation in JS outside of something like WASM. Once the BVH to switched to use that data started being stored in bitmasked fields when possible the memory allocation was significantly improved and now it's a fairly tight and optimal representation. Array buffers also are significantly faster to transfer between web workers and can be easily serialized and reloaded into the app so that was a big win, as well. I'm not sure if there are any intrinsic performance benefits to them over using basic arrays and objects but there could be. I just haven't tested that.

Thank you for the explanation! Looks like I need to start looking into Array Buffers / Typed Arrays sooner than later.

3. For each piece in the new split sets of triangles use something like the ray-crossing algorithm determine whether or not a vertex is inside or outside the other mesh and include or discard it depending on whether union, difference, or intersection are being used.

When you say ray-crossing algorithm do you mean raycasting? or is it something different? If it's something that can help the CSG operations I will start looking into it. How will you classify triangles that were split and have 1 or 2 vertices "touching" the other object?

[gkjohnson]

I'm still not 100% following the logic behind bvhcast and shapecast, I will get there :)

Shapecast is for arbitrary geometrically-defined shapes while "bvhcast" is for running fast traversals between two geometries with BVHs. The latter is likely closer to what's needed for this work.

This basically means keeping "history" of the objects (octrees) so when one piece changes there is no need to recalculate the rest of the pieces, correct?

Kind of -- not a history like undo but definitely a cache of data for each intermediate object throughout the hierarchy that are required to regenerate the next step if something changes so that data (octree, etc) doesn't need to be regenerated. Of course this is where memory footprint becomes important, as well.

When you say ray-crossing algorithm do you mean raycasting? or is it something different?

There's a wiki article on the raycrossing algorithm here. It's basically just a way to determine whether a point is within a closed shape by counting the number of faces a ray hits in a direction. With well-formed geometry it should be enough to just to check the first face hits normal direction, though (which may be what you're doing already).

How will you classify triangles that were split and have 1 or 2 vertices "touching" the other object?

Can you elaborate on what you mean by this? As in two vertices of object B lying on the surface of object A? It seems like a corner case where you'd just not consider the triangles intersecting because there'd be no clipping result.

[giladdarshan]

There's a wiki article on the raycrossing algorithm here. It's basically just a way to determine whether a point is within a closed shape by counting the number of faces a ray hits in a direction. With well-formed geometry it should be enough to just to check the first face hits normal direction, though (which may be what you're doing already).

Right, currently for triangles (from object A) that do not intersect with triangles from object B I do raycasting test of the triangle's vertices against object B and then check if the first face hits the normal direction. If all 3 vertices are intersecting it classifies the triangle as inside, if none are intersecting it will be outside, if less than 3 but not 0 it will attempt to do classification by plane intersection with all triangles in object B (since it failed the triangle intersection test, it cannot filter intersecting triangles).

let insideCount = 0;
for (let i = 0; i < currentPolygon.vertices.length; i++) {
    let vertex = currentPolygon.vertices[i];
    let point = vertex.pos;
    let direction = vertex.normal.clone().normalize();

    _raycaster1.set(point, direction);
    let intersects = _raycaster1.intersectObject(targetOctree.mesh);
    if (intersects.length) {
        if (direction.dot(intersects[0].face.normal) > 0) {
            insideCount++;
        }
    }

}
if (insideCount >= currentPolygon.vertices.length) {
    currentPolygon.setState("inside");
}
else if (insideCount > 0) {
    currentPolygon.checkAll = true;
    polygonStack.push(currentPolygon);
}
else {
    currentPolygon.setState("outside");
}

Do you recommend using the raycrossing algorithm which is mentioned in the wiki article instead? Unless I misunderstood, it should be something simple like this?

let raycaster = new THREE.Raycaster(point, direction);
let intersects = raycaster.intersectObject(mesh1);
console.log("Intersects Object?", (intersects.length % 2) == 0 ? "outside" : "inside");

Can you elaborate on what you mean by this? As in two vertices of object B lying on the surface of object A? It seems like a corner case where you'd just not consider the triangles intersecting because there'd be no clipping result.

Yes, one or two vertices of a triangle from object A lying on the surface of object B, it is rather common since these triangles are usually the result of splitting intersecting triangles. For example, view from the top of the red triangle from object A is intersecting with (triangles from) object B (yellowish box/cube):

After splitting the red triangle with object B's intersecting triangle's plane, we get 3 triangles (A, B & C):

Triangle A should be classified as inside and has 1 vertex inside object B and 2 vertices on the surface / boundary. Triangle B should be classified as outside and has 1 vertex outside of object B and 2 vertices on the surface / boundary. Triangle C should be classified as outside and has 2 vertex outside of object B and 1 vertices on the surface / boundary.

The issue I'm hitting with raycasting (which can be very likely due to my lack of knowledge / experience with raycasting) is that depends on the direction, results may vary. Lets take the following object, point and direction for example:

let cubeGeometry = new THREE.BoxGeometry(10, 10, 10);
let cubeMaterial = new THREE.MeshBasicMaterial({color: 0x0000ff});
let mesh1 = new THREE.Mesh(cubeGeometry, cubeMaterial);
// mesh1.geometry.computeBoundsTree();
let point1 = new THREE.Vector3(0, 0, 5);
let point2 = new THREE.Vector3(0, 0, -5);
let direction1 = new THREE.Vector3(0, 1, 0);
let direction2 = new THREE.Vector3(0, 0, 1);

If I test raycast mesh1 with points 1 & 2 using direction1, both return as inside. If I test point1 with direction2:

• direction.dot(intersects[0].face.normal) > 0 return as inside
• (intersects.length % 2) == 0 ? "outside" : "inside" return as outside

If I test point2 with direction2:

• direction.dot(intersects[0].face.normal) > 0 return as outside
• (intersects.length % 2) == 0 ? "outside" : "inside" return as outside

[gkjohnson]

Do you recommend using the raycrossing algorithm which is mentioned in the wiki article instead? Unless I misunderstood, it should be something simple like this?

Nope! Checking the face orientation is good enough (and potentially faster) with well formed geometry like you have. Ray crossing is just a common algorithm for checking within bounds like lasso shapes, etc with no face orientation. I've used it in this lasso demo.

The issue I'm hitting with raycasting (which can be very likely due to my lack of knowledge / experience with raycasting) is that depends on the direction, results may vary.

From your code it looks like you're picking positions on the surface of the box to test with which can be ambiguous and numerically unstable. If you want to test whether a triangle is inside or outside a geometry I would pick the center of the triangle as a test point.

[giladdarshan]

Nope! Checking the face orientation is good enough (and potentially faster) with well formed geometry like you have. Ray crossing is just a common algorithm for checking within bounds like lasso shapes, etc with no face orientation. I've used it in this lasso demo.

From your code it looks like you're picking positions on the surface of the box to test with which can be ambiguous and numerically unstable. If you want to test whether a triangle is inside or outside a geometry I would pick the center of the triangle as a test point.

Ok, makes sense to only test the triangle's center of post-split triangles, this looks better?

let _v1 = new THREE.Vector3();
let triangleNormal = new THREE.Vector3();
let point = currentPolygon.triangle.getMidpoint(_v1);
currentPolygon.triangle.getNormal(triangleNormal);
// let direction = currentPolygon.vertices[0].normal.clone().normalize();
let direction = triangleNormal.normalize(); 

let _raycaster1 = new THREE.Raycaster(point, direction);
let intersects = _raycaster1.intersectObject(targetOctree.mesh);
let inside = false;
if (intersects.length) {
    if (direction.dot(intersects[0].face.normal) > 0) {
        inside = true;
    }
}
currentPolygon.setState(inside ? "inside" : "outside");

I thought about using the triangle's center in the beginning, I think I moved away because of lack of confidence in the triangle intersection algorithm I'm using (and the issue I mentioned). Going to give it another try.

[gkjohnson]

Yup!

let direction = triangleNormal.normalize();

There's not necessarily a need for this to be the face normal, though. You could use a consistent direction like 0, 0, 1

[giladdarshan]

There's not necessarily a need for this to be the face normal, though. You could use a consistent direction like 0, 0, 1

Thank you, this question has bothered me for a while. Alright, going to start working on some modifications based on what we discussed here and see how it goes.

On a side note, I went over the raycast flow in three-mesh-bvh and I have a couple of questions. From what I understand, the raycast function (called from acceleratedRaycast as bvh.raycast) recursively checks if the node is leaf or not. If it is a leaf node, it will call the intersectTris function and continue the raycast flow. (this part I'm ok with) If it is not a leaf node, it will check if the ray intersects the left / right node's bounding boxes and if they intersect, call the raycast function using the left / right node details.

1. If it's checking the node's bounding boxes before continuing to the raycast flow, it means raycast is not being performed on ALL triangles in the object like in regular threejs raycast, correct?
2. If question 1 is correct, how is it I'm getting the same number of intersection results as if I'm checking against all triangles in the object?

[gkjohnson]

If question 1 is correct, how is it I'm getting the same number of intersection results as if I'm checking against all triangles in the object?

Question 1 is correct. The bounds are there to cull triangles that are definitely not in the path of the ray. If the ray does not intersect the parent bounding box then it definitely does not intersect the child boxes or triangles so there's no use in checking them.

[giladdarshan]

If I take the hollowed / nested cube from earlier (8mm cube subtracted from 10mm cube) -

And use the following point & direction for the raycasting:

let point = new THREE.Vector3(0, 0, 4.5);
let direction = new THREE.Vector3(0, 0, 1);

If I use threejs's raycast or three-mesh-bvh's raycast I get 2 results in the intersection array. If I use the raycast function I'm working on for OctreeCSG and the following function to get the triangles to test against, I get 1 result in the intersection array, if instead I get all triangles in the octree I get 2 results like threejs and three-mesh-bvh.

getRayPolygons(ray, polygons = []) {
    if (this.polygons.length > 0) {
        for (let i = 0; i < this.polygons.length; i++) {
            if (polygons.indexOf(this.polygons[i]) === -1) {
                polygons.push(this.polygons[i]);
            }
        }
    }
    for (let i = 0; i < this.subTrees.length; i++) {
        if (ray.intersectsBox(this.subTrees[i].box)) {
            this.subTrees[i].getRayPolygons(ray, polygons);
        }
    }

    return polygons;
}

If I test the point (0, 0, -4.5), I get 6 results with threejs and three-mesh-bvh, with OctreeCSG and the above function I get 4 results, if I change to get all triangles I get 6 results.

Am I overthinking this? the discrepancy when I use the function above vs. getting all triangles makes me think that there can be time when I don't get the right results.

[giladdarshan]

Ok I'm not overthinking it, I don't get results when using the following example and the mentioned getRayPolygons, but when I change to get all the triangles in the octree (instead of only those in subtrees that the ray intersects with their bounding box) it works fine same as when I use threejs and three-mesh-bvh raycasts.

EDIT - I think I figured it out, when I build the octree, I expand the box of the subtrees by the polygons inside it but not the parent octrees so they ended up being smaller than the subtree and failed the ray.intersectsBox check, I now added a function to recursively go to the parent of the subtree and expand it as well. Now I get the relevant triangles and pass the raycast test.

Will continue testing.

[gkjohnson]

Without looking deeply at your code it can be hard to tell what might be going on. These types of spatial data structures can be difficult to debug and the algorithms are sensitive - it took awhile and a lot of tests to ensure three-mesh-bvh was returning the exact same results as three.js. Occasionally a change is made something breaks again, as well, despite my best efforts 😅

[giladdarshan]

Yep.

I'm still trying to find the best way to handle the inconsistent raycast intersection result when the point is on the boundary / surface of the object, do you have a suggestion?

Is there a downside to looping over the intersection results to check if the point is inside?

let _v1 = new THREE.Vector3();
let triangleNormal = new THREE.Vector3();
let point = currentPolygon.triangle.getMidpoint(_v1);
let direction = currentPolygon.plane.normal.clone().normalize();

let _raycaster1 = new THREE.Raycaster(point, direction);
let intersects = _raycaster1.intersectObject(targetOctree.mesh);
let inside = false;
if (intersects.length) {
    for (let i = 0; i < intersects.length; i++)
    if (direction.dot(intersects[i].face.normal) > 0) {
        inside = true;
        break;
    }
}
currentPolygon.setState(inside ? "inside" : "outside");

[gkjohnson]

there's likely not a great way to test a point on the surface of a mesh like that due to numerical precision issues so you might need to have a fallback in that case. And testing multiple intersections past the first does not seem like a valid approach.

[giladdarshan]

Hmm, what about looping over the 3 directions (0, 0, 1), (0, 1, 0) and (1, 0, 0)? In theory if a point is really outside it should be outside no matter the direction, right?

let _v1 = new THREE.Vector3();
let point = currentPolygon.triangle.getMidpoint(_v1);
let _raycaster1 = new THREE.Raycaster();
let inside = false;
let raycastDirections = [
    new THREE.Vector3(0, 0, 1),
    new THREE.Vector3(0, 1, 0),
    new THREE.Vector3(1, 0, 0)
];
_raycaster1.ray.origin.copy(point);
for (let i = 0; i < 3; i++) {
    _raycaster1.ray.direction.copy(raycastDirections[i]);
    let intersects = _raycaster1.intersectObject(targetOctree.mesh);
    if (intersects.length) {
        if (raycastDirections[i].dot(intersects[0].face.normal) > 0) {
            inside = true;
            break;
        }
    }
}
currentPolygon.setState(inside ? "inside" : "outside");

Did some tests and so far it works, will continue testing in the meantime.

[giladdarshan]

Another issue I came across with raycasting (unrelated to the above), for the below cone (CylinderGeometry) as target and triangle, raycasting shows that the midpoint of the triangle (4.440892098500626e-16, 1, -4.333333333333333) is inside the target cone:

let geometry = new THREE.CylinderGeometry( 5, 0, 10, 6 );
let material = new THREE.MeshBasicMaterial({color: 0x0000ff});
let mesh1 = new THREE.Mesh(geometry, material);
mesh1.material.side = THREE.DoubleSide;
// mesh1.geometry.computeBoundsTree();
let triangle = new THREE.Triangle(new THREE.Vector3(-3.464101791381836, 1, -5), new THREE.Vector3(3.4641017913818377, 1, -5), new THREE.Vector3(-4.440892098500626e-16, 1, -2.999999999999999));
let _v1 = new THREE.Vector3();
let point = triangle.getMidpoint(_v1); // (4.440892098500626e-16, 1, -4.333333333333333)
let direction = new THREE.Vector3(0, 0, 1);
let _raycaster1 = new THREE.Raycaster(point, direction);
let intersects = _raycaster1.intersectObject(mesh1);
if (intersects.length) {
    if (direction.dot(intersects[0].face.normal) > 0) {
        console.log("Point is inside the object", direction.dot(intersects[0].face.normal));
    }
    else {
        console.log("Point is outside the object", direction.dot(intersects[0].face.normal));
    }
}

The yellow circle is the point, blue is the cone -

There are two ways I found so far that works around the problem:

1. Change the radiusBottom from 0 to 0.01 (or any number that is not 0), while it does fix the problem, I would rather have CSG be able to work with cones.
2. Round the point from 16+ decimal points (point.x = 4.440892098500626e-16 which is 0.00000000000000004440892098500626) to 15 decimal points or lower.

What are your thoughts on this? I'm thinking maybe rounding the point's x,y,z for raycasting to prevent such small numbers?

function pointRounding(point) {
    point.x = +point.x.toFixed(15);
    point.y = +point.y.toFixed(15);
    point.z = +point.z.toFixed(15);
    return point;
}

EDIT - Seems like the point rounding workaround doesn't solve all the issues with cones, the below shows as inside as well. What are your thoughts on dealing with cones?

let geometry = new THREE.CylinderGeometry(0, 5, 10, 6);
let triangle = new THREE.Triangle(new THREE.Vector3(-2.5980763435363765, -1, -1.5000000000000002), new THREE.Vector3(4.440892098500626e-16, -1, -3), new THREE.Vector3(-2.5980763435363774, -1, -4.500000000000001));
let direction = new THREE.Vector3(1, 0, 0);

[giladdarshan]

Due to all the issues with raycast's precision I'm starting to think it might not be the right solution for determining if a point is inside or outside the target object for CSG operations since the results must be correct at all times or the result mesh / object won't be correct with extra triangles on the inside or missing triangles.

I've started looking into the Winding Numbers algorithm and adapted some C++ code I found online, while it does work (so far), it has obvious performance impact (due to the need to calculate the number using all triangles from target object with every point we need to check).

// Winding Number algorithm adapted from https://github.com/grame-cncm/faust/blob/master-dev/tools/physicalModeling/mesh2faust/vega/libraries/windingNumber/windingNumber.cpp
const EPSILON = 1e-5;
const _wV1 = new THREE.Vector3();
const _wV2 = new THREE.Vector3();
const _wV3 = new THREE.Vector3();
const _wP = new THREE.Vector3();
const _wP_EPS_ARR = [
    new THREE.Vector3(EPSILON, 0, 0),
    new THREE.Vector3(0, EPSILON, 0),
    new THREE.Vector3(0, 0, EPSILON),
    new THREE.Vector3(-EPSILON, 0, 0),
    new THREE.Vector3(0, -EPSILON, 0),
    new THREE.Vector3(0, 0, -EPSILON)
];
const _wP_EPS_ARR_COUNT = _wP_EPS_ARR.length;
const _matrix3 = new THREE.Matrix3();
const wNPI = 4 * Math.PI;

function calcWindingNumber(polygons, point) {
    let wN = 0;
    for (let i = 0; i < polygons.length; i++) {
        let pT = polygons[i].triangle;
        _wV1.subVectors(pT.a, point);
        _wV2.subVectors(pT.b, point);
        _wV3.subVectors(pT.c, point);
        let lenA = _wV1.length();
        let lenB = _wV2.length();
        let lenC = _wV3.length();
        _matrix3.set(_wV1.x, _wV1.y, _wV1.z, _wV2.x, _wV2.y, _wV2.z, _wV3.x, _wV3.y, _wV3.z);
        let omega = 2 * Math.atan2(_matrix3.determinant(), (lenA * lenB * lenC + _wV1.dot(_wV2) * lenC + _wV2.dot(_wV3) * lenA + _wV3.dot(_wV1) * lenB));
        wN += omega;
    }
    wN = Math.round(wN / wNPI);
    return wN;
}

function polyInside_WindingNumber(targetPolygons, polygon, boundingBox) {
    let result = false;
    let point = polygon.triangle.midPoint;
    if (!boundingBox.containsPoint(point)) {
        return false;
    }
    _wP.copy(point);
    let wN = calcWindingNumber(targetPolygons, _wP);
    if (wN === 0) {
        if (polygon.coplanar) {
            for (let j = 0; j < _wP_EPS_ARR_COUNT; j++) {
                _wP.copy(point).add(_wP_EPS_ARR[j]);
                wN = calcWindingNumber(targetPolygons, _wP);
                if (wN !== 0) {
                    result = true;
                    break;
                }
            }
        }
    }
    else {
        result = true;
    }
    return result;
}

I'm going to continue looking for better implementation since according to the wiki page there should be an algorithm with performance similar to raycast.

If there are any ideas how to approach the raycasting issues, I will be more than happy to hear, I'm out of ideas and worried by the correctness issues :sweat:

[gkjohnson]

I'm having trouble following exactly what the issues are but it's inevitable that there will be precision issues in some cases such as a point lying on the surface or a face or possibly raycasting between faces. It's something you have to build your code around being resilient to. In cases where something like a triangle center point is guaranteed to be on the surface (such as two geometries that meet at a coplanar surface) - that's a case you'll have to handle in other ways.

[giladdarshan]

Hey @gkjohnson, sorry for the spam and thank you again for the assistance! :)

It's something you have to build your code around being resilient to. In cases where something like a triangle center point is guaranteed to be on the surface (such as two geometries that meet at a coplanar surface) - that's a case you'll have to handle in other ways.

Understood, but the concern I have with raycasting is that these precision issues and their fixes are starting to conflict.

I'm having trouble following exactly what the issues are

I will summarize the issues to bullet points without examples/code to keep it short:

1. Triangle mid point on the surface sometimes being detected as inside and sometimes as outside depends on which surface the point is (which side of a cube) and which direction is being used. <- This one is ok-ish, since I know if the triangle is coplanar or not, I can do additional tests by adding an epsilon to the point and retest.
2. Triangle mid point is outside of the object (not on surface/boundary) and is detected as inside when it is supposed to be outside (noticed when the target object is a Cylinder with radius bottom of 0).
3. Triangle mid point is inside of the object (not on surface/boundary) and is detected as outside when it is supposed to be inside (noticed when the target object is a Icosahedron with radius 5 and detail 10).
4. With #⁠2 and #⁠3 conflicting, how can I guarantee the inside/outside classification will be correct?

Please let me know if you want code examples for any of the above, I tried to keep it short in this comment.

Edit - While testing I also test with threejs's builtin raycast and three-mesh-bvh's raycast in addition to the raycasting function I'm working on in OctreeCSG, the issues mentioned above applies to all.

[gkjohnson]

Can you recreate issues 2 and 3 in a js fiddle using just three.js' raycasting? This sounds like a logic issue.

[giladdarshan]

Issue 2 - https://jsfiddle.net/giladdarshan/mao1prs4/ Issue 3 - https://jsfiddle.net/giladdarshan/9r1te8pk/ I've added a yellow sphere to show where the point is, it's not playing any part in the calculations.

[gkjohnson]

This is a precision issue. If you change the ray origin to ( 0, -5, -4.333333333333333 ) rather than ( 4.440892098500626e-16, -5, -4.333333333333333 ) in the issue 2 example you'll see it works - there's not a whole lot to do about this, I don't think, without doing deep dives into three.js raycasting. Usually these are cases that are rare to hit, though. I'd have to think more about how it could be addressed more effectively.

[giladdarshan]

I thought so too, so I added the below function:

function pointRounding(point, num) {
    point.x = +point.x.toFixed(num);
    point.y = +point.y.toFixed(num);
    point.z = +point.z.toFixed(num);
    return point;
}

But it doesn't solve the issue completely, another example for issue 2 this time with rounding down the points - https://jsfiddle.net/giladdarshan/srbdxyhz/ In the example I rounded the point's x,y,z down to integers ((-1.732050895690918, -1, -3) to (-2, -1, -3)) and the raycasting result is still inside instead of outside.

Edit - Not sure if it matters, I've updated the examples and added three-mesh-bvh, just uncomment the mesh1.geometry.computeBoundsTree(); line.

[giladdarshan]

Hey @gkjohnson, I started playing around with Interactive CSG, I'm just starting but here is an example:

https://user-images.githubusercontent.com/75711568/171998134-2d58eb6d-b2f6-4dc3-9e5e-76a6b11f8241.mp4

If you want to test it, you can run the following code in 3dd.dev

OctreeCSG.useWindingNumber = false;
OctreeCSG.rayIntersectTriangleType = "MollerTrumbore";
// OctreeCSG.rayIntersectTriangleType = "regular";

function calcTransformMatrix(original, target) {
    let newMatrix = new THREE.Matrix4();
    newMatrix.elements[12] = target.elements[12] - original.elements[12];
    newMatrix.elements[13] = target.elements[13] - original.elements[13];
    newMatrix.elements[14] = target.elements[14] - original.elements[14];
    return newMatrix;
}
import('threeModules/controls/TransformControls.js').then((module) => {
    const { TransformControls } = module;
    const transformControl = new TransformControls(camera, renderer.domElement);
    transformControl.addEventListener('change', render);
    transformControl.addEventListener('dragging-changed', function (event) {
        orbitControl.enabled = !event.value;
    });

    let baseMaterial = new THREE.MeshPhongMaterial({ color: 0x0000ff });
    let baseCylinderGeometry = new THREE.CylinderGeometry(40, 40, 100, 12);
    let cylinderMesh1 = new THREE.Mesh(baseCylinderGeometry, baseMaterial.clone());

    cylinderMesh1.material.color.set(0x371d04);
    let windowMesh = createWindow();
    windowMesh.geometry.clearGroups();
    windowMesh.material.color.set(0x342f36);
    let windowDimensions = getDimensions(windowMesh.geometry);

    let hollowWindowGeometry = new THREE.BoxGeometry(windowDimensions.width - 0.1, windowDimensions.height - 0.1, windowDimensions.depth - 0.1);
    let hollowWindowMesh = new THREE.Mesh(hollowWindowGeometry, baseMaterial.clone());

    windowMesh.position.set(60, 0, 35);
    hollowWindowMesh.position.copy(windowMesh.position);

    let towerOctree = OctreeCSG.fromMesh(cylinderMesh1);
    let windowOctree = OctreeCSG.fromMesh(windowMesh);
    let hollowWindowOctree = OctreeCSG.fromMesh(hollowWindowMesh);
    let originalMatrix = windowMesh.matrix.clone();

    towerOctree.setPolygonIndex(0);
    windowOctree.setPolygonIndex(1);
    let towerWithWindowMesh = OctreeCSG.operation({
        op: "union",
        material: [cylinderMesh1.material, windowMesh.material],
        objA: {
            op: "subtract",
            objA: towerOctree.clone(),
            objB: hollowWindowOctree.clone()
        },
        objB: windowOctree.clone()
    });

    windowMesh.visible = false;
    scene.add(towerWithWindowMesh);
    scene.add(windowMesh);

    transformControl.attach(windowMesh);
    scene.add(transformControl);

    render();
    function render(event) {
        if (event && event.target) {
            if (event.target.dragging) {
                windowMesh.updateMatrix();
                let newTransformedMatrix = calcTransformMatrix(originalMatrix, windowMesh.matrix);
                windowOctree.applyMatrix(newTransformedMatrix);
                hollowWindowOctree.applyMatrix(newTransformedMatrix);
                originalMatrix.copy(windowMesh.matrix);
                towerWithWindowMesh.geometry.dispose();
                scene.remove(towerWithWindowMesh);
                towerWithWindowMesh = undefined;
                towerWithWindowMesh = OctreeCSG.operation({
                    op: "union",
                    material: [cylinderMesh1.material, windowMesh.material],
                    objA: {
                        op: "subtract",
                        objA: towerOctree.clone(),
                        objB: hollowWindowOctree.clone()
                    },
                    objB: windowOctree.clone()
                });
                scene.add(towerWithWindowMesh);
            }
        }
        renderer.render(scene, camera);
    }
});

function getDimensions(geometry) {
    geometry = geometry.isMesh ? geometry.geometry : geometry;
    geometry.computeBoundingBox();
    let dimensions = new THREE.Vector3().subVectors(geometry.boundingBox.max, geometry.boundingBox.min);
    return {
        width: dimensions.x,
        height: dimensions.y,
        depth: dimensions.z
    };
}

function createWindow() {
    let baseMaterial = new THREE.MeshPhongMaterial({ color: 0x0000ff });
    let backPlateGeometry = new THREE.BoxGeometry(20, 40, 2);
    let verticalWallGeometry = new THREE.BoxGeometry(2, 40, 20);
    let horizontalWallGeometry = new THREE.BoxGeometry(20, 2, 20);

    let backPlateMesh = new THREE.Mesh(backPlateGeometry, baseMaterial.clone());
    let verticalWallMesh = new THREE.Mesh(verticalWallGeometry, baseMaterial.clone());
    let verticalWallMesh2 = new THREE.Mesh(verticalWallGeometry.clone(), baseMaterial.clone());
    let verticalWallMesh3 = new THREE.Mesh(verticalWallGeometry.clone(), baseMaterial.clone());
    let horizontalWallMesh = new THREE.Mesh(horizontalWallGeometry, baseMaterial.clone());
    let horizontalWallMesh2 = new THREE.Mesh(horizontalWallGeometry.clone(), baseMaterial.clone());
    let horizontalWallMesh3 = new THREE.Mesh(horizontalWallGeometry.clone(), baseMaterial.clone());
    let horizontalWallMesh4 = new THREE.Mesh(horizontalWallGeometry.clone(), baseMaterial.clone());

    backPlateMesh.material.color.set(0x0ca0ad);
    verticalWallMesh.material.color.set(0x856614);
    verticalWallMesh2.material.color.set(0x856614);
    verticalWallMesh3.material.color.set(0x856614);
    horizontalWallMesh.material.color.set(0x856614);
    horizontalWallMesh2.material.color.set(0x856614);
    horizontalWallMesh3.material.color.set(0x856614);
    horizontalWallMesh4.material.color.set(0x856614);

    backPlateMesh.position.z = -9;
    verticalWallMesh.position.x = -9;
    verticalWallMesh3.position.x = 9;
    horizontalWallMesh.position.y = -19;
    horizontalWallMesh2.position.y = -9;
    horizontalWallMesh3.position.y = 9;
    horizontalWallMesh4.position.y = 19;
    let result = OctreeCSG.operation({
        op: "union",
        material: verticalWallMesh.material,
        objA: backPlateMesh,
        objB: {
            op: "union",
            objA: {
                op: "union",
                objA: verticalWallMesh,
                objB: verticalWallMesh2,
            },
            objB: {
                op: "union",
                objA: {
                    op: "union",
                    objA: verticalWallMesh3,
                    objB: horizontalWallMesh,
                },
                objB: {
                    op: "union",
                    objA: horizontalWallMesh2,
                    objB: {
                        op: "union",
                        objA: horizontalWallMesh3,
                        objB: horizontalWallMesh4,
                    },
                },
            },
        }
    });
    return result;
}

[gkjohnson]

Nice work! When I have a bit more time I might do a few experiments around real time CSG with webworkers - this has all piqued my interest again.

[giladdarshan]

Thank you, I started experimenting with web workers when I implemented the Winding Number algorithm and it works but I need to get more familiarized with web workers and promises to make it more efficient.

A couple of questions about web workers if I may:

1. What is the best practice when it comes to the number of web workers, spin up thousands of web workers at the same time or limit to the number of cores (using navigator.hardwareConcurrency)?
2. If you do limit the number of web workers, how do you efficiently wait for a web worker to become available?

I assume you will be experimenting real time CSG with three-mesh-bvh? I'm curious to see how you will overcome the precision issues with Three.js's built-in ray.intersectTriangle algorithm, I replaced it with Möller–Trumbore intersection algorithm for ray-triangle intersection tests in OctreeCSG and it works much better precision-wise.

[gkjohnson]

I'll deal with intersection issues when I get to them - haven't thought through it all yet. Regarding your webworker questions:

1 - I'm not sure about a best practice. I wouldn't spin up thousands of workers, though. I would start up to the hardwareConcurrency count or maybe a couple more and keep them around.

2 - In the past I've created a queue of tasks and pulled from a pool of workers and wait until one has finished to start the next task.

[gkjohnson]

Here's a little peak at some of the interactive CSG I've done so far. Still no workers but it's just a start. Once I'm further along I'll make the project public, as well. This is definitely a fun problem to work on!

I'm curious to see how you will overcome the precision issues with Three.js's built-in ray.intersectTriangle algorithm

For the moment I've addressed this issue by jittering the ray origin slightly and taking several samples and assuming the dominant results are correct. It would probably be best to use a more regular pattern for consistency and it's definitely not as robust as an improved triangle intersection function like you're using. At some point I think I'm going to move to a half edge data structure that can be traversed to avoid do so many raycasts.

[giladdarshan]

Awesome! I'm looking forward to see the code behind it.

If you want to test Möller–Trumbore intersection algorithm, here is a quick implementation that replaces the built-in one:

function rayIntersectsTriangle(ray, triangle, target = new THREE.Vector3()) {
    let edge1 = new THREE.Vector3();
    let edge2 = new THREE.Vector3();
    let h = new THREE.Vector3();
    let s = new THREE.Vector3();
    let q = new THREE.Vector3();

    const EPSILON = 0.0000001;
    edge1.subVectors(triangle.b, triangle.a);
    edge2.subVectors(triangle.c, triangle.a);
    h.crossVectors(ray.direction, edge2);
    let a = edge1.dot(h);
    if (a > -EPSILON && a < EPSILON) {
        return null; // Ray is parallel to the triangle
    }
    let f = 1 / a;
    s.subVectors(ray.origin, triangle.a);
    let u = f * s.dot(h);
    if (u < 0 || u > 1) {
        return null;
    }
    q.crossVectors(s, edge1);
    let v = f * ray.direction.dot(q);
    if (v < 0 || u + v > 1) {
        return null;
    }
    // Check where intersection is
    let t = f * edge2.dot(q);
    if (t > EPSILON) {
        return target.copy(ray.direction).multiplyScalar(t).add(ray.origin);
    }
    return null;
}
THREE.Ray.prototype.intersectTriangle = function(a, b, c, backfaceCulling, target) {
    return rayIntersectsTriangle(this, {a, b, c}, target);
}

[gkjohnson]

Thanks! I may give that a try once I'm further along - code is now available here if you want take a look!

https://github.com/gkjohnson/three-bvh-csg

[giladdarshan]

Nice, good job! I started going over the code, will take a deeper dive probably next week when I'm done with some small VR project for work.