This project contains the following files (right-click files you'd like to download):
acceleration.cppteapotdata.h//[header]
// Example of Acceleration Structures for Ray-Tracing (BBox, BVH and Grid)
//[/header]
//[compile]
// Download the acceleration.cpp and teapotdata.h file to a folder.
// Open a shell/terminal, and run the following command where the files is saved:
//
// clang++ -std=c++14 -o acceleration acceleration.cpp -O3 -DACCEL_BBOX
//
// clang++ -std=c++14 -o acceleration acceleration.cpp -O3 -DACCEL_BVH
//
// clang++ -std=c++14 -o acceleration acceleration.cpp -O3 -DACCEL_GRID
//
// You can use c++ if you don't use clang++
//
// Run with: ./acceleration. Open the file ./image.png in Photoshop or any program
// reading PPM files.
//[/compile]
//[ignore]
// Copyright (C) 2016 www.scratchapixel.com
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//[/ignore]
#include <atomic>
#include <memory>
#include <cassert>
#include <vector>
#include <iostream>
#include <fstream>
#include <limits>
#include <cmath>
#include <chrono>
#include <queue>
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
const float kEpsilon = 1e-8;
const float kInfinity = std::numeric_limits<float>::max();
std::atomic<uint32_t> numPrimaryRays(0);
std::atomic<uint32_t> numRayTriangleTests(0);
std::atomic<uint32_t> numRayTriangleIntersections(0);
std::atomic<uint32_t> numRayBBoxTests(0);
std::atomic<uint32_t> numRayBoundingVolumeTests(0);
template<typename T>
class Vec3
{
public:
Vec3() : x(0), y(0), z(0) {}
Vec3(T xx) : x(xx), y(xx), z(xx) {}
Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}
Vec3 operator * (const T& r) const { return Vec3(x * r, y * r, z * r); }
Vec3 operator + (const Vec3& v) const { return Vec3(x + v.x, y + v.y, z + v.z); }
Vec3 operator - (const Vec3& v) const { return Vec3(x - v.x, y - v.y, z - v.z); }
template<typename U>
Vec3 operator / (const Vec3<U>& v) const { return Vec3(x / v.x, y / v.y, z / v.z); }
friend Vec3 operator / (const T r, const Vec3& v)
{ return Vec3(r / v.x, r / v.y, r / v.z); }
const T& operator [] (size_t i) const { return (&x)[i]; }
T& operator [] (size_t i) { return (&x)[i]; }
T length2() const{ return x * x + y * y + z * z; }
friend Vec3 operator * (const float&r, const Vec3& v)
{ return Vec3(v.x * r, v.y * r, v.z * r); }
friend std::ostream& operator << (std::ostream& os, const Vec3<T>& v)
{ os << v.x << " " << v.y << " " << v.z << std::endl; return os; }
T x, y, z;
};
template<typename T>
Vec3<T> cross(const Vec3<T>& a, const Vec3<T>& b)
{
return Vec3<T>(a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x);
}
template<typename T>
T dot(const Vec3<T>& va, const Vec3<T>& vb)
{ return va.x * vb.x + va.y * vb.y + va.z * vb.z; }
template<typename T>
void normalize(Vec3<T>& vec)
{
T len2 = vec.length2();
if (len2 > 0) {
T invLen = 1 / sqrt(len2);
vec.x *= invLen, vec.y *= invLen, vec.z *= invLen;
}
}
template<typename T>
class Matrix44
{
public:
Matrix44() { /* ... define identity matrix ... */ }
Matrix44(T m00, T m01, T m02, T m03,
T m10, T m11, T m12, T m13,
T m20, T m21, T m22, T m23,
T m30, T m31, T m32, T m33)
{
m[0][0] = m00; m[0][1] = m01; m[0][2] = m02; m[0][3] = m03;
m[1][0] = m10; m[1][1] = m11; m[1][2] = m12; m[1][3] = m13;
m[2][0] = m20; m[2][1] = m21; m[2][2] = m22; m[2][3] = m23;
m[3][0] = m30; m[3][1] = m31; m[3][2] = m32; m[3][3] = m33;
}
Matrix44 inverse() const { Matrix44 matInv = *this; return matInv; }
T* operator [] (size_t i) { return &m[i][0]; }
const T* operator [] (size_t i) const { return &m[i][0]; }
T m[4][4] = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
};
template<typename T>
void matVecMult(const Matrix44<T>& m, Vec3<T>& v)
{
Vec3<T> vt;
vt.x = v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0],
vt.y = v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1],
vt.z = v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2];
v = vt;
}
template<typename T>
void matPointMult(const Matrix44<T>& m, Vec3<T>& p)
{
Vec3<T> pt;
pt.x = p.x * m[0][0] + p.y * m[1][0] + p.z * m[2][0] + m[3][0];
pt.y = p.x * m[0][1] + p.y * m[1][1] + p.z * m[2][1] + m[3][1];
pt.z = p.x * m[0][2] + p.y * m[1][2] + p.z * m[2][2] + m[3][2];
T w = p.x * m[0][3] + p.y * m[1][3] + p.z * m[2][3] + m[3][3];
if (w != 1) {
pt.x /= w, pt.y /= w, pt.z /= w;
}
p = pt;
}
using Vec3f = Vec3<float>;
using Vec3b = Vec3<bool>;
using Vec3i = Vec3<int32_t>;
using Vec3ui = Vec3<uint32_t>;
using Matrix44f = Matrix44<float>;
template<typename T = float>
class BBox
{
public:
BBox() {}
BBox(Vec3<T> min_, Vec3<T> max_)
{
bounds[0] = min_;
bounds[1] = max_;
}
BBox& extendBy(const Vec3<T>& p)
{
if (p.x < bounds[0].x) bounds[0].x = p.x;
if (p.y < bounds[0].y) bounds[0].y = p.y;
if (p.z < bounds[0].z) bounds[0].z = p.z;
if (p.x > bounds[1].x) bounds[1].x = p.x;
if (p.y > bounds[1].y) bounds[1].y = p.y;
if (p.z > bounds[1].z) bounds[1].z = p.z;
return *this;
}
/*inline */ Vec3<T> centroid() const { return (bounds[0] + bounds[1]) * 0.5; }
Vec3<T>& operator [] (bool i) { return bounds[i]; }
const Vec3<T> operator [] (bool i) const { return bounds[i]; }
bool intersect(const Vec3<T>&, const Vec3<T>&, const Vec3b&, float&) const;
Vec3<T> bounds[2] = { kInfinity, -kInfinity };
};
template<typename T>
bool BBox<T>::intersect(const Vec3<T>& orig, const Vec3<T>& invDir, const Vec3b& sign, float& tHit) const
{
numRayBBoxTests++;
float tmin, tmax, tymin, tymax, tzmin, tzmax;
tmin = (bounds[sign[0] ].x - orig.x) * invDir.x;
tmax = (bounds[1 - sign[0]].x - orig.x) * invDir.x;
tymin = (bounds[sign[1] ].y - orig.y) * invDir.y;
tymax = (bounds[1 - sign[1]].y - orig.y) * invDir.y;
if ((tmin > tymax) || (tymin > tmax))
return false;
if (tymin > tmin)
tmin = tymin;
if (tymax < tmax)
tmax = tymax;
tzmin = (bounds[sign[2] ].z - orig.z) * invDir.z;
tzmax = (bounds[1 - sign[2]].z - orig.z) * invDir.z;
if ((tmin > tzmax) || (tzmin > tmax))
return false;
if (tzmin > tmin)
tmin = tzmin;
if (tzmax < tmax)
tmax = tzmax;
tHit = tmin;
return true;
}
template<typename T> inline T clamp(const T &v, const T &lo, const T &hi)
{ return std::max(lo, std::min(v, hi)); }
class GeomPrimitive
{
public:
GeomPrimitive(Matrix44f objectToWorld_) :
objectToWorld(objectToWorld_), worldToObject(objectToWorld.inverse())
{}
virtual ~GeomPrimitive() {}
virtual bool intersect(const Vec3f&, const Vec3f&, float &) const = 0;
Matrix44f objectToWorld;
Matrix44f worldToObject;
BBox<> bbox;
uint32_t test;
};
class Mesh : public GeomPrimitive
{
public:
Mesh(
uint32_t numPolygons,
const std::vector<uint32_t>& polygonNumVertsArray,
const std::vector<uint32_t>& polygonIndicesInVertexPool,
std::vector<Vec3f> vertexPool_,
Matrix44f objectToWorld_ = Matrix44f()) :
GeomPrimitive(objectToWorld_),
vertexPool(vertexPool_)
{
// pass by value (move constructor shouldn't even be called here ?
for (uint32_t i = 0; i < vertexPool.size(); ++i) {
matPointMult(objectToWorld, vertexPool[i]);
bbox.extendBy(vertexPool[i]);
}
// compute total number of triangles
for (uint32_t i = 0; i < numPolygons; ++i) {
assert(polygonNumVertsArray[i] >= 3);
numTriangles += polygonNumVertsArray[i] - 2;
}
// create array to store the triangle indices in the vertex pool
triangleIndicesInVertexPool.reserve(numTriangles * 3);
mailbox.reserve(numTriangles); // should all be initialized with 0
// for each face
for (uint32_t i = 0, offset = 0, currTriangleIndex = 0; i < numPolygons; ++i) {
// for each triangle in the face
for (uint32_t j = 0; j < polygonNumVertsArray[i] - 2; ++j) {
triangleIndicesInVertexPool[currTriangleIndex ] = polygonIndicesInVertexPool[offset ];
triangleIndicesInVertexPool[currTriangleIndex + 1] = polygonIndicesInVertexPool[offset + j + 1];
triangleIndicesInVertexPool[currTriangleIndex + 2] = polygonIndicesInVertexPool[offset + j + 2];
currTriangleIndex += 3;
}
offset += polygonNumVertsArray[i];
}
}
bool intersect(const Vec3f&, const Vec3f&, float &t) const;
uint32_t numTriangles = { 0 };
std::vector<uint32_t> triangleIndicesInVertexPool;
std::vector<Vec3f> vertexPool;
// [comment]
// Mailboxes are used by the Grid acceleration structure
// [/comment]
mutable std::vector<uint32_t> mailbox;
};
// use Moller-Trumbor method
bool rayTriangleIntersect(
const Vec3f& orig, const Vec3f &dir,
const Vec3f& v0, const Vec3f& v1, const Vec3f& v2,
float &t, float &u, float &v)
{
numRayTriangleTests++;
Vec3f v0v1 = v1 - v0;
Vec3f v0v2 = v2 - v0;
Vec3f pvec = cross(dir, v0v2);
float det = dot(v0v1, pvec);
// ray and triangle are parallel if det is close to 0
if (fabs(det) < kEpsilon) return false;
float invDet = 1 / det;
Vec3f tvec = orig - v0;
u = dot(tvec, pvec) * invDet;
if (u < 0 || u > 1) return false;
Vec3f qvec = cross(tvec, v0v1);
v = dot(dir, qvec) * invDet;
if (v < 0 || u + v > 1) return false;
t = dot(v0v2, qvec) * invDet;
numRayTriangleIntersections++;
return true;
}
bool Mesh::intersect(const Vec3f& rayOrig, const Vec3f& rayDir, float &tNear) const
{
// naive approach, loop over all triangles in the mesh and return true if one
// of the triangles at least is intersected
float t, u, v;
uint32_t intersectedTriIndex;
bool intersected = false;
// tNear should be set inifnity first time this function is called and it
// will get eventually smaller as the ray intersects geometry
for (uint32_t i = 0; i < numTriangles; ++i) {
if (rayTriangleIntersect(rayOrig, rayDir,
vertexPool[triangleIndicesInVertexPool[i * 3]],
vertexPool[triangleIndicesInVertexPool[i * 3 + 1]],
vertexPool[triangleIndicesInVertexPool[i * 3 + 2]], t, u, v) && t < tNear)
{
tNear = t;
intersectedTriIndex = i;
intersected = true;
}
}
return intersected;
}
#include "teapotdata.h"
Vec3f evalBezierCurve(const Vec3f *P, const float &t)
{
float b0 = (1 - t) * (1 - t) * (1 - t);
float b1 = 3 * t * (1 - t) * (1 - t);
float b2 = 3 * t * t * (1 - t);
float b3 = t * t * t;
return P[0] * b0 + P[1] * b1 + P[2] * b2 + P[3] * b3;
}
Vec3f evalBezierPatch(const Vec3f *controlPoints, const float &u, const float &v)
{
Vec3f uCurve[4];
for (size_t i = 0; i < 4; ++i)
uCurve[i] = evalBezierCurve(controlPoints + 4 * i, u);
return evalBezierCurve(uCurve, v);
}
Vec3f derivBezier(const Vec3f *P, const float &t)
{
return -3 * (1 - t) * (1 - t) * P[0] +
(3 * (1 - t) * (1 - t) - 6 * t * (1 - t)) * P[1] +
(6 * t * (1 - t) - 3 * t * t) * P[2] +
3 * t * t * P[3];
}
Vec3f dUBezier(const Vec3f* controlPoints, float u, float v)
{
Vec3f P[4];
Vec3f vCurve[4];
for (size_t i = 0; i < 4; ++i) {
P[0] = controlPoints[i];
P[1] = controlPoints[4 + i];
P[2] = controlPoints[8 + i];
P[3] = controlPoints[12 + i];
vCurve[i] = evalBezierCurve(P, v);
}
return derivBezier(vCurve, u);
}
Vec3f dVBezier(const Vec3f* controlPoints, float u, float v)
{
Vec3f uCurve[4];
for (size_t i = 0; i < 4; ++i) {
uCurve[i] = evalBezierCurve(controlPoints + 4 * i, u);
}
return derivBezier(uCurve, v);
}
std::vector<std::unique_ptr<const Mesh>> createUtahTeapot()
{
Matrix44f rotate90(1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1);
std::vector<std::unique_ptr<const Mesh>> meshes;
uint32_t width = 8, height = 8;
uint32_t numPolygons = width * height;
std::vector<uint32_t> polyNumVertsArray(numPolygons, 4);
std::vector<uint32_t> polyIndicesInVertPool(numPolygons * 4);
// set indices
for (uint32_t y = 0, offset = 0; y < height; ++y) {
for (uint32_t x = 0; x < width; ++x, offset += 4) {
// counter-clockwise to get the normal pointing in the right direction
polyIndicesInVertPool[offset ] = (width + 1) * y + x;
polyIndicesInVertPool[offset + 3] = (width + 1) * y + x + 1;
polyIndicesInVertPool[offset + 2] = (width + 1) * (y + 1) + x + 1;
polyIndicesInVertPool[offset + 1] = (width + 1) * (y + 1) + x;
}
}
Vec3f controlPoints[16];
for (uint32_t i = 0; i < kTeapotNumPatches; ++i) {
std::vector<Vec3f> vertPool((width + 1) * (height + 1));
for (uint32_t j = 0; j < 16; ++j) {
controlPoints[j].x = teapotVertices[teapotPatches[i][j] - 1][0],
controlPoints[j].y = teapotVertices[teapotPatches[i][j] - 1][1],
controlPoints[j].z = teapotVertices[teapotPatches[i][j] - 1][2];
}
for (uint32_t y = 0, currVertIndex = 0; y <= height; ++y) {
float v = y / (float)height;
for (uint32_t x = 0; x <= width; ++x, ++currVertIndex) {
float u = x / (float)width;
vertPool[currVertIndex] = evalBezierPatch(controlPoints, u, v);
matVecMult(rotate90, vertPool[currVertIndex]);
Vec3f dU = dUBezier(controlPoints, u, v);
Vec3f dV = dVBezier(controlPoints, u, v);
Vec3f N = cross(dU, dV);
}
}
meshes.emplace_back(new Mesh(numPolygons, polyNumVertsArray, polyIndicesInVertPool, vertPool));
}
return meshes;
}
void makeScene(std::vector<std::unique_ptr<const Mesh>>& meshes)
{
meshes = std::move(createUtahTeapot());
}
template<typename T>
inline
T degToRad(const T& angle) { return angle / 180.f * M_PI; }
struct Options
{
float fov = { 90 };
uint32_t width = { 640 };
uint32_t height = { 480 };
Matrix44f cameraToWorld, worldToCamera;
};
// [comment]
// The most basic acceleration class (the parent class of all the other acceleration structures)
// could have a *pure* virtual intersect() method but instead we decided in this implementation
// to have it supporting the basic ray-mesh intersection routine.
// [/comment]
class AccelerationStructure
{
public:
// [comment]
// We transfer owner ship of the mesh list to the acceleration structure. This makes
// more sense from a functional/structure stand point because the objects/meshes themselves
// should be destroyed/deleted when the acceleration structure is being deleted
// Ideally this means the render function() itself should be bounded (in terms of lifespan)
// to the lifespan of the acceleration structure (aka we should wrap the accel struc instance
// and the render method() within the same object, so that when this object is deleted,
// the render function can't be called anymore.
// [/comment]
AccelerationStructure(std::vector<std::unique_ptr<const Mesh>>& m) : meshes(std::move(m)) {}
virtual ~AccelerationStructure() {}
virtual bool intersect(const Vec3f& orig, const Vec3f& dir, const uint32_t& rayId, float& tHit) const
{
// [comment]
// Because we don't want to change the content of the mesh itself, just get a point to it so
// it's safer to make it const (which doesn't mean we can't change its assignment just that
// we can't do something like intersectedMesh->color = blue. You would get something like:
// "read-only variable is not assignable" error message at compile time)
// [/comment]
const Mesh* intersectedMesh = nullptr;
float t = kInfinity;
for (const auto& mesh: meshes) {
if (mesh->intersect(orig, dir, t) && t < tHit) {
intersectedMesh = mesh.get();
tHit = t;
}
}
return (intersectedMesh != nullptr);
}
protected:
const std::vector<std::unique_ptr<const Mesh>> meshes;
};
// [comment]
// Implementation of the ray-bbox method. If the ray intersects the bbox of a mesh then
// we test if the ray intersects the mesh contained by the bbox itself.
// [/comment]
class BBoxAcceleration : public AccelerationStructure
{
public:
BBoxAcceleration(std::vector<std::unique_ptr<const Mesh>>& m) : AccelerationStructure(m) {}
// [comment]
// Implement the ray-bbox acceleration method. The method consist of intersecting the
// ray against the bbox of the mesh first, and if the ray inteesects the boudning box
// then test if the ray intersects the mesh itsefl. It is obvious that the ray can't
// intersect the mesh if it doesn't intersect its boudning volume (a box in this case)
// [/comment]
virtual bool intersect(const Vec3f& orig, const Vec3f& dir, const uint32_t& rayId, float& tHit) const
{
const Mesh* intersectedMesh = nullptr;
const Vec3f invDir = 1 / dir;
const Vec3b sign(dir.x < 0, dir.y < 0, dir.z < 0);
float t = kInfinity;
for (const auto& mesh : meshes) {
// If you intersect the box
if (mesh->bbox.intersect(orig, invDir, sign, t)) {
// Then test if the ray intersects the mesh and if does then first check
// if the intersection distance is the nearest and if we pass that test as well
// then update tNear variable with t and keep a pointer to the intersected mesh
if (mesh->intersect(orig, dir, t) && t < tHit) {
tHit = t;
intersectedMesh = mesh.get();
}
}
}
// Return true if the variable intersectedMesh is not null, false otherwise
return (intersectedMesh != nullptr);
}
};
// [comment]
// Implementation of the Bounding Volume Hieratchy (BVH) acceleration structure
// [/comment]
class BVH : public AccelerationStructure
{
static const uint8_t kNumPlaneSetNormals = 7;
static const Vec3f planeSetNormals[kNumPlaneSetNormals];
struct Extents
{
Extents()
{
for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i)
d[i][0] = kInfinity, d[i][1] = -kInfinity;
}
void extendBy(const Extents& e)
{
for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) {
if (e.d[i][0] < d[i][0]) d[i][0] = e.d[i][0];
if (e.d[i][1] > d[i][1]) d[i][1] = e.d[i][1];
}
}
/* inline */
Vec3f centroid() const
{
return Vec3f(
d[0][0] + d[0][1] * 0.5,
d[1][0] + d[1][1] * 0.5,
d[2][0] + d[2][1] * 0.5);
}
bool intersect(const float*, const float*, float&, float&, uint8_t&) const;
float d[kNumPlaneSetNormals][2];
const Mesh* mesh;
};
struct Octree
{
Octree(const Extents& sceneExtents)
{
float xDiff = sceneExtents.d[0][1] - sceneExtents.d[0][0];
float yDiff = sceneExtents.d[1][1] - sceneExtents.d[1][0];
float zDiff = sceneExtents.d[2][1] - sceneExtents.d[2][0];
float maxDiff = std::max(xDiff, std::max(yDiff, zDiff));
Vec3f minPlusMax(
sceneExtents.d[0][0] + sceneExtents.d[0][1],
sceneExtents.d[1][0] + sceneExtents.d[1][1],
sceneExtents.d[2][0] + sceneExtents.d[2][1]);
bbox[0] = (minPlusMax - maxDiff) * 0.5;
bbox[1] = (minPlusMax + maxDiff) * 0.5;
root = new OctreeNode;
}
~Octree() { deleteOctreeNode(root); }
void insert(const Extents* extents) { insert(root, extents, bbox, 0); }
void build() { build(root, bbox); };
struct OctreeNode
{
OctreeNode* child[8] = { nullptr };
std::vector<const Extents *> nodeExtentsList; // pointer to the objects extents
Extents nodeExtents; // extents of the octree node itself
bool isLeaf = true;
};
struct QueueElement
{
const OctreeNode *node; // octree node held by this element in the queue
float t; // distance from the ray origin to the extents of the node
QueueElement(const OctreeNode *n, float tn) : node(n), t(tn) {}
// priority_queue behaves like a min-heap
friend bool operator < (const QueueElement &a, const QueueElement &b) { return a.t > b.t; }
};
OctreeNode* root = nullptr; // make unique son don't have to manage deallocation
BBox<> bbox;
private:
void deleteOctreeNode(OctreeNode*& node)
{
for (uint8_t i = 0; i < 8; i++) {
if (node->child[i] != nullptr) {
deleteOctreeNode(node->child[i]);
}
}
delete node;
}
void insert(OctreeNode*& node, const Extents* extents, const BBox<>& bbox, uint32_t depth)
{
if (node->isLeaf) {
if (node->nodeExtentsList.size() == 0 || depth == 16) {
node->nodeExtentsList.push_back(extents);
}
else {
node->isLeaf = false;
// Re-insert extents held by this node
while (node->nodeExtentsList.size()) {
insert(node, node->nodeExtentsList.back(), bbox, depth);
node->nodeExtentsList.pop_back();
}
// Insert new extent
insert(node, extents, bbox, depth);
}
}
else {
// Need to compute in which child of the current node this extents should
// be inserted into
Vec3f extentsCentroid = extents->centroid();
Vec3f nodeCentroid = (bbox[0] + bbox[1]) * 0.5;
BBox<> childBBox;
uint8_t childIndex = 0;
// x-axis
if (extentsCentroid.x > nodeCentroid.x) {
childIndex = 4;
childBBox[0].x = nodeCentroid.x;
childBBox[1].x = bbox[1].x;
}
else {
childBBox[0].x = bbox[0].x;
childBBox[1].x = nodeCentroid.x;
}
// y-axis
if (extentsCentroid.y > nodeCentroid.y) {
childIndex += 2;
childBBox[0].y = nodeCentroid.y;
childBBox[1].y = bbox[1].y;
}
else {
childBBox[0].y = bbox[0].y;
childBBox[1].y = nodeCentroid.y;
}
// z-axis
if (extentsCentroid.z > nodeCentroid.z) {
childIndex += 1;
childBBox[0].z = nodeCentroid.z;
childBBox[1].z = bbox[1].z;
}
else {
childBBox[0].z = bbox[0].z;
childBBox[1].z = nodeCentroid.z;
}
// Create the child node if it doesn't exsit yet and then insert the extents in it
if (node->child[childIndex] == nullptr)
node->child[childIndex] = new OctreeNode;
insert(node->child[childIndex], extents, childBBox, depth + 1);
}
}
void build(OctreeNode*& node, const BBox<>& bbox)
{
if (node->isLeaf) {
for (const auto& e: node->nodeExtentsList) {
node->nodeExtents.extendBy(*e);
}
}
else {
for (uint8_t i = 0; i < 8; ++i) {
if (node->child[i]) {
BBox<> childBBox;
Vec3f centroid = bbox.centroid();
// x-axis
childBBox[0].x = (i & 4) ? centroid.x : bbox[0].x;
childBBox[1].x = (i & 4) ? bbox[1].x : centroid.x;
// y-axis
childBBox[0].y = (i & 2) ? centroid.y : bbox[0].y;
childBBox[1].y = (i & 2) ? bbox[1].y : centroid.y;
// z-axis
childBBox[0].z = (i & 1) ? centroid.z : bbox[0].z;
childBBox[1].z = (i & 1) ? bbox[1].z : centroid.z;
// Inspect child
build(node->child[i], childBBox);
// Expand extents with extents of child
node->nodeExtents.extendBy(node->child[i]->nodeExtents);
}
}
}
}
};
std::vector<Extents> extentsList;
Octree* octree = nullptr;
public:
BVH(std::vector<std::unique_ptr<const Mesh>>& m);
bool intersect(const Vec3f&, const Vec3f&, const uint32_t&, float&) const;
~BVH() { delete octree; }
};
const Vec3f BVH::planeSetNormals[BVH::kNumPlaneSetNormals] = {
Vec3f(1, 0, 0),
Vec3f(0, 1, 0),
Vec3f(0, 0, 1),
Vec3f( sqrtf(3) / 3.f, sqrtf(3) / 3.f, sqrtf(3) / 3.f),
Vec3f(-sqrtf(3) / 3.f, sqrtf(3) / 3.f, sqrtf(3) / 3.f),
Vec3f(-sqrtf(3) / 3.f, -sqrtf(3) / 3.f, sqrtf(3) / 3.f),
Vec3f( sqrtf(3) / 3.f, -sqrtf(3) / 3.f, sqrtf(3) / 3.f)
};
BVH::BVH(std::vector<std::unique_ptr<const Mesh>>& m) : AccelerationStructure(m)
{
Extents sceneExtents; // that's the extent of the entire scene which we need to compute for the octree
extentsList.reserve(meshes.size());
for (uint32_t i = 0; i < meshes.size(); ++i) {
for (uint8_t j = 0; j < kNumPlaneSetNormals; ++j) {
for (const auto vtx : meshes[i]->vertexPool) {
float d = dot(planeSetNormals[j], vtx);
// set dNEar and dFar
if (d < extentsList[i].d[j][0]) extentsList[i].d[j][0] = d;
if (d > extentsList[i].d[j][1]) extentsList[i].d[j][1] = d;
}
}
sceneExtents.extendBy(extentsList[i]); // expand the scene extent of this object's extent
extentsList[i].mesh = meshes[i].get(); // the extent itself needs to keep a pointer to the object its holds
}
// Now that we have the extent of the scene we can start building our octree
// Using C++ make_unique function here but you don't need to, just to learn something...
octree = new Octree(sceneExtents);
for (uint32_t i = 0; i < meshes.size(); ++i) {
octree->insert(&extentsList[i]);
}
// Build from bottom up
octree->build();
}
bool BVH::Extents::intersect(
const float* precomputedNumerator,
const float* precomputedDenominator,
float& tNear, // tn and tf in this method need to be contained
float& tFar, // within the range [tNear:tFar]
uint8_t& planeIndex) const
{
numRayBoundingVolumeTests++;
for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) {
float tNearExtents = (d[i][0] - precomputedNumerator[i]) / precomputedDenominator[i];
float tFarExtents = (d[i][1] - precomputedNumerator[i]) / precomputedDenominator[i];
if (precomputedDenominator[i] < 0) std::swap(tNearExtents, tFarExtents);
if (tNearExtents > tNear) tNear = tNearExtents, planeIndex = i;
if (tFarExtents < tFar) tFar = tFarExtents;
if (tNear > tFar) return false;
}
return true;
}
bool BVH::intersect(const Vec3f& orig, const Vec3f& dir, const uint32_t& rayId, float& tHit) const
{
tHit = kInfinity;
const Mesh* intersectedMesh = nullptr;
float precomputedNumerator[BVH::kNumPlaneSetNormals];
float precomputedDenominator[BVH::kNumPlaneSetNormals];
for (uint8_t i = 0; i < kNumPlaneSetNormals; ++i) {
precomputedNumerator[i] = dot(planeSetNormals[i], orig);
precomputedDenominator[i] = dot(planeSetNormals[i], dir);
}
/*
tNear = kInfinity; // set
for (uint32_t i = 0; i < meshes.size(); ++i) {
numRayVolumeTests++;
float tn = -kInfinity, tf = kInfinity;
uint8_t planeIndex;
if (extents[i].intersect(precomputedNumerator, precomputedDenominator, tn, tf, planeIndex)) {
if (tn < tNear) {
intersectedMesh = meshes[i].get();
tNear = tn;
// normal = planeSetNormals[planeIndex];
}
}
}
*/
uint8_t planeIndex;
float tNear = 0, tFar = kInfinity; // tNear, tFar for the intersected extents
if (!octree->root->nodeExtents.intersect(precomputedNumerator, precomputedDenominator, tNear, tFar, planeIndex) || tFar < 0)
return false;
tHit = tFar;
std::priority_queue<BVH::Octree::QueueElement> queue;
queue.push(BVH::Octree::QueueElement(octree->root, 0));
while (!queue.empty() && queue.top().t < tHit) {
const Octree::OctreeNode *node = queue.top().node;
queue.pop();
if (node->isLeaf) {
for (const auto& e: node->nodeExtentsList) {
float t = kInfinity;
if (e->mesh->intersect(orig, dir, t) && t < tHit) {
tHit = t;
intersectedMesh = e->mesh;
}
}
}
else {
for (uint8_t i = 0; i < 8; ++i) {
if (node->child[i] != nullptr) {
float tNearChild = 0, tFarChild = tFar;
if (node->child[i]->nodeExtents.intersect(precomputedNumerator, precomputedDenominator, tNearChild, tFarChild, planeIndex)) {
float t = (tNearChild < 0 && tFarChild >= 0) ? tFarChild : tNearChild;
queue.push(BVH::Octree::QueueElement(node->child[i], t));
}
}
}
}
}
return (intersectedMesh != nullptr);
}
// [comment]
// Implementation of the Grid acceleration structure
// [/comment]
class Grid : public AccelerationStructure
{
struct Cell
{
Cell() {}
struct TriangleDesc
{
TriangleDesc(const Mesh* m, const uint32_t &t) : mesh(m), tri(t) {}
const Mesh* mesh;
uint32_t tri;
};
void insert(const Mesh* mesh, uint32_t t)
{ triangles.push_back(Grid::Cell::TriangleDesc(mesh, t)); }
bool intersect(const Vec3f&, const Vec3f&, const uint32_t&, float&, const Mesh*&) const;
std::vector<TriangleDesc> triangles;
};
public:
Grid(std::vector<std::unique_ptr<const Mesh>>& m);
~Grid()
{
for (uint32_t i = 0; i < resolution[0] * resolution[1] * resolution[2]; ++i)
if (cells[i] != NULL) delete cells[i];
delete [] cells;
}
bool intersect(const Vec3f&, const Vec3f&, const uint32_t&, float&) const;
Cell **cells;
BBox<> bbox;
Vec3<uint32_t> resolution;
Vec3f cellDimension;
};
Grid::Grid(std::vector<std::unique_ptr<const Mesh>>& m) : AccelerationStructure(m)
{
uint32_t totalNumTriangles = 0;
for (const auto& m : meshes) {
bbox.extendBy(m->bbox[0]);
bbox.extendBy(m->bbox[1]);
totalNumTriangles += m->numTriangles;
}
// Create the grid
Vec3f size = bbox[1] - bbox[0];
float cubeRoot = std::powf(totalNumTriangles / (size.x * size.y * size.z), 1. / 3.f);
for (uint8_t i = 0; i < 3; ++i) {
resolution[i] = std::floor(size[i] * cubeRoot);
if (resolution[i] < 1) resolution[i] = 1;
if (resolution[i] > 128) resolution[i] = 128;
}
cellDimension = size / resolution;
// [comment]
// Allocate memory - note that we don't create the cells yet at this
// point but just an array of pointers to cell. We will create the cells
// dynamically later when we are sure to insert something in them
// [/comment]
uint32_t numCells = resolution.x * resolution.y * resolution.z;
cells = new Grid::Cell* [numCells];
memset(cells, 0x0, sizeof(Grid::Grid*) * numCells);
for (const auto& m : meshes) {
for (uint32_t i = 0, off = 0; i < m->numTriangles; ++i, off += 3) {
Vec3f min(kInfinity), max(-kInfinity);
const Vec3f& v0 = m->vertexPool[m->triangleIndicesInVertexPool[off]];
const Vec3f& v1 = m->vertexPool[m->triangleIndicesInVertexPool[off + 1]];
const Vec3f& v2 = m->vertexPool[m->triangleIndicesInVertexPool[off + 2]];
for (uint8_t j = 0; j < 3; ++j) {
if (v0[j] < min[j]) min[j] = v0[j];
if (v1[j] < min[j]) min[j] = v1[j];
if (v2[j] < min[j]) min[j] = v2[j];
if (v0[j] > max[j]) max[j] = v0[j];
if (v1[j] > max[j]) max[j] = v1[j];
if (v2[j] > max[j]) max[j] = v2[j];
}
// Convert to cell coordinates
min = (min - bbox[0]) / cellDimension;
max = (max - bbox[0]) / cellDimension;
uint32_t zmin = clamp<uint32_t>(std::floor(min[2]), 0, resolution[2] - 1);
uint32_t zmax = clamp<uint32_t>(std::floor(max[2]), 0, resolution[2] - 1);
uint32_t ymin = clamp<uint32_t>(std::floor(min[1]), 0, resolution[1] - 1);
uint32_t ymax = clamp<uint32_t>(std::floor(max[1]), 0, resolution[1] - 1);
uint32_t xmin = clamp<uint32_t>(std::floor(min[0]), 0, resolution[0] - 1);
uint32_t xmax = clamp<uint32_t>(std::floor(max[0]), 0, resolution[0] - 1);
// Loop over the cells the triangle overlaps and insert
for (uint32_t z = zmin; z <= zmax; ++z) {
for (uint32_t y = ymin; y <= ymax; ++y) {
for (uint32_t x = xmin; x <= xmax; ++x) {
uint32_t index = z * resolution[0] * resolution[1] + y * resolution[0] + x;
if (cells[index] == NULL) cells[index] = new Grid::Cell;
cells[index]->insert(m.get(), i);
}
}
}
}
}
}
bool Grid::Cell::intersect(
const Vec3f& orig, const Vec3f& dir, const uint32_t& rayId,
float& tHit, const Mesh*& intersectedMesh) const
{
float uhit, vhit;
for (uint32_t i = 0; i < triangles.size(); ++i) {
// [comment]
// Be sure that rayId is never 0 - because all mailbox values
// in the array are initialized with 0 too
// [/comment]
if (rayId != triangles[i].mesh->mailbox[triangles[i].tri]) {
triangles[i].mesh->mailbox[triangles[i].tri] = rayId;
const Mesh *mesh = triangles[i].mesh;
uint32_t j = triangles[i].tri * 3;
const Vec3f &v0 = mesh->vertexPool[mesh->triangleIndicesInVertexPool[j ]];
const Vec3f &v1 = mesh->vertexPool[mesh->triangleIndicesInVertexPool[j + 1]];
const Vec3f &v2 = mesh->vertexPool[mesh->triangleIndicesInVertexPool[j + 2]];
float t, u, v;
if (rayTriangleIntersect(orig, dir, v0, v1, v2, t, u, v)) {
if (t < tHit) {
tHit = t;
uhit = u;
vhit = v;
intersectedMesh = triangles[i].mesh;
}
}
}
}
return (intersectedMesh != nullptr);
}
bool Grid::intersect(const Vec3f& orig, const Vec3f& dir, const uint32_t& rayId, float& tHit) const
{
const Vec3f invDir = 1 / dir;
const Vec3b sign(dir.x < 0, dir.y < 0, dir.z < 0);
float tHitBox;
if (!bbox.intersect(orig, invDir, sign, tHitBox)) return false;
// initialization step
Vec3i exit, step, cell;
Vec3f deltaT, nextCrossingT;
for (uint8_t i = 0; i < 3; ++i) {
// convert ray starting point to cell coordinates
float rayOrigCell = ((orig[i] + dir[i] * tHitBox) - bbox[0][i]);
cell[i] = clamp<uint32_t>(std::floor(rayOrigCell / cellDimension[i]), 0, resolution[i] - 1);
if (dir[i] < 0) {
deltaT[i] = -cellDimension[i] * invDir[i];
nextCrossingT[i] = tHitBox + (cell[i] * cellDimension[i] - rayOrigCell) * invDir[i];
exit[i] = -1;
step[i] = -1;
}
else {
deltaT[i] = cellDimension[i] * invDir[i];
nextCrossingT[i] = tHitBox + ((cell[i] + 1) * cellDimension[i] - rayOrigCell) * invDir[i];
exit[i] = resolution[i];
step[i] = 1;
}
}
// Walk through each cell of the grid and test for an intersection if
// current cell contains geometry
const Mesh* intersectedMesh = nullptr;
while (1) {
uint32_t o = cell[2] * resolution[0] * resolution[1] + cell[1] * resolution[0] + cell[0];
if (cells[o] != nullptr) {
cells[o]->intersect(orig, dir, rayId, tHit, intersectedMesh);
//if (intersectedMesh != nullptr) { ray.color = cells[o]->color; }
}
uint8_t k =
((nextCrossingT[0] < nextCrossingT[1]) << 2) +
((nextCrossingT[0] < nextCrossingT[2]) << 1) +
((nextCrossingT[1] < nextCrossingT[2]));
static const uint8_t map[8] = {2, 1, 2, 1, 2, 2, 0, 0};
uint8_t axis = map[k];
if (tHit < nextCrossingT[axis]) break;
cell[axis] += step[axis];
if (cell[axis] == exit[axis]) break;
nextCrossingT[axis] += deltaT[axis];
}
return (intersectedMesh != nullptr);
}
// [comment]
// Main Render() function. Loop over each pixel in the image and trace a primary
// ray starting from the camera origin and passing through the current pixel. If they
// ray intersects geometry in the scene return some color (the color of the object)
// otherwise nothing (black or the color of the background)
// [/comment]
void render(const std::unique_ptr<AccelerationStructure>& accel, const Options& options)
{
std::unique_ptr<Vec3f []> buffer(new Vec3f[options.width * options.height]);
Vec3f orig(0,0,5);
matPointMult(options.cameraToWorld, orig);
float scale = std::tan(degToRad<float>(options.fov * 0.5));
float imageAspectRatio = options.width / static_cast<float>(options.height);
assert(imageAspectRatio > 1);
uint32_t rayId = 1; // Start at 1 not 0!! (see Grid code and mailboxing)
for (uint32_t j = 0; j < options.height; ++j) {
for (uint32_t i = 0; i < options.width; ++i) {
Vec3f dir((2 * (i + 0.5f) / options.width - 1) * scale * imageAspectRatio,
(1 - 2 * (j + 0.5) / options.height) * scale,
-1);
matVecMult(options.cameraToWorld, dir);
normalize(dir);
numPrimaryRays++;
float tHit = kInfinity;
buffer[j * options.width + i] = (accel->intersect(orig, dir, rayId++, tHit)) ? Vec3f(1) : Vec3f(0);
}
}
// store to PPM file
std::ofstream ofs;
ofs.open("image.ppm");
ofs << "P6\n" << options.width << " " << options.height << "\n255\n";
for (uint32_t i = 0; i < options.width * options.height; ++i) {
Vec3<uint8_t> pixRgb;
pixRgb.x = static_cast<uint8_t>(255 * std::max(0.f, std::min(1.f, buffer[i].x)));
pixRgb.y = static_cast<uint8_t>(255 * std::max(0.f, std::min(1.f, buffer[i].y)));
pixRgb.z = static_cast<uint8_t>(255 * std::max(0.f, std::min(1.f, buffer[i].z)));
ofs << pixRgb.x << pixRgb.y << pixRgb.z;
}
ofs.close();
}
void exportMesh(const std::vector<std::unique_ptr<const Mesh>>& meshes)
{
std::ofstream f;
f.open("mesh.obj");
uint32_t k = 0, off = 0;
for (const auto& mesh : meshes) {
f << "g default" << std::endl;
for (uint32_t i = 0; i < mesh->vertexPool.size(); ++i) {
f << "v " << mesh->vertexPool[i].x << " " << mesh->vertexPool[i].y << " " << mesh->vertexPool[i].z << std::endl;
}
f << "g mesh" << k++ << std::endl;
for (uint32_t i = 0; i < mesh->numTriangles; ++i) {
f << "f " << mesh->triangleIndicesInVertexPool[i * 3] + 1 + off << " " <<
mesh->triangleIndicesInVertexPool[i * 3 + 1] + 1 + off << " " <<
mesh->triangleIndicesInVertexPool[i * 3 + 2] + 1 + off << std::endl;
}
off += mesh->vertexPool.size();
}
f.close();
}
int main(int argc, char **argv)
{
std::vector<std::unique_ptr<const Mesh>> meshes;
makeScene(meshes);
// [comment]
// Create the acceleration structure
// [/comment]
#if defined(ACCEL_BBOX)
std::unique_ptr<AccelerationStructure> accel(new BBoxAcceleration(meshes));
#elif defined(ACCEL_BVH)
std::unique_ptr<AccelerationStructure> accel(new BVH(meshes));
#elif defined(ACCEL_GRID)
std::unique_ptr<AccelerationStructure> accel(new Grid(meshes));
#else
std::unique_ptr<AccelerationStructure> accel(new AccelerationStructure(meshes));
#endif
//exportMesh(meshes);
uint32_t numTriangles = 0;
for (const auto& mesh : meshes) {
numTriangles += mesh->numTriangles;
}
Options options;
using Time = std::chrono::high_resolution_clock;
using fsec = std::chrono::duration<float>;
auto t0 = Time::now();
render(accel, options);
auto t1 = Time::now();
fsec fs = t1 - t0;
std::cout << "Render time | " << fs.count() << " sec" << std::endl;
std::cout << "Total number of triangles | " << numTriangles << std::endl;
std::cout << "Total number of primary rays | " << numPrimaryRays << std::endl;
std::cout << "Total number of ray-bbox tests | " << numRayBBoxTests << std::endl;
std::cout << "Total number of ray-boundvolume tests | " << numRayBoundingVolumeTests << std::endl;
std::cout << "Total number of ray-triangles tests | " << numRayTriangleTests << std::endl;
std::cout << "Total number of ray-triangles intersections | " << numRayTriangleIntersections << std::endl;
return 0;
}