# Efficient N-Dimensional AABB Subdivision Algorithm

This is supposed to subdivide a n-dimensional aabb (axis aligned bounding box) into 2ⁿ equal pieces as if sliced by n planes that are perpendicular to only one axis and parallel to the rest of the axes at a time.

The code:

CGL_bool CGL_aabb_subdivide_ndim(CGL_int n, CGL_float* aabb_min, CGL_float* aabb_max, CGL_float* aabbs_min_out, CGL_float* aabbs_max_out)
{
static CGL_float aabb_mid_pos[1024], aabb_axis_deltas[1024]; // I don't think we'll ever need more than 1024 dimensions
CGL_float nth_dim_val_holder = 0.0f, *pm_a[2];
CGL_int result_count = 1 << n, nth_dim_val_holder_index = 0;
for (CGL_int i = 0; i < n; i++)
{
aabb_mid_pos[i] = (aabb_min[i] + aabb_max[i]) * 0.5f;
aabb_axis_deltas[i] = aabb_mid_pos[i] - aabb_min[i];
}

for (CGL_int j = 0; j < n; j++)
{
pm_a[0] = &aabb_min[j]; pm_a[1] = &aabb_mid_pos[j];
nth_dim_val_holder = *pm_a[nth_dim_val_holder_index = 0];
for (CGL_int i = 0; i < result_count; i++)
{
aabbs_min_out[i * n + (n - j - 1)] = nth_dim_val_holder;
if ((i - 1) % (1 << j) == 0) nth_dim_val_holder = *pm_a[nth_dim_val_holder_index = (nth_dim_val_holder_index + 1) % 2];
}
}

for (CGL_int i = 0; i < result_count; i++)
{
for (CGL_int j = 0; j < n; j++)
{
aabbs_max_out[i * n + j] = aabbs_min_out[i * n + j] + aabb_axis_deltas[j];
}
}

return CGL_TRUE;
}


Here CGL_float, CGL_int, etc are just standard C types renamed with typedefs.

