BVH Tree Construction - Компиляция дает случайные ошибки

Question

Большое спасибо за помощь в дополнении. Я пытаюсь построить дерево BVH с поверхностной областью Heuristi c, но каждый раз, когда я компилирую свой код, он выдает мне случайные ошибки, такие как:

"Место чтения нарушения доступа"

"Ошибка проверки времени выполнения # 2 - стек вокруг переменной 'x' поврежден."

"Переполнение стека"

Ошибки происходят в BVH :: buildSAH ( ) функция. И я пытался найти решение на весь день, бессмысленно. Может ли это быть что-то из функции std :: partition или из отправки переменных с указателями в рекурсию?

Я читаю из книги «Физически обоснованный рендеринг: от теории к реализации Мэтта Фарра, Грега Хамфриза»

Это работает для 2 примитивов в области, но это тривиально ... Если вы хотите клонировать: https://github.com/vkaytsanov/MortonCode-BVH-KD Мой 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);

}

и 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;

}

Answer 1

Вероятно, было бы целесообразно сначала очистить ваш github. Это значит обновить материал до недавнего стандарта c ++. Кажется, что вы можете использовать c ++ 17, так что используйте его. Также, пожалуйста, посмотрите на некоторые имена. Например, «узлы» используются как переменная-член, а также имя параметра, это путаница. Пожалуйста, также инициализируйте соответствующие (все) переменные-члены.

Теперь кажется, что код в buildSAH переопределяет память. Кажется, что он может записать в конец массива ведра.

Answer 2

Исправлена ​​ошибка при запуске на CLion, где показывалось, что он выходит за границы массива buckets с целым числом, по-видимому, Visual Studio плохо описывает ошибки в рекурсивных функциях