Search code examples
c++treeraytracing

BVH Tree Construction - Compiling gives Random mistakes


Much thanks for the help in additionally. I'm trying to build a BVH Tree with Surface Area Heuristic, but everytime I compile my code it gives me random errors like:

"Access violation reading location"

"Run-Time Check Failure #2 - Stack around the variable 'x' was corrupted."

"Stack overflow "

The errors happen in the BVH::buildSAH() function. And I have tried to find a solution for the whole day, meaningless. Could it be something from the std::partition function or from sending variables with pointers to a recursion?

I'm reading from the book "Physically Based Rendering: From Theory to Implementation By Matt Pharr, Greg Humphreys"

It works for 2 primitives in the area, but thats trivial... If you would like to clone: https://github.com/vkaytsanov/MortonCode-BVH-KD My BVH.hpp:

#include <vector>
#include <cassert>
#include <algorithm>
#include "memory.hpp"
#include "Screen.hpp"
#include "Point3D.hpp"
#include "BoundBox.hpp"


#pragma once
enum Axis{
    X, Y, Z
};
struct MortonPrimitive{
    int primitiveIndex;
    uint32_t mortonCode;
};

struct BVHPrimitiveInfo {
    BVHPrimitiveInfo() {}
    BVHPrimitiveInfo(int primitiveNumber, const BoundBox& box) : primitiveNumber(primitiveNumber), box(box),
        centroid(Point3D(box.pMin.x* 0.5f + box.pMax.x * 0.5f, box.pMin.y* 0.5f + box.pMax.y * 0.5f, box.pMin.z* 0.5f + box.pMax.z * 0.5f)) {}

    int primitiveNumber;
    BoundBox box;
    Point3D centroid;
};

struct BVHNode {
    void InitLeaf(int first, int n, const BoundBox& b) {
        firstPrimOffset = first;
        nPrimitives = n;
        box = b;
        children[0] = children[1] = nullptr;
    }
    void InitInterior(int axis, BVHNode* c0, BVHNode* c1) {
        assert(c0 != NULL || c1 != NULL);
        children[0] = c0;
        children[1] = c1;
        this->box = Union(c0->box, c1->box);
        splitAxis = axis;
        nPrimitives = 0;
    }
    BoundBox box;
    BVHNode* children[2];
    int splitAxis, firstPrimOffset, nPrimitives;
};

struct LinearBVHNode {
    BoundBox bounds;
    union {
        int primitivesOffset;   // leaf
        int secondChildOffset;  // interior
    };
    uint16_t nPrimitives;  // 0 -> interior node
    uint8_t axis;          // interior node: xyz
    uint8_t pad[1];        // ensure 32 byte total size
};


struct BVHLittleTree {
    int startIndex;
    int numPrimitives;
    BVHNode* nodes;
    
};

struct BVH {
    BVH(std::vector<std::shared_ptr<Primitive>> p) : primitives(std::move(p)) {
        std::vector<BVHPrimitiveInfo> BVHPrimitives;
        BVHPrimitives.reserve(primitives.size());
        for (int i = 0; i < primitives.size(); i++) {
            BVHPrimitives.push_back({ i, primitives[i]->box });
        }
        MemoryArena arena(1024 * 1024);
        int totalNodes = 0;
        std::vector<std::shared_ptr<Primitive>> orderedPrimitives;
        orderedPrimitives.reserve(primitives.size());

        BVHNode* root;
        root = HLBVHBuild(arena, BVHPrimitives, &totalNodes, orderedPrimitives);
        primitives.swap(orderedPrimitives);
        BVHPrimitives.resize(0);
        printf("BVH created with %d nodes for %d "
            "primitives (%.4f MB), arena allocated %.2f MB\n",
            (int)totalNodes, (int)primitives.size(),
            float(totalNodes * sizeof(LinearBVHNode)) /
            (1024.f * 1024.f),
            float(arena.TotalAllocated()) /
            (1024.f * 1024.f));
        assert(root != NULL);
        nodes = AllocAligned<LinearBVHNode>(totalNodes);
        int offset = 0;
        flattenBVHTree(root, &offset);

    }
    ~BVH() { FreeAligned(nodes); }

