If I have a 3d grid of cubes, how can I efficiently find all the cubes that intersect a sphere of radius r and position p?

Question

Imagine you have a 3D space between 0 and d meters in the x,y and z directions (it's a cube of length d). now let's assume I want to break this up into a 3D array of smaller cubes. We want n cubes in each direction. This means each cube has a length, width and height of d/n (you can assume that d/n divides perfectly). This can be held in the following array:

Cube* cubes[n][n][n];

Now let us assume we can have a sphere of radius r and position p, where p can be anywhere inside the 3d space, and r has no constraints. What is the most efficient way of finding and returning all the cubes that intersect with this sphere?

Currently, I am looping through all the cubes in the array and using the Pythagorean theorem to check the distance between each cube and p, which I wrote the pseudo code for below.

Cube* getIntersectingCubes(r,p) { 
    Cube* cubes[n][n][n];
    Cube* result[];
    for(int i=0;i<n;i++) {
        for(int j=0;j<n;j++) {
            for(int k=0;k<n;k++) {
                if(distanceBetween(cubes[i][j][k],p) < r) {
                    result.add(cubes[i][j][k]);
                }
            }
        }
    }
    return result;
}

However, there is definitely a better way to do this. How could I return the same output but check as few cubes as possible?

Answer 1

If you are still having trouble working out a coordinate system and your algorithm, start with a diagram:

Essentially what you have is a 3D grid of cubes and your points p can be chosen any way you like so long as they describe a cube completely. What you want to determine is for which cubes does the radius r intersect any part of the cube. You can handle that with three vectors. You simply treat the cube origin as a vector. You know the lengths of the sides of each cube (d / n) you can use to form a second vector (call it the extent of p or p_extent) which you can add to the origin vector to determine the farthest point in the cube. So p_extent would equal d/ni + d/nj + d/nk.

You can create a struct to help in handling the coordinates as a vector, e.g.:

typedef struct {
    double x, y, z;
} vect_pos;

With your two vectors, you can simply add them together (e.g. p + p_extent) to form a vector describing the farthest point in the cube from the origin. (you want the direction of the vector inside your cube to have the same sign direction as your origin vector (so you can simply multiply it with a unit vector in that direction). For example is your point is -3i -1j and 2k, you can multiple by the p_extent vector by -1i -1j 1k. (picture the cube origin vector in each of the 8-quadrants and the values that x, y, z would have in each.

To add two vector_pos together, you can simply use vector addition, e.g.

/* vector addition */
vect_pos vect_add (const vect_pos *a, const vect_pos *b)
{
    vect_pos c = { .x = 0. };
    
    c.x = a->x + b->x;
    c.y = a->y + b->y;
    c.z = a->z + b->z;
    
    return c;
}

The last wrinkle that choosing the coordinate at a corner of the cube will require you solve is that the origin vector will not always point to the corner closest to the origin (picture in the diagram above where p would point if you rotated the origin counter-clockwise by 90 deg about the axis. (your alternative is to choose the center of the cube as its origin, but then you have twice the number of calculations).

You can handle the closest corner issue by a simple check and shifting the corner for the origin by a cube side length for purposes of the calculation of whether r is between the closest and farthest corners (the cube is still the same cube) With your cube origin and the farthest corner described, you can simply use the vector-distance formula to calculate the distance each is from the origin and if r is between the two, the sphere intersects that cube.

To choose the closest point to the origin before calculating the distances you can use:

/* point in cube closest to origin */
vect_pos closest_pt (const vect_pos *p, const int len)
{
    vect_pos c = {  p->x >= 0 ? p->x : p->x + len,
                    p->y >= 0 ? p->y : p->y + len,
                    p->z >= 0 ? p->z : p->z + len  };
    
    return c;
}

(choosing your cube at origin 0, 0, 0 with positive sides allows you to simply add a cube side length in the case that coordinate is negative to choose the closest point on the cube as the origin)

and for the distance calculations, a simple square-root of the sum of the squares:

/* vector distance */
double vect_dist (const vect_pos *p)
{
    return sqrt (p->x * p->x + p->y * p->y + p->z * p->z);
}

To check if r falls between the closest and farthest distance within each cube you can put it altogether with a short function that handles choosing the extent vector in the same direction as the origin, adds the vectors and then compares the distances of each to r, e.g.

/* does r intersect cube at postition p */
int intersect (vect_pos *p, const double r, const int len)
{
    vect_pos c = closest_pt (p, len),                   /* closest pt to origin */
            p_extent = { .x = c.x >= 0 ? len : -len,    /* length, unit vector */
                         .y = c.y >= 0 ? len : -len,    /* pointing outward */
                         .z = c.z >= 0 ? len : -len },
             p_farthest = vect_add (&c, &p_extent);     /* farthest point in cube */
    double   nearest  = vect_dist (&c),                 /* distance to nearest */
             farthest = vect_dist (&p_farthest);        /* distance to farthest */
    
    return nearest <= r && r <= farthest;               /* return radius in between */
}

That's basically it. You can add a short function to output the points to as a convenience, e.g.

void prn_pos (const vect_pos *p)
{
    printf ("(% d, % d, % d)\n",
            (int)p->x, (int)p->y, (int)p->z);
}

and then add a main() that allows entry of r, d and n and outputs the results:

int main (int argc, char **argv) {
    
    double  r = argc > 1 ? strtod (argv[1], NULL) : .9;     /* r */
    int     d = argc > 2 ? strtol (argv[2], NULL, 0) : 1,   /* d */
            n = argc > 3 ? strtol (argv[3], NULL, 0) : d,   /* n */
            dn = 0,             /* d/n */
            nc = 0,             /* number of cubes in each direction */
            total = 0,          /* total cubes in search grid */
            intersecting = 0;   /* number of cubes that intersect sphere */
    
    if (r < 0 || d < 0) {   /* validate length inputs */
        fputs ("error: negtive length in input of r or d.\n", stderr);
        return 1;
    }
    if (!n || n < 0) {      /* validate n not zero or negative */
        fputs ("error: n - divide by zero or negative.\n", stderr);
        return 1;
    }
    if (d < r) {            /* ensure d >= r */
        int min_d = ceil (r);
        n *= min_d / d;
        d = min_d;
    }
    if (d % n) {            /* validate d equally divisible by n */
        fputs ("error: d not equally divisible by n.\n", stderr);
        return 1;
    }
    
    dn = d / n;                     /* cube side length */
    nc = ceil (r / dn) + 1;         /* limit of cubes needed per-direction */
    total = nc * nc * nc * 8;       /* total number of cubes in search grid */
    
    printf ("\nOut of %d cubes, the following intersect sphere of radius %.2f\n\n",
            total, r);
    
    for (int i = -d; i <= d; i += dn)               /* loop -d to d */
        for (int j = -d; j <= d; j += dn)
            for (int k = -d; k <= d; k += dn) {
                vect_pos p = { i, j, k };           /* vector to cube origin */
                if (intersect (&p, r, dn)) {        /* does it intersect sphere? */
                    prn_pos (&p);                   /* output cube origin */
                    intersecting++;                 /* increment count */
                }
            }
    
    printf ("\nintersecting cubes: %d\n", intersecting);
}

Add the headers:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

and then compile linking against the math library -lm and you have a working program to test. If no arguments for r, d or n are provided, the code uses .90, 1, and 1 as defaults. Note the loop limits are at minimum from -(next int above r) to (next int above r) in each dimension or as given by d whichever is greater and not shorter than r.

Example Use/Output

Default case with r = .9 and d = 1 and n = 1:

$ ./bin/cube-sphere-intersect

Out of 64 cubes, the following intersect sphere of radius 0.90

(-1, -1, -1)
(-1, -1,  0)
(-1,  0, -1)
(-1,  0,  0)
( 0, -1, -1)
( 0, -1,  0)
( 0,  0, -1)
( 0,  0,  0)

intersecting cubes: 8

If you draw the cubes out on a grid, those are the 8 cubes surrounding the origin.

Choosing r = 5.5, d = 6 and n = 1 would give the same result other than the sides of the cubes being of different length:

$ ./bin/cube-sphere-intersect 5.5 6 1

Out of 64 cubes, the following intersect sphere of radius 5.50

(-6, -6, -6)
(-6, -6,  0)
(-6,  0, -6)
(-6,  0,  0)
( 0, -6, -6)
( 0, -6,  0)
( 0,  0, -6)
( 0,  0,  0)

intersecting cubes: 8

If you enter an r length of exacly 1, then this is where your choice of whether to consider r equal to the min or max length will come into play. There are arguments for saying "touching" a cube doesn't really mean it intersects as the point will technically lie on a tangent that distance out, but there are equally good arguments that say that intersects (the choice of whether you want to use <= or < in the comparison is left to you)

For example if you enter anything between 1.0 < sqrt(2) (1.414...) you will find 32 cube intersect, e.g.

$ ./bin/cube-sphere-intersect 1.1

Out of 216 cubes, the following intersect sphere of radius 1.10

(-2, -1, -1)
(-2, -1,  0)
(-2,  0, -1)
(-2,  0,  0)
(-1, -2, -1)
(-1, -2,  0)
(-1, -1, -2)
(-1, -1, -1)
(-1, -1,  0)
(-1, -1,  1)
(-1,  0, -2)
(-1,  0, -1)
(-1,  0,  0)
(-1,  0,  1)
(-1,  1, -1)
(-1,  1,  0)
( 0, -2, -1)
( 0, -2,  0)
( 0, -1, -2)
( 0, -1, -1)
( 0, -1,  0)
( 0, -1,  1)
( 0,  0, -2)
( 0,  0, -1)
( 0,  0,  0)
( 0,  0,  1)
( 0,  1, -1)
( 0,  1,  0)
( 1, -1, -1)
( 1, -1,  0)
( 1,  0, -1)
( 1,  0,  0)

intersecting cubes: 32

Then for radius between sqrt(2) <= sqrt(3) you would find 56 cube intersect and so on.

This is just one way to implement it. If you are curious, go rewrite it to use the center of each cube as the origin of the cube. You don't have to use a struct, you can just operate on each x, y, z directly, but the struct helps keep the logic clear. Look things over and let me know if you have further questions.