八叉树剖分和按可见性遍历

Tree.cpp

#include "stdafx.h"
#include "Tree.h"
#include <libmorton/morton.h>

Tree::Tree() // sharing face, edge, vertex
	:order{ { 0, 1, 2, 4, 3, 5, 6, 7 },{ 1, 0, 3, 5, 2, 4, 7, 6 },
	{ 2, 0, 3, 6, 1, 4, 7, 5 },{ 3, 1, 2, 7, 0, 5, 6, 4 },
	{ 4, 0, 5, 6, 1, 2, 7, 3 },{ 5, 1, 4, 7, 0, 3, 6, 2 },
	{ 6, 2, 4, 7, 0, 3, 5, 1 },{ 7, 3, 5, 6, 1, 2, 4, 0 } }
{
}


Tree::~Tree()
{
}

void Tree::splat(const glm::vec3* lower, float stride, const FaceVec& faceVec, FaceArray& faceArray)
{
	// test each face against all axis-aligned cubes
	const auto& v = vertex;
	const glm::vec3* p = lower;
	const glm::vec3 deltap(stride);
	for (glm::uvec3 f : faceVec)
	{
		// Fast Parallel Surface and Solid Voxelization on GPUs
		const glm::vec3& v0 = v[f[0]].first, v1 = v[f[1]].first, v2 = v[f[2]].first;
		const glm::vec3 l = glm::min(v0, glm::min(v1, v2)), u = glm::max(v0, glm::max(v1, v2)); // AABB
		const glm::vec3 e0 = v1 - v0, e1 = v2 - v1, e2 = v0 - v2, n = glm::cross(e0, e1);
		// setup data for 3D plane overlap test
		const glm::vec3 c = (1.f - glm::step(0.f, -n)) * deltap; // critical point
		const float d0 = glm::dot(n, c - v0), d1 = glm::dot(n, deltap - c - v0);
		// setup data for 2D projection overlap test
		const glm::vec3 s = glm::step(0.f, n) * 2.f - 1.f; // 1 for 0 <= n
														   // setup data for edge functions
		glm::vec2 n_e0_xy = glm::vec2(-e0.y, e0.x) * s.z;
		glm::vec2 n_e1_xy = glm::vec2(-e1.y, e1.x) * s.z;
		glm::vec2 n_e2_xy = glm::vec2(-e2.y, e2.x) * s.z;
		float d_e0_xy = -glm::dot(n_e0_xy, glm::vec2(v0.xy)) + glm::max(0.f, deltap.x * n_e0_xy.x) + glm::max(0.f, deltap.y * n_e0_xy.y);
		float d_e1_xy = -glm::dot(n_e1_xy, glm::vec2(v1.xy)) + glm::max(0.f, deltap.x * n_e1_xy.x) + glm::max(0.f, deltap.y * n_e1_xy.y);
		float d_e2_xy = -glm::dot(n_e2_xy, glm::vec2(v2.xy)) + glm::max(0.f, deltap.x * n_e2_xy.x) + glm::max(0.f, deltap.y * n_e2_xy.y);
		glm::vec2 n_e0_zx = glm::vec2(-e0.x, e0.z) * s.y;
		glm::vec2 n_e1_zx = glm::vec2(-e1.x, e1.z) * s.y;
		glm::vec2 n_e2_zx = glm::vec2(-e2.x, e2.z) * s.y;
		float d_e0_zx = -glm::dot(n_e0_zx, glm::vec2(v0.zx)) + glm::max(0.f, deltap.z * n_e0_zx.x) + glm::max(0.f, deltap.x * n_e0_zx.y);
		float d_e1_zx = -glm::dot(n_e1_zx, glm::vec2(v1.zx)) + glm::max(0.f, deltap.z * n_e1_zx.x) + glm::max(0.f, deltap.x * n_e1_zx.y);
		float d_e2_zx = -glm::dot(n_e2_zx, glm::vec2(v2.zx)) + glm::max(0.f, deltap.z * n_e2_zx.x) + glm::max(0.f, deltap.x * n_e2_zx.y);
		glm::vec2 n_e0_yz = glm::vec2(-e0.z, e0.y) * s.x;
		glm::vec2 n_e1_yz = glm::vec2(-e1.z, e1.y) * s.x;
		glm::vec2 n_e2_yz = glm::vec2(-e2.z, e2.y) * s.x;
		float d_e0_yz = -glm::dot(n_e0_yz, glm::vec2(v0.yz)) + glm::max(0.f, deltap.y * n_e0_yz.x) + glm::max(0.f, deltap.z * n_e0_yz.y);
		float d_e1_yz = -glm::dot(n_e1_yz, glm::vec2(v1.yz)) + glm::max(0.f, deltap.y * n_e1_yz.x) + glm::max(0.f, deltap.z * n_e1_yz.y);
		float d_e2_yz = -glm::dot(n_e2_yz, glm::vec2(v2.yz)) + glm::max(0.f, deltap.y * n_e2_yz.x) + glm::max(0.f, deltap.z * n_e2_yz.y);

		for (int i = 0; i < numChildren; i++) {
			if (glm::any(glm::max(glm::lessThan(p[i] + deltap, l), glm::greaterThan(p[i], u)))) // AABB-AABB test
				continue; // no intersection if separated along an axis

			const float np = glm::dot(n, p[i]);
			if ((np + d0) * (np + d1) > 0) // 3D plane overlap test
				continue;

			if ((glm::dot(n_e0_xy, glm::vec2(p[i].xy)) + d_e0_xy >= 0) &&
				(glm::dot(n_e1_xy, glm::vec2(p[i].xy)) + d_e1_xy >= 0) &&
				(glm::dot(n_e2_xy, glm::vec2(p[i].xy)) + d_e2_xy >= 0) &&
				(glm::dot(n_e0_zx, glm::vec2(p[i].zx)) + d_e0_zx >= 0) &&
				(glm::dot(n_e1_zx, glm::vec2(p[i].zx)) + d_e1_zx >= 0) &&
				(glm::dot(n_e2_zx, glm::vec2(p[i].zx)) + d_e2_zx >= 0) &&
				(glm::dot(n_e0_yz, glm::vec2(p[i].yz)) + d_e0_yz >= 0) &&
				(glm::dot(n_e1_yz, glm::vec2(p[i].yz)) + d_e1_yz >= 0) &&
				(glm::dot(n_e2_yz, glm::vec2(p[i].yz)) + d_e2_yz >= 0)) {
				faceArray[i].push_back(f);
			}
		}
	}
}

size_t Tree::build()
{
	// depth-first traversal of [-1, 1]^3 using morton code as stack
	InnerNode root;
	build(1, 0, faceBuffer, root); // initially all the face as an outer node
	innerBuffer.push_back(root);
	outerBuffer.emplace_back(faceBuffer.size());
	mutation.reserve(outerBuffer.size() - 1);

	return innerBuffer.size();
}

void Tree::build(int depth, uint_fast64_t morton, FaceVec& faceVec, InnerNode& innerNode)
{
	InnerArray innerArray;
	FaceArray faceArray;

	// depth-1 for [0, 1]^3 -> [0, 2]^3, displace for [0, 2]^3 -> [-1, 1]^3
	const float stride = glm::pow(0.5f, depth - 1), displace = -1.f;
	// generate lower corners for all children's AABB
	glm::vec3 lower[numChildren];
	for (int i = 0; i < numChildren; i++)
	{
		uint_fast32_t x, y, z;
		morton3D_64_decode(morton + i, x, y, z);
		lower[i] = glm::vec3(x, y, z) * stride + displace;
		//glm::vec3 upper = lower[i] + stride;
		//printf("%f %f %f %f %f %f\n", lower[i].x, lower[i].y, lower[i].z, upper.x, upper.y, upper.z);
	}

	splat(lower, stride, faceVec, faceArray);
	FaceVec().swap(faceVec); // deallocating as it has been splatted into the faceArray

	// test all the triangle vector if it requires further distribution
	for (int i = 0; (i < numChildren) && (depth < maxDepth); i++) {
		if (faceArray[i].size() > batchSize) {
			build(depth + 1, (morton + i) << dimension, faceArray[i], innerArray[i]);
		}
	}

	// allocate buffer while backtracing such that siblings can be stored together
	for (int i = 0; i < numChildren; i++)
	{
		if (!faceArray[i].empty() && innerArray[i].empty()) // outer node, i.e. face vector
		{
			if (innerNode.outerBase == -1) {
				innerNode.outerBase = outerBuffer.size();
			}
			innerNode.childrenOffset[i] = -char(outerBuffer.size() - innerNode.outerBase + 1);

			//glm::mat4 xform(stride * 0.5f);
			//xform[3] = glm::vec4(lower[i] + stride * 0.5f, 1.f);
			//glm::vec4 v = xform * glm::vec4(-1.f, -1.f, -1.f, 1.f);

			outerBuffer.emplace_back(faceBuffer.size(), glm::vec4(lower[i] + stride * 0.5f, stride * 0.5f));
			faceBuffer.insert(std::end(faceBuffer), std::begin(faceArray[i]), std::end(faceArray[i]));
		}
		else if (faceArray[i].empty() && !innerArray[i].empty()) // inner node
		{
			if (innerNode.innerBase == -1) {
				innerNode.innerBase = innerBuffer.size();
			}
			innerNode.childrenOffset[i] = char(innerBuffer.size() - innerNode.innerBase + 1);
			innerBuffer.push_back(innerArray[i]);
		}
		else {
			std::logic_error("Illegal nodes during construction!\n");
		}
	}
}

int Tree::dump() const
{
	std::ofstream file("dump.obj");
#if _DEBUG
	for (auto v : vertex) {
		file << "v " << v.first.x << " " << v.first.y << " " << v.first.z << " "
			<< v.second.x << " " << v.second.y << " " << v.second.z << std::endl;
}
	file << std::endl;
#endif

	return dump(file, innerBuffer.back());
}

int Tree::dump(std::ofstream& file, const InnerNode& innerNode) const
{
	int depth = 1;
	for (int i = 0; i < numChildren; i++)
	{
		const char offset = innerNode.childrenOffset[i];
		if (offset != 0)
		{
			if (offset > 0) {
				depth = std::max(depth, 1 + dump(file, innerBuffer[innerNode.innerBase + offset - 1]));
			}
			else
			{
				const size_t j = innerNode.outerBase - offset - 1;
				size_t first = outerBuffer[j].first, last = outerBuffer[j + 1].first;
#if _DEBUG
				for (; first < last; first++) {
					const auto& f = faceBuffer[first];
					file << "f " << (f.x + 1) << " " << (f.y + 1) << " " << (f.z + 1) << std::endl;
				}

				file << std::endl;
#endif
			}
		}
	}

	return depth;
}

void Tree::traverse(glm::vec3 viewpoint)
{
	mutation.clear();
	traverse(viewpoint, 1, 0, innerBuffer.back());
}

void Tree::traverse(glm::vec3 viewpoint, int depth, uint_fast64_t morton, const InnerNode& innerNode)
{
	const float stride = glm::pow(0.5f, depth - 1);
	auto toVertex = [stride](uint_fast64_t morton) {
		const float displace = -1.f;
		uint_fast32_t x, y, z;
		morton3D_64_decode(morton, x, y, z);
		return glm::vec3(x, y, z) * stride + displace;
	};
	// the last corner is just the center of the parent cell
	const int o = octant(viewpoint, toVertex(morton + 7));

	for (int i = 0; i < numChildren; i++)
	{
		const int j = order[o][i]; // mutated order of children
		const char offset = innerNode.childrenOffset[j];
		if (offset != 0)
		{
			if (offset > 0) {
				traverse(viewpoint, depth + 1, (morton + j) << dimension, innerBuffer[innerNode.innerBase + offset - 1]);
			}
			else {
				mutation.push_back(innerNode.outerBase - offset - 1); // record view-dependent order
			}
		}
	}
}

Tree.h

#pragma once

static const int numChildren = 8, dimension = 3;
struct InnerNode {
	InnerNode() : innerBase(-1), outerBase(-1) { memset(childrenOffset, 0, sizeof(childrenOffset)); }
	char childrenOffset[numChildren];
	size_t innerBase, outerBase;
	bool empty() const { return (innerBase == -1) && (outerBase == -1); }
};
struct OuterNode
{
	OuterNode(size_t first, glm::vec4 xform = glm::vec4(0.f, 0.f, 0.f, 1.f))
		: first(first), xform(xform) {}
	size_t first; // index into the expanded face vector
	glm::vec4 xform; // bounding cube with respect to [-1,1]^3
	glm::mat4 toMat() { // transformation matrix for [-1, 1]^3 -> cube of outer node
		glm::mat4 mat(xform.w); // scale
		mat[3] = glm::vec4(xform.xyz, 1.f); // translate
		return mat;
	}
};

using InnerArray = std::array<InnerNode, numChildren>;
using FaceVec = std::vector<glm::uvec3>;
using FaceArray = std::array<FaceVec, numChildren>;

class Tree
{
public:
	Tree();
	virtual ~Tree();

	size_t build();
	// dump all the face into a wavefront obj
	int dump() const;

	// visit outer node front-to-back by depth-first traversal
	void traverse(glm::vec3 viewpoint);

	friend class Blader;
	friend class Dicer;
	friend class Piler;
protected:
	std::vector<std::pair<glm::vec3, glm::vec3>> vertex;
	// index into the above vector
	FaceVec faceBuffer;

	static const int maxDepth = 20, batchSize = 64512;
	std::vector<InnerNode> innerBuffer;
	std::vector<OuterNode> outerBuffer;

	// distribute an outer node into children according to their bounding box
	void splat(const glm::vec3* lower, float stride, const FaceVec& faceVec, FaceArray& faceArray);
	// construct an inner node according to the face
	void build(int depth, uint_fast64_t morton, FaceVec& faceVec, InnerNode& innerNode);
	// write out an inner node
	int dump(std::ofstream& file, const InnerNode& innerNode) const;

	void traverse(glm::vec3 viewpoint, int depth, uint_fast64_t morton, const InnerNode& innerNode);
	// determine the octant in which the viewpoint lies
	int octant(glm::vec3 viewpoint, glm::vec3 origin) const {
		glm::vec3 mask = glm::step(0.f, viewpoint - origin);
		return (int(mask.z) * 0x4) | (int(mask.y) * 0x2) | (int(mask.x) * 0x1);
	}
	const std::array<char, numChildren> order[numChildren];
	std::vector<size_t> mutation; // mutated order of outer nodes
};