    BVHNode* build(std::vector<MortonPrimitive>&, std::vector<Primitive>&);
    BVHNode* HLBVHBuild(MemoryArena& arena, const std::vector<BVHPrimitiveInfo>& BVHPrimitives, int* totalNodes, std::vector<std::shared_ptr<Primitive>>& orderedPrims);
    BVHNode* emit(BVHNode*& nodes, const std::vector<BVHPrimitiveInfo>& BVHPrimitives, MortonPrimitive* mortonPrimitives, std::vector<std::shared_ptr<Primitive>>&, int, int*, int*, int);
    BVHNode* buildSAH(MemoryArena& arena, std::vector<BVHNode*>& treeRoots, int start, int end, int* total) const;
    int flattenBVHTree(BVHNode*, int*);

    std::vector<std::shared_ptr<Primitive>> primitives;
    LinearBVHNode* nodes = nullptr;
    int maxPrimsInNode = 1;
};



inline uint32_t LeftShift3(uint32_t x) {
    if (x == (1 << 10)) --x;
    x = (x | (x << 16)) & 0b00000011000000000000000011111111;
    x = (x | (x << 8)) & 0b00000011000000001111000000001111;
    x = (x | (x << 4)) & 0b00000011000011000011000011000011;
    x = (x | (x << 2)) & 0b00001001001001001001001001001001;
    return x;
}

uint32_t EncodeMorton3(const Point3D& p) {
    return (LeftShift3(p.z) << 2) |
           (LeftShift3(p.y) << 1) |
           (LeftShift3(p.x) << 0);
}


short bitValue(uint32_t& number, uint32_t& mask) {
    return number & mask ? 1 : 0;
}


static void radixSort(std::vector<MortonPrimitive>* v)
{
    std::vector<MortonPrimitive> tempVector(v->size());
    const int bitsPerPass = 6;
    const int nBits = 30;
    static_assert((nBits % bitsPerPass) == 0,
        "Radix sort bitsPerPass must evenly divide nBits");
    const int nPasses = nBits / bitsPerPass;

    for (int pass = 0; pass < nPasses; ++pass) {
        // Perform one pass of radix sort, sorting _bitsPerPass_ bits
        int lowBit = pass * bitsPerPass;

        // Set in and out vector pointers for radix sort pass
        std::vector<MortonPrimitive>& in = (pass & 1) ? tempVector : *v;
        std::vector<MortonPrimitive>& out = (pass & 1) ? *v : tempVector;

        // Count number of zero bits in array for current radix sort bit
        const int nBuckets = 1 << bitsPerPass;
        int bucketCount[nBuckets] = { 0 };
        const int bitMask = (1 << bitsPerPass) - 1;
        for (const MortonPrimitive& mp : in) {
            int bucket = (mp.mortonCode >> lowBit) & bitMask;
            ++bucketCount[bucket];
        }

        // Compute starting index in output array for each bucket
        int outIndex[nBuckets];
        outIndex[0] = 0;
        for (int i = 1; i < nBuckets; ++i)
            outIndex[i] = outIndex[i - 1] + bucketCount[i - 1];

        // Store sorted values in output array
        for (const MortonPrimitive& mp : in) {
            int bucket = (mp.mortonCode >> lowBit) & bitMask;
            out[outIndex[bucket]++] = mp;
        }
    }
    // Copy final result from _tempVector_, if needed
    if (nPasses & 1) std::swap(*v, tempVector);
}

//BVHNode* BVH::build(std::vector<MortonPrimitive>& mortonPrimitives, std::vector<Primitive>& prims) {
//  
//  
//}

struct BucketInfo {
    int count = 0;
    BoundBox bounds;
};

BVHNode* BVH::HLBVHBuild(MemoryArena& arena, const std::vector<BVHPrimitiveInfo>& BVHPrimitives, int* totalNodes, std::vector<std::shared_ptr<Primitive>>& orderedPrims)  {
    BoundBox box;
    for (const BVHPrimitiveInfo& pi : BVHPrimitives) {
        box = box.Union(box, pi.centroid); // maybe it should be UNION @TODO
    }

    std::vector<MortonPrimitive> mortonPrims(BVHPrimitives.size());
    for (int i = 0; i < BVHPrimitives.size(); i++) {
        const int mortonBits = 10;
        const int mortonScale = 1 << mortonBits;

        mortonPrims[i].primitiveIndex = BVHPrimitives[i].primitiveNumber;
        Point3D p = box.offset(BVHPrimitives[i].centroid);
        p.x = p.x * mortonScale;
        p.y = p.y * mortonScale;
        p.z = p.z * mortonScale;
        mortonPrims[i].mortonCode = EncodeMorton3(p);
    }

    radixSort(&mortonPrims);

    //for (MortonPrimitive mp : mortonPrims) {
    //  std::cout << mp.primitiveIndex << " " << mp.mortonCode << std::endl;
    //}
    std::vector<BVHLittleTree> treesToBuild;

    uint32_t mask = 0b00111111111111000000000000000000; // first 12 bits describe the position of the primitive
    for (int start = 0, end = 1; end <= (int)mortonPrims.size(); end++) {
        if (end == mortonPrims.size() || ((mortonPrims[start].mortonCode & mask) != (mortonPrims[end].mortonCode & mask))) {
            int n = end - start;
            int maxNodes = 2 * n;
            BVHNode* nodes = arena.Alloc<BVHNode>(maxNodes, false);
            treesToBuild.push_back({ start, n, nodes });
            start = end;
        }
    }

    int orderedPrimsOffset = 0;
    orderedPrims.resize(primitives.size());
    int nodesCreated = 0;
    int firstBitIndex = 29 - 12;

    for (int i = 0; i < treesToBuild.size(); i++) {
        treesToBuild[i].nodes = BVH::emit(treesToBuild[i].nodes, BVHPrimitives, &mortonPrims[treesToBuild[i].startIndex], orderedPrims, treesToBuild[i].numPrimitives, &nodesCreated, &orderedPrimsOffset, firstBitIndex);
        *totalNodes += nodesCreated;
    }
    totalNodes += nodesCreated;
    std::vector<BVHNode*> finishedTrees;
    finishedTrees.reserve(treesToBuild.size());
    for (BVHLittleTree& tr : treesToBuild) {
        finishedTrees.emplace_back(tr.nodes);
    }
    return buildSAH(arena, finishedTrees, 0, finishedTrees.size(), totalNodes);

}

BVHNode* BVH::emit(BVHNode*& nodes, const std::vector<BVHPrimitiveInfo>& BVHPrimitive, MortonPrimitive* mortonPrimitives, std::vector<std::shared_ptr<Primitive>>& orderedPrimitives, int primitivesCount, int* totalNodes, int* orderedPrimsOffset, int bitIndex) {
    if (bitIndex == -1 || primitivesCount < maxPrimsInNode) {
        (*totalNodes)++;
        BVHNode* tmp = nodes++;
        BoundBox box;
        int firstPrimOffset = *orderedPrimsOffset;
        for (int i = 0; i < primitivesCount; i++) {
            int index = mortonPrimitives[i].primitiveIndex;
            orderedPrimitives[firstPrimOffset + i] = primitives[index];
            box = box.Union(box, BVHPrimitive[index].box);
        }
        tmp->InitLeaf(0, primitivesCount, box);
        return tmp;

    }
    else {
        int mask = 1 << bitIndex;
        if ((mortonPrimitives[0].mortonCode & mask) == (mortonPrimitives[primitivesCount - 1].mortonCode & mask)){ // Next tree if nothing to split for this bit
            return emit(nodes, BVHPrimitive, mortonPrimitives, orderedPrimitives, primitivesCount, totalNodes, orderedPrimsOffset, bitIndex - 1);
        }
        int start = 0;
        int end = primitivesCount - 1;
        while (start + 1 != end) {
            int mid = (end - start) / 2 + start; // (start-end)/2
            if ((mortonPrimitives[start].mortonCode & mask) == (mortonPrimitives[mid].mortonCode & mask)) {
                start = mid;
            }
            else {
                end = mid;
            }
        }
        int split = end;
        (*totalNodes)++;
        BVHNode* tmp = nodes++;
        BVHNode* lbvh[2];
        lbvh[0] = emit(nodes, BVHPrimitive, mortonPrimitives, orderedPrimitives, split, totalNodes, orderedPrimsOffset, bitIndex-1);
        lbvh[1] = emit(nodes, BVHPrimitive, &mortonPrimitives[split], orderedPrimitives, primitivesCount - split, totalNodes, orderedPrimsOffset, bitIndex - 1);
        int axis = bitIndex % 3;
        tmp->InitInterior(axis, lbvh[0], lbvh[1]);
        return tmp;
    }
}

BVHNode* BVH::buildSAH(MemoryArena& arena, std::vector<BVHNode*>& treeRoots, int start, int end, int* total) const {
    
    int nodesCount = end - start;
    if (nodesCount == 1) {
        return treeRoots[start];
    }
    assert(nodesCount > 1);
    (*total)++;
    
    BVHNode* node = arena.Alloc<BVHNode>();
    BoundBox box;
    for (int i = start; i < end; i++) {
        box = Union(box, treeRoots[i]->box);
    }
    BoundBox centroidBox;
    for (int i = start; i < end; i++) {
        Point3D centroid = Point3D((treeRoots[i]->box.pMin.x + treeRoots[i]->box.pMax.x) * 0.5f, (treeRoots[i]->box.pMin.y + treeRoots[i]->box.pMax.y) * 0.5f, (treeRoots[i]->box.pMin.z + treeRoots[i]->box.pMax.z) * 0.5f);
        centroidBox = Union(centroidBox, centroid);
    }
    const int dimension = centroidBox.MaximumExtent();
    const int nBuckets = 12;
    struct Buckets {
        int count = 0;
        BoundBox box;
    };
    Buckets buckets[nBuckets];
    for (int i = start; i < end; i++) {
        float centroid = (treeRoots[i]->box.pMin[dimension] * 0.5f + treeRoots[i]->box.pMax[dimension] * 0.5f) ;
        int b = nBuckets * ((centroid - centroidBox.pMin[dimension]) / (centroidBox.pMax[dimension] - centroidBox.pMin[dimension]));
        if (b == nBuckets) b = nBuckets - 1;
        //assert(b < nBuckets);
        buckets[b].count++;
        buckets[b].box = Union(buckets[b].box, treeRoots[i]->box);
    }

    float cost[nBuckets - 1];
    for (int i = 0; i < nBuckets - 1; i++) {
        BoundBox b0, b1;
        int count0 = 0, count1 = 0;
        for (int j = 0; j <= i; j++) {
            b0 = Union(b0, buckets[j].box);
            count0 += buckets[j].count; 
        }
        for (int j = i+1; j < nBuckets; j++) {
            b1 = Union(b1, buckets[j].box);
            count1 += buckets[j].count;
        }

        cost[i] = (.125f + (count0 * b0.surfaceArea() + count1 * b1.surfaceArea())) / box.surfaceArea();
    }

    double minCost = cost[0];
    int minCostSplitBucket = 0;
    for (int i = 1; i < nBuckets - 1; ++i) {
        if (cost[i] < minCost) {
            minCost = cost[i];
            minCostSplitBucket = i;
        }
    }

    BVHNode** pmid = std::partition(&treeRoots[start], &treeRoots[end - 1] + 1, [=](const BVHNode* node) {
            double centroid = (node->box.pMin[dimension]*0.5f + node->box.pMax[dimension] * 0.5f) ;
            int b = nBuckets * ((centroid - centroidBox.pMin[dimension]) / (centroidBox.pMax[dimension] - centroidBox.pMin[dimension]));
            if (b == nBuckets) b = nBuckets - 1;
            return b <= minCostSplitBucket;
        });
    assert(pmid != NULL);
    //std::cout << pmid << "  " << &treeRoots[0];
    int mid = pmid - &treeRoots[0];
    //std::cout << start << " " << mid << std::endl;
    //std::cout << mid << " " << end << std::endl;
    std::cout << dimension << std::endl;
    //assert(dimension < 3);

    node->InitInterior(dimension, this->buildSAH(arena, treeRoots, start, mid, total), this->buildSAH(arena, treeRoots, mid, end, total));
    return node;
}

