Skip to content

Instantly share code, notes, and snippets.

@jamornsriwasansak

jamornsriwasansak/sbvh.cpp Secret

Last active Jun 20, 2021
Embed
What would you like to do?
split bounding volume hierarchy
#pragma once
#include <vector>
#include "common\reflectcuts.h"
#include "common\shape.h"
#include "common\accel.h"
#include "math\math.h"
#include "math\aabb.h"
#include "accels\linearbvh.h"
class SplitBvh : public Accel
{
public:
const float AlphaSplit = 1e-5f;
SplitBvh() {};
bool intersect(Intersection * isectPtr, const Ray & r) const override;
bool intersectP(const Ray & r) const override;
void build(const std::vector<shared_ptr<Shape>> & hitablePtrs) override;
inline Aabb computeBbox() const override { return this->bound; }
friend class Qbvh4;
private:
bool intersectBranchSimple(Intersection * isectPtr, const Ray & r, const size_t index) const;
bool intersectBranch(Intersection * isectPtr, const Ray & r, const size_t index) const;
bool intersectBranchP(const Ray & r, const size_t index) const;
int buildBranch(std::vector<std::pair<size_t, Aabb>> &hitablePairs, const Aabb & bound);
Aabb bound;
std::vector<shared_ptr<Shape>> _mShapePtrs;
std::vector<LinearBvhNode> _mBvhNodes;
};
#include "sbvh.h"
#include "common\intersection.h"
bool SplitBvh::intersect(Intersection * isectPtr, const Ray & r) const
{
isectPtr->t = r.tmax;
if (!bound.intersectP(r))
{
return false;
}
return this->intersectBranchSimple(isectPtr, r, 0);
}
bool SplitBvh::intersectP(const Ray & r) const
{
if (!bound.intersectP(r))
{
return false;
}
return this->intersectBranchP(r, 0);
}
bool SplitBvh::intersectBranchSimple(Intersection * isectPtr, const Ray & r, const size_t index) const
{
const LinearBvhNode & currentNode = this->_mBvhNodes[index];
bool isIntersect = false;
for (size_t i = 0;i < 2;i++)
{
size_t childIndex = i;
if (currentNode.child[childIndex] != 0 && currentNode.bbox[i].intersectP(r))
{
// check if leaf
if (currentNode.child[childIndex] < 0)
{
int32_t shapeIndex = -currentNode.child[childIndex] - 1;
bool isHit = this->_mShapePtrs[shapeIndex]->intersect(isectPtr, Ray(r.origin, r.direction, r.tmin, isectPtr->t));
isIntersect = isHit || isIntersect;
}
else if (currentNode.child[childIndex] > 0)
{
bool isHit = intersectBranchSimple(isectPtr, Ray(r.origin, r.direction, r.tmin, isectPtr->t), currentNode.child[childIndex]);
isIntersect = isHit || isIntersect;
}
}
}
return isIntersect;
}
bool SplitBvh::intersectBranch(Intersection * isectPtr, const Ray & r, const size_t index) const
{
const LinearBvhNode & currentNode = this->_mBvhNodes[index];
bool isIntersect = false;
Float tNear[2];
bool isHit[2];
size_t childIndices[2];
childIndices[0] = 0;
childIndices[1] = 1;
for (size_t i = 0;i < 2;i++)
{
Float tFar;
if (currentNode.bbox[i].intersect(&tNear[i], &tFar, r))
{
tNear[i] = r.tmax;
isHit[i] = true;
}
else
{
isHit[i] = false;
}
}
if (tNear[0] > tNear[1])
{
std::swap(isHit[0], isHit[1]);
std::swap(childIndices[0], childIndices[1]);
}
for (size_t i = 0;i < 2;i++)
{
size_t childIndex = childIndices[i];
if (isHit[i])
{
// check if leaf
if (currentNode.child[childIndex] < 0)
{
int32_t shapeIndex = -currentNode.child[childIndex] - 1;
bool isHit = this->_mShapePtrs[shapeIndex]->intersect(isectPtr, Ray(r.origin, r.direction, r.tmin, isectPtr->t));
isIntersect = isHit || isIntersect;
}
else if (currentNode.child[childIndex] > 0)
{
bool isHit = intersectBranch(isectPtr, Ray(r.origin, r.direction, r.tmin, isectPtr->t), currentNode.child[childIndex]);
isIntersect = isHit || isIntersect;
}
}
}
return isIntersect;
}
bool SplitBvh::intersectBranchP(const Ray & r, const size_t index) const
{
const LinearBvhNode & currentNode = this->_mBvhNodes[index];
for (size_t i = 0;i < 2;i++)
{
if (currentNode.bbox[i].intersectP(r))
{
if (currentNode.child[i] < 0)
{
int32_t shapeIndex = -currentNode.child[i] - 1;
if (this->_mShapePtrs[shapeIndex]->intersectP(r))
{
return true;
}
}
else if (currentNode.child[i] > 0)
{
if (intersectBranchP(r, currentNode.child[i]))
{
return true;
}
}
}
}
return false;
}
void SplitBvh::build(const std::vector<shared_ptr<Shape>>& shapePtrs)
{
// refine all shape
std::vector<shared_ptr<Shape>> refinedShapePtrs;
for (size_t i = 0;i < shapePtrs.size();i++)
{
if (shapePtrs[i]->canIntersect())
{
refinedShapePtrs.push_back(shapePtrs[i]);
}
else
{
shapePtrs[i]->refine(&refinedShapePtrs);
}
}
// construct shape pairs
size_t numShapes = refinedShapePtrs.size();
this->_mShapePtrs = refinedShapePtrs;
std::vector<std::pair<size_t, Aabb>> shapePairs(refinedShapePtrs.size()); // shape are identified with index
Aabb sumBound;
for (size_t i = 0;i < refinedShapePtrs.size();i++)
{
Aabb bound = refinedShapePtrs[i]->computeBbox();
sumBound = Aabb::Union(sumBound, bound);
shapePairs[i] = std::pair<size_t, Aabb>(i, bound);
}
this->bound = sumBound;
this->buildBranch(shapePairs, sumBound);
}
int SplitBvh::buildBranch(std::vector<std::pair<size_t, Aabb>> &shapePairs, const Aabb & bound)
{
std::vector<shared_ptr<Shape>> & shapePtrs = this->_mShapePtrs;
std::vector<LinearBvhNode> & linearBvhNodes = this->_mBvhNodes;
LinearBvhNode result;
const size_t numShapes = shapePairs.size();
const size_t maxExtentDim = bound.maxExtent();
// handle a very rare special case where numShapes = 0 / 1
if (numShapes == 0)
{
result.child[0] = 0;
result.child[1] = 0;
size_t currentNodeIndex = this->_mBvhNodes.size();
this->_mBvhNodes.push_back(result);
return static_cast<int32_t>(currentNodeIndex + 1);
}
else if (numShapes == 1)
{
result.bbox[0] = shapePairs[0].second;
result.child[0] = -(int32_t(shapePairs[0].first) + 1);
result.child[1] = 0; // means it doesn't have child
size_t currentNodeIndex = this->_mBvhNodes.size();
this->_mBvhNodes.push_back(result);
return static_cast<int32_t>(currentNodeIndex + 1);
}
else if (numShapes == 2)
{
result.bbox[0] = shapePairs[0].second;
result.child[0] = -(int32_t(shapePairs[0].first) + 1);
result.bbox[1] = shapePairs[1].second;
result.child[1] = -(int32_t(shapePairs[1].first) + 1);;
size_t currentNodeIndex = this->_mBvhNodes.size();
this->_mBvhNodes.push_back(result);
return static_cast<int32_t>(currentNodeIndex + 1);
}
// find object split
struct ObjectSplitInfo
{
Float cost = std::numeric_limits<Float>::infinity();
size_t index = 0;
Aabb leftBound;
Aabb rightBound;
} minOs; // min object split info
{
std::sort(shapePairs.begin(), shapePairs.end(),
[&] (const std::pair<size_t, Aabb> & a, const std::pair<size_t, Aabb> & b) {
return a.second.pMin[maxExtentDim] < b.second.pMin[maxExtentDim];
});
// precompute area
std::vector<Float> leftArea = std::vector<Float>(numShapes),
rightArea = std::vector<Float>(numShapes);
{
Aabb leftBBoxTmp;
for (size_t i = 0;i < numShapes;i++)
{
leftBBoxTmp = Aabb::Union(leftBBoxTmp, shapePairs[i].second);
leftArea[i] = leftBBoxTmp.surfaceArea();
}
Aabb rightBBoxTmp;
for (int i = static_cast<int>(numShapes - 1);i >= 0;i--)
{
rightBBoxTmp = Aabb::Union(rightBBoxTmp, shapePairs[i].second);
rightArea[i] = rightBBoxTmp.surfaceArea();
}
}
// find index where it has minSahCost
for (size_t i = 0;i < numShapes - 1;i++)
{
Float sahCost = (i + 1) * leftArea[i] + (numShapes - i - 1) * rightArea[i + 1];
if (sahCost < minOs.cost)
{
minOs.cost = sahCost;
minOs.index = i + 1;
}
}
std::sort(shapePairs.begin(), shapePairs.end(),
[&] (const std::pair<size_t, Aabb> & a, const std::pair<size_t, Aabb> & b) {
return a.second.pMin[maxExtentDim] < b.second.pMin[maxExtentDim];
});
for (size_t i = 0;i < minOs.index;i++)
{
minOs.leftBound = Aabb::Union(minOs.leftBound, shapePairs[i].second);
}
for (size_t i = minOs.index;i < shapePairs.size();i++)
{
minOs.rightBound = Aabb::Union(minOs.rightBound, shapePairs[i].second);
}
} // end object split computing
// compute lambda
//const Float Alpha = 1;
Float lambda = Aabb::Intersect(minOs.leftBound, minOs.rightBound).surfaceArea() / bound.surfaceArea();
bool shouldTrySplit = lambda > SplitBvh::AlphaSplit // SBVH threshold
&& (bound.surfaceArea() >= 0.000001) // prevent bad mesh
&& (numShapes > 128); // prevent doing spatial when number of nodes is not too many
struct SpatialSplitInfo
{
Float cost = std::numeric_limits<Float>::infinity();
Float split = 0;
uint8_t dim = 4;
size_t leftCount = 0;
size_t rightCount = 0;
Aabb leftBound;
Aabb rightBound;
} minSs; // min spatial split
const size_t numBins = 8;
struct BinInfo
{
size_t numEnter = 0;
size_t numExit = 0;
Aabb bound;
};
// find spatial split
//lambda = 0.0f;
if (shouldTrySplit)
{
for (uint8_t dim = 0;dim < 3;dim++)
{
// compute chopped binned
BinInfo bins[numBins];
// compute all split plane and its bound
// (splitplane 0) (binbound 0) (splitplane 1) (binbound 1) ...
Float splitPlanes[numBins + 1];
Aabb binBounds[numBins];
// compute all split planes
for (size_t ibin = 0;ibin <= numBins;ibin++)
{
splitPlanes[ibin] = (Float(ibin) / Float(numBins)) * (bound.pMax[dim] - bound.pMin[dim]) + bound.pMin[dim];
}
for (size_t ibin = 0;ibin < numBins;ibin++)
{
binBounds[ibin] = bound;
binBounds[ibin].pMin[dim] = splitPlanes[ibin];
binBounds[ibin].pMax[dim] = splitPlanes[ibin + 1];
}
for (const auto shapePair : shapePairs)
{
Shape &t = *shapePtrs[shapePair.first];
size_t entryBin = 0, exitBin = 0;
bool isClip = false;
for (size_t ibin = 0;ibin < numBins;ibin++)
{
Aabb choppedBbox;
Float leftPlane = std::max(shapePair.second.pMin[dim], splitPlanes[ibin]);
Float rightPlane = std::min(shapePair.second.pMax[dim], splitPlanes[ibin + 1]);
if (leftPlane <= rightPlane && t.clipAabb(&choppedBbox, leftPlane, rightPlane, dim))
{
choppedBbox = Aabb::Intersect(choppedBbox, binBounds[ibin]);
bins[ibin].bound = Aabb::Union(bins[ibin].bound, choppedBbox);
if (!isClip)
{
entryBin = ibin;
isClip = true;
}
exitBin = ibin;
}
}
if (isClip)
{
bins[entryBin].numEnter++;
bins[exitBin].numExit++;
}
}
// grow bound from both side
Aabb leftBboxes[numBins];
Aabb rightBboxes[numBins];
size_t leftCounts[numBins];
size_t rightCounts[numBins];
{
size_t leftCount = 0;
Aabb leftBbox;
for (size_t ibin = 0;ibin < numBins;ibin++)
{
leftBbox = Aabb::Union(leftBbox, bins[ibin].bound);
leftBboxes[ibin] = leftBbox;
leftCount += bins[ibin].numEnter;
leftCounts[ibin] = leftCount;
}
size_t rightCount = 0;
Aabb rightBbox;
for (int ibin = numBins - 1;ibin >= 0;ibin--)
{
rightBbox = Aabb::Union(rightBbox, bins[ibin].bound);
rightBboxes[ibin] = rightBbox;
rightCount += bins[ibin].numExit;
rightCounts[ibin] = rightCount;
}
}
for (size_t i = 0;i < numBins - 1;i++)
{
Float sahCost = leftCounts[i] * leftBboxes[i].surfaceArea() + rightCounts[i + 1] * rightBboxes[i + 1].surfaceArea();
if (sahCost < minSs.cost)
{
minSs.cost = sahCost;
minSs.split = splitPlanes[i + 1];
minSs.dim = dim;
minSs.leftCount = leftCounts[i];
minSs.rightCount = rightCounts[i + 1];
minSs.leftBound = leftBboxes[i];
minSs.rightBound = rightBboxes[i + 1];
}
}
} // end each dimension for spatial split
} // end spatial split computing
std::vector<std::pair<size_t, Aabb>> leftPairs;
std::vector<std::pair<size_t, Aabb>> rightPairs;
Aabb leftBound;
Aabb rightBound;
bool useObjectSplitOnly = minOs.cost < minSs.cost ||
minSs.leftCount == 0 ||
minSs.rightCount == 0 ||
!shouldTrySplit;
if (!useObjectSplitOnly)
{
Float boundLeftPlane = bound.pMin[minSs.dim];
Float boundMidPlane = minSs.split;
Float boundRightPlane = bound.pMax[minSs.dim];
// partition to left side
const auto leftEnd = std::partition(shapePairs.begin(), shapePairs.end(), [&] (const std::pair<size_t, Aabb> & bound) {
return bound.second.pMax[minSs.dim] < minSs.split;
});
leftPairs = std::vector<std::pair<size_t, Aabb>>(shapePairs.begin(), leftEnd);
for (std::vector<std::pair<size_t, Aabb>>::const_iterator it = shapePairs.begin();it != leftEnd;it++)
{
leftBound = Aabb::Union(leftBound, it->second);
}
// partition to right side
const auto rightBegin = std::partition(leftEnd, shapePairs.end(), [&] (const std::pair<size_t, Aabb> & bound) {
return bound.second.pMin[minSs.dim] < minSs.split;
});
rightPairs = std::vector<std::pair<size_t, Aabb>>(rightBegin, shapePairs.end());
for (std::vector<std::pair<size_t, Aabb>>::const_iterator it = rightBegin;it != shapePairs.end();it++)
{
rightBound = Aabb::Union(rightBound, it->second);
}
size_t leftCount = minSs.leftCount;
size_t rightCount = minSs.rightCount;
Aabb estimatedLeftBound = minSs.leftBound;
Aabb estimatedRightBound = minSs.rightBound;
// decide the rest whether we do split or not
size_t numUndecided = rightBegin - leftEnd;
for (std::vector<std::pair<size_t, Aabb>>::const_iterator it = leftEnd;it != rightBegin;it++)
{
// compute bound delta
const Aabb & boundDelta = it->second;
Float costSplit = leftCount * estimatedLeftBound.surfaceArea() + rightCount * estimatedRightBound.surfaceArea();
Float cost1 = leftCount * Aabb::Union(estimatedLeftBound, boundDelta).surfaceArea() + (rightCount - 1) * estimatedRightBound.surfaceArea();
Float cost2 = (leftCount - 1) * estimatedLeftBound.surfaceArea() + rightCount * Aabb::Union(estimatedRightBound, boundDelta).surfaceArea();
if (cost1 <= costSplit)
{
leftPairs.push_back(*it);
estimatedLeftBound = Aabb::Union(estimatedLeftBound, boundDelta);
leftBound = Aabb::Union(leftBound, boundDelta);
rightCount--;
}
else if (cost2 <= costSplit)
{
rightPairs.push_back(*it);
estimatedRightBound = Aabb::Union(estimatedRightBound, boundDelta);
rightBound = Aabb::Union(rightBound, boundDelta);
leftCount--;
}
else
{
// split current shape in to leftbounddelta and rightbounddelta
Aabb leftBoundDelta;
{
Float leftPlane = std::max(it->second.pMin[minSs.dim], boundLeftPlane);
Float midPlane = std::min(it->second.pMax[minSs.dim], boundMidPlane);
bool isClip = shapePtrs[it->first]->clipAabb(&leftBoundDelta, leftPlane, midPlane, minSs.dim);
leftBoundDelta = Aabb::Intersect(leftBoundDelta, it->second);
}
Aabb rightBoundDelta;
{
Float midPlane = std::max(it->second.pMin[minSs.dim], boundMidPlane);
Float rightPlane = std::min(it->second.pMax[minSs.dim], boundRightPlane);
bool isClip = shapePtrs[it->first]->clipAabb(&rightBoundDelta, midPlane, rightPlane, minSs.dim);
rightBoundDelta = Aabb::Intersect(rightBoundDelta, it->second);
}
leftPairs.push_back(std::pair<size_t, Aabb>(it->first, leftBoundDelta));
rightPairs.push_back(std::pair<size_t, Aabb>(it->first, rightBoundDelta));
leftBound = Aabb::Union(leftBound, leftBoundDelta);
rightBound = Aabb::Union(rightBound, rightBoundDelta);
}
}
}
if (useObjectSplitOnly || leftPairs.size() == 0 || rightPairs.size() == 0)
{
leftPairs.clear();
rightPairs.clear();
leftPairs = std::vector<std::pair<size_t, Aabb>>(shapePairs.begin(), shapePairs.begin() + minOs.index);
leftBound = minOs.leftBound;
rightPairs = std::vector<std::pair<size_t, Aabb>>(shapePairs.begin() + minOs.index, shapePairs.end());
rightBound = minOs.rightBound;
}
result.bbox[0] = leftBound;
result.bbox[1] = rightBound;
size_t currentNodeIndex = this->_mBvhNodes.size();
this->_mBvhNodes.push_back(result);
size_t resultIndex = currentNodeIndex + 1;
if (leftPairs.size() == 1)
{
this->_mBvhNodes[currentNodeIndex].child[0] = -static_cast<int32_t>(leftPairs[0].first + 1);
}
else
{
this->_mBvhNodes[currentNodeIndex].child[0] = static_cast<int32_t>(resultIndex);
resultIndex = this->buildBranch(leftPairs, leftBound);
}
if (rightPairs.size() == 1)
{
this->_mBvhNodes[currentNodeIndex].child[1] = -static_cast<int32_t>(rightPairs[0].first + 1);
}
else
{
this->_mBvhNodes[currentNodeIndex].child[1] = static_cast<int32_t>(resultIndex);
resultIndex = this->buildBranch(rightPairs, rightBound);
}
return static_cast<int32_t>(resultIndex);
}
#pragma once
#include <vector>
#include "common\reflectcuts.h"
#include "common\shape.h"
#include "common\accel.h"
#include "math\math.h"
#include "math\aabb.h"
#include "accels\linearbvh.h"
class SplitBvh : public Accel
{
public:
const float AlphaSplit = 1e-5f;
SplitBvh() {};
bool intersect(Intersection * isectPtr, const Ray & r) const override;
bool intersectP(const Ray & r) const override;
void build(const std::vector<shared_ptr<Shape>> & hitablePtrs) override;
inline Aabb computeBbox() const override { return this->bound; }
friend class Qbvh4;
private:
bool intersectBranchSimple(Intersection * isectPtr, const Ray & r, const size_t index) const;
bool intersectBranch(Intersection * isectPtr, const Ray & r, const size_t index) const;
bool intersectBranchP(const Ray & r, const size_t index) const;
int buildBranch(std::vector<std::pair<size_t, Aabb>> &hitablePairs, const Aabb & bound);
Aabb bound;
std::vector<shared_ptr<Shape>> _mShapePtrs;
std::vector<LinearBvhNode> _mBvhNodes;
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment