Introduction

This implementation focuses on physically-based rendering through radiometric principles. This article explains core radiometry concepts, presents key code implementations with references, and details the microfacet shading model. While the model produces correct results, significant noise artifacts appear in the path tracing output that remain unresolved at this stage.

Radiometric Fundamentals

The following provides an intuitive understanding of radiometric quantities:

• Radiant Energy Q [J = Joule]: The total energy of electromagnetic radiation

Radiant energy represents pure energy. Consider a comparison with displacement: starting at 0m and ending at 100m tells us nothing useful by itself—we need velocity to describe motion meaningfully. Similarly, a 10W lamp might emit 1000J over 10 hours while a 100W lamp emits the same 1000J in 1 hour. We intuitively sense the 10W lamp is dimmer, leading to the concept of power.

• Radiant Flux (Power) Φ = dQ/dt [W=Watt] [lm=lumen]: Energy per unit time

Using a particle metaphor, radiant flux describes how many energy "particles" are emitted per second. Direction vectors represent emission directions, with vector magnitude representing particle count. Note that in physics, vectors have only magnitude and direction—no fixed origin.

• Radiant Intensity I(ω) = dΦ/dω: Power per unit solid angle
• Solid Angle Ω = A/r²: Ratio of subtended area on a sphere to radius squared

For a point source emitting particles in all directions, radiant intensity counts the particles passing through a conic solid angle per second. The differential notation indicates limiting processes: as area approaches an infinitesimally small point, the solid angle becomes an infinitesimal cone.

• Irradiance
  E(x) = dΦ(x)/dA
  [W/m²] [lm/m² = lux]

Irradiance measures incoming power per unit surface area from all directions at a point, derived through limiting processes.

• Radiance
  L(p, ω) = d²Φ(p,ω)/(dω dA cosθ)
  L(p, ω) = dE(p)/(dω cosθ)
  L(p, ω) cosθ dω = dE(p)

Radiance characterizes light rays during propagation (between surface interactions). When light strikes a surface, energy loss depends on the angle between the surface normal and ray direction—Lambert's Cosine Law. A ray carrying 100 units of power might deliver only 50 units due to angle effects, but the ray itself maintains 100 units. When receiving irradiance, the inverse cosine correction restores the full radiance value.

Implementation Details

Path Tracing Pseudocode

shade(p, wo)
    // Sample light source uniformly
    sample light at position x (pdf_light = 1 / A)
    cast ray from p to x
    if ray is unobstructed:
        L_dir = L_i * f_r * cos_theta * cos_theta_x / |x-p|^2 / pdf_light
    
    L_indir = 0.0
    // Russian roulette termination
    if random() < P_RR:
        sample hemisphere uniformly toward wi (pdf_hemi = 1 / 2pi)
        trace ray r(p, wi)
        if r hits non-emissive object at q:
            L_indir = shade(q, wi) * f_r * cos_theta / pdf_hemi / P_RR
    
    return L_dir + L_indir

Uniform Random Sampling in Triangles

void Triangle::Sample(Intersection& position, float& probability) {
    float u = std::sqrt(random_generator::next_float());
    float v = random_generator::next_float();
    position.coords = v0 * (1.0f - u) + v1 * (u * (1.0f - v)) + v2 * (u * v);
    position.normal = this->normal;
    probability = 1.0f / surface_area;
}

Reference: sIGGRAPH 2002 Tutorial Notes

Multithreaded Rendering

int samples_per_pixel = 64;
std::vector<std::thread> worker_threads;
std::vector<Vector3f> pixel_samples(samples_per_pixel);

// Launch threads for each ray sample
Ray primary_ray(eye_position, ray_direction);
for (int idx = 0; idx < samples_per_pixel; ++idx) {
    worker_threads.emplace_back(
        std::thread(&Renderer::RenderThread, this,
                   std::ref(pixel_samples[idx]),
                   std::ref(scene),
                   primary_ray)
    );
}

// Synchronize all threads
for (auto& t : worker_threads) {
    t.join();
}

// Accumulate samples with proper averaging
for (int idx = 0; idx < samples_per_pixel; ++idx) {
    framebuffer[pixel_index] += pixel_samples[idx] / samples_per_pixel;
}
worker_threads.clear();

BRDF Evaluation

Using Lambertian diffuse material with two key assumptions:

1. Incident radiance is uniform across all directions
2. Outgoing radiance is uniform across all output directions

When ρ = 1 (albedo), energy conservation holds: one incoming ray carrying energy x spreads uniformly across all directions, but infinitely many incoming rays each carrying energy x combine to produce x in any given outgoing direction.

Vector3f Material::eval(const Vector3f& incident_dir,
                        const Vector3f& outgoing_dir,
                        const Vector3f& surface_normal) {
    switch (material_type) {
        case DIFFUSE:
        {
            float cos_angle = dot(surface_normal, outgoing_dir);
            if (cos_angle > 0.0f) {
                // Transform incident direction to outward convention
                return microfacetBRDF(-incident_dir, outgoing_dir,
                                     surface_normal, 0.5f, 0.2f,
                                     diffuse_color, 0.5f);
            }
            return Vector3f(0.0f);
        }
    }
}

Hemisphere Sampling

Convert polar coordinates (θ, φ) to Cartesian coordinates using trigonometric functions, then transform to world space. The resulting wi vector is a unit vector. For arbitrary normal orientations, check alignment: if the sample direction opposes the normal, invert the vector.

Vector3f Material::sample(const Vector3f& incident_dir,
                          const Vector3f& surface_normal) {
    switch (material_type) {
        case DIFFUSE:
        {
            // Uniform hemisphere sampling via cosine-weighted method
            float u1 = random_generator::next_float();
            float u2 = random_generator::next_float();
            float cos_theta = 1.0f - 2.0f * u1;
            float sin_theta = std::sqrt(std::max(0.0f, 1.0f - cos_theta * cos_theta));
            float phi = 2.0f * M_PI * u2;
            Vector3f local_sample(sin_theta * std::cos(phi),
                                 sin_theta * std::sin(phi),
                                 cos_theta);
            return transform_to_world(local_sample, surface_normal);
        }
    }
}

Microfacet Shading Model

The GGX distribution function D(m) quantifies what fraction of microfacets orient parallel to the half-vector m, ranging from 0 to 1. The Smith masking-shadowing function G combines two遮挡 components from both view and light directions.

GGX Distribution

Two equivalent formulations exist, where θm is the angle between the half-vector and surface normal:

DGGX(m) = α² / (π((n·m)²(α² - 1) + 1)²)

DGGX(m) = α² / (π cos⁴θm (α² + tan²θm)²)

Where:

n·m = cos θm

tan θm = sin θm / cos θm

// Walter 2007 Microfacet Model
// Parameters:
//   incident: incoming light direction
//   view: view direction
//   macro_normal: surface normal
//   metallic: dielectric(0.0f) to metal(1.0f) factor
//   roughness: GGX distribution parameter
//   reflectance: controls F0 for Fresnel (e.g., 0.04 for glass)
//   base_color: diffuse contribution and metal F0
Vector3f microfacetBRDF(Vector3f Light, Vector3f View, Vector3f N,
                        float metallic, float roughness,
                        Vector3f baseColor, float reflectance);

Vector3f fresnelSchlick(float cos_angle, Vector3f F0);
float GGXDistribution(float NdotH, float roughness);
float SmithMasking(float NdotV, float NdotL, float roughness);
float SmithGGXSchlick(float NdotV, float roughness);

Vector3f Material::eval(const Vector3f& incident_dir,
                        const Vector3f& outgoing_dir,
                        const Vector3f& surface_normal) {
    switch (material_type) {
        case DIFFUSE:
        {
            float cos_angle = dot(surface_normal, outgoing_dir);
            if (cos_angle > 0.0f) {
                return microfacetBRDF(-incident_dir, outgoing_dir,
                                     surface_normal, 0.5f, 0.2f,
                                     diffuse_color, 0.5f);
            }
            return Vector3f(0.0f);
        }
    }
}

Vector3f Material::microfacetBRDF(Vector3f Light, Vector3f View,
                                  Vector3f N, float metallic,
                                  float roughness, Vector3f baseColor,
                                  float reflectance) {
    // Half-angle vector
    Vector3f H = normalize(Light + View);

    float NdotL = saturate(dot(N, Light));
    float NdotV = saturate(dot(N, View));
    float NdotH = saturate(dot(N, H));
    float VdotH = saturate(dot(View, H));

    Vector3f F0(0.16f * (reflectance * reflectance));
    F0 = mix(F0, baseColor, metallic);

    Vector3f F = fresnelSchlick(VdotH, F0);
    float D = GGXDistribution(NdotH, roughness);
    float G = SmithMasking(NdotV, NdotL, roughness);

    Vector3f specular = (F * D * G) / (4.0f * NdotV * NdotL);

    // Diffuse component (energy-conserving)
    Vector3f diffuse_color = baseColor;
    diffuse_color = diffuse_color * (Vector3f(1.0f) - F);
    diffuse_color = diffuse_color * (1.0f - metallic);
    Vector3f diffuse = diffuse_color / M_PI;

    return diffuse + specular;
}

Vector3f Material::fresnelSchlick(float cos_angle, Vector3f F0) {
    return F0 + (Vector3f(1.0f) - F0) * std::pow(1.0f - cos_angle, 5.0f);
}

float Material::GGXDistribution(float NdotH, float roughness) {
    float alpha = roughness * roughness;
    float alpha_squared = alpha * alpha;
    float NdotH_squared = NdotH * NdotH;
    float denominator = NdotH_squared * (alpha_squared - 1.0f) + 1.0f;
    return alpha_squared / (M_PI * denominator * denominator);
}

float Material::SmithMasking(float NdotV, float NdotL, float roughness) {
    return SmithGGXSchlick(NdotV, roughness) *
           SmithGGXSchlick(NdotL, roughness);
}

float Material::SmithGGXSchlick(float NdotV, float roughness) {
    float alpha = roughness * roughness;
    float k = alpha * 0.5f;
    return NdotV / (NdotV * (1.0f - k) + k);
}

References:

• Mitsuba2 Documentation: Rough Conductor Material
• Walter et al., "Microfacet Models for Refraction through Rough Surfaces" (EGSR 2007)
• Youtube: Microfacet BRDF - Theory and Implementation of Basic PBR Materials