Code Sample from Path Tracing assignment in CMU 15-462 Computer Graphics, creating a Bounding Volume Hierarchy. By William Ozeas.

BVH.cpp

// The following functions build a Bounding Volume Hierarchy to speed up rendering
// using path tracing. We speed up collision checking by checking with the constructed 
// bounding volumes instead of primitives. We want to construct these volumes in a way 
// that minimizes a surface area heuristic. The hierarchy is organized into a tree 
// structure with leaves of max_leaf_size and constructed recursively by splitting 
// each volume 8 different ways along each axis and measuring which scores the highest 
// using the heuristic. 

// This is an excerpt of the path tracing assignment in CMU 15-462, Computer
// Graphics. The function header, BVH & BBox class, and first four lines of code were 
// provided as starter code. Code written by William Ozeas in November 2021.
// Redistribution is not permitted.

// Begin Starter Code
template<typename Primitive>
void BVH<Primitive>::build(std::vector<Primitive>&& prims, size_t max_leaf_size) {

    
    nodes.clear();
    primitives = std::move(prims);
    BBox box;
    for(const Primitive& prim : primitives) box.enclose(prim.bbox());
    // End Starter Code

    if(primitives.size() <= max_leaf_size) {
        root_idx = new_node(box, 0, primitives.size());
    } else {
        root_idx = build_recursive(max_leaf_size, box, 0, primitives.size(), 0);
    }
}

//returns index of created bounding volume
template<typename Primitive>
size_t BVH<Primitive>::build_recursive(size_t max_leaf_size, BBox box, size_t start, size_t end, int lvl) {

    //create a leaf with the start and end points instead of start and size
    auto make_leaf = [&](BBox box, size_t start1, size_t end1) {
        return new_node(box, start1, end1 - start1);
    };

    //compute the surface area heuristic
    auto compute_SAH = [&](BBox B0, int primCount0, BBox B1, int primCount1, BBox parent) {
        return (B0.surface_area()/parent.surface_area())*(float)primCount0 + (B1.surface_area()/parent.surface_area())*(float)primCount1;
    };

    const int bucketCount = 8;
    BBox bestB0, bestB1;
    float bestCoord = -1.0f;
    float bestSAH = -1.0f;
    int bestAxis = -1; //axes: x = 0, y = 1, z = 2
    size_t bestPrimCount0 = 0;
    size_t bestPrimCount1 = 0;
    
    //along each axis, one at a time, split the bounding box into several pieces. Enclose resulting partitioned
    //primitives into a new bounding box and measure heuristic, remembering the minimal split.
    for(int axis = 0; axis < 3; axis++) {
        
        float bucketSize = (box.max[axis] - box.min[axis])/((float)bucketCount);
        float coord = box.min[axis];

        for(int bucketIdx = 0; bucketIdx < bucketCount; bucketIdx++) {
            coord += bucketSize;

            //partition the primitives along coordinate coord
            auto part = std::partition(std::next(primitives.begin(), start), std::next(primitives.begin(), end),
                [&](Primitive& centroid) {
                    return centroid.bbox().center()[axis] <= coord;
                }
            );

            BBox B0, B1;
            int primCount0 = 0;
            int primCount1 = 0;

            for(auto prim = std::next(primitives.begin(), start); prim < part; prim++) {
                B0.enclose(prim->bbox());
                primCount0++;
            }

            for(auto prim = part; prim < std::next(primitives.begin(), end); prim++) {
                B1.enclose(prim->bbox());
                primCount1++;
            }

            float SAH = compute_SAH(B0, primCount0, B1, primCount1, box);
            if(bestSAH < 0 || SAH < bestSAH) {
                bestSAH = SAH;
                bestCoord = coord;
                bestAxis = axis;
                bestB0 = B0;
                bestB1 = B1;
                bestPrimCount0 = primCount0;
                bestPrimCount1 = primCount1;
            }
        } 
    }

    //re-partition according to best measured heuristic
    std::partition(std::next(primitives.begin(), start), std::next(primitives.begin(), end),
        [&](Primitive& centroid) {
            return centroid.bbox().center()[bestAxis] <= bestCoord;
        }
    ); 
    
    size_t l, r; //child bounding volumes
    if(bestPrimCount1 == 0 || bestPrimCount0 == 0) {
        return make_leaf(box, start, end);
    }

    if(bestPrimCount0 <= max_leaf_size) {
        l = make_leaf(bestB0, start, start + bestPrimCount0);
    } else {
        l = build_recursive(max_leaf_size, bestB0, start, start + bestPrimCount0, lvl + 1);
    }
    if(bestPrimCount1 <= max_leaf_size) {
        r = make_leaf(bestB1, start + bestPrimCount0, end);
    } else {
        r = build_recursive(max_leaf_size, bestB1, start + bestPrimCount0, end, lvl + 1);
    }

    return new_node(box, start, end - start, l, r);
}