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.