Efficient N-Dimensional AABB Subdivision Algorithm

Question

This is supposed to subdivide a n-dimensional aabb (axis aligned bounding box) into 2n 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 (2n * 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)

Answer 1

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 */
 }

• Jaysmito Mukherjee: 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
• Jaysmito Mukherjee: 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
• Toby Speight: 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).
• Jaysmito Mukherjee: 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.
• Toby Speight: 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++.

Answer 2

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];
 }
 }

• greybeard: sizeof(size_t) * NBBY - 1 - logb(sizeof(CGL_float))?