CGAL: efficiently finding all triangles in a volume (bounding box or ball query)

Question

I want to efficiently find all triangles (e.g. facets of a polyhedral mesh) that are contained in or intersected by a volume (e.g. AABB or sphere query). I want to do this, if possible and reasonable, with CGALs built-in features.

My current solution is simple brute force: for all facets in mesh, check if distance of facet to ball center is smaller than the ball radius. I used a fuzzy sphere query of the KD-tree on the vertices before, but it was not accurate enough for my application. I need full sphere-triangle intersections.

The CGAL AABB tree (https://doc.cgal.org/latest/AABB_tree/index.html) seems like the best datastructure but I need the 3D linear Kernel which has no intersection test for triangles and any kind of volume (https://doc.cgal.org/latest/Kernel_23/group__intersection__linear__grp.html). The CGAL KD tree doesn't work because it can only contain points.

My ideas are:

1. Add a custom do_intersect() for Triangle_3 and Sphere_3 that the AABB tree can use. Is that even possible? This will probably require a custom squared_distance() as well.
2. Convert the volume query into two area queries by projecting the triangles on e.g. the XY and YZ planes. Can the AABB tree even handle 2D search? There is no intersection for circle and 2D triangle but I could use the intersection of Iso_rectangle_2 and Triangle_2 as good first guess.
3. Somehow traverse the internals of the AABB tree and stop before I find a node that doesn't contain my query anymore.
4. Modify the closest_point_and_primitive() method, give an optional max_distance parameter and return all primitives within that distance. This would require me to mess with the CGAL code.
5. Write my own datastructure. This will take some time to get it right.

Did I miss anything? Which solution would be the least effort?

Answer 1

Using some ideas from sloriot's answer, I was able to solve my problem.

Since the documentation of do_intersect() doesn't show an overload for Sphere_3 and Triangle_3, it's not surprising that the AABB tree doesn't support queries with Sphere_3.

However, the AABB tree supports queries with BBox_3, which is not mentioned in the do_intersect() docu.

Calling all_intersected_primitives() with the Bbox_3 returns all primitives of the AABB tree that are contained or intersection by the Bbox_3. This is a first good guess to get the actual intersections with the sphere.

With this optimization, I could reduce the query time from 20 ms (simple brute force) to 3 ms on a mesh with 100k faces where one query returns roughly 10 faces.

Here is the relevant code:

double r = CGAL::sqrt(patch_radius_squared);
CGAL::Bbox_3 query(
    sphere_center.x() - r, 
    sphere_center.y() - r, 
    sphere_center.z() - r, 
    sphere_center.x() + r, 
    sphere_center.y() + r, 
    sphere_center.z() + r);

std::list<Tree::Primitive_id> primitives;
tree.all_intersected_primitives(query, std::back_inserter(primitives));

std::vector<Triangle_3> intersected_facets;
for (const Tree::Primitive_id& p : primitives)
{
    // intersection with bb gives only a good first guess -> check actual intersection
    if (CGAL::squared_distance(*p, sphere_center) <= patch_radius_squared)
    {
        intersected_facets.push_back(*p);
    }
}

Answer 2

For the set of intersecting triangles you can use:

std::vector<Primitive> primitives;
tree.all_intersected_primitives(sphere, std::back_inserter(primitives));
//or
tree.all_intersected_primitives(tree2, std::back_inserter(primitives));

with tree2 being an AABB-tree of triangles of your bounding volume.

This will gives you the elements that are intersected by the boundary of your volume. Then you can use Surface_mesh properties to set a bool to true for all the faces that are intersected. Then create a new property on edges and set it to true if you set the property of one of the incident face was set to true.

Then calling connected_components() you'll be able to split the mesh into part inside and outside the bounding volume (ignoring the intersected part).

To finish, you need to pick one point per connected component and look if it is inside or outside of the bounding volume. It is straight forward for the sphere and you can use Side_of_triangle_mesh for the mesh case (without duplicating the AABB-tree to save memory and time).

If your bounding volume is a bbox you can do like for the sphere.