aabb_min is the min point of the n-dim box (it's an array of n elements) Similar for aabb_max. aabb_min_out and aabb_max_out are pre-allocated float arrays of size (2ⁿ * n) to store the result aabb min and max values for all the resulting aabb.

I tested this using:

CGL_float test_for_n_dim(CGL_int n)
{
CGL_float* aabb_minn = (CGL_float*)malloc(sizeof(CGL_float) * n); if (!aabb_minn) return;
CGL_float* aabb_maxn = (CGL_float*)malloc(sizeof(CGL_float) * n); if (!aabb_maxn) return;
for (CGL_int i = 0; i < n; i++) { aabb_minn[i] = 0; aabb_maxn[i] = 1; }
CGL_float* aabb_out_minn = (CGL_float*)malloc(sizeof(CGL_float) * n * (1 << n)); if (!aabb_out_minn) return;
CGL_float* aabb_out_maxn = (CGL_float*)malloc(sizeof(CGL_float) * n * (1 << n)); if (!aabb_out_maxn) return;
CGL_float time_before = CGL_utils_get_time();
CGL_aabb_subdivide_ndim(n, aabb_minn, aabb_maxn, aabb_out_minn, aabb_out_maxn);
CGL_float time_after = CGL_utils_get_time();
if (n <= 3)
{
CGL_info("For %dD", n);
for (CGL_int i = 0; i < (1 << n); i++)
{
CGL_printf_green("{(");
for (CGL_int j = 0; j < n; j++)
{
CGL_printf_green("%f", aabb_out_minn[i * n + j]);
if (j != n - 1) CGL_printf_green(", ");
}
CGL_printf_green(") -> (");
for (CGL_int j = 0; j < n; j++)
{
CGL_printf_green("%f", aabb_out_maxn[i * n + j]);
if (j != n - 1) CGL_printf_green(", ");
}
CGL_printf_green(")}\n");
}
}
free(aabb_minn); free(aabb_maxn);free(aabb_out_minn); free(aabb_out_maxn);
return time_after - time_before;
}


Results:

[INTERNAL] [Tue May  9 22:04:28 202] : Started Logger Session
[INFO] [Tue May  9 22:04:28 202] : Testing NDIM AABB Subdivision
[INFO] [Tue May  9 22:04:28 202] : For 1D
{(0.000000) -> (0.500000)}
{(0.500000) -> (1.000000)}
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000002 (1D)
[INFO] [Tue May  9 22:04:28 202] : For 2D
{(0.000000, 0.000000) -> (0.500000, 0.500000)}
{(0.000000, 0.500000) -> (0.500000, 1.000000)}
{(0.500000, 0.000000) -> (1.000000, 0.500000)}
{(0.500000, 0.500000) -> (1.000000, 1.000000)}
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000000 (2D)
[INFO] [Tue May  9 22:04:28 202] : For 3D
{(0.000000, 0.000000, 0.000000) -> (0.500000, 0.500000, 0.500000)}
{(0.000000, 0.000000, 0.500000) -> (0.500000, 0.500000, 1.000000)}
{(0.500000, 0.500000, 0.000000) -> (1.000000, 1.000000, 0.500000)}
{(0.500000, 0.500000, 0.500000) -> (1.000000, 1.000000, 1.000000)}
{(0.500000, 0.000000, 0.000000) -> (1.000000, 0.500000, 0.500000)}
{(0.500000, 0.000000, 0.500000) -> (1.000000, 0.500000, 1.000000)}
{(0.000000, 0.500000, 0.000000) -> (0.500000, 1.000000, 0.500000)}
{(0.000000, 0.500000, 0.500000) -> (0.500000, 1.000000, 1.000000)}
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000000 (3D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000001 (4D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000001 (5D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000001 (6D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000002 (7D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000005 (8D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000010 (9D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000030 (10D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000091 (11D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000207 (12D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.000394 (13D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.001144 (14D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.001890 (15D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.004921 (16D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.009553 (17D)
[INFO] [Tue May  9 22:04:28 202] : Time taken: 0.092833 (18D)
[INFO] [Tue May  9 22:04:29 202] : Time taken: 0.201626 (19D)
[INFO] [Tue May  9 22:04:29 202] : Time taken: 0.657279 (20D)
[INFO] [Tue May  9 22:04:30 202] : Time taken: 1.104113 (21D)
[INFO] [Tue May  9 22:04:32 202] : Time taken: 1.859703 (22D)
[INFO] [Tue May  9 22:04:36 202] : Time taken: 3.357090 (23D)
[INFO] [Tue May  9 22:04:43 202] : Time taken: 7.311569 (24D)
[INFO] [Tue May  9 22:05:09 202] : Time taken: 25.260529 (25D)


Please don't hide important code off the right-hand side of the screen like this:

    CGL_float* aabb_minn = (CGL_float*)malloc(sizeof(CGL_float) * n); if (!aabb_minn) return;
CGL_float* aabb_maxn = (CGL_float*)malloc(sizeof(CGL_float) * n); if (!aabb_maxn) return;


One line per statement, please! And make the return value consistent with the function's declaration - getting away with this is a clear sign that you didn't enable enough warning options when you compiled.

There's no need to cast the void* which malloc() returns. Also, it's generally clearer to use the assigned-to pointer in the sizeof expression rather than its type (this matters more when the declaration is far away):

    CGL_float* aabb_minn = malloc((sizeof *aabb_min) * n);


Did you see the bug in these lines? When the allocation for aabb_maxn fails, we leak the memory pointed to by aabb_minn. That's why we don't hide the clean-up code off the right-hand margin!

We really need something more like

    CGL_float *const aabb_minn = malloc((sizeof *aabb_minn) * n);
CGL_float *const aabb_maxn = malloc((sizeof *aabb_maxn) * n);
if (!aabb_minn || !aabb_maxn) {
free(aabb_minn);
free(aabb_maxn);
return 0.0/0.0;  /* NaN */
}

• The cast is so that this can be used from c++ too. Also i was wondering if there are any possible optimization in th algorithm itself Commented May 10, 2023 at 8:04
• Also i didn't bother checking if malloc failed properly as I was just testing the algorithm with some thow away code and that's not anything important Commented May 10, 2023 at 8:14
• In C++, prefer not to use malloc() - use a smart pointer instead (so you never need to explicitly free()). It's usually a mistake to share source code between languages like that (only do so for declarations that are accessed from both languages). Commented May 10, 2023 at 8:21
• This code is originally for c. But i want it to be able to compile with c++ too so. I guess a extra typecast doesn't hurt much. Commented May 10, 2023 at 8:31
• Still a bad idea. If you need it for programs in both languages, write it in C, and use extern "C" to link it into your C++ program. IMO, that's better than writing code that's simultaneously poor C and poor C++. Commented May 10, 2023 at 8:54

In The code, there is no indication what CGL_aabb_subdivide_ndim() is supposed to do beyond function and parameter names.
This makes it impossible to tell happenstance from requirement.
(Example: Is the order of result coordinates of any significance?)

// I don't think we'll ever need more than 1024 dimensions
So #define MAX_DIMENSIONS 1024, code a check and document behaviour in case n < 0, n = 0, n = 1, MAX_DIMENSIONS < n - as well as CGL_aabb_subdivide_ndim() being non-reentrant.
(1024 dimensions would mean about $$\2^{1037}\$$ bytes output.)

You assign 1 << n (not 1LL << n) to a variable of circumscribed type CGL_int.
That will not work from n = #bits in int - 1 up.

In the "main loop" assigning values to aabbs_min_out elements, why take and dereference addresses of values instead of just manipulating the values?

    for (int dim = 0 ; dim < n ; dim++) {
CGL_float
lower[] = { aabb_min[dim], aabb_min[dim] + aabb_max[dim]) * 0.5f },
upper[] = { lower[1], aabb_max[dim] };
for (int box = 0 ; box < result_count ; box++) {
aabbs_min_out[box * n + (n - dim - 1)] = lower[(box >> dim) & 1];
aabbs_max_out[box * n + (n - dim - 1)] = upper[(box >> dim) & 1];
}
}


Walking aabbs_min_out in steps of n isn't cache friendly; populating aabbs_max_out in this same loop doubles the pressure, worse if the memory system can't handle two address streams concurrently.

Your last loop handles coordinates for a single box in order of dimensions, while your main loop assigned mins in reverse order.
I can't imagine that to work - test with different sizes and offsets for each dimension.

Maybe better

    static const int MAX_DIMENSIONS = sizeof(size_t) * NBBY - sizeof(CGL_float);
if (n <= 0 || MAX_DIMENSIONS < n)
handle()
CGL_float aabb_mid_pos[MAX_DIMENSIONS];
for (int dim = 0 ; dim < n ; dim++)
aabb_mid_pos[dim] = (aabb_min[dim] + aabb_max[dim]) * 0.5f;
CGL_float
*lower[] = { aabb_min, aabb_mid_pos },
*upper[] = { aabb_mid_pos, aabb_max };
for (size_t box = 0, coordinate_index = 0 ; box < result_count ; box++)
for (int dim = 0, bit = box ; dim < n
; dim++, coordinate_index++, bit >>= 1) {
aabbs_min_out[coordinate_index] = lower[bit & 1][dim];
aabbs_max_out[coordinate_index] = upper[bit & 1][dim];
}
}

• sizeof(size_t) * NBBY - 1 - logb(sizeof(CGL_float))? Commented May 11, 2023 at 9:54