Raytracer

I started from Peter Shirley's Ray Tracing in One Weekend series, and the camera, the scatter() material interface, the Lambertian, metal and dielectric materials and the sampling helpers follow it closely. I extended it with planes and triangles, explicit sampling of point and area lights, normal maps, JSON scene files and an interactive progressive viewer.

In November 2026 I came back for version 2: a BVH, OpenMP rendering, Russian roulette, a WebAssembly build, and a careful pass over the sampling and shading maths, which turned up five bugs. This page covers what was wrong, why it showed up in the image, and how I checked each fix.

What it does now

• Iterative path tracing with next-event estimation: one shadow-ray sample per light at each diffuse hit, with the Lambertian BRDF albedo/π, and emission counted only where light sampling didn't already account for it.
• Russian roulette after 3 bounces, with the survival probability set from the path throughput and capped at 0.95.
• Materials: Lambertian (cosine-weighted bounce), textured Lambertian with a tangent-space normal map, metal with fuzz, dielectric (Snell refraction, total internal reflection, Schlick Fresnel), emissive.
• Spheres, infinite planes, Möller-Trumbore triangles and Wavefront OBJ meshes, smooth-shaded from vertex normals when the scene asks for it.
• A bounding volume hierarchy built with a 16-bin surface area heuristic, stored as one flat node array and traversed iteratively, near child first.
• OpenMP over image rows, with a thread_local PCG32 generator reseeded per pixel, so the image doesn't depend on the thread count.
• A headless command-line renderer (PNG or PPM, progressive snapshots), the GLFW viewer, and a wasm32-wasi build. The browser version on my home page is that same C++, one instance per Web Worker.

Fixes

Every before/after pair below is the same scene rendered from the commit before and the commit after the fix, with the headless renderer compiled to WebAssembly. Old commits only had the viewer, so they were rendered through a small driver that calls their own scene loader, camera and renderer.

The random number range was frozen by a static

Fixed in 351cc24

A function-local static is initialised once, the first time control passes through it. The distribution object stores its range, so it kept the min and max of the first call forever. The first call in a render is random_float() for pixel jitter, so the range was frozen at [0, 1), and every later random_float(-1, 1) silently returned [0, 1). random_int had the same problem.

float random_float(float min, float max) {
    static std::uniform_real_distribution<float> distribution(min, max);
    static std::mt19937 generator;
    return distribution(generator);
}

Three things are built from random_float(-1, 1) by rejection sampling: points in the unit sphere, points in the unit disk and unit vectors. All three came from one octant or quadrant. The Lambertian bounce is n + random_unit_vector(), which is cosine-weighted only if the unit vector is uniform over the whole sphere, so indirect light came from one side. Metal fuzz spread over an eighth of the ball, and the depth-of-field lens offset covered a quarter of the aperture.

Now only the engine is persistent and the range is applied on every call. Using 24 random bits means every value is exactly representable as a float, so the result never rounds up to max.

// One engine per thread. Only the engine is persistent: the [min, max)
// range is applied on every call (a static distribution object would keep
// the range of the first call).
thread_local Pcg32 tls_rng;
...
// Uniform float in [min, max). The range is applied on every call.
float random_float(float min, float max) {
    // 24 random bits -> exactly representable float in [0, 1)
    float unit = static_cast<float>(tls_rng.next() >> 8) * (1.0f / 16777216.0f);
    return min + (max - min) * unit;
}

The TBN matrix was transposed

Fixed in 3c6446a

GLM matrices are column-major: the nine-float constructor fills column 0 first. Typing T, B and N in component by component built a matrix whose rows were T, B and N, which is the world-to-tangent transform. Even an unperturbed normal (0, 0, 1) came out as (T.z, B.z, N.z) instead of N. The fix is glm::mat3(rec.tangent, rec.bitangent, rec.normal), three column vectors.

Area lights were offset by half their size

Fixed in 8d651dd

The light stores a centre and two edge vectors but sampled center + s·u + t·v with s and t in [0, 1), so the light sat (u + v) / 2 away from where the scene file put it. Sampling s and t in [-0.5, 0.5) fixes it, and the pdf doesn't change because the area is the same.

Direct light: the missing 1/π, and only for diffuse surfaces

Fixed in 704fd0c

Next-event estimation estimates the direct part of the rendering equation with one light sample: L ≈ f · L_e · cos θ / p. For a Lambertian surface f = albedo / π. The old code computed albedo · L · cos / p, which is π times too bright. It also added that term for mirrors, glass and fuzzy metal. Their BRDF is a delta (or, for fuzzy metal, one this code can't evaluate), so a random light sample contributes nothing and their light has to come from the scattered ray. The term made glass and chrome look partly matte.

    const float inv_pi = 1.0f / glm::pi<float>();
    const Color brdf = albedo * inv_pi;
    Color direct(0.0f);

    for (const auto& light : scene.lights()) {
        auto sample = light->sample(rec.point);
        if (!(sample.pdf > 0.0f)) continue;

        float cos_theta = glm::dot(shading_normal, sample.direction);
        // Skip lights behind the shading normal or behind the actual surface
        if (cos_theta <= 0.0f || glm::dot(rec.geometric_normal, sample.direction) <= 0.0f) continue;

        // Shadow ray
        core::Ray shadow_ray(rec.point, sample.direction);
        geometry::HitRecord shadow_rec;
        if (scene.hit(shadow_ray, 0.001f, sample.distance - 0.001f, shadow_rec)) continue;

        direct += brdf * sample.radiance * cos_theta / sample.pdf;
    }

Light counted twice

Fixed in f615451

At every diffuse hit the integrator sampled the lights and also followed a bounce ray, adding the emission of whatever it hit. If an emitter is both scene geometry and a sampled light, it gets counted by both. None of the shipped scenes had a light with visible geometry, so this one didn't change any pixels, but it would have as soon as one did. The integrator is now a loop with a path throughput and a flag: after a diffuse vertex that sampled the lights, emission from light-sampled emitters is skipped; after a camera ray, mirror or glass it is counted. Multiple importance sampling would weight both strategies instead of switching one off, but it needs a pdf for every material; the flag is the textbook first step and stays unbiased.

        const materials::Material& material = *rec.material;
        if (count_light_emission || !material.emission_sampled_by_lights())
            radiance += throughput * material.emit();
        ...
        if (material.diffuse_brdf(rec, albedo, shading_normal)) {
            radiance += throughput * direct_lighting(rec, albedo, shading_normal, scene);
            count_light_emission = false;
        } else {
            count_light_emission = true;
        }

        throughput *= attenuation;
        ray = scattered;

        if (rr_start_ > 0 && depth + 1 >= rr_start_) {
            float p = glm::min(glm::max(throughput.r, glm::max(throughput.g, throughput.b)), 0.95f);
            if (core::random_float() >= p)
                break;
            throughput /= p;
        }

Russian roulette

Added in ab74a56

After 3 bounces the path continues with probability p = min(max(throughput.r, .g, .b), 0.95), and if it continues the throughput is divided by p. If X is the rest of the path's contribution, the estimator is X / p with probability p and 0 otherwise, so its expected value is still X. Dim paths stop early, and the cap keeps bright paths from running forever; the depth limit of 50 stays as a hard stop.

The BVH

Added in c1ddbc5

The old BVHNode::hit() returned false, no bounding boxes were computed and nothing used the class, so every ray tested every object. Now every bounded primitive returns an AABB (infinite planes stay in a short list that's tested linearly), and the tree is built top-down. For each node, the primitive centroids go into 16 bins per axis, and each of the 15 boundaries is scored with the surface area heuristic:

cost = 1 + (A_L · N_L + A_R · N_R) / A_node

A_L / A_node is the chance that a ray hitting the node also hits the left box, so the cost is the expected work for one ray. The cheapest split wins, and a node becomes a leaf (at most 4 primitives) when splitting isn't cheaper than testing everything in it.

            for (int i = 0; i < kBins - 1; ++i) {
                if (left_count[i] == 0 || right_count[i] == 0) continue;
                float cost = 1.0f + (left_area[i] * left_count[i] + right_area[i] * right_count[i]) /
                                        node_area;
                if (cost < best_cost) {
                    best_cost = cost;
                    best_axis = a;
                    best_split = i;
                }
            }

The nodes live in one flat std::vector, and a node's two children are adjacent so it stores one index. Traversal uses a fixed stack, clips each slab test against the closest hit so far, and visits the child nearer along the ray first, so far boxes are often culled.

    while (sp > 0) {
        const Node& node = nodes_[stack[--sp]];
        if (!slab_hit(node.bmin, node.bmax, origin, inv_dir, t_min, closest)) continue;

        if (node.count > 0) {
            for (uint32_t i = 0; i < node.count; ++i) {
                if (prims_[node.first + i]->hit(ray, t_min, closest, temp)) {
                    hit_anything = true;
                    closest = temp.t;
                    rec = temp;
                }
            }
        } else {
            // Push the far child first so the near child is popped next; a
            // near hit then shrinks 'closest' and lets the far box be culled.
            uint32_t near_child = node.first, far_child = node.first + 1;
            if (dir_neg[node.axis]) std::swap(near_child, far_child);
            stack[sp++] = far_child;
            stack[sp++] = near_child;
        }
    }

To check it, tests/bvh_test.cpp builds random scenes (triangles of mixed sizes, spheres, axis-aligned triangles, a plane) and fires 200,000 rays, including rays that start inside the scene and shadow rays with a finite t_max. BVH and brute force must agree on hit or miss, distance and which primitive was hit, for both SAH and midpoint builds. They do: 0 mismatches.

OpenMP and per-thread random numbers

Added in 42267be

CMake found OpenMP and added the flags, but no loop used it. The row loops now run with schedule(dynamic, 1), because rows through glass and metal take longer than sky rows. The single shared mt19937 would have been a data race, so the engine is a thread_local PCG32 reseeded for every pixel from a hash of (x, y, pass, seed). Every pixel gets the same random numbers whichever thread renders it, so the image is the same for any thread count.

#if defined(_OPENMP)
#pragma omp parallel for schedule(dynamic, 1)
#endif
    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            core::seed_rng_for_pixel(x, y, 0, seed_);

All renders and timings on this page come from the single-threaded WebAssembly build, so I haven't measured the OpenMP speed-up on this machine yet. The browser version gets its parallelism differently: each Web Worker runs its own instance on interleaved rows, and because pixels are seeded per pixel and pass, the result is bit-identical to the command-line renderer.

Ray-triangle intersection

Triangles use the Möller-Trumbore test. It solves for the hit distance t and the barycentric coordinates u and v directly with two cross products and a few dot products, without first intersecting the triangle's plane. Each early return rejects a case: the ray parallel to the triangle, u or v outside the triangle, or t outside the valid range. Smooth shading reuses the same u and v to interpolate vertex normals.

bool Triangle::hit(const core::Ray& ray, float t_min, float t_max, HitRecord& rec) const {
    // Möller-Trumbore algorithm
    Vec3 edge1 = v1_ - v0_;
    Vec3 edge2 = v2_ - v0_;
    Vec3 h = glm::cross(ray.direction(), edge2);
    float a = glm::dot(edge1, h);

    // Ray is parallel to triangle
    if (a > -1e-8f && a < 1e-8f) {
        return false;
    }

    float f = 1.0f / a;
    Vec3 s = ray.origin() - v0_;
    float u = f * glm::dot(s, h);

    if (u < 0.0f || u > 1.0f) {
        return false;
    }

    Vec3 q = glm::cross(s, edge1);
    float v = f * glm::dot(ray.direction(), q);

    if (v < 0.0f || u + v > 1.0f) {
        return false;
    }

    float t = f * glm::dot(edge2, q);

    if (t < t_min || t > t_max) {
        return false;
    }

Smooth shading for meshes

Added in 2e10902

Every triangle used to shade with its face normal, so the 17,280-triangle knot looked faceted: in a metal, neighbouring faces reflect visibly different directions and the facets read as stripes. A triangle hit now carries two normals. The geometric normal, the normalised cross product of two edges, is the real surface, so it still decides which side of the glass a ray is on and whether a direction is above the surface. The shading normal, which the materials and the BRDF use, is interpolated from the three vertex normals with the same barycentric weights Möller-Trumbore already returns: n = normalize((1 − u − v)·n0 + u·n1 + v·n2).

Vertex normals come from the OBJ's vn lines when it has them. Otherwise each vertex gets the sum of the unnormalised cross products of the faces around it; a cross product's length is twice the face's area, so big faces count for more than slivers. Metal absorbs a reflected ray that points below either normal, and direct light skips lights below the geometric normal, so a tilted shading normal near a silhouette can't send light through the mesh. Smooth shading is off unless the scene sets "smooth": true, and tests/smooth_normals_test.cpp checks the area weighting, the barycentric interpolation on a tessellated sphere and the inside/outside orientation.

Still open

• No multiple importance sampling, so caustics and small emitters that aren't registered lights are noisy at low sample counts.
• Mesh::hit is a stub; OBJ files load as individual triangles instead.
• Tangent frames aren't consistent between primitives: normal maps on planes have their green channel flipped.
• ImageTexture clamps UVs instead of wrapping, and its raw-data constructor allocates with new[] while the destructor calls stbi_image_free.
• The GLFW viewer ignores the scene file's camera, and it wasn't rebuilt on the machine used for the fixes.
• No tone mapping beyond gamma 2.2 and clamping. The shipped scenes were tuned with the π-too-bright light and now render darker.

How I built this

I wrote the original tracer while working through the book series. I led the 2026 review and fixes, and used Claude Code as a coding assistant for them; the renders and tests on this page are how I checked each one.

Back to all projects