Octree for Voxel Outlines of a Triangular Mesh

by Martin Cavarga

Problem statement:
Find uniform voxels of size centered at positions which intersect a given mesh geometry (triangle soup) .

If we are able to find volume cells intersecting a given mesh within reasonable time, we can generate a volumetric representation of the object, which leads to a point-cloud approximation of the distance field. We know that cells intersecting at least one mesh triangle are close to a small number of mesh primitives and evaluating the distance values at any point inside the cell returns a value close to zero.

Introducing: Octree

Recall that in my last post, I introduced a method BoxIntersectsAPrimitive which searches through an AABBTree generated from a mesh geometry until it finds a leaf AABBNode and then iterates through all of its primitives (e.g.: triangles) to find out whether at least one of them intersects the examined box volume. If it finds one, the test returns true. This means that any box volume can be tested for intersection with mesh primitives.

What if we start with a uniform box – a cube ? For now it can be of any axis-aligned position and scale in space, as long as it contains all of the mesh primitives. Now subdivide into 8 subcells (splitting every edge in half) if intersects the mesh. Clearly if contains all mesh triangles the subdivision condition is satisfied.

We continue this process for each subcell recursively until where is the minimum cell size (resolution), or until MAX_OCTREE_DEPTH is reached.

We have implemented a bounding volume hierarchy based on a tree with up to 8 children for each node – an Octree. When leaf conditions are met the node can be processed for further means. For example: we might want to compute an exact value of distance from the leaf cell centroid to the closest mesh triangle and store it for further use.

The overall construction algorithm can be written in pseudocode as:

Note that the distances to primitives are first compared as squared (because sqrt operation slows things down) and then we take the sqrt of the result.

Octree to Grid Transformation and the Initial Cube Cell

The advantage of octrees is that they can subdivide the surrounding space into self-similar cells. Starting with a bounding cube we can then produce cells that can be thought of as a subset of a regular grid. As a result, we can easily write values from the octree leaf cells to a scalar grid, in particular by using a transformation from real space to “grid index” space:

where is a leaf cell cube, the grid cell size, and are minimum coordinates of the scalar grid’s bounding box.

The above transformation can be used, for example, in a method of Grid object which utilizes an existing Octree:

void Grid::ApplyLeafValuesFromOctree(Octree* octree)
{
	std::vector<Box3*> boxBuffer = {};
	std::vector<double> valueBuffer = {};
	octree->GetLeafBoxesAndValues(&boxBuffer, &valueBuffer);

	size_t NLeaves = boxBuffer.size();
	uint Nx = m_Nx, Ny = m_Ny;
	
	double gMinX = m_BBox.min.x, gMinY = m_BBox.min.y, gMinZ = m_BBox.min.z;
	
	uint ix, iy, iz, gridPos;

	for (uint i = 0; i < NLeaves; i++) {
		// transform from real space to grid index space
		ix = (uint)std::floor((0.5 * (boxBuffer[i]->min.x + boxBuffer[i]->max.x) - gMinX) / m_CellSize);
		iy = (uint)std::floor((0.5 * (boxBuffer[i]->min.y + boxBuffer[i]->max.y) - gMinY) / m_CellSize);
		iz = (uint)std::floor((0.5 * (boxBuffer[i]->min.z + boxBuffer[i]->max.z) - gMinZ) / m_CellSize);

		gridPos = Nx * Ny * iz + Nx * iy + ix;
		m_Field[gridPos] = valueBuffer[i];
		m_FrozenCells[gridPos] = true; // freeze voxel mesh outline
	}
}

The Octree::GetLeafBoxesAndValues method is relatively simple. Similarly to AABB intersection queries in my previous post, it uses a stack to sort through Octree nodes and fills the respective buffers with leaf data (boxes and values).

The simplest way to build the initial cube box is to take the maximal size of the geometry’s bounding box: and expand the bounding box in the remaining directions, so that it forms a cube containing the geometry. Subdividing such volume, however, produces unsatisfactory results with shifted cell positions. The shift is proportional to the minumum cell size .

Note that the “leaf-cube-to-index-space” transformation relies on the fact that we use a regular grid with cube voxels. On top of this, we need to properly position the octree root box’s min and max points, so that:

1. The Octree‘s bounding cube really is a cube
2. The bounding cube is properly centered – i.e.: the true center of is near the center of the mesh bounding box, but coincides with the global grid with cell size .
3. the half size of satisfies: where is the expected Octree depth. This ensures that centroids of and its further subdivisions up to depth coincide with the global grid.

This can be done in the constructor as follows:

Octree::Octree(const AABBTree& aabbTree, const Box3& bbox, double minCellSize) {

        Vector3 boxCenter = bbox.getCenter();

        Vector3 gridCenterMin = Vector3(std::floor(boxCenter.x / minCellSize - 0.5), std::floor(boxCenter.y / minCellSize - 0.5), std::floor(boxCenter.z / minCellSize - 0.5));
        Vector3 gridCenterMax = Vector3(std::ceil(boxCenter.x / minCellSize - 0.5), std::ceil(boxCenter.y / minCellSize - 0.5), std::ceil(boxCenter.z / minCellSize - 0.5));
        Vector3 newGridCenter = (gridCenterMin + gridCenterMax) * minCellSize / 2.0;

        Vector3 boxSize = bbox.getSize();
        double maxDim = std::max({ boxSize.x, boxSize.y, boxSize.z });
        m_Depth = std::floor(log2(maxDim / minCellSize));
        double boxHalfDim = pow(2, m_Depth) * minCellSize;

        m_CubeBox = Box3(newGridCenter.clone().subScalar(boxHalfDim), newGridCenter.clone().addScalar(boxHalfDim));

        m_MinCellSize = minCellSize;
        m_AABBTree = aabbTree;

        m_Root = new OctreeNode(this, m_CubeBox);
    }

Final Remarks

The resulting Grid with proper bounding box position and size, and with suitable values written into the outline voxels (which are “frozen”, so they’re not affected by the upcoming computation) serves as an initial condition for an Eikonal solver I will elaborate on in my next post. The average complexities of the preceding steps – constructing an AABBTree and an Octree – are and respectively, where is the mesh complexity (number of primitives) and target scalar grid resolution.

My implementation of the full Octree construction takes approx. 0.1 seconds for grid resolution , making it mere 8% of the total computation time (AABBTree construction takes about 4% of total time). The upcoming steps for completing the computation of the signed distance field will take way more time because they need to process significantly more data (the rest of the scalar grid).

1. 

   very nice publish, i definitely love this website, keep on it

2. 

   Oooo, velmi dobre ! Dňa pi 27. 8. 2021, 16:50 Martin Cavarga napísal(a): > mshgrid posted: ” A Possible Resolution…