In order to represent spheres with octal trees, I need to determine overlaps between spheres and axis-aligned cubes. More specifically, there are three cases :
I'm asking here to know if there are any simpler ways (meaning lesser complexity / fewer operations)
Thanks to an answer to this previous question, I found this algorithm to detect any overlap. However, it is not able to distinguish between the cube being fully inside the sphere, or just partially.
According to suggestions given here, I'm first checking for total overlap by checking whether the furthest corner of the cube is inside the sphere. If it isn't, then we use Jim Arvo's algorithm to determine if there is any overlap at all.
So here is what I came up with : (using c23)
#include <OpenGL/OpenGL.h>
#include <math.h>
/* A 3D point with type-punning to access coordinates by name or index */
typedef union point_3d {
struct {
GLdouble x, y, z;
} coord;
GLdouble t[ 3 ];
} point_3d;
/* An axis-aligned cube */
typedef struct cube {
/* corner with smallest coordinates */
point_3d p_min;
/* side length */
GLdouble len;
} cube_t;
typedef struct sphere {
unsigned int radius;
point_3d center;
} sphere_t;
/* Clearly defined values for our result */
typedef enum overlap {
no_overlap = -1,
partial_overlap = 0,
total_overlap = 1,
} overlap_t;
static inline GLdouble sq( GLdouble val ) {
return val * val;
}
static bool point_in_sphere( point_3d p, struct sphere s ) {
GLdouble dx = p.t[ 0 ] - s.center.t[ 0 ];
GLdouble dy = p.t[ 1 ] - s.center.t[ 1 ];
GLdouble dz = p.t[ 2 ] - s.center.t[ 2 ];
return sq( s.radius ) >= sq( dx ) + sq( dy ) + sq( dz );
}
overlap_t cube_overlaps_sphere( cube_t c, sphere_t s ) {
point_3d p_min = c.p_min;
point_3d p_max = { .t= {
p_min.t[0] + c.taille,
p_min.t[1] + c.taille,
p_min.t[2] + c.taille,
} };
point_3d c_centre = { .t = {
p_min.t[0] + c.taille / 2,
p_min.t[1] + c.taille / 2,
p_min.t[2] + c.taille / 2,
} };
point_3d v = { .t = {
c_centre.t[ 0 ] - s.centre.t[ 0 ],
c_centre.t[ 1 ] - s.centre.t[ 1 ],
c_centre.t[ 2 ] - s.centre.t[ 2 ],}
};
point_3d far_corner = { .t = {
v.t[ 0 ] < 0 ? p_min.t[ 0 ] : p_max.t[ 0 ],
v.t[ 1 ] < 0 ? p_min.t[ 1 ] : p_max.t[ 1 ],
v.t[ 2 ] < 0 ? p_min.t[ 2 ] : p_max.t[ 2 ],
} };
if ( point_in_sphere( far_corner, s ) ) {
return total_overlap;
}
/* Algorithm by Jim Arvo */
GLdouble dmin = 0;
for ( unsigned int i = 0; i < 3; i++ ) {
if ( s.center.t[ i ] < p_min.t[ i ] ) {
dmin += pow( s.center.t[ i ] - p_min.t[ i ], 2 );
} else if ( s.center.t[ i ] > p_max.t[ i ] ) {
dmin += pow( s.center.t[ i ] - p_max.t[ i ], 2 );
}
}
if ( dmin <= pow( s.radius, 2 ) ) {
return partial_overlap;
} else {
return no_overlap;
}
}
The number of operations required to compute your overlap parameter scales linearly in the number of spacial dimensions. For any particular numbers of dimensions, it is O(1). The strategy that yields the fewest operations for some sphere, cube pairs will be different than that for some other pairs, so best strategy overall for a given application will probably be the one that is best for the largest number of cases tested. For your particular application, I anticipate that that will be the partial overlap case, because it is near the surface of the sphere that you will be branching most deeply.
The Avro algorithm you have incorporated computes the shortest distance from the center of the sphere to any part of the solid cube (corner, edge, face, or interior), and compares that to the sphere's radius to determine whether there is any overlap. It is effective and efficient. I don't see a way to improve on that for the job it does.
Avro does not distinguish between partial and complete overlap of the cube, whereas you need to do. Courtesy of @EricPostpischil, however, we have the observation that if the corner of the cube most distant from the center of the sphere is in or on the surface of the sphere, then the cube is completely encompassed, otherwise not, and also that we can determine which corner that is from the vector from sphere center to cube center. It turns out that a variation on that computation can be structured analogously enough to the Avro algorithm to combine the two in the same loop. That might look like so:
overlap_t compute_overlap(cube_t cube, sphere_t sphere) {
GLdouble d2min = 0; // minimum sq distance from the center of the sphere to any part of the cube
GLdouble d2max = 0; // maximum sq distance from the center of the sphere to [a corner of] the cube
for (int i = 0; i < 3; i++) {
GLdouble delta = sphere.center[i] - cube.pmin[i];
if (delta < cube.len / 2) {
// the sphere center is on the decreasing side of the cube center
d2max += sq(delta - cube.len);
if (delta < 0) {
// the sphere center is on the decreasing side of the whole cube
d2min += sq(delta);
}
} else {
// the sphere center is on the increasing side of the cube center
d2max += sq(delta);
if (delta > cube.len) {
// the sphere center is on the increasing side of the whole cube
d2min += sq(delta - cube.len);
}
}
}
if (sq(sphere.radius) <= d2min) {
return no_overlap;
} else if (sq(sphere.radius) < d2max) {
return partial_overlap;
} else {
return total_overlap;
}
}
I can't prove that that is optimal for your case, and especially not that there are no minor rearrangements that your particular compiler's optimizer might find more pleasing. I do, however, think that modulo such minor micro-optimizations, you're going to have trouble doing better for your case.