32x32x32 Units Octree With 15-Bit Voxels: Second Edition

Question

A while back I made an octree system that stores 15-bit voxel in a 32x32x32 voxel octree: 32x32x32 units octree, that supports 15-bit data units

I have since improved upon this system. The new octree structure is largely the same except for a key difference, being that that index nodes are now relative to the branch (referred to as set in the code) they are in as opposed to being relative to the start of the entire structure. The variables have been given more descriptive names.

Here is a description of the new system:

• Voxels are stored as 15-bit values in an octree comprising a 32^3 region of voxels.
• Each node is a 16-bit object that is either a 15-bit value or a 15-bit index to another set of 8 nodes, as determined by the 16th bit. This index is relative to the start of the set the node is in. The relative indices mean that any given branch, including all its descendants, are independent of where they are in the tree, vastly reducing the number of indices to update when inserting or deleting items in the center of the tree.
• Maximum of 5 layers (6 including the root node (Octree.base)), so for nodes in the 5th-divided layer, the 16th bit no longer affects structure of the octree and is used for a different purpose, specifying if the object is either a 15-bit value or yet another 15-bit index to a list of 64-bit objects (Octree.data), allowing for more complicated voxels that require extra data to be stored in the tree, without inflating the size of all voxels.
• The idea is to have a recursive hierarchical structure, where any given node can represent the 8 nodes in the next divided layer, so if any given set of 8 nodes are identical, their value can just be stored once in the previous divided layer. This has potential memory saving and performance benefits as any 8 identical nodes do not need to be individually stored or handled.
• For example, if a region contains only 0s, the 0 is stored only once as the root node, instead of being stored 32^3 = 32768 times individually as would be the case with an array. And checking that this region is all 0s would entail one operation, checking that the root node is 0, which is much faster than checking if 32768 objects are 0 individually.
• When the octree is reallocated, the new capacity is 2^Octree.set_alloc plus all powers of 8 less than 2^Octree.set_alloc where Octree.set_alloc is set to the smallest solution where the resulting capacity is not smaller than the new size. Since this scheme naturally leads up to the maximum possible size of the octree (4681 sets of 8 nodes) where all nodes are fully divided, it seems less arbitrary, (to me at least) than simply doubling each time and capping the allocation size at the maximum size.
• The purpose of the node_dup variables are vectorization, to assign and compare the values 4 nodes (8 bytes) at at time.
• The purpose of the extra fields in the NodeN unions is again vectorization, so that in the future I can add other code that can handle 2 or 4 nodes at a time.

The code:

#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <assert.h>
#include <stdio.h>
typedef union Node2 {
    uint16_t x[2];
    uint32_t z;
} Node2;
typedef union Node4 {
    uint16_t x[4];
    Node2 z[2];
    uint64_t y;
} Node4;
typedef union Node8 {
    uint16_t x[8];
    Node2 z[4];
    Node4 y[2];
} Node8;
typedef union Ptr {
    union Ptr *p;
    uint64_t u;
} Ptr;
typedef struct Octree {
    Ptr data;
    uint8_t data_alloc, set_alloc;
    uint16_t data_size, set_size, base;
    Node8 set[];
} Octree;
static_assert(sizeof(Node8)==sizeof(uint16_t[8]), "sizeof(Node8)!=sizeof(uint16_t[8])");
static_assert(sizeof(Ptr)==sizeof(uint64_t), "sizeof(Ptr)!=sizeof(uint64_t)");
static_assert(sizeof(Octree)==sizeof(Node8), "sizeof(Octree)!=sizeof(Node8)");
typedef struct State {
    unsigned level, offset[5];
} State;
static unsigned octree_index(const unsigned x, const unsigned z, const unsigned y, const unsigned level) {
    return (x>>level&1)|(z>>level&1)<<1|(y>>level&1)<<2;
}
unsigned octree_get(Octree *const octree, State *const state, const unsigned x, const unsigned z, const unsigned y) {
    unsigned node = octree->base, set = 0;
    state->level = 5;
    while (node&0x8000 && state->level) {
        --state->level;
        node = ((uint16_t *)octree->set)[state->offset[state->level] = (set += node&0x7FFF)<<3|octree_index(x, z, y, state->level)];
    }
    return node;
}
static size_t octree_alloc(const unsigned alloc) { // Calculate the size from an allocation step
    const size_t power = 1<<alloc;
    return (power-1&4681)|power>>1;
}
Octree *octree_set(Octree *octree, const unsigned x, const unsigned z, const unsigned y, const unsigned new) {
    State state;
    const unsigned node = octree_get(octree, &state, x, z, y);
    if (node != new) {
        if (state.level) {
            const unsigned prev_size = octree->set_size;
            octree->set_size += state.level;
            { // Reallocate the tree
                size_t alloc_size;
                unsigned alloc = octree->set_alloc;
                while ((alloc_size = octree_alloc(alloc)) < octree->set_size)
                    ++alloc;
                if (octree->set_alloc != alloc) {
                    if (!(octree = realloc(octree, alloc_size*sizeof(Node8)+sizeof(Octree))))
                        abort();
                    octree->set_alloc = alloc;
                }
            }
            Node8 *set_ptr = octree->set;
            if (state.level == 5)
                octree->base = 0x8000;
            else {
                unsigned set = 0;
                unsigned level = state.level;
                do { // This loop does 2 things: Find the insertion point for new branches, and offsets all the relevant indices in all the ancestry
                    unsigned offset = state.offset[level];
                    while (++offset&7) {
                        uint16_t *const node_ptr = (uint16_t *)octree->set+offset;
                        if (*node_ptr&0x8000) {
                            if (!set)
                                set = (offset>>3)+(*node_ptr&0x7FFF);
                            *node_ptr += state.level;
                        }
                    }
                } while (++level != 5);
                const unsigned offset = state.offset[state.level];
                ((uint16_t *)octree->set)[offset] = (set // If no set was found, it means there is nothing to shift over, and the new data should be placed at the end of the tree
                    ? (set_ptr += set, memmove(set_ptr+state.level, set_ptr, (prev_size-set)*sizeof(Node8)), set)
                    : (set_ptr += prev_size, prev_size)
                )-(offset>>3)|0x8000;
            }
            uint16_t *node_ptr;
            Node4 node_dup;
            node_dup.x[1] = node_dup.x[0] = node;
            node_dup.z[1] = node_dup.z[0];
            for (;;) {
                set_ptr->y[1] = set_ptr->y[0] = node_dup;
                node_ptr = (uint16_t *)set_ptr+octree_index(x, z, y, --state.level);
                if (state.level) {
                    *node_ptr = 0x8001;
                    ++set_ptr;
                    continue;
                }
                break;
            }
            *node_ptr = new;
        } else {
            unsigned set;
            Node4 node_dup;
            node_dup.x[1] = node_dup.x[0] = new;
            node_dup.z[1] = node_dup.z[0];
            for (;;) {
                const unsigned offset = state.offset[state.level], next_set = offset>>3;
                ((uint16_t *)octree->set)[offset] = new;
                Node8 *const set_ptr = octree->set+next_set;
                if (set_ptr->y[0].y == node_dup.y && set_ptr->y[1].y == node_dup.y) {
                    set = next_set;
                    if (++state.level == 5) { // If level is 5, then all of the data is in the entire octree is the same, and all the sets can be deleted right off the bat, however this is extremely rare, so is including these lines worth it instead of just skipping to the reallocation part?
                        octree->set_alloc = 0;
                        octree->set_size = 0;
                        octree->base = new;
                        Octree *const prev_octree = realloc(octree, sizeof(Octree));
                        return prev_octree ? prev_octree : octree;
                    }
                    continue;
                }
                break;
            }
            if (state.level) {
                Node8 *set_ptr = octree->set+set;
                memmove(set_ptr, set_ptr+state.level, (octree->set_size-set)*sizeof(Node8));
                unsigned level = state.level;
                do {
                    uint_fast16_t offset = state.offset[level];
                    while (++offset&7) {
                        uint16_t *const node_ptr = (uint16_t *)octree->set+offset;
                        if (*node_ptr&0x8000) {
                            *node_ptr -= state.level;
                        }
                    }
                } while (++level != 5);
                octree->set_size -= state.level;
                size_t alloc_size;
                unsigned alloc;
                { // Calculate new allocation size, starting at the current size and then stopping once the newest try is less than the actual octree size
                    size_t next_alloc_size;
                    unsigned next_alloc = octree->set_alloc;
                    while (alloc = next_alloc, (next_alloc_size = octree_alloc(--next_alloc)) >= octree->set_size) {
                        alloc_size = next_alloc_size;
                    }
                }
                if (octree->set_alloc != alloc) {
                    Octree *const prev_octree = realloc(octree, alloc_size*sizeof(Node8)+sizeof(Octree));
                    if (prev_octree)
                        octree = prev_octree;
                    octree->set_alloc = alloc;
                }
            }
        }
    }
    return octree;
}

Some concerns:

• Stability: This code very complicated and, while it has passed some brief testing that should invoke all code paths, such as inserting nodes before and after other nodes, to see if the parent indices are updated properly, etc. I am still not 100 percent confident that all of this code is 100 percent correct.
• Duplicate code: The code to update the indices in the parents after making new space for new branches is almost the same as the code to update the indices in the parents after filling in a gap after deleting a branch (if all nodes in a branch are identical, the branch is deleted and the node referencing it is simply replaced with the value), except for the if (!set) set = (offset>>3)+(*node_ptr&0x7FFF); and += parts. The reason the first line is necessary in the first part but not the second is the correct insertion point for the new branches needs to be located, and it just so happens that the first parent index that needs to be updated is a reference to the location where the new branches need to be inserted. Related: How can I make code that is both DRY and fast where intermediate values in a calculation may or may not be needed?
• Portability: As usual, I try to make my code portable within reason (yes, uintN_t types are optional and unavailable on some systems) but other than that, I try to avoid compiler extensions or anything not outlined in the standard.
• Uninitialized variable warnings: Some variables such as set and alloc_size in the second branch will in fact never be used while they are not initialized, but the compiler is not so sure. Can the logic be improved so that these warnings will disappear without the cop-out of initializing the variables to silence the warnings?
• Is the special case that all of the nodes are the same and all the sets can be deleted immediately being handled by its own few lines of code worth it?

I would like an overall review, but with some emphasis on the concerns listed above.

Here is the code I used to test the system (no review required):

void dump(Octree *octree) {
    printf("%d\n", octree->base);
    for (int i = 0; i < octree->set_size; ++i) {
        for (int j = 0; j < 8; ++j) {
            printf("%d ", octree->set[i].x[j]);
        }
        putchar('\n');
    }
}
#include <stdio.h>
size_t malloc_size(const void *);
int main(void) {
    Octree *octree = memset(malloc(sizeof(Octree)), 0, sizeof(Octree));
    octree->base = 0;
    for (;;) {
        char buf[64], *p = buf;
        fgets(buf, 64, stdin);
        switch (*(p++)) {
            case 's':
                octree = octree_set(octree, strtol(p, &p, 0), strtol(p, &p, 0), strtol(p, &p, 0), strtol(p, NULL, 0));
                break;
            case 'g':
                {State state; printf("%d %d\n", octree_get(octree, &state, strtol(p, &p, 0), strtol(p, &p, 0), strtol(p, &p, 0)), state.level);}
                break;
            case 'm':
                printf("%d %zu\n", octree->set_size, malloc_size(octree));
                break;
            case 'd':
                dump(octree);
                break;
            case 'q':
                free(octree);
                return 0;
        }
    }
}

Here is a test that proves that deletion and insertion work (this is command line I/O while using the test code):

s 31 31 31 1
s 15 15 15 2
s 3 3 3 4
s 1 1 1 5
d
32768
32769 0 0 0 0 0 0 32777 
32769 0 0 0 0 0 0 32773 
32769 0 0 0 0 0 0 0 
32769 0 0 0 0 0 0 32770 
0 0 0 0 0 0 0 5 
0 0 0 0 0 0 0 4 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 2 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 1 
s 1 1 1
d
32768
32769 0 0 0 0 0 0 32776 
32769 0 0 0 0 0 0 32772 
32769 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 4 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 2 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 1 
s 1 1 1 5
d
32768
32769 0 0 0 0 0 0 32777 
32769 0 0 0 0 0 0 32773 
32769 0 0 0 0 0 0 0 
32769 0 0 0 0 0 0 32770 
0 0 0 0 0 0 0 5 
0 0 0 0 0 0 0 4 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 2 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 1 
s 31 31 31
d
32768
32769 0 0 0 0 0 0 0 
32769 0 0 0 0 0 0 32773 
32769 0 0 0 0 0 0 0 
32769 0 0 0 0 0 0 32770 
0 0 0 0 0 0 0 5 
0 0 0 0 0 0 0 4 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 32769 
0 0 0 0 0 0 0 2 
g 31 31 31
0 4

And even when the X Y and Z coordinates differ significantly:

s 1 15 31 7888
s 30 16 2 886
d
32768
0 0 0 32769 32773 0 0 0 
0 32769 0 0 0 0 0 0 
0 32769 0 0 0 0 0 0 
0 0 0 0 0 32769 0 0 
886 0 0 0 0 0 0 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 0 7888 
g 1 15 31
7888 0
g 30 16 2
886 0
g 0 15 31
0 0
s 24 8 18 70
d
32768
0 0 0 32769 32773 32777 0 0 
0 32769 0 0 0 0 0 0 
0 32769 0 0 0 0 0 0 
0 0 0 0 0 32769 0 0 
886 0 0 0 0 0 0 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 0 7888 
0 0 0 32769 0 0 0 0 
32769 0 0 0 0 0 0 0 
0 0 0 0 32769 0 0 0 
70 0 0 0 0 0 0 0 
s 30 16 2
d
32768
0 0 0 0 32769 32773 0 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 32769 0 
0 0 0 0 0 0 0 7888 
0 0 0 32769 0 0 0 0 
32769 0 0 0 0 0 0 0 
0 0 0 0 32769 0 0 0 
70 0 0 0 0 0 0 0 
g 30 16 2
0 4
g 24 8 18
70 0
g 1 15 31
7888 0

It is really supposed to behave like a 3-dimentional array. Set a value at certain coordinates and retrieve the same value if you query those same coordinates. The difference being that in the octree, any power-of-2 aligned region of all the same type are only stored as a single entry.

Answer 1

In the time between now and when I asked the question, there are a few improvements I have made to the code.

Magic numbers

I have decided to use these instead of using the magic numbers in the code:

#define FLAG_MASK 0x8000
#define DATA_MASK 0x7FFF
#define MAX_LEVEL 5
#define node_to_set(ptr) (Node8 *)((uintptr_t)(ptr)&-sizeof(Node8))

Duplicate ancestry updating code

It turns out there are are 2 ways to obtain the correct bounds for the second memmove(). One is to find the start of the sets to delete as the destination and add state.level for the source (the start of the sets to delete is the last set to be all the same type, or the one encountered before the first set with not all nodes of the same type).

The other is to find the start of the next sets after the deleted region as the source, and subtract state.level to get the destination. The latter means no previous set needs to be recorded in the first loop in the delete branch, but also means that the 'first index that needs to be updated' needs to be found via the loop that updates the indices.

As such, the duplication can be resolved by creating this function:

static unsigned octree_offset(Octree *const octree, const State *const state, const int delta) {
    unsigned set = 0, level = state->level;
    do {
        unsigned offset = state->offset[level];
        while (++offset&7) {
            uint16_t *const node_ptr = (uint16_t *)octree->set+offset;
            if (*node_ptr&FLAG_MASK) {
                if (!set)
                    set = (offset>>3)+(*node_ptr&DATA_MASK);
                *node_ptr += delta;
            }
        }
    } while (++level != MAX_LEVEL);
    if (set) {
        Node8 *const set_ptr = octree->set+set;
        memmove(set_ptr+delta, set_ptr, (octree->set_size-set)*sizeof(Node8));
        return set;
    }
    return octree->set_size;
}

The entire else branch containing the first do while in the old code was replaced with:

else {
    const unsigned
        offset = state.offset[state.level],
        set = octree_offset(octree, &state, state.level);
    ((uint16_t *)octree->set)[offset] = set-(offset>>3)|FLAG_MASK;
    set_ptr += set;
}
octree->set_size = new_size;

However, this requires the old octree->set_size so without the prev_size variable, I created a new_size variable and simply set octree->set_size to new_size afterwards.

const unsigned new_size = octree->set_size+state.level;

This also means the condition in the reallocation while () needs to be updated too:

while ((alloc_size = octree_alloc(alloc)) < new_size)

At this point, the prev_size variable is not needed and was removed entirely.

Duplicate reallocation code/special case

I was worried about duplicating the realloc call and subsequent check in the 'special case' that all the represented nodes are identical and thus no sets need to be allocated at all. While breaking the loop instead of invoking the second case would result in some number of redundant tests and unnecessary calculation, I have decided to simply goto the relevant code in the case of the special case.

if (++state.level == MAX_LEVEL) {
    *octree = (struct Octree){.base = new};
    alloc_size = 0;
    goto dealloc;
}

The declaration of alloc_size had to be moved to the start of the main else branch.

if (octree->set_alloc != alloc) {
    octree->set_alloc = alloc;
dealloc:
    {
        Octree *const prev_octree = realloc(octree, alloc_size*sizeof(Node8)+sizeof(Octree));
        if (prev_octree)
            octree = prev_octree;
    }
}

Subtle bug

I discovered a subtle subtle bug, that would be a big vulnerability if this code had to be secure against attacks (which it currently does not). But vulnerabilities are still bad practice, so I went ahead and fixed it. If the set function is called with a new where new&FLAG_MASK is true, with identical values for all 8 nodes in a bottom level branch, they will be merged by the set function resulting in a non-bottom level node with its index flag bit set and an arbitrary data, which is interpreted as an arbitrary array index if the flag bit is set in a non-bottom level branch. This could potentially have catastrophic consequences. Values bigger than 32767 are only supported at the bottom level where the flag bit is irrelevant, but I failed to check that this bit was not set before attempting to merge. This was fixed here:

uint16_t *node_ptr = (uint16_t *)octree->set+*state.offset;
if (new&FLAG_MASK)
    *node_ptr = new;
else {
    // All the code that was originally in the main `else` branch
}

Index of octree->base

Currently all the indices contained in nodes are relative to the start of the set they are in. This means that the set referenced by any given node can be obtained using only a pointer to the node. (Unlike with absolute indices, which also require a pointer to the start of the entire octree). All, that is, except octree->base which will always refer to the same place, but to stay consistent with the others, should be able to be handled the same way, i.e. node_to_set(node_ptr)+(*node_ptr&DATA_MASK) should work with node_ptr == &octree->base. To make this work, when the base node is split (and is assigned to an index) the index should be 0x8001 instead of 0x8000:

if (state.level == 5)
    octree->base = 0x8001;
else // ...

But doing this breaks the octree_get() function, which already treats octree->base as being in the first set even though conceptually it is not. The fix to this is easy too, which is to initialize the set variable in the octree_get() function to -1 instead of 0 (and make it signed).

Although this capability is not used in any of the provided code, it becomes useful down the road.

Answer 2

Make your code more readable

Same as in the review of your previous version: the code is very compact, making it hard to read. Use whitespace to separate functions and to demarcate the various parts within a function. Also ensure you are consistent, particularly with whitespace around operators.

Also try to add more names to things. For example, what does that constant 4681 in octree_alloc() mean? Just make it a named constant, where the name can explain what it is. You can also give names to blocks of code by putting them into separate functions. This would also greatly help to reduce the size and number of indentation levels of octree_set().

Avoid cramming so much into single statements. Instead of:

node = ((uint16_t *)octree->set)[state->offset[state->level] = (set += node&0x7FFF)<<3|octree_index(x, z, y, state->level)];

Write:

static const index_mask = 0x7FFF;
uint16_t *octree_sets = (uint16_t *)octree->set;

set += node & index_mask;
unsigned child_node_offset = set << 3 | octree_index(x, z, y, state->level);
state->offset[state->level] = child_node_offset;
node = octree_sets[child_node_offset];

This already makes things a lot clearer; I don't have to entangle what gets set anymore, some constants are now explained, I don't have to scroll my editor window far to the right to see the whole line, and so on.

Readable code is automatically perceived to be less complex, it is more maintainable, and it is easier to spot and fix bugs in it. Code density does not improve performance in any way.

Add a test suite

You are worried about stability, which is good. While a manual code review can help with that, you can also increase the confidence in your code's stability by writing a proper test suite that excercises your code. Combine this with static analysis (using Clang's and GCC's analyzers, tools like cppcheck), runtime analysis (using sanitizers, tools like Valgrind), and code coverage analysis to ensure your test suite actually tests all your code.

Avoiding code duplication

Remember that compilers optimize code for you. They are extremely good at finding useless code and removing it altogether. So if you need the same thing twice, but in the first case you also want to update some variable that you don't need in the second case, reuse that code anyway, and just ignore the unused variable (you might need to cast it to void, annotate it with [[maybe_unused]] or do something similar to avoid compiler warnings).

Trust that the compiler will optimize this for you. If you don't have that trust yet, consider using a tool like the Compiler Explorer to see what kind of assembly output compilers actually generate, and experiment with the different optimization levels.

Uninitialized variable warnings

It depends on the compiler, but they are quite good at evaluating all the possible code paths. If they warn about uninitialized variables, I would take that seriously. Are you absolutely sure that there is no way that next_alloc_size >= octree->set_size could be false in the first iteration of that while-loop?

You could try to reorganize the code such that the uninitialized variable situation does not occur, without having to force a dummy initialization which might itself hide errors. Maybe you could even calculate alloc without having to use a loop?

Handling special cases

It's hard to tell whether handling special cases is worth it; that depends on how it is used. Is it common that all nodes will be set to the same value? If this happens very infrequently, you could omit this case. Shrinking a memory allocation should be a very fast operation, and it might free up memory for other things. On the other hand, if the octree needs to be grown later back to its original size, then this might suddenly be much more expensive.

You could take a hint from how C++'s std::vector is implemented, and keep track both of the amount of memory actually used, and how much memory is reserved. Then you can decide to never shrink, and only call realloc() if you need more memory than you have reserved. This avoids potentially wasteful calls to realloc() in general. You can even add an octree_reserve() function to reserve memory up front, in case the caller knows roughly how much memory is going to be needed.