split bounding volume hierarchy

sbvh.cpp

#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);
}

sbvh.hpp

#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;
};