//////////////////////////////////////////////////////////////////////
// Provides the framework for a raytracer.
////////////////////////////////////////////////////////////////////////

#include "raytrace.h"

#define STB_IMAGE_IMPLEMENTATION
#define STBI_FAILURE_USERMSG
#include "stb_image.h"

// A good quality *thread-safe* Mersenne Twister random number generator.
#include <random>
std::random_device device;
std::mt19937_64 RNGen(device());
std::uniform_real_distribution<> myrandom(0.0, 1.0);
// Call myrandom(RNGen) to get a uniformly distributed random number in [0,1].

void Scene::Finit()
{
    this->accelStruc = new AccelerationBvh(this->shapes);
}

void Scene::triangleMesh(MeshData* mesh)
{
    // Use the currentMat from the scene file (set by brdf command)
    // This allows the scene file to control mesh appearance
    Material* meshMat = currentMat;

    // Convert mesh data into individual Triangle shapes and add to scene
    for (const auto& tri : mesh->triangles) {
        // tri is an ivec3 containing three vertex indices
        vec3 v0 = mesh->vertices[tri[0]].pnt;
        vec3 v1 = mesh->vertices[tri[1]].pnt;
        vec3 v2 = mesh->vertices[tri[2]].pnt;

        // Create a Triangle shape and add it to the scene with the scene file's material
        this->shapes.push_back(new Triangle(v0, v1, v2, meshMat));
    }
}

Texture::Texture(const std::string& bpath) : id(0)
{
    // Replace backslashes with forward slashes -- Good for Linux, and maybe Windows?
    std::string path = bpath;
    std::string bs = "\\";
    std::string fs = "/";
    while (path.find(bs) != std::string::npos) {
        path.replace(path.find(bs), 1, fs);
    }

    // Does the file exist?
    std::ifstream find_it(path.c_str());
    if (find_it.fail()) {
        std::cerr << "Texture file not found: " << path << std::endl;
        exit(-1);
    }
    else {
        // Read image, and check for success
        stbi_set_flip_vertically_on_load(true);
        image = stbi_load(path.c_str(), &width, &height, &depth, 4);
        printf("%d %d %d %s\n", depth, width, height, path.c_str());
        if (!image) {
            printf("\nRead error on file %s:\n  %s\n\n", path.c_str(), stbi_failure_reason());
            exit(-1);
        }
    }
}

quat Orientation(int i,
    const std::vector<std::string>& strings,
    const std::vector<float>& f)
{
    quat q(1, 0, 0, 0); // Unit quaternion
    while (i < strings.size()) {
        std::string c = strings[i++];
        if (c == "x")
            q *= angleAxis(f[i++] * Radians, Xaxis());
        else if (c == "y")
            q *= angleAxis(f[i++] * Radians, Yaxis());
        else if (c == "z")
            q *= angleAxis(f[i++] * Radians, Zaxis());
        else if (c == "q") {
            q *= quat(f[i + 0], f[i + 1], f[i + 2], f[i + 3]);
            i += 4;
        }
        else if (c == "a") {
            q *= angleAxis(f[i + 0] * Radians, normalize(vec3(f[i + 1], f[i + 2], f[i + 3])));
            i += 4;
        }
    }
    return q;
}

void Scene::Command(const std::vector<std::string>& strings,
    const std::vector<float>& f)
{
    if (strings.size() == 0) return;
    std::string c = strings[0];

    if (c == "screen") {
        // syntax: screen width height        
        width = int(f[1]);
        height = int(f[2]);
    }

    else if (c == "camera") {
        // syntax: camera x y z   ry   <orientation spec>
        // Eye position (x,y,z),  view orientation (qw qx qy qz),  frustum height ratio ry        
        this->camera = new Camera(vec3(f[1], f[2], f[3]), Orientation(5, strings, f), f[4]);
    }


    else if (c == "ambient") {
        // syntax: ambient r g b
        // Sets the ambient color.  Note: This parameter is temporary.
        // It will be ignored once your raytracer becomes capable of
        // accurately *calculating* the true ambient light.        ; 
        this->ambientLight = vec3(f[1], f[2], f[3]);
    }

    else if (c == "brdf") {
        // syntax: brdf  r g b   r g b  alpha
        // later:  brdf  r g b   r g b  alpha  r g b ior
        // First rgb is Diffuse reflection, second is specular reflection.
        // third is beer's law transmission followed by index of refraction.
        // Creates a Material instance to be picked up by successive shapes
        currentMat = new Material(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), f[7]);
    }

    else if (c == "light") {
        // syntax: light  x y z   
        // The position of the light source
        // Light color defaults to white     
        currentMat = new Light(vec3(f[1], f[2], f[3]));
    }

    else if (c == "sphere") {
        // syntax: sphere x y z   r
        // Creates a Shape instance for a sphere defined by a center and radius
        this->shapes.push_back(new Sphere(vec3(f[1], f[2], f[3]), f[4], currentMat));

        if (currentMat->isLight()) {
            this->light = this->shapes.back(); // Store reference to light shape for sampling
            lightPos = vec3(f[1], f[2], f[3]);
        }
    }

    else if (c == "box") {
        // syntax: box bx by bz   dx dy dz
        // Creates a Shape instance for a box defined by a corner point and diagonal vector
        this->shapes.push_back(new Box(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), currentMat));
    }

    else if (c == "cylinder") {
        // syntax: cylinder bx by bz   ax ay az  r
        // Creates a Shape instance for a cylinder defined by a base point, axis vector, and radius
        this->shapes.push_back(new Cylinder(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), f[7], currentMat));
    }


    else if (c == "mesh") {
        // syntax: mesh   filename   tx ty tz   s   <orientation>
        // Creates many Shape instances (one per triangle) by reading
        // model(s) from filename. All triangles are rotated by a
        // quaternion (qw qx qy qz), uniformly scaled by s, and
        // translated by (tx ty tz) .
        mat4 modelTr = translate(vec3(f[2], f[3], f[4]))
            * scale(vec3(f[5], f[5], f[5]))
            * toMat4(Orientation(6, strings, f));
        ReadAssimpFile(strings[1], modelTr);
    }


    else {
        fprintf(stderr, "\n*********************************************\n");
        fprintf(stderr, "* Unknown command: %s\n", c.c_str());
        fprintf(stderr, "*********************************************\n\n");
    }
}


glm::vec3 Scene::SampleBRDF(const glm::vec3& normal, const glm::vec3& outDir, Material* mat)
{
	float e = myrandom(RNGen); // Uniform random number for sampling    
    
    float s = glm::length(mat->Kd) + glm::length(mat->Ks);
	float pd = glm::length(mat->Kd) / s;

	float e1 = myrandom(RNGen);
	float e2 = myrandom(RNGen) * 2.0f * PI;

    if (e < pd) {
		return SampleLobe(normal, sqrt(e1), e2); // Sample diffuse lobe
    }
    else {

        float alpha = mat->alpha;
        float cosTheta = pow(e1, 1.0f / (alpha + 2.0f));

        glm::vec3 m = SampleLobe(normal, cosTheta, e2);
        return 2.0f * dot(outDir, m) * m - outDir;
    }

}

float Scene::pdfBrdf(const glm::vec3& normal, const glm::vec3& inDir, const glm::vec3& outDir, Material* mat)
{
 
    // Normalize at the very start
    vec3 N = normalize(normal);
    vec3 L = normalize(inDir);
    vec3 V = normalize(outDir);
    vec3 m = normalize(L + V);

    float s = glm::length(mat->Kd) + glm::length(mat->Ks);
    if (s <= 0.0f) return 0.0f; // Guard against black materials

    float pd = glm::length(mat->Kd) / s;
    float pr = glm::length(mat->Ks) / s;

    float lDotN = dot(L, N);
    if (lDotN <= 0.0f) return 0.0f;

    // Diffuse PDF: cos(theta) / PI
    float Pd = lDotN / PI;

    // Specular PDF: (D * mDotN) / (4 * vDotM)
    float alpha = mat->alpha;
    float dm = DistributionPhong(N, m, alpha);
    float vDotM = dot(V, m);

    float Pr = (vDotM > 1e-3f) ? (dm * dot(m, N)) / (4.0f * vDotM) : 0.0f;

    return pd * Pd + pr * Pr;

}

float Scene::DistributionPhong(glm::vec3 normal, glm::vec3 m, float alpha)
{
    float mDotN = dot(m, normal);    
    //float mDotN = glm::clamp(glm::dot(m, normal), 0.0f, 1.0f);
	return xPlus(mDotN) * ( (alpha + 2) / (2.0f * PI) ) * pow(mDotN, alpha);
}

float Scene::GeometryPhong(glm::vec3 normal, glm::vec3 v, glm::vec3 m, float alpha) 
{
    float vDotM = dot(v, m);
	float vDotN = dot(v, normal);

    if( vDotM<= 0.0f || vDotN <= 0.0f) {
        return 0.0f;
    }

    if (vDotN >= 0.99999f) {
        return 1.0f;
    }
    
	float tanThetaV = sqrt(1 - vDotN * vDotN) / vDotN; // sin(thetaV) / cos(thetaV)

	if (tanThetaV == 0.0f) {
        return 1.0f;
    }
    
	float a = sqrt(alpha * 0.5f + 1.0f) / tanThetaV;

    if (a >= 1.6f) {
        return 1.0f;
	}

	return (3.535f * a + 2.181f * a * a) / (1.0f + 2.276f * a + 2.577f * a * a);

}

float Scene::GeometrySmith(glm::vec3 normal, glm::vec3 inDir, glm::vec3 outDir, float alpha)
{    
	
    glm::vec3 m = normalize(inDir + outDir);
    
    float G1_in = GeometryPhong(normal, inDir, m, alpha);
    
    fl

    fl
    
    float G1_in = GeometryPhong(normal, inDir, m, alpha);
    
    glm::vec3 m = normalize(inDir + outDir);
	
{    
float Scene::GeometrySmith(glm::vec3 normal, glm::vec3 inDir, glm::vec3 outDir, float alpha)

}

	return (3.535f * a + 2.181f * a * a) / (1.0f + 2.276f * a + 2.577f * a * a);

	}
        return 1.0f;
    if (a >= 1.6f) {

	float a = sqrt(alpha * 0.5f + 1.0f) / tanThetaV;
    
    }
        return 1.0f;
	if (tanThetaV == 0.0f) {

	float tanThetaV = sqrt(1 - vDotN * vDotN) / vDotN; // sin(thetaV) / cos(thetaV)
    
    }
        return 1.0f;
    if (vDotN >= 0.99999f) {

    }
        return 0.0f;
    if( vDotM<= 0.0f || vDotN <= 0.0f) {

	float vDotN = dot(v, normal);
    float vDotM = dot(v, m);
{
float Scene::GeometryPhong(glm::vec3 normal, glm::vec3 v, glm::vec3 m, float alpha) 

}
	return xPlus(mDotN) * ( (alpha + 2) / (2.0f * PI) ) * pow(mDotN, alpha);
    //float mDotN = glm::clamp(glm::dot(m, normal), 0.0f, 1.0f);
    float mDotN = dot(m, normal);    
{
float Scene::DistributionPhong(glm::vec3 normal, glm::vec3 m, float alpha)

}

    return pd * Pd + pr * Pr;

    float Pr = (vDotM > 1e-3f) ? (dm * dot(m, N)) / (4.0f * vDotM) : 0.0f;

    float vDotM = dot(V, m);
    float dm = DistributionPhong(N, m, alpha);
    float alpha = mat->alpha;
    // Specular PDF: (D * mDotN) / (4 * vDotM)

    float Pd = lDotN / PI;
    // Diffuse PDF: cos(theta) / PI

    if (lDotN <= 0.0f) return 0.0f;
    float lDotN = dot(L, N);

    float pr = glm::length(mat->Ks) / s;
    float pd = glm::length(mat->Kd) / s;

    if (s <= 0.0f) return 0.0f; // Guard against black materials
    float s = glm::length(mat->Kd) + glm::length(mat->Ks);

    vec3 m = normalize(L + V);
    vec3 V = normalize(outDir);
    vec3 L = normalize(inDir);
    vec3 N = normalize(normal);
    // Normalize at the very start
 
{
float Scene::pdfBrdf(const glm::vec3& normal, const glm::vec3& inDir, const glm::vec3& outDir, Material* mat)

}

    }
        return 2.0f * dot(outDir, m) * m - outDir;
        glm::vec3 m = SampleLobe(normal, cosTheta, e2);

        float cosTheta = pow(e1, 1.0f / (alpha + 2.0f));
        float alpha = mat->alpha;

    else {
    }
		return SampleLobe(normal, sqrt(e1), e2); // Sample diffuse lobe
    if (e < pd) {

	float e2 = myrandom(RNGen) * 2.0f * PI;
	float e1 = myrandom(RNGen);

	float pd = glm::length(mat->Kd) / s;
    float s = glm::length(mat->Kd) + glm::length(mat->Ks);
    
	float e = myrandom(RNGen); // Uniform random number for sampling    
{
glm::vec3 Scene::SampleBRDF(const glm::vec3& normal, const glm::vec3& outDir, Material* mat)


}
    }
        fprintf(stderr, "*********************************************\n\n");
        fprintf(stderr, "* Unknown command: %s\n", c.c_str());
        fprintf(stderr, "\n*********************************************\n");
    else {


    }
        ReadAssimpFile(strings[1], modelTr);
            * toMat4(Orientation(6, strings, f));
            * scale(vec3(f[5], f[5], f[5]))
        mat4 modelTr = translate(vec3(f[2], f[3], f[4]))
        // translated by (tx ty tz) .
        // quaternion (qw qx qy qz), uniformly scaled by s, and
        // model(s) from filename. All triangles are rotated by a
        // Creates many Shape instances (one per triangle) by reading
        // syntax: mesh   filename   tx ty tz   s   <orientation>
    else if (c == "mesh") {


    }
        this->shapes.push_back(new Cylinder(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), f[7], currentMat));
        // Creates a Shape instance for a cylinder defined by a base point, axis vector, and radius
        // syntax: cylinder bx by bz   ax ay az  r
    else if (c == "cylinder") {

    }
        this->shapes.push_back(new Box(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), currentMat));
        // Creates a Shape instance for a box defined by a corner point and diagonal vector
        // syntax: box bx by bz   dx dy dz
    else if (c == "box") {

    }
        }
            lightPos = vec3(f[1], f[2], f[3]);
            this->light = this->shapes.back(); // Store reference to light shape for sampling
        if (currentMat->isLight()) {

        this->shapes.push_back(new Sphere(vec3(f[1], f[2], f[3]), f[4], currentMat));
        // Creates a Shape instance for a sphere defined by a center and radius
        // syntax: sphere x y z   r
    else if (c == "sphere") {

    }
        currentMat = new Light(vec3(f[1], f[2], f[3]));
        // Light color defaults to white     
        // The position of the light source
        // syntax: light  x y z   
    else if (c == "light") {

    }
        currentMat = new Material(vec3(f[1], f[2], f[3]), vec3(f[4], f[5], f[6]), f[7]);
        // Creates a Material instance to be picked up by successive shapes
        // third is beer's law transmission followed by index of refraction.
        // First rgb is Diffuse reflection, second is specular reflection.
        // later:  brdf  r g b   r g b  alpha  r g b ior
        // syntax: brdf  r g b   r g b  alpha
    else if (c == "brdf") {

    }
        this->ambientLight = vec3(f[1], f[2], f[3]);
        // accurately *calculating* the true ambient light.        ; 
        // It will be ignored once your raytracer becomes capable of
        // Sets the ambient color.  Note: This parameter is temporary.
        // syntax: ambient r g b
    else if (c == "ambient") {


    }
        this->camera = new Camera(vec3(f[1], f[2], f[3]), Orientation(5, strings, f), f[4]);
        // Eye position (x,y,z),  view orientation (qw qx qy qz),  frustum height ratio ry        
        // syntax: camera x y z   ry   <orientation spec>
    else if (c == "camera") {

    }
        height = int(f[2]);
        width = int(f[1]);
        // syntax: screen width height        
    if (c == "screen") {

    std::string c = strings[0];
    if (strings.size() == 0) return;
{
    const std::vector<float>& f)
void Scene::Command(const std::vector<std::string>& strings,

}
    return q;
    }
        }
            i += 4;
            q *= angleAxis(f[i + 0] * Radians, normalize(vec3(f[i + 1], f[i + 2], f[i + 3])));
        else if (c == "a") {
        }
            i += 4;
            q *= quat(f[i + 0], f[i + 1], f[i + 2], f[i + 3]);
        else if (c == "q") {
            q *= angleAxis(f[i++] * Radians, Zaxis());
        else if (c == "z")
            q *= angleAxis(f[i++] * Radians, Yaxis());
        else if (c == "y")
            q *= angleAxis(f[i++] * Radians, Xaxis());
        if (c == "x")
        std::string c = strings[i++];
    while (i < strings.size()) {
    quat q(1, 0, 0, 0); // Unit quaternion
{
    const std::vector<float>& f)
    const std::vector<std::string>& strings,
quat Orientation(int i,

}
    }
        }
            exit(-1);
            printf("\nRead error on file %s:\n  %s\n\n", path.c_str(), stbi_failure_reason());
        if (!image) {
        printf("%d %d %d %s\n", depth, width, height, path.c_str());
        image = stbi_load(path.c_str(), &width, &height, &depth, 4);
        stbi_set_flip_vertically_on_load(true);
        // Read image, and check for success
    else {
    }
        exit(-1);
        std::cerr << "Texture file not found: " << path << std::endl;
    if (find_it.fail()) {
    std::ifstream find_it(path.c_str());
    // Does the file exist?

    }
        path.replace(path.find(bs), 1, fs);
    while (path.find(bs) != std::string::npos) {
    std::string fs = "/";
    std::string bs = "\\";
    std::string path = bpath;
    // Replace backslashes with forward slashes -- Good for Linux, and maybe Windows?
{
Texture::Texture(const std::string& bpath) : id(0)

}
    }
        this->shapes.push_back(new Triangle(v0, v1, v2, meshMat));
        // Create a Triangle shape and add it to the scene with the scene file's material

        vec3 v2 = mesh->vertices[tri[2]].pnt;
        vec3 v1 = mesh->vertices[tri[1]].pnt;
        vec3 v0 = mesh->vertices[tri[0]].pnt;
        // tri is an ivec3 containing three vertex indices
    for (const auto& tri : mesh->triangles) {
    // Convert mesh data into individual Triangle shapes and add to scene

    Material* meshMat = currentMat;
    // This allows the scene file to control mesh appearance
    // Use the currentMat from the scene file (set by brdf command)
{
void Scene::triangleMesh(MeshData* mesh)

}
    this->accelStruc = new AccelerationBvh(this->shapes);
{
void Scene::Finit()

// Call myrandom(RNGen) to get a uniformly distributed random number in [0,1].
std::uniform_real_distribution<> myrandom(0.0, 1.0);
std::mt19937_64 RNGen(device());
std::random_device device;
#include <random>
// A good quality *thread-safe* Mersenne Twister random number generator.

#include "stb_image.h"
#define STBI_FAILURE_USERMSG
#define STB_IMAGE_IMPLEMENTATION

#include "raytrace.h"

////////////////////////////////////////////////////////////////////////
// Provides the framework for a raytracer.
//////////////////////////////////////////////////////////////////////

Project 3

Reflections

Replacing the diffuse-only BRDF with a Cook-Torrance microfacet model (Phong NDF, Smith shadowing, Schlick Fresnel), and fixed-weight blending with the balance heuristic.

Course: CS 500, Ray Tracing
Submitted: 14 March 2026
Status: Complete

Project overview

This project extends the explicit path tracer from Project 2 to support physically based materials via a microfacet BRDF model. Project 2 established the core path tracing infrastructure (implicit and explicit sampling, Russian Roulette termination, and fixed-weight blending) but evaluated surfaces using only a Lambertian (purely diffuse) BRDF. Project 3 replaces that with a full Cook-Torrance microfacet model, which renders glossy and reflective surfaces physically accurately.

The framework introduces a toggle (ENABLE_BRDF) that switches between the legacy diffuse-only mode and the new full BRDF evaluation, keeping all other path tracing infrastructure identical. The linear blending of weights is also upgraded from the fixed 0.5f split of Project 2 to the Balance Heuristic (discussed as part of next project), which further reduces variance when the BRDF lobe and explicit light sampling strategies differ in quality.

Test scenes

Three scenes of increasing complexity validate the implementation, from a single reflective sphere on a table to a full scene with high-polygon meshes and multiple inter-reflecting spheres.

BRDF model overview

The implemented BRDF is a microfacet specular model combined with a Lambertian diffuse term, following the Cook-Torrance formulation. The full BRDF evaluated in EvalScattering is:

f(L, V) = Kd/π * (N·L) + D(m) * G(L, V) * F(L, m) / (4 * (N·V))

Where:

D: Phong Normal Distribution Function (NDF), which models the statistical distribution of microfacet normals.
G: Smith geometry / shadowing-masking term, which accounts for mutual occlusion of microfacets.
F: Fresnel reflectance (Schlick approximation), which models the fraction of light reflected vs. refracted at a surface.
m: half-vector, normalize(L + V), the microfacet normal that reflects L toward V.
α (alpha): roughness parameter, which controls the width of the specular lobe.

Phong normal distribution function (NDF)

The NDF D(m) describes what fraction of microfacets are oriented with their normal aligned to the half-vector m. A Phong-based NDF is used:

float Scene::DistributionPhong(glm::vec3 normal, glm::vec3 m, float alpha)
{
    float mDotN = dot(m, normal);
    return xPlus(mDotN) * ((alpha + 2) / (2.0f * PI)) * pow(mDotN, alpha);
}

As alpha increases, the distribution becomes tighter and more mirror-like. At low alpha values, the distribution is broad and the material appears diffuse-dominant; at high values, it produces sharp, concentrated specular highlights.

Geometry / shadowing-masking term

The geometry term G(L, V) accounts for self-shadowing and masking between microfacets. The implementation uses a Smith factorization, computing an independent G1 term for both the incident and outgoing directions and multiplying them together:

float Scene::GeometrySmith(glm::vec3 normal, glm::vec3 inDir, glm::vec3 outDir, float alpha)
{
    glm::vec3 m = normalize(inDir + outDir);
    float G1_in  = GeometryPhong(normal, inDir,  m, alpha);
    float G1_out = GeometryPhong(normal, outDir, m, alpha);
    return G1_in * G1_out;
}

Each G1 factor is evaluated using the Phong-specific approximation:

float Scene::GeometryPhong(glm::vec3 normal, glm::vec3 v, glm::vec3 m, float alpha)
{
    float a = sqrt(alpha * 0.5f + 1.0f) / tanThetaV;
    // Rational approximation for a < 1.6
    return (3.535f * a + 2.181f * a * a) / (1.0f + 2.276f * a + 2.577f * a * a);
}

Rays that are nearly grazing (small N·V) produce small G values, correctly attenuating the specular contribution to avoid energy non-conservation at extreme viewing angles.

Fresnel reflectance: Schlick approximation

The Fresnel term models the angle-dependent ratio of reflected light. At normal incidence, a surface reflects a fraction determined by Ks (the specular color). As the viewing angle approaches grazing, reflectance tends toward 1 (total reflection). The Schlick approximation is used:

glm::vec3 Scene::FresnelSchlick(glm::vec3 inDir, glm::vec3 m, Material* mat)
{
    vec3 Ks = mat->Ks;
    float d = dot(inDir, m);
    return Ks + (vec3(1.0f) - Ks) * pow((1.0f - abs(d)), 5.0f);
}

Ks acts as the base reflectance at normal incidence. For metals, Ks is a tinted color reflecting the metal’s wavelength-dependent reflectance. For dielectrics, it is a low-value grey. The pow(..., 5) term gives the physically plausible rapid rise in reflectance at grazing angles.

Full BRDF evaluation: `EvalScattering`

The full BRDF is assembled from the D, G, and F terms in EvalScattering:

vec3 Scene::EvalScattering(glm::vec3 normal, glm::vec3 inDir, glm::vec3 outDir, Material* mat)
{
    float D = DistributionPhong(N, m, alpha);
    float G = GeometrySmith(N, L, V, alpha);
    vec3  F = FresnelSchlick(L, m, mat);

    float denom = 4.0f * glm::max(outDotN, 0.0001f);

    vec3 diffuse  = mat->Kd / PI;
    vec3 specular = (D * G * F) / denom;

    return inDotN * (diffuse) + specular;
}

The denominator clamp (glm::max(outDotN, 0.0001f)) prevents numerical blow-up at extreme grazing configurations. The final result combines the diffuse and specular lobes, multiplied by the cosine factor inDotN = N·L to account for the projected surface area facing the light.

BRDF sampling: `SampleBRDF`

With a mixed diffuse-specular BRDF, the ray direction at each bounce must be sampled proportionally from both lobes. The sampling probability is determined by the relative magnitudes of Kd and Ks:

glm::vec3 Scene::SampleBRDF(const glm::vec3& normal, const glm::vec3& outDir, Material* mat)
{
    float s  = glm::length(mat->Kd) + glm::length(mat->Ks);
    float pd = glm::length(mat->Kd) / s;

    float e  = myrandom(RNGen);
    float e1 = myrandom(RNGen);
    float e2 = myrandom(RNGen) * 2.0f * PI;

    if (e < pd) {
        return SampleLobe(normal, sqrt(e1), e2);  // Cosine-weighted diffuse
    } else {
        float cosTheta = pow(e1, 1.0f / (alpha + 2.0f));
        glm::vec3 m = SampleLobe(normal, cosTheta, e2);
        return 2.0f * dot(outDir, m) * m - outDir;  // Phong specular reflection
    }
}

With probability pd, a cosine-weighted diffuse direction is chosen. Otherwise, a specular microfacet normal m is sampled from the Phong NDF, and the outgoing direction is computed as the specular reflection of outDir about m. Sampling in proportion to Kd and Ks keeps either lobe from being over- or under-sampled regardless of material composition.

BRDF PDF: `pdfBrdf`

To correctly weight the Monte Carlo estimate, the PDF must account for both the diffuse and specular sampling strategies:

float Scene::pdfBrdf(const glm::vec3& normal, const glm::vec3& inDir, const glm::vec3& outDir, Material* mat)
{
    float pd = glm::length(mat->Kd) / s;
    float pr = glm::length(mat->Ks) / s;

    // Diffuse PDF: cos(theta) / PI
    float Pd = lDotN / PI;

    // Specular PDF: (D * mDotN) / (4 * vDotM)
    float dm = DistributionPhong(N, m, alpha);
    float Pr = (dm * dot(m, N)) / (4.0f * vDotM);

    return pd * Pd + pr * Pr;
}

The mixed PDF is a weighted sum of both lobe PDFs using the same pd / pr probabilities as the sampler, which keeps the Monte Carlo weight f / p unbiased across all path configurations.

Equal-weight sampling

When both explicit shadow rays and implicit BRDF-sampled rays are active simultaneously, both strategies can contribute light from the same source. To combine them without double-counting, the project baseline assigns each strategy an equal fixed weight of 0.5f:

// Explicit shadow ray contribution
C += 0.5f * W * (f / p) * hit.shape->mat->Kd;

// Implicit BRDF ray that happens to hit the light
C += W * nextHit.shape->mat->Kd * 0.5f;

This is simple, unbiased, and sufficient for purely diffuse scenes where both strategies contribute roughly equally at every surface point. On some specular and glossy surfaces, though, the two strategies are far from equal in quality: the BRDF sampler concentrates rays tightly around the specular lobe, while the light sampler draws directions uniformly toward the light regardless of the surface’s reflective peak. Assigning each a flat 0.5f weight lets the weaker strategy inject disproportionate variance into the estimate, which shows up as bright noise speckles and fireflies concentrated near reflective surfaces.

Multiple importance sampling: the Balance Heuristic

To resolve fireflies introduced by equal weighting on specular surfaces, the fixed 0.5f weights are replaced with the Balance Heuristic. The root cause of fireflies is a BRDF-sampled ray that accidentally strikes a bright light source from a direction with very low BRDF probability: the Monte Carlo weight f/p becomes large because p is small, spiking that pixel far above its neighbors.

For the explicit shadow ray connection:

float p_brdf = pdfBrdf(P.normal, inDir, outDir, P.shape->mat);
float weight = p / (p + p_brdf);  // p = light PDF
C += weight * W * (f / p) * hit.shape->mat->Kd;

For the implicit BRDF-sampled ray that hits the light:

float p_brdf  = pdfBrdf(N, inDir, outDir, P.shape->mat) * RUSSIAN_ROULETTE_PROB;
float p_light = PdfLight(nextHit) / GeometryFactor(P, nextHit);
float weight  = p_brdf / (p_brdf + p_light);
C += W * nextHit.shape->mat->Kd * weight;

The Balance Heuristic suppresses this by incorporating the light sampling PDF into the denominator alongside the BRDF PDF, acknowledging that the explicit strategy could have reached that same direction with much higher probability:

weight_brdf = p_brdf / (p_brdf + p_light)

This larger denominator pulls the weight down for exactly the rays that cause fireflies: those sampled by the BRDF at low probability but easily reachable by the light sampler. On a diffuse surface where p_brdf and p_light are comparable in magnitude, the heuristic recovers weights close to 0.5f, matching the baseline behavior. Fireflies on specular surfaces are suppressed with no additional bias introduced into the estimator.

Comparison: equal-weight vs. Balance Heuristic

Passes	Equal-Weight Sampling	Balance Heuristic (MIS)
256	Click to zoom	Click to zoom
512	Click to zoom	Click to zoom
2048	Click to zoom	Click to zoom
4096	Click to zoom	Click to zoom

Test scenes

Three test scenes of increasing complexity validate the BRDF implementation. Each scene is rendered at pass counts growing in powers of 2 from 8 to 4096, which shows how quickly the specular model converges under the Balance Heuristic MIS.

Scene 1: reflective table with a single sphere

The baseline test scene has a rectangular table with a reflective top surface and metallic edges, with a single highly reflective sphere at its center. Colored diffuse walls (green at the back, blue on the side, neutral on the floor) surround the scene and contribute indirect bounce lighting. A spherical area light positioned directly above the sphere provides the primary illumination, and the table legs are thin white cylinders.

The metallic table edges pick up green and blue tints from the surrounding walls, which is colored indirect transport. The sphere, being nearly mirror-like, prominently reflects the overhead light source as a sharp highlight and captures the colored environment across its surface. A faint reflection of the sphere itself is also visible on the table top, which confirms specular inter-reflection between the two surfaces.

Test images

256 Passes	512 Passes
Click to zoom	Click to zoom
2048 Passes	4096 Passes
Click to zoom	Click to zoom

Scene 2: table with sphere, rabbit, and dwarf

The second scene extends Scene 1 by placing two high-polygon mesh objects on the table surface: a Stanford Rabbit on the right and an armored Dwarf character on the left. Both meshes use a diffuse-dominant white material with minimal specular response, which contrasts with the highly reflective table and sphere. The colored-wall environment and overhead light from Scene 1 are unchanged.

The rabbit and dwarf accumulate soft colored indirect illumination from the surrounding walls, with green and blue tints visible across their matte surfaces. Their geometry casts shadows onto the reflective table top, producing visible occlusion patterns in the table’s specular reflection. The sphere’s curved surface also reflects both figures, though with significant angular distortion due to its curvature: the rabbit’s ear in particular is identifiable in the reflection. Faint specular reflections of both meshes are also visible on the table surface, consistent with the inter-reflection behavior established in Scene 1.

Test images

256 Passes	512 Passes
Click to zoom	Click to zoom
2048 Passes	4096 Passes
Click to zoom	Click to zoom

Scene 3: table with meshes and multiple reflective spheres

The third scene is the most complex, retaining the rabbit and dwarf from Scene 2 while replacing the single central sphere with two reflective spheres positioned across the table surface. Each sphere reflects the other, the mesh objects, and the colored environment simultaneously, forming a network of multi-bounce specular paths that places the highest demand on the path tracer.

Specular-to-specular transport between the two spheres requires multiple bounces to resolve correctly. The rabbit and dwarf, now flanked by spheres on both sides, accumulate colored specular bounce light from the adjacent sphere surfaces. The rabbit’s silhouette and the dwarf’s form are both identifiable as distorted reflections on the second sphere. A close inspection of the first sphere reveals the reflection of the second sphere nested within it, and within that reflection a further recursive specular highlight, which confirms that the renderer correctly propagates specular-to-specular paths across multiple bounces.

Test images

256 Passes	512 Passes
Click to zoom	Click to zoom
2048 Passes	4096 Passes
Click to zoom	Click to zoom

Summary

Microfacet BRDF: the renderer now evaluates a full Cook-Torrance specular model (Phong NDF + Smith G + Schlick F) alongside Lambertian diffuse, which covers a wide range of real-world surface appearances.
Mixed sampling: SampleBRDF probabilistically samples from either the diffuse or specular lobe proportional to |Kd| and |Ks|, and pdfBrdf returns the corresponding mixed PDF, keeping the Monte Carlo estimator unbiased.
MIS, equal weights to Balance Heuristic: combining explicit and implicit light contributions with a fixed 0.5f weight sometimes produces fireflies near specular surfaces. The Balance Heuristic replaces it and weights each strategy in proportion to its PDF at the sampled direction, which suppresses high-frequency noise on reflective geometry while preserving unbiasedness.
Framework extensibility: the clean separation of EvalScattering, pdfBrdf, and SampleBRDF into per-material, direction-aware overloads makes it straightforward to incorporate additional material types (e.g., glass / transmission, anisotropic microfacets) in future projects.

Project overview

Test scenes

BRDF model overview

Phong normal distribution function (NDF)

Geometry / shadowing-masking term

Fresnel reflectance: Schlick approximation

Full BRDF evaluation: EvalScattering

BRDF sampling: SampleBRDF

BRDF PDF: pdfBrdf

Equal-weight sampling

Multiple importance sampling: the Balance Heuristic

Comparison: equal-weight vs. Balance Heuristic

Test scenes

Scene 1: reflective table with a single sphere

Test images

Scene 2: table with sphere, rabbit, and dwarf

Test images

Scene 3: table with meshes and multiple reflective spheres

Test images

Summary

Full BRDF evaluation: `EvalScattering`

BRDF sampling: `SampleBRDF`

BRDF PDF: `pdfBrdf`