int BVH::flattenBVHTree(BVHNode* node, int* offset) {
    LinearBVHNode* linearNode = &nodes[*offset];
    linearNode->bounds = node->box;
    int myOffset = (*offset)++;
    if (node->nPrimitives > 0) {
        linearNode->primitivesOffset = node->firstPrimOffset;
        linearNode->nPrimitives = node->nPrimitives;
    }
    else {
        // Create interior flattened BVH node
        linearNode->axis = node->splitAxis;
        linearNode->nPrimitives = 0;
        flattenBVHTree(node->children[0], offset);
        linearNode->secondChildOffset = flattenBVHTree(node->children[1], offset);
    }
    return myOffset;
}

My Point3D.hpp

#include <cstdint>
#pragma once
struct Point3D {
    float x;
    float y;
    float z;
    Point3D(uint32_t, uint32_t, uint32_t);
    Point3D();

    int operator[](int);
    int operator[](int) const;
    Point3D operator+(int);
    Point3D operator-(int);
    Point3D operator-(Point3D&);
};

Point3D::Point3D() {
    x = 0;
    y = 0;
    z = 0;
}

Point3D::Point3D(uint32_t x, uint32_t y, uint32_t z) {
    this->x = x;
    this->y = y;
    this->z = z;
}

bool operator<(Point3D a, Point3D b) {
    uint32_t xSquare = a.x * a.x;
    uint32_t ySquare = a.y * a.y;
    uint32_t zSquare = a.z * a.z;

    uint32_t x2Square = b.x * b.x;
    uint32_t y2Square = b.y * b.y;
    uint32_t z2Square = b.z * b.z;

    int64_t sum = std::sqrt(xSquare + ySquare + z2Square) - std::sqrt(x2Square + y2Square + z2Square);
    return sum < 0 ||
        sum == 0 && xSquare < x2Square ||
        sum == 0 && xSquare == x2Square && ySquare < y2Square ||
        sum == 0 && xSquare == x2Square && ySquare == y2Square && zSquare < z2Square;
}

bool operator>(Point3D a, Point3D b) {
    uint32_t xSquare = a.x * a.x;
    uint32_t ySquare = a.y * a.y;
    uint32_t zSquare = a.z * a.z;

    uint32_t x2Square = b.x * b.x;
    uint32_t y2Square = b.y * b.y;
    uint32_t z2Square = b.z * b.z;

    int32_t sum = std::sqrt(xSquare + ySquare + z2Square) - std::sqrt(x2Square + y2Square + z2Square);
    return sum > 0 ||
        sum == 0 && xSquare > x2Square ||
        sum == 0 && xSquare == x2Square && ySquare > y2Square ||
        sum == 0 && xSquare == x2Square && ySquare == y2Square && zSquare > z2Square;
}

int Point3D::operator[](int i) {
    if (i == 0) return x;
    if (i == 1) return y;
    return z;
}

Point3D Point3D::operator+(int i) {
    this->x += i;
    this->y += i;
    this->z += i;
    return *this;
}

Point3D Point3D::operator-(const int i) {
    this->x -= i;
    this->y -= i;
    this->z -= i;
    return *this;
}

Point3D Point3D::operator-(Point3D& p) {
    this->x -= p.x;
    this->y -= p.y;
    this->z -= p.z;
    return *this;
}

int Point3D::operator[](const int i) const {
    if (i == 0) return x;
    if (i == 1) return y;
    return z;
}

My BoundBox.hpp

#include "Point3D.hpp"
#include "Vector3D.hpp"
#pragma once


struct BoundBox { 
    Point3D pMin;
    Point3D pMax;

    BoundBox(Point3D);
    BoundBox(Point3D, Point3D);
    BoundBox();
    void setBounds(BoundBox);
    void Union(BoundBox);

    BoundBox Union(BoundBox&, Point3D&);
    BoundBox Union(BoundBox, BoundBox);
    BoundBox unite(BoundBox, BoundBox);
    BoundBox unite(BoundBox);

    const Point3D offset(const Point3D&);
    Point3D diagonal();
    const int MaximumExtent();
    float surfaceArea();
};


BoundBox::BoundBox() {
    float minNum = 0;
    pMin = Point3D(800, 600, 300);
    pMax = Point3D(minNum, minNum, minNum);
}
BoundBox::BoundBox(Point3D p){
    pMin = p;
    pMax = p;
}

BoundBox::BoundBox(Point3D p1, Point3D p2) {
    pMin = Point3D(std::min(p1.x, p2.x), std::min(p1.y, p2.y), std::min(p1.z, p2.z));
    pMax = Point3D(std::max(p1.x, p2.x), std::max(p1.y, p2.y), std::max(p1.z, p2.z));
        
}

BoundBox BoundBox::Union(BoundBox& box, Point3D& p) {
    BoundBox newBox;
    newBox.pMin = Point3D(std::min(box.pMin.x, p.x), std::min(box.pMin.y, p.y), std::min(box.pMin.z, p.z));
    newBox.pMax = Point3D(std::max(box.pMax.x, p.x), std::max(box.pMax.y, p.y), std::max(box.pMax.z, p.z));
    return newBox;
}

BoundBox BoundBox::Union(BoundBox box1, BoundBox box2) {
    BoundBox newBox;
    newBox.pMin = std::min(box1.pMin, box2.pMin);
    newBox.pMax = std::max(box1.pMax, box2.pMax);
    return newBox;
}

BoundBox Union(BoundBox box1, BoundBox box2) {
    BoundBox newBox;
    newBox.pMin = std::min(box1.pMin, box2.pMin);
    newBox.pMax = std::max(box1.pMax, box2.pMax);
    return newBox;
}

BoundBox BoundBox::unite(BoundBox b1, BoundBox b2) {
    bool x = (b1.pMax.x >= b2.pMin.x) && (b1.pMin.x <= b2.pMax.x);
    bool y = (b1.pMax.y >= b2.pMin.y) && (b1.pMin.y <= b2.pMax.y);
    bool z = (b1.pMax.z >= b2.pMin.z) && (b1.pMin.z <= b2.pMax.z);
    if (x && y && z) {
        return Union(b1, b2);
    }
}

BoundBox BoundBox::unite(BoundBox b2) {
    bool x = (this->pMax.x >= b2.pMin.x) && (this->pMin.x <= b2.pMax.x);
    bool y = (this->pMax.y >= b2.pMin.y) && (this->pMin.y <= b2.pMax.y);
    bool z = (this->pMax.z >= b2.pMin.z) && (this->pMin.z <= b2.pMax.z);
    if (x && y && z) {
        return Union(*this, b2);
    }
    else return *this;
}

const int BoundBox::MaximumExtent() {
    Point3D d = Point3D(this->pMax.x - this->pMin.x, this->pMax.y - this->pMin.y, this->pMax.z - this->pMin.z); // diagonal
    if (d.x > d.y && d.x > d.z) {
        return 0;
    }
    else if (d.y > d.z) {
        return 1;
    }
    else {
        return 2;
    }
}

float BoundBox::surfaceArea() {
    Point3D d = Point3D(this->pMax.x - this->pMin.x, this->pMax.y - this->pMin.y, this->pMax.z - this->pMin.z); // diagonal
    return 2 * (d.x * d.y + d.x * d.z + d.y * d.z);
}

const Point3D BoundBox::offset(const Point3D& p) {
    Point3D o = Point3D(p.x - pMin.x, p.y - pMin.y, p.z - pMin.z);
    
    if (pMax.x > pMin.x) o.x /= pMax.x - pMin.x;
    if (pMax.y > pMin.y) o.y /= pMax.y - pMin.y;
    if (pMax.z > pMin.z) o.z /= pMax.z - pMin.z;
    return o;
}

My memory.hpp

#include <list>
#include <cstddef>
#include <algorithm>
#include <malloc.h>
#include <stdlib.h> 
#pragma once
#define ARENA_ALLOC(arena, Type) new ((arena).Alloc(sizeof(Type))) Type
void* AllocAligned(size_t size);
template <typename T>
T* AllocAligned(size_t count) {
    return (T*)AllocAligned(count * sizeof(T));
}

void FreeAligned(void*);
class
#ifdef PBRT_HAVE_ALIGNAS
    alignas(PBRT_L1_CACHE_LINE_SIZE)
#endif // PBRT_HAVE_ALIGNAS
    MemoryArena {
public:
    // MemoryArena Public Methods
    MemoryArena(size_t blockSize = 262144) : blockSize(blockSize) {}
    ~MemoryArena() {
        FreeAligned(currentBlock);
        for (auto& block : usedBlocks) FreeAligned(block.second);
        for (auto& block : availableBlocks) FreeAligned(block.second);
    }
    void* Alloc(size_t nBytes) {
        // Round up _nBytes_ to minimum machine alignment
#if __GNUC__ == 4 && __GNUC_MINOR__ < 9
        // gcc bug: max_align_t wasn't in std:: until 4.9.0
        const int align = alignof(::max_align_t);
#elif !defined(PBRT_HAVE_ALIGNOF)
        const int align = 16;
#else
        const int align = alignof(std::max_align_t);
#endif
#ifdef PBRT_HAVE_CONSTEXPR
        static_assert(IsPowerOf2(align), "Minimum alignment not a power of two");
#endif
        nBytes = (nBytes + align - 1) & ~(align - 1);
        if (currentBlockPos + nBytes > currentAllocSize) {
            // Add current block to _usedBlocks_ list
            if (currentBlock) {
                usedBlocks.push_back(
                    std::make_pair(currentAllocSize, currentBlock));
                currentBlock = nullptr;
                currentAllocSize = 0;
            }

            // Get new block of memory for _MemoryArena_

            // Try to get memory block from _availableBlocks_
            for (auto iter = availableBlocks.begin();
                iter != availableBlocks.end(); ++iter) {
                if (iter->first >= nBytes) {
                    currentAllocSize = iter->first;
                    currentBlock = iter->second;
                    availableBlocks.erase(iter);
                    break;
                }
            }
            if (!currentBlock) {
                currentAllocSize = std::max(nBytes, blockSize);
                currentBlock = AllocAligned<uint8_t>(currentAllocSize);
            }
            currentBlockPos = 0;
        }
        void* ret = currentBlock + currentBlockPos;
        currentBlockPos += nBytes;
        return ret;
    }
    template <typename T>
    T* Alloc(size_t n = 1, bool runConstructor = true) {
        T* ret = (T*)Alloc(n * sizeof(T));
        if (runConstructor)
            for (size_t i = 0; i < n; ++i) new (&ret[i]) T();
        return ret;
    }
    void Reset() {
        currentBlockPos = 0;
        availableBlocks.splice(availableBlocks.begin(), usedBlocks);
    }
    size_t TotalAllocated() const {
        size_t total = currentAllocSize;
        for (const auto& alloc : usedBlocks) total += alloc.first;
        for (const auto& alloc : availableBlocks) total += alloc.first;
        return total;
    }

private:
    MemoryArena(const MemoryArena&) = delete;
    MemoryArena & operator=(const MemoryArena&) = delete;
    // MemoryArena Private Data
    const size_t blockSize;
    size_t currentBlockPos = 0, currentAllocSize = 0;
    uint8_t * currentBlock = nullptr;
    std::list<std::pair<size_t, uint8_t*>> usedBlocks, availableBlocks;
};

template <typename T, int logBlockSize>
class BlockedArray {
public:
    // BlockedArray Public Methods
    BlockedArray(int uRes, int vRes, const T* d = nullptr)
        : uRes(uRes), vRes(vRes), uBlocks(RoundUp(uRes) >> logBlockSize) {
        int nAlloc = RoundUp(uRes) * RoundUp(vRes);
        data = AllocAligned<T>(nAlloc);
        for (int i = 0; i < nAlloc; ++i) new (&data[i]) T();
        if (d)
            for (int v = 0; v < vRes; ++v)
                for (int u = 0; u < uRes; ++u) (*this)(u, v) = d[v * uRes + u];
    }
    const int BlockSize() const { return 1 << logBlockSize; }
    int RoundUp(int x) const {
        return (x + BlockSize() - 1) & ~(BlockSize() - 1);
    }
    int uSize() const { return uRes; }
    int vSize() const { return vRes; }
    ~BlockedArray() {
        for (int i = 0; i < uRes * vRes; ++i) data[i].~T();
        FreeAligned(data);
    }
    int Block(int a) const { return a >> logBlockSize; }
    int Offset(int a) const { return (a & (BlockSize() - 1)); }
    T& operator()(int u, int v) {
        int bu = Block(u), bv = Block(v);
        int ou = Offset(u), ov = Offset(v);
        int offset = BlockSize() * BlockSize() * (uBlocks * bv + bu);
        offset += BlockSize() * ov + ou;
        return data[offset];
    }
    const T & operator()(int u, int v) const {
        int bu = Block(u), bv = Block(v);
        int ou = Offset(u), ov = Offset(v);
        int offset = BlockSize() * BlockSize() * (uBlocks * bv + bu);
        offset += BlockSize() * ov + ou;
        return data[offset];
    }
    void GetLinearArray(T * a) const {
        for (int v = 0; v < vRes; ++v)
            for (int u = 0; u < uRes; ++u) * a++ = (*this)(u, v);
    }

private:
    // BlockedArray Private Data
    T * data;
    const int uRes, vRes, uBlocks;
};

void* AllocAligned(size_t size) {
    return _aligned_malloc(size, 32);
}

void FreeAligned(void* ptr) {
    if (!ptr) return;
    _aligned_free(ptr);

}

and My Source.cpp

#include <iostream>
#include <vector>
#include <chrono>

#include "Point3D.hpp"
#include "Screen.hpp"
#include "BVH.hpp"


#define N 150

int main(){
    auto startTime = std::chrono::high_resolution_clock::now();

    Screen* screen = new Screen(800, 600, 300);
    screen->generatePoints(N);
    
    
    
    //for (MortonPrimitive m : mortonPrims) {
    //  std::cout << m.mortonCode << std::endl;
    //}
    
    std::vector<std::shared_ptr<Primitive>> primitives;
    primitives.reserve(N);
    for (int i = 0; i < N; i++) {
        primitives.emplace_back(screen->castPointToPrimitive(i));
    }

    BVH test(primitives);
    auto endTime = std::chrono::high_resolution_clock::now();

    std::cout << "Time spent: " << std::chrono::duration_cast<std::chrono::milliseconds>(endTime - startTime).count() << "ms\n";
    getchar();
    delete screen;

}

Solution

  • Probably it would be wise to first cleanup your github. This mean update stuff to the recent c++ standard. It seems that you can use c++17 so use it. Also please look at some names. For example 'nodes' is used as member variable as well as parameter name, this is confusion. Please also initialize relevant (all) member variables.

    Now it seems that the code in buildSAH override memory. It seems that it it can write over the end of buckets array.