Voxel in an octree - a c++ implementation

Introduction

A voxel is the 3D counterpart of a 2D pixel, it represents a value on a regular grid in 3D space [from Wikipedia], see below an illustration of a regular voxel grid.

An octree is a tree structure in which each internal node has exactly eight children. A node that has no children is called a leaf node. Octrees are most often used to partition a three dimensional space by recursively subdividing it into eight octants [wikipedia], see below.

Typically, each voxel is represented as the leaf node of the octree, and doesn’t have their position (their coordinates) explicitly encoded along with their values. Instead, the position of a voxel is inferred based upon its position relative to other voxels.

In this post, I’ll show a design of voxel index that is independent on which level this particular voxel is in. The implementation is from Simon Fuhrmann’s implementation of Floating Scale Surface Reconstruction. Thanks the original authors for sharing the awesome code, the git repo.

Implementation

Each voxel is indexed by a 64-bit ID, the 63rd index is the unused bit (As convention, the rightmost bit is the 0th bit). Let’s first start with the simplest case, see the node below with 8 children, the index of the eight leaf nodes (voxels) from the first to the eighth octant are: 111, 110, 100, 101, 011, 010, 000, 001. Each voxel is indexed in the order of zyx. In the current case, the level of the root node (the parent node) is 0 and the level of each leaf node is 1. So for the 64-bit voxel index, 21 bits are used to represent each coordinate of a voxel, and the maximum level of the voxel is 21, the maximum level of a node is 20. So if the index of the voxel is stored in the lth bit, the first (l-1) bits are used to store the coordinates of the direct parent node. For example, if the 21st bit stores the dimension index of the octant that the voxel is in, then the first 20 bits store the path to the parent node.

There are generally two methods that interest us:

Get the index of a certain voxel given the level of voxel, the path to the parent node and the octant.

Get the position of a certain voxel given the size of the root node, the position of the center of the root node, and the index of this certain voxel.

Index of Voxel

The index is irrelevant to the level that it’s in, i.e. the voxel shared by nodes from different levels has fixed index. Since 000 can mean either octant 7 or no node, in order to get the correct parent node index for the voxel, it’s important to have both node path and node level.

We use a temporary 64-bit variable index to store the index of a specific level l, and a uint64_t array variale to store the index for the x, y, z dimensions seperately, the x, y and z coordinates for each level are stored in element 0, 1, and 2 with higher level in higher bit.

uint64_t index = (path >> (3 * l)) & 7;
for (int i = 0; i < 3; --i)
{
    coords[i] = coords[i] << 1; // move the index for the previous level up 1 bit
    coords[i] = (coords[i] & (1 << i)) >> i;
}

After we get the index of the direct parent node, we need to retrieve the index of the voxel according to the octant. First we get the index for each dimension and store it in the 0 bit of the cooresponding array variable element. Then we move this bit up 20-level, where level is the level that parent node is in. The reason why it should 20, not 21 is because there are 21-level vacant bits, therefore the left most one is 20-lth.

for (int i = 0; i < 3; ++i)
{
    coords[i] += (octant >> i) & 1; // don't think + is necessary
    coords[i] = coords[i] << (20 - level);
}

Position of Voxel

Let’s start with the simplest case, see above the node with 8 children, and use the 2nd octant as an example, whose index is 110. The maximum index for each dimension is 1 = 2^0. So if the index for one dimension is 1, then center - size / 2 + 1 / dim_max * size = center + size / 2 is the coordinate for that dimension, otherwise, it’s center - size / 2 + 0 / dim_max * size = center - size / 2. I haven’t thought of an intuitive explanation for multiple levels yet.

math::Vec3d
VoxelIndex::compute_distance(math::Vec3d const& center, double size) const
{
    double const dim_max = static_cast<double>(1 << 20);
    double const fx = static_cast<double>(get_offset_x()) / dim_max;
    double const fy = static_cast<double>(get_offset_y()) / dim_max;
    double const fz = static_cast<double>(get_offset_z()) / dim_max;

return center - size / 2 + math::Vec3d(fx, fy, fz) * size;
}

int32_t
VoxelIndex::get_offset_x(void) const
{
    return static_cast<int32_t>(this->index & 0x1fffff);
}