Performance of intersection of sphere with cuboid

[krober10nd]

Hello CGAL developers,

I have a psuedo-code posted below that seems to perform okay to detect whether a sphere intersects one of two axis-aligned cuboids. In my case, I have many spheres and one only two cuboids. The balls are circumballs of a 3d Delaunay triangulation. I need to repeat this calculation many times and it becomes performance critical for me to get it as fast possible.

I'd like to know and would be quite grateful if someone could comment on whether the below psuedo-code is the most performant way of accomplishing this task?

Should I rely on some other component of the CGAL library instead? Is there perhaps some better way to take advantage of the kernels/arithmetic type available in CGAL? My coordinates are Cartesian in double precision format.

Thank you in advance for any advice or feedback you can provide,

• Psuedo-code starts here

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Iso_cuboid_3.h>
#include <CGAL/Sphere_3.h>
#include <CGAL/intersections.h>

#include <assert.h>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point3;
typedef K::Sphere_3 Sphere;
typedef K::Iso_cuboid_3 Cuboid;

// build the two cuboids 
  Cuboid cube1 = Cuboid(Point3(llc[0 * 3], llc[0 * 3 + 1], llc[0 * 3 + 2]),
                Point3(urc[0 * 3], urc[0 * 3 + 1], urc[0 * 3 + 2]));
  Cuboid cube2 = Cuboid(Point3(llc[1 * 3], llc[1 * 3 + 1], llc[1 * 3 + 2]),
                Point3(urc[1 * 3], urc[1 * 3 + 1], urc[1 * 3 + 2]));

// loop through a lot of cells
for (std::size_t ie = 0; ie < candidate_cells.size(); ++ie) {

      int cell = cells[ie];

      // Indices of element into points
      int nm1 = faces[cell * 4];
      int nm2 = faces[cell * 4 + 1];
      int nm3 = faces[cell * 4 + 2];
      int nm4 = faces[cell * 4 + 3];

      // Coordinates of each vertex of element
      Point3 pnm1 =
          Point3(points[nm1 * 3], points[nm1 * 3 + 1], points[nm1 * 3 + 2]);
      Point3 pnm2 =
          Point3(points[nm2 * 3], points[nm2 * 3 + 1], points[nm2 * 3 + 2]);
      Point3 pnm3 =
          Point3(points[nm3 * 3], points[nm3 * 3 + 1], points[nm3 * 3 + 2]);
      Point3 pnm4 =
          Point3(points[nm4 * 3], points[nm4 * 3 + 1], points[nm4 * 3 + 2]);

      if (volume(pnm1, pnm2, pnm3, pnm4) == 0) {
        continue;
      }

      bool isCoplanar = CGAL::coplanar(pnm1, pnm2, pnm3, pnm4);
      if (isCoplanar) {
        continue;
      }

      // Calculate the circumball of cell
      Point3 cc = CGAL::circumcenter(pnm1, pnm2, pnm3, pnm4);
      double sqr_radius = CGAL::squared_radius(pnm1, pnm2, pnm3, pnm4);
      Sphere sphere = Sphere(cc, sqr_radius, CGAL::CLOCKWISE);

      // Does this circumball intersect with box above or box below?
      bool intersects;
      intersects = CGAL::do_intersect(sphere, cube1);
      if (intersects){
          result[iv] = 0;
          below += 1.0;
          break;
      }
     intersects = CGAL::do_intersect(sphere, cube2);
     if (intersects){
          result[iv] = 1;
          above += 1.0;
      }

    }
  }

[sloriot]

if you have many spheres vs only 2 cubes, I would simply add a filter test using bboxes of the cubes + the sphere before doing the test.
Something like CGAL::do_overlap(sphere.bbox(), cube1.bbox()) (storing the bboxes into local object to avoid rebuiding them again and again of course).

[krober10nd]

Thanks @sloriot so first building all the geometric objects and doing a bbox filter test to eliminate the potential unnecessary tests?

I suppose this filter test is extremely cheap in comparison to the sphere cube intersection test?

Do you have any recommendations for how to store many sphere objects? Can I pack them into a std::array for instance?

[sloriot]

I'd do it on-the-fly while in the loop, even before construction each sphere. Only the bbox for the 2 cubes can be outside of the loop.

[krober10nd]

this worked great and sped things up. thank you for the commentary and discussion!