Raytracer (CPU)
During my first year of study at BUas, I made this CPU whited style raytracer. The raytracer was written in C++ over 8 weeks.
This project is one of the first to introduce me to graphics development, writing a math library, profiling and doing optimizations. This project also made me realize how fast computers are these days and what they are capable of.
The project also inspired many other raytracing projects such as Luminous and the Rust-Tracer.
In the next section, I will be going more in-depth on some of the more interesting bits of this project.
The tech behind making pretty pictures
Raytracing
Raytracing is a rendering technique where we simulate how light behaves. The technique itself is relatively old but has recently been popularized by the specialized raytracing cores on Graphics Cards (for example RTX).
With raytracing you start by shooting rays from the camera, these rays could then intersect geometry. If this is the case, you check the material properties and decide to reflect/refract or shade the object. For the shading of the objects, you shoot rays towards light sources to see if they are lit or in shade.
Show Code
renderer.cpp
Vector3 Renderer::Trace(Ray* ray, int lastPrimitiveIndex, int count)
{
if (count > _MAX_RAY_BOUNCE)
{
return {1.f,1.f,1.f};
}
int pIndex = -1;
FindIntersection(ray, pIndex);
Vector3 colour = { 0.04f,0.627f,0.804f }; //Sky Colour.
if (ray->t > 0 && ray->t < 100 && pIndex > -1)
{
Vector3 point = ray->origin + (ray->t * ray->direction);
Primitive* primitive = _scene.GetPrimitive(pIndex);
Material* mat = primitive->GetMaterial();
Vector3 normal = primitive->GetNormal(point);
unsigned int flags = mat->GetFlags();
if (flags & (unsigned int)MaterialFlags::Reflective)
{
// Is half reflective?
if (flags & (unsigned int)MaterialFlags::Opaque)
{
Vector3 baseColour = mat->GetBaseColour(point);
Vector3 opaqueColour = Illuminate(point, normal, ray->direction, mat->GetBaseColour(point));
Vector3 reflectiveColour = Reflect(point, normal, ray->direction, pIndex, count);
float reflectiveness = mat->GetReflectiveness();
return (1.f - reflectiveness) * (opaqueColour)+reflectiveness * reflectiveColour * baseColour;
}
return Reflect(point, normal, ray->direction,pIndex, count);
}
if (flags & (unsigned int)MaterialFlags::Opaque)
{
Vector3 baseColour = mat->GetBaseColour(point);
return Illuminate(point, normal, ray->direction, baseColour);
}
if (flags & (unsigned int)MaterialFlags::Dielectric)
{
Vector3 baseColour = mat->GetBaseColour(point);
float n = 1.0f;
float nt = mat->GetRefractionIndex();
float fresnel = Fresnel(normal, ray->direction, n, nt);
if (normal.Dot(ray->direction) > 0)
{
n = 1.f/n;
nt = 1.f/nt;
normal = -1.f * normal;
}
// Still in glass
Vector3 absorb = { 1,1,1 };
if (pIndex == lastPrimitiveIndex)
{
// Beers law.
Vector3 col = Vector3(1, 1, 1) - baseColour;
absorb.x = expf(-col.x * ray->t);
absorb.y = expf(-col.y * ray->t);
absorb.z = expf(-col.z * ray->t);
}
return (fresnel * Reflect(point, normal, ray->direction, pIndex, count) + (1 - fresnel) * Refract(point, normal, ray->direction, n, nt, pIndex, count)) * absorb;
}
// No flags = unlit.
return mat->GetBaseColour(point);
}
return colour;
}
bool Renderer::TraceShadowRay(Ray* ray, float maxDistance)
{
return _scene.CheckIntersectionFast(ray, 0.f, maxDistance - 0.0001f);
}
#pragma warning(disable:4100)
Vector3 Renderer::Illuminate(Vector3 point, Vector3 normal, Vector3 rayDirection, Vector3 baseColour)
{
Vector3 colour = { 0,0,0 };
#ifdef _USE_AREA_LIGHTS
for (int i = 0; i < MAX_AREA_LIGHTS; i++)
{
Quad* areaLight = static_cast<Quad*>(_scene.GetAreaLight(i));
// SAMPLE 4 points.
for (int j = 0; j < _NUM_AREA_SAMPLES; j++)
{
Ray lightRay;
Vector3 randPoint = areaLight->GetRandomPosition();
Vector3 lightDir = (randPoint - point);
float lightDistance = lightDir.Magnitude();
lightRay.origin = point + (lightDir * 0.0001f);
lightRay.direction = lightDir.Normalized();
bool castShadow = TraceShadowRay(&lightRay, lightDistance);
if (!castShadow) {
float intensity = clamp(normal.Dot(lightDir), 0.f, 1.f);
float inv = 1.f / (lightDistance * 2);
Vector3 c = baseColour * ((intensity * areaLight->GetMaterial()->GetBaseColour(randPoint)) * inv); // Difuse Colour
colour += c;
}
else
{
colour += Vector3(0.04f, 0.627f, 0.804f) * (0.7f / 4);
}
}
}
constexpr int samplesDividedByLights = (MAX_AREA_LIGHTS * _NUM_AREA_SAMPLES);
colour *= 1.f / samplesDividedByLights;
#else
for (int i = 0; i < _scene.GetLightCount(); i++)
{
Light* lamp = _scene.GetLight(i);
Vector3 lightDir = (lamp->position - point);
float lightDistance = lightDir.Magnitude();
lightDir = lightDir.Normalized();
Ray lightRay;
lightRay.origin = point + (lightDir * 0.0001f);
lightRay.direction = lightDir;
bool castShadow = TraceShadowRay(&lightRay, lightDistance);
if (!castShadow) {
float intensity = clamp(normal.Dot(lightDir), 0.f, 1.f);
Vector3 viewDirection = -1.f * rayDirection;
Vector3 halfVector = (lightDir + viewDirection).Normalized();
float halfDot = normal.Dot(halfVector);
float inv = 1.f / (lightDistance * 2);
Vector3 c = baseColour * ((intensity * lamp->colour) * inv); // Difuse Colour
if (halfDot > 0) {
float specIntensity = powf(clamp(halfDot, 0.f, 1.f), 50.f);
c += (specIntensity * lamp->colour) * inv; // Specular color
}
colour += c;
}
}
#endif
return colour;
}
Vector3 Renderer::Reflect(Vector3 point, Vector3 normal, Vector3 rayDirection, int lastPrimitiveIndex, int count)
{
Vector3 reflection = (rayDirection - 2.f * normal * rayDirection.Dot(normal)).Normalized();
Ray reflectionRay;
reflectionRay.origin = point + (reflection * 0.0001f);
reflectionRay.direction = reflection;
return Trace(&reflectionRay,lastPrimitiveIndex, ++count);
}
float Renderer::Fresnel(Vector3 normal, Vector3 rayDirection, float n, float nt)
{
if (normal.Dot(rayDirection) < 0) { rayDirection = -1.f * rayDirection; }
//// Reflectance at normal incidence
float reflectance = powf((n - nt) / (n + nt), 2.f);
//// Schlick approximation
return reflectance + (1.f - reflectance) * (powf(1.f - normal.Dot(rayDirection), 5));
}
Vector3 Renderer::Refract(Vector3 point, Vector3 normal, Vector3 rayDirection, float n, float nt, int lastPrimitiveIndex, int count)
{
Vector3 a = n * (rayDirection - normal * rayDirection.Dot(normal)) / nt;
float power = powf(rayDirection.Dot(normal), 2.f);
float sqr = sqrtf(1.f - (((n * n) * 1.f - power) / (nt * nt)));
if (sqr >= 0) {
Vector3 b = normal * sqr;
Vector3 refractDir = a - b;
Ray refraction;
refraction.direction = refractDir.Normalized();
refraction.origin = point + refractDir * 0.0001f;
return Trace(&refraction,lastPrimitiveIndex, ++count);
}
return { 1,0,1 };
}Math Library
For the project, another requirement was that the raytracer must be built using your own Math Library. This was a great opportunity to brush up on my math skills as well as understand the math involved in making these projects.
Show Code
vector3.h
#pragma once
#include "FloatOperations.h"
#include "Vector2.h"
#include "nmmintrin.h"
// Disables nameless union warning, allowing the use of v.x for example.
#pragma warning (disable:4201)
namespace MathLib
{
class Vector3
{
public:
union
{
float f[3];
struct
{
float x;
float y;
float z;
};
};
public:
Vector3() : x(0.f), y(0.f), z(0.f) {}
Vector3(float x, float y, float z) : x(x), y(y), z(z) {}
Vector3(Vector2 a, float z = 0) : x(a.x), y(a.y), z(z) {}
Vector3(float v) : x(v), y(v), z(v) {}
inline float Dot(const Vector3& vector) const
{
return x * vector.x + y * vector.y + z * vector.z;
}
inline Vector3 Cross(const Vector3& v) const
{
return Vector3(y * v.z - v.y * z, z * v.x - v.z * x, x * v.y - v.x * y);
}
inline float Magnitude() const
{
return sqrtf(x * x + y * y + z * z);
}
inline float SqrMagnitude() const
{
return x * x + y * y + z * z;
}
void Normalize()
{
float l = Magnitude();
if (l > 1)
{
float inverse = 1.f / l;
x *= inverse;
y *= inverse;
z *= inverse;
}
}
Vector3 Normalized()
{
Vector3 out = *this;
out.Normalize();
return out;
}
#pragma region operators
Vector3 operator+(const Vector3& v) const
{
return Vector3(x + v.x, y + v.y, z + v.z);
}
Vector3 operator-(const Vector3& v) const
{
return Vector3(x - v.x, y - v.y, z - v.z);
}
Vector3 operator*(float value) const
{
return Vector3(x * value, y * value, z * value);
}
Vector3 operator/(float value) const
{
return Vector3(x / value, y / value, z / value);
}
Vector3 operator*(const Vector3& v) const
{
return Vector3(x * v.x, y * v.y, z * v.z);
}
void operator*=(const Vector3& v)
{
x *= v.x;
y *= v.y;
z *= v.z;
}
void operator +=(const Vector3& v)
{
x += v.x;
y += v.y;
z += v.z;
}
void operator -=(const Vector3& v)
{
x -= v.x;
y -= v.y;
z -= v.z;
}
void operator *=(const float value)
{
x *= value;
y *= value;
z *= value;
}
void operator /=(const float value)
{
x /= value;
y /= value;
z /= value;
}
bool operator==(const Vector3& other) const
{
return x == other.x && y == other.y && z == other.z;
}
bool operator!=(const Vector3& other) const
{
return !(*this == other);
}
#pragma endregion
};
/// Multiplication with Rhs Vector
inline Vector3 operator*(float val, const Vector3& rhs)
{
return Vector3(rhs.x * val, rhs.y * val, rhs.z * val);
}
}matrix3x3.h
#pragma once
#include "Vector3.h"
namespace MathLib
{
class Matrix3x3
{
public:
float f[9];
public:
#pragma region
// Get value from Matrix in column major order.
float& operator()(unsigned int row, unsigned int col)
{
return f[row + 3 * col];
}
static Matrix3x3 Identity()
{
return
{
1,0,0,
0,1,0,
0,0,1
};
}
static Matrix3x3 RotationX(const float& angle)
{
return
{
1,0,0,
0,cos(angle),-sin(angle),
0,sin(angle),cos(angle)
};
}
static Matrix3x3 RotationY(const float& angle)
{
return
{
cos(angle), 0, sin(angle),
0,1,0,
-sin(angle),0,cos(angle)
};
}
static Matrix3x3 RotationZ(const float& angle)
{
return
{
cos(angle),-sin(angle),0,
sin(angle),cos(angle),0,
0,0,1
};
}
Matrix3x3 operator+(const Matrix3x3& m) const
{
Matrix3x3 result;
for (int i = 0; i < 9; i++)
{
result.f[i] = f[i] + m.f[i];
}
return result;
}
Matrix3x3 operator-(const Matrix3x3& m) const
{
Matrix3x3 result;
for (int i = 0; i < 9; i++)
{
result.f[i] = f[i] - m.f[i];
}
return result;
}
// Scalar multiplication
Matrix3x3 operator*(const float& s) const
{
Matrix3x3 result;
for (int i = 0; i < 9; i++)
{
result.f[i] = f[i] * s;
}
return result;
}
Vector3 operator*(const Vector3& v) const
{
Vector3 result;
result.x = f[0] * v.x + f[3] * v.y + f[6] * v.z;
result.y = f[1] * v.x + f[4] * v.y + f[7] * v.z;
result.z = f[2] * v.x + f[5] * v.y + f[8] * v.z;
return result;
}
// 0 3 6
// 1 4 7
// 2 5 8
Matrix3x3 operator*(const Matrix3x3& m) const
{
Matrix3x3 result;
result.f[0] = f[0] * m.f[0] + f[3] * m.f[1] + f[6] * m.f[2];
result.f[1] = f[1] * m.f[0] + f[4] * m.f[1] + f[7] * m.f[2];
result.f[2] = f[2] * m.f[0] + f[5] * m.f[1] + f[8] * m.f[2];
result.f[3] = f[0] * m.f[3] + f[3] * m.f[4] + f[6] * m.f[5];
result.f[4] = f[1] * m.f[3] + f[4] * m.f[4] + f[7] * m.f[5];
result.f[5] = f[2] * m.f[3] + f[5] * m.f[4] + f[8] * m.f[5];
result.f[6] = f[0] * m.f[6] + f[3] * m.f[7] + f[6] * m.f[8];
result.f[7] = f[1] * m.f[6] + f[4] * m.f[7] + f[7] * m.f[8];
result.f[8] = f[2] * m.f[6] + f[5] * m.f[7] + f[8] * m.f[8];
return result;
}
};
}
Bounding Volume Hierarchy
To make the raytracer somewhat real-time, optimization is required, One way we can optimize is by minimizing the required intersection checks we have to do for an individual ray.
Instead of looping over all objects that our ray can hit, we can divide the scene into volumes, covering a single or a group of objects. We can then create a tree-like structure by combining multiple of these volumes in another volume. When we want to see if our ray hits an object, we iteratively check the volumes to see if the ray could hit the volumes, and thus potentially the child volumes. Eventually, we either find a volume that holds some objects with a high chance of being intersected by the ray, or we don’t and we hit nothing, in which case we can render a sky.
Show Code
bvh.cpp
#include "headers\BVH.h"
#include "headers\primitives\AABB.h"
#include <iostream>
void BVH::CheckIntersection(Ray* ray, int& primitiveIndex)
{
float distance = INFINITY;
Traverse(ray, &_pool[0], primitiveIndex, distance);
ray->t = distance;
}
void BVH::Traverse(Ray* ray, BVHNode* node, int& primitiveIndex, float& t)
{
if (!(node->bounds.Intersect(ray)))
{
return;
};
if (ray->t > t) { return; }
if (node->count < MAX_BVH_NODE_PRIMITIVES && node->bounds.GetVolume() < MAX_BVH_NODE_VOLUME)
{
IntersectPrimitives(ray, node, primitiveIndex, t);
}
else
{
Traverse(ray, node->left, primitiveIndex,t);
Traverse(ray, node->right, primitiveIndex,t);
}
}
void BVH::IntersectPrimitives(Ray* ray, BVHNode* node, int& primitiveIndex, float& t)
{
ray->t = t;
for (int i = node->leftFirst; i < node->leftFirst + node->count; i++)
{
if (_primitives[_indices[i]]->Intersect(ray))
{
if (ray->t < t)
{
primitiveIndex = _indices[i];
t = ray->t;
}
}
}
}
void BVH::ConstructBVH(Primitive** primitives, const unsigned int primitiveCount)
{
_primitives = primitives;
for (unsigned int i = 0; i < primitiveCount; i++)
{
_indices[i] = i;
}
//Setup debug colours
for (unsigned int i = 0; i < (primitiveCount * 2) - 1; i++)
{
float r = (static_cast<float>(rand()) / static_cast<float>(RAND_MAX / (+1)));
float g = (static_cast<float>(rand()) / static_cast<float>(RAND_MAX / (+1)));
float b = (static_cast<float>(rand()) / static_cast<float>(RAND_MAX / (+1)));
_debugColours[i] = Vector3(r,g,b);
}
_pool = new BVHNode[(primitiveCount * 2) - 1];
_rootNode = &_pool[0];
_poolPtr = 1;
_rootNode->leftFirst = 0;
_rootNode->count = primitiveCount;
CalculateBounds(_rootNode);
Subdivide(_rootNode);
}
void BVH::Subdivide(BVHNode* node)
{
if (node->count < MAX_BVH_NODE_PRIMITIVES && node->bounds.GetVolume() < MAX_BVH_NODE_VOLUME)
{
std::cout << "NODE <LEAF>" << std::endl;
for (int i = node->leftFirst; i < node->leftFirst + node->count; i++)
{
std::cout << "NODE CONTAINS: " << i << " - " << _indices[i] << std::endl;;
}
std::cout << "NODE <LEAF> END" << std::endl;
std::cout << std::endl;
return;
}
int newCount = node->count / 2;
BVHNode* left = &_pool[_poolPtr++];
left->leftFirst = node->leftFirst;
left->count = newCount;
BVHNode* right = &_pool[_poolPtr++];
right->leftFirst = left->leftFirst + newCount;
right->count = right->leftFirst + newCount;
node->left = left;
node->right = right;
SplitPlane split = GetSplitPlane(node);
Partition(node->leftFirst, node->count, split.value, split.direction);
std::cout << "SORT" << std::endl;
std::cout << "PIVOT WAS: " << split.value << std::endl;
for (int i = node->leftFirst; i < node->leftFirst + node->count; i++)
{
std::cout << "INDEX: " << i << " - " << _indices[i] << " VAL: " << EvaluatePrimitive(_indices[i],split.direction) << std::endl;;
}
std::cout << "SORT END" << std::endl;
int c = -1;
for (int i = node->leftFirst; i < node->leftFirst + node->count; i++)
{
c++;
float debug = EvaluatePrimitive(_indices[i], split.direction);
if (debug >= split.value)
{
break;
}
}
left->leftFirst = node->leftFirst;
left->count = c;
right->leftFirst = left->leftFirst + c;
right->count = (node->count - c);
CalculateBounds(left);
CalculateBounds(right);
Subdivide(left);
Subdivide(right);
}
Vector3 BVH::DrawBVH(Ray* r)
{
Vector3 col = { 0,0,0 };
for (int i = 0; i < (MAX_PRIMITIVES * 2) - 1; i++)
{
BVHNode* n = &_pool[i];
if (n->count > 0 && n->count <= MAX_BVH_NODE_PRIMITIVES && n->bounds.GetVolume() < MAX_BVH_NODE_VOLUME)
{
if (n->bounds.Intersect(r))
{
col = _debugColours[i];
}
}
}
return col;
}
void BVH::Partition(int first, int count, float pivot, Direction dir)
{
int i = first;
for (int j = first; j < first + count; j++)
{
if (EvaluatePrimitive(_indices[j], dir) < pivot)
{
// Swap I with J
int temp = _indices[i];
_indices[i] = _indices[j];
_indices[j] = temp;
i++;
}
}
}
SplitPlane BVH::GetSplitPlane(BVHNode* node)
{
AABB bounds = node->bounds;
Vector3 min = bounds.GetMin();
Vector3 max = bounds.GetMax();
Vector3 lenght;
lenght.x = max.x - min.x;
lenght.y = max.y - min.y;
lenght.z = max.z - min.z;
if (lenght.x >= lenght.y && lenght.x >= lenght.z)
{
return SplitPlane((min.x + max.x) * 0.5f, Direction::X);
}
if (lenght.y >= lenght.x && lenght.y >= lenght.z)
{
return SplitPlane((min.y + max.y) * 0.5f, Direction::Y);
}
if (lenght.z >= lenght.x && lenght.z >= lenght.y)
{
return SplitPlane((min.z + max.z) * 0.5f, Direction::Z);
}
// If all axis are the same lenght, split on X axis.
return SplitPlane(lenght.x * 0.5f, Direction::X);
}
inline float BVH::EvaluatePrimitive(int index, Direction dir) const
{
return EvaluateVector(((Primitive*)_primitives[index])->GetPosition(), dir);
}
inline float BVH::EvaluateVector(Vector3 vector, Direction dir) const
{
switch (dir)
{
case Direction::X:
return vector.x;
case Direction::Y:
return vector.y;
case Direction::Z:
return vector.z;
default:
return -1.f;
}
}
void BVH::CalculateBounds(BVHNode* node)
{
Vector3 _min = Vector3(INFINITY,INFINITY,INFINITY);
Vector3 _max = Vector3(-INFINITY, -INFINITY, -INFINITY);
for (int i = node->leftFirst; i < node->leftFirst + node->count; i++)
{
Vector3 min = ((Primitive*)_primitives[_indices[i]])->GetMin();
Vector3 max = ((Primitive*)_primitives[_indices[i]])->GetMax();
// Compare and set min.
for (int j = 0; j < 3; j++)
{
if (min.f[j] < _min.f[j]) { _min.f[j] = min.f[j]; }
}
// Compare and set max.
for (int j = 0; j < 3; j++)
{
if (max.f[j] > _max.f[j]) { _max.f[j] = max.f[j]; }
}
}
node->bounds.SetMin(_min);
node->bounds.SetMax(_max);
}