Ray Tracing & Global Illumination

A ray is a half-line: a point of origin o and a unit direction d, parameterized as r(t) = o + td for t ≥ 0. Every pixel in the final image corresponds to at least one ray cast from the camera through that pixel's position on the virtual film plane.

Camera model. The standard pinhole camera places the eye at a point eye, looks toward a lookat target, and defines an up vector. From these you compute an orthonormal basis (u, v, w) for camera space:

• w = normalize(eye − lookat) (points away from the scene, backward along the view)
• u = normalize(up × w) (right)
• v = w × u (up, recomputed for orthogonality)

The image plane sits at distance d (the focal length) in front of the eye. For a film of width W and height H pixels, pixel (i, j) maps to the 3D point:

film_point = eye − d*w + (2i/W − 1)*half_w*u + (2j/H − 1)*half_h*v

where half_w and half_h are the half-extents of the film in world units (related to the field of view). The ray for pixel (i, j) has origin eye and direction normalize(film_point − eye).

Why this parameterization? The scalar t measures distance along the ray. Finding an intersection reduces to solving for t. The closest intersection is the one with the smallest positive t.

import numpy as np

def make_camera_rays(eye, lookat, up, fov_y, W, H):
    w = (eye - lookat); w /= np.linalg.norm(w)
    u = np.cross(up, w);  u /= np.linalg.norm(u)
    v = np.cross(w, u)
    half_h = np.tan(np.radians(fov_y / 2))
    half_w = half_h * W / H
    rays = []
    for j in range(H):
        for i in range(W):
            # normalized device coords in [-1,1]
            s = (2*(i+0.5)/W - 1) * half_w
            t = (1 - 2*(j+0.5)/H) * half_h
            direction = -w + s*u + t*v
            direction /= np.linalg.norm(direction)
            rays.append((eye, direction))
    return rays

The core computational primitive of ray tracing is the ray-object intersection test: given a ray r(t) = o + td, find the smallest t > t_min such that r(t) lies on the surface of some object.

Sphere intersection. A sphere of center c and radius r satisfies |p − c|² = r². Substituting p = o + td:

|(o + td) − c|² = r²

Let oc = o − c. Expanding:

t²(d·d) + 2t(d·oc) + (oc·oc − r²) = 0

This is a quadratic in t: at² + bt + c = 0. The discriminant Δ = b² − 4ac tells you: Δ < 0 means no intersection (ray misses), Δ = 0 means tangent, Δ > 0 means two intersections (entry and exit). The smaller positive root is the nearest hit.

Triangle intersection (Möller–Trumbore algorithm). Almost all real scene geometry is triangle meshes. A point on a triangle with vertices V0, V1, V2 is expressed in barycentric coordinates: p = (1−u−v)V0 + uV1 + vV2, with u ≥ 0, v ≥ 0, u+v ≤ 1. Setting o + td = p gives a 3×3 linear system solved by Cramer's rule in ~20 FLOPs.

Axis-aligned bounding box (AABB) intersection. For spatial acceleration, you need fast slab tests. An AABB is the intersection of three pairs of parallel planes (slabs). For each axis k, compute t_min_k = (box_min_k − o_k)/d_k and t_max_k = (box_max_k − o_k)/d_k. The ray intersects the box iff max(t_min_x, t_min_y, t_min_z) ≤ min(t_max_x, t_max_y, t_max_z) and this interval overlaps [t_near, t_far].

Surface normals. At the intersection point, you need the outward unit normal n for lighting. For a sphere: n = (hit − c) / r. For a triangle: n = normalize((V1−V0) × (V2−V0)).

def intersect_sphere(ray_o, ray_d, center, radius, t_min=1e-4, t_max=1e9):
    oc = ray_o - center
    a = np.dot(ray_d, ray_d)
    b = 2.0 * np.dot(oc, ray_d)
    c = np.dot(oc, oc) - radius**2
    disc = b*b - 4*a*c
    if disc < 0:
        return None  # miss
    sqrt_disc = np.sqrt(disc)
    t = (-b - sqrt_disc) / (2*a)  # near root
    if t < t_min or t > t_max:
        t = (-b + sqrt_disc) / (2*a)  # far root
        if t < t_min or t > t_max:
            return None
    hit_point = ray_o + t * ray_d
    normal = (hit_point - center) / radius
    return t, hit_point, normal

Once you find the nearest intersection, you need to determine the color of the surface at that point. Shading models specify how light interacts with a material.

Phong shading. The classical local illumination model decomposes reflectance into three terms:

1. Ambient: a constant k_a * I_a — a hack to approximate indirect light.
2. Diffuse (Lambertian): k_d * I * max(0, n·l) — light scattered uniformly in all directions; intensity depends only on the angle between the surface normal n and the light direction l.
3. Specular: k_s * I * max(0, r·v)^α — a glossy highlight. r = 2(n·l)n − l is the perfect mirror direction; v is the view direction; α controls shininess.

Physically based shading. Modern renderers use the BRDF (Bidirectional Reflectance Distribution Function), a function f_r(ω_i, ω_o) that gives the ratio of outgoing radiance to incoming irradiance for any pair of directions. A valid BRDF must be non-negative, satisfy Helmholtz reciprocity (f_r(ω_i, ω_o) = f_r(ω_o, ω_i)), and conserve energy (the integral of f_r * cos(θ) over the hemisphere ≤ 1).

The Cook-Torrance microfacet BRDF models surfaces as collections of tiny mirror facets:

f_r = (D * F * G) / (4 * (n·ω_i) * (n·ω_o))

where D is the normal distribution function (how many facets point in the half-vector direction), F is the Fresnel term (fraction of light reflected vs. refracted), and G is the geometric attenuation (self-shadowing between facets).

The key insight: shading is ultimately about evaluating the rendering equation (covered later). Local models like Phong are cheap approximations. Global illumination methods evaluate it correctly.

def phong_shade(hit_pt, normal, view_dir, light_pos, light_color,
               ka, kd, ks, alpha):
    l = light_pos - hit_pt
    dist = np.linalg.norm(l); l /= dist
    # diffuse
    diff = kd * light_color * max(0.0, np.dot(normal, l))
    # specular
    r = 2*np.dot(normal, l)*normal - l
    spec = ks * light_color * max(0.0, np.dot(r, view_dir))**alpha
    ambient = ka * light_color
    return ambient + diff + spec

Turner Whitted's 1980 paper "An Improved Illumination Model for Shaded Display" introduced recursive ray tracing. The key insight: at an intersection point, you can spawn secondary rays to compute the effects of reflection and refraction, and recursively evaluate them the same way.

Algorithm:

function trace(ray, depth):
    if depth > MAX_DEPTH: return BLACK
    (t, hit, normal, material) = nearest_intersection(ray)
    if no hit: return background_color

    color = local_shading(hit, normal, material)

    if material.is_reflective:
        reflect_dir = reflect(ray.dir, normal)
        color += material.k_r * trace(Ray(hit, reflect_dir), depth+1)

    if material.is_refractive:
        refract_dir = refract(ray.dir, normal, material.ior)
        color += material.k_t * trace(Ray(hit, refract_dir), depth+1)

    return color

Reflection. The perfect mirror reflection direction is: r = d − 2(d·n)n, where d is the incoming direction and n is the surface normal. This is just the vector reflection formula.

Refraction (Snell's law). When light crosses a boundary between media of refractive indices n1 and n2, its direction bends. Snell's law: n1 sin(θ_i) = n2 sin(θ_t). Vectorized:

Let η = n1/n2. The refracted direction is: t = ηd + (η cos(θ_i) − cos(θ_t))n

where cos(θ_i) = −(d·n) (assuming n points away from the surface) and cos(θ_t) = sqrt(1 − η²(1 − cos²(θ_i))). If the expression under the sqrt is negative, total internal reflection occurs — no refraction, only reflection.

Shadow rays. For each intersection, cast a shadow ray toward each light source. If anything intersects the shadow ray before it reaches the light, the point is in shadow and diffuse/specular terms are skipped. This is the cleanest, most elegant way to compute hard shadows — it requires no shadow maps, no depth buffer tricks.

Depth limit. The recursion must be bounded (MAX_DEPTH ~ 5-10 is typical). Rays that contribute negligibly to the final image can be terminated early using the Russian roulette technique.

def reflect(d, n):
    return d - 2*np.dot(d, n)*n

def refract(d, n, ni_over_nt):
    uv = d / np.linalg.norm(d)
    dt = np.dot(uv, n)
    disc = 1.0 - ni_over_nt**2 * (1 - dt**2)
    if disc > 0:
        return ni_over_nt*(uv - n*dt) - n*np.sqrt(disc)
    return None  # total internal reflection

def trace(ray_o, ray_d, scene, depth=0, MAX_DEPTH=5):
    if depth >= MAX_DEPTH:
        return np.zeros(3)
    hit = nearest_intersection(ray_o, ray_d, scene)
    if hit is None:
        return scene.background
    t, pt, n, mat = hit
    color = local_shade(pt, n, mat, scene)
    if mat.k_r > 0:
        r_dir = reflect(ray_d, n)
        color += mat.k_r * trace(pt + 1e-4*r_dir, r_dir, scene, depth+1)
    if mat.k_t > 0:
        t_dir = refract(ray_d, n, mat.ior)
        if t_dir is not None:
            color += mat.k_t * trace(pt + 1e-4*t_dir, t_dir, scene, depth+1)
    return color

Whitted ray tracing handles mirror reflection and refraction beautifully but fails for soft shadows, diffuse inter-reflections, color bleeding, and glossy materials. The reason: it doesn't account for light arriving from arbitrary directions — only from point lights and perfect mirrors.

In 1986, Jim Kajiya published "The Rendering Equation," which places all of rendering on a rigorous physical foundation. The rendering equation states that the outgoing radiance L_o at point x in direction ω_o is:

L_o(x, ω_o) = L_e(x, ω_o) + ∫_{Ω} f_r(x, ω_i, ω_o) L_i(x, ω_i) cos(θ_i) dω_i

where:

• L_e is the emitted radiance (the surface is a light source)
• f_r is the BRDF of the material at x
• L_i(x, ω_i) is the incoming radiance from direction ω_i — which itself satisfies the rendering equation recursively
• cos(θ_i) = n·ω_i is the geometric foreshortening factor
• The integral is over the hemisphere Ω of incoming directions

Why this is hard. L_i(x, ω_i) is itself the outgoing radiance from the surface visible from x in direction −ω_i, so the equation is recursive/self-referential — it's a Fredholm integral equation of the second kind. There is no general closed-form solution. All practical rendering algorithms are methods for approximating this integral.

Physical basis. Radiance L (Watts per square meter per steradian) is conserved along rays in a vacuum — so L_i(x, ω_i) = L_o(x', −ω_i) where x' is the first surface hit when marching from x in direction ω_i. This is why ray tracing naturally evaluates the rendering equation: firing a ray in direction ω_i from x computes L_i(x, ω_i).

Special cases:

• Point light, mirror BRDF → Whitted ray tracing
• Lambertian BRDF + area lights → radiosity (finite element solution on diffuse surfaces)
• Full BRDF, stochastic sampling → path tracing

# Rendering equation in pseudocode:
# L_o(x, wo) = L_e(x, wo)
#            + integral over hemisphere of:
#                 f_r(x, wi, wo) * L_i(x, wi) * cos(theta_i) dwi
#
# In a path tracer, we ESTIMATE this integral by sampling:
# wi ~ p(wi)  (some probability density over the hemisphere)
# estimate = f_r(x, wi, wo) * L_i(x, wi) * cos(theta_i) / p(wi)
#
# L_i(x, wi) is computed recursively by tracing a ray from x in direction wi

Path tracing is the Monte Carlo method applied to the rendering equation. Instead of computing the hemisphere integral analytically, we sample it stochastically.

Monte Carlo integration. To estimate ∫ f(x) dx over a domain, sample N points x_i from a distribution with PDF p(x), and compute:

(1/N) Σ f(x_i) / p(x_i)

By the law of large numbers, this converges to the true integral as N → ∞. The variance scales as 1/N, so standard deviation scales as 1/√N — doubling quality requires 4× more samples.

The path tracer. Each pixel gets many samples. For each sample, trace a path through the scene:

1. Cast a primary ray from the camera.
2. At each intersection, sample a new direction ω_i from the BRDF (or uniformly from the hemisphere).
3. Trace the new ray. Multiply the accumulated path weight by f_r * cos(θ) / p(ω_i).
4. Continue until the ray escapes the scene or is terminated by Russian roulette.
5. If the path hits a light source, accumulate its emission weighted by the path weight.

The result is an unbiased estimator of the rendering equation — it converges to the physically correct answer. This is what makes path tracing the ground truth algorithm.

Importance sampling. Naive uniform sampling over the hemisphere is inefficient. BRDF importance sampling instead samples directions proportional to the BRDF (or BRDF * cos(θ)), concentrating samples where they contribute most. For a Lambertian surface, cosine-weighted hemisphere sampling has PDF p(ω) = cos(θ)/π, which exactly cancels the cos(θ) factor, giving:

estimate = f_r * L_i (no denominator, minimal variance)

Next event estimation (direct light sampling). At each bounce, additionally sample a point on each light source and compute its contribution directly (with a visibility test). This dramatically reduces variance for scenes with small lights.

Russian roulette. To terminate paths without bias, at each bounce keep the path with probability q (typically proportional to the surface albedo) and boost the weight by 1/q if kept. Paths that are terminated contribute zero. In expectation, the estimator is unchanged, but computation is saved on low-contribution paths.

import random, math

def cosine_sample_hemisphere(n):
    """Sample direction proportional to cos(theta) around normal n."""
    r1, r2 = random.random(), random.random()
    phi = 2 * math.pi * r1
    cos_t = math.sqrt(1 - r2)
    sin_t = math.sqrt(r2)  # sin(theta) = sqrt(1 - cos^2)
    # local coords
    wi_local = (sin_t*math.cos(phi), sin_t*math.sin(phi), cos_t)
    return local_to_world(wi_local, n)  # rotate to world frame

def path_trace(ray_o, ray_d, scene, depth=0, MAX_DEPTH=10):
    hit = nearest_intersection(ray_o, ray_d, scene)
    if hit is None:
        return scene.env_light(ray_d)
    t, pt, n, mat = hit
    # Emission
    L = mat.emission
    # Russian roulette
    q = max(mat.albedo)  # survival prob
    if depth > 3 and random.random() > q:
        return L
    boost = 1.0 / q if depth > 3 else 1.0
    # Sample new direction
    wi = cosine_sample_hemisphere(n)
    # PDF for cosine sampling: cos(theta)/pi
    # f_r for Lambertian: albedo/pi
    # weight = f_r * cos(theta) / pdf = albedo (cos cancels!)
    weight = mat.albedo
    L += boost * weight * path_trace(pt + 1e-4*wi, wi, scene, depth+1)
    return L

A naive ray tracer tests every ray against every object — O(n) per ray, O(n * W * H * samples) per image. With millions of triangles and millions of rays, this is completely intractable. Acceleration structures reduce intersection to O(log n) per ray.

Bounding Volume Hierarchy (BVH). A BVH is a binary tree where each node stores an AABB that tightly bounds all primitives in its subtree. Leaves contain a small number of actual primitives.

Traversal: at each node, test the ray against the node's AABB. If the ray misses the AABB, skip the entire subtree (O(1) cull). If it hits, recurse into both children. Only leaf nodes require actual primitive intersections.

Construction: sort primitives by their centroid along the longest axis, split at the median (or use the Surface Area Heuristic — SAH — which minimizes expected traversal cost: split at the position minimizing SA(left)*N_left + SA(right)*N_right, where SA is surface area and N is primitive count).

kd-tree. Partitions 3D space with axis-aligned splitting planes. Efficient for static scenes; BVH is generally preferred for dynamic scenes and GPU implementation.

GPU ray tracing (RTX). Modern NVIDIA RTX GPUs have dedicated BVH traversal hardware (RT cores) that performs ray-AABB and ray-triangle tests in fixed-function silicon. The BVH is built by the driver, traversal and leaf intersection happen in hardware, and the user supplies shaders (closest hit, any hit, miss) for the hit and miss events. This achieves real-time performance for complex scenes.

Practical impact. With a BVH, a scene with 10M triangles requires ~23 AABB tests per ray on average — down from 10M. This is the difference between 1 hour per frame and 1/10 of a second.

class BVHNode:
    def __init__(self, primitives):
        self.bbox = AABB.union(p.bbox for p in primitives)
        if len(primitives) <= 4:
            self.left = self.right = None
            self.prims = primitives
        else:
            axis = self.bbox.longest_axis()
            primitives.sort(key=lambda p: p.centroid[axis])
            mid = len(primitives) // 2
            self.left  = BVHNode(primitives[:mid])
            self.right = BVHNode(primitives[mid:])
            self.prims = []

    def intersect(self, ray):
        if not self.bbox.hit(ray):
            return None  # cull entire subtree
        if self.prims:  # leaf
            return min_hit(p.intersect(ray) for p in self.prims)
        hit_l = self.left.intersect(ray)
        hit_r = self.right.intersect(ray)
        return closer_hit(hit_l, hit_r)

A pixel is a small square, not a point. If you shoot exactly one ray through the center of each pixel, you get aliasing — jagged edges on diagonal lines and moire patterns in detailed geometry. This is the classic sampling theorem problem in disguise.

The pixel as an integral. The true pixel value is the integral of radiance over the pixel's area (and lens aperture, and time if motion blur is needed). We estimate this by sampling multiple sub-pixel positions and averaging:

I(pixel) = (1/n) Σ L(ray_from_sample_k)

Supersampling. The simplest approach: shoot n rays per pixel, jittered (not purely random, not purely grid) for better variance properties. Jittered sampling divides the pixel into an n×n grid and places one sample randomly within each cell — this is stratified sampling, which reduces variance compared to pure random (Monte Carlo) sampling.

Fourier perspective. Aliasing occurs when the signal (scene) has frequencies higher than the Nyquist limit for the sample rate. Regular supersampling converts aliasing into noise (high-frequency, less objectionable). Importance-weighted sampling (e.g., more samples near edges) is the adaptive equivalent.

Depth of field. A pinhole camera has infinite depth of field. A physical lens focuses light from a specific focal plane; objects outside it are blurry. In a ray tracer, simulate this by sampling a random point on the lens aperture for each ray, adjusting the ray direction so that the focal plane is unchanged. Multiple rays per pixel average into the correct bokeh.

Motion blur. Sample the time dimension: give each ray a random time in [t0, t1] and move scene geometry accordingly. Objects in motion will be rendered blurry along their trajectory. This comes for free in a path tracer once you have stochastic sampling.

def render_pixel_aa(eye, i, j, W, H, scene, n_samples=64):
    color = np.zeros(3)
    # n_samples stratified samples in the pixel
    n = int(np.sqrt(n_samples))
    for si in range(n):
        for sj in range(n):
            # jitter within stratified cell
            u_off = (si + random.random()) / n
            v_off = (sj + random.random()) / n
            # pixel center + sub-pixel offset
            s = (i + u_off) / W
            t = (j + v_off) / H
            ray_d = pixel_to_ray_dir(s, t, eye, scene.camera)
            color += path_trace(eye, ray_d, scene)
    return color / n_samples

When light hits a surface, how much is reflected vs. refracted (transmitted)? The answer is given by the Fresnel equations, derived from Maxwell's equations of electromagnetism.

Fresnel reflectance. For unpolarized light at angle θ_i (from the normal) hitting a surface with refractive indices n1 (incoming) and n2 (outgoing):

R_s = ((n1 cos θ_i − n2 cos θ_t) / (n1 cos θ_i + n2 cos θ_t))²

R_p = ((n1 cos θ_t − n2 cos θ_i) / (n1 cos θ_t + n2 cos θ_i))²

R = (R_s + R_p) / 2

T = 1 − R (energy conservation)

Schlick's approximation. Computing exact Fresnel is expensive. Christophe Schlick (1994) proposed a fast approximation:

R(θ) ≈ R_0 + (1 − R_0)(1 − cos θ)^5

where R_0 = ((n1−n2)/(n1+n2))². This is the standard in real-time and path-tracing implementations.

Physical meaning. At normal incidence (θ = 0), only a small fraction reflects (e.g., ~4% for glass-air). At grazing incidence (θ → 90°), all light reflects — that's why wet roads are mirror-like when viewed at a glance. Metals have complex refractive indices with a large imaginary part, giving high reflectance at all angles and characteristic tinted reflections (gold, copper).

Implementation. In a path tracer, at each glass surface, use the Fresnel term to stochastically choose: with probability R, follow the reflected ray; with probability T=1−R, follow the refracted ray. This IS unbiased importance sampling of the glass BRDF.

def schlick(cosine, n1, n2):
    r0 = ((n1 - n2) / (n1 + n2))**2
    return r0 + (1 - r0) * (1 - cosine)**5

def scatter_glass(ray_d, normal, ior):
    """Return (new_ray_dir, attenuation) for a glass surface."""
    # Determine if inside or outside the surface
    if np.dot(ray_d, normal) < 0:
        n1, n2, n = 1.0, ior, normal
    else:
        n1, n2, n = ior, 1.0, -normal
    cosine = -np.dot(ray_d, n)
    R = schlick(cosine, n1, n2)
    refracted = refract(ray_d, n, n1/n2)
    if refracted is None or random.random() < R:
        return reflect(ray_d, n), np.array([1.0,1.0,1.0])  # reflect
    else:
        return refracted, np.array([1.0,1.0,1.0])  # transmit

To reason correctly about the rendering equation, you need precise physical quantities. Radiometry is the science of measuring electromagnetic radiation.

Radiant energy Q [J]: total energy in photons.

Radiant flux (power) Φ [W = J/s]: energy per unit time.

Irradiance E [W/m²]: flux arriving per unit area. The cos(θ) factor in the rendering equation arises because a beam of light striking a surface at angle θ is spread over a larger area by factor 1/cos(θ).

Radiance L [W/(m²·sr)]: flux per unit projected area per unit solid angle. This is the fundamental quantity in rendering — it's what the camera measures, what lights emit, and what rays carry. Crucially, radiance is conserved along a ray in a vacuum (the intensity of a beam doesn't change with distance, unlike irradiance).

Solid angle [sr]: the 2D angle generalized to 3D. The solid angle subtended by a cone is its area on the unit sphere. The full sphere is 4π sr. A differential solid angle is dω = sin(θ) dθ dφ.

Connection to the rendering equation: The rendering equation is written in radiance. The integral ∫ L cos(θ) dω converts incoming radiance to irradiance — then the BRDF converts irradiance to outgoing radiance. This dimensional chain is exact.

BRDF units. The BRDF f_r has units 1/sr, ensuring that f_r * L_i * cos(θ) * dω has units of radiance. A Lambertian BRDF f_r = ρ/π (where ρ ∈ [0,1] is albedo) is the canonical example: integrating over the hemisphere gives outgoing radiance = ρ * E (albedo times irradiance — all diffuse reflection).

The central challenge in path tracing is variance: with few samples, images are noisy. Variance reduction techniques reduce noise without bias — they still converge to the correct answer, just faster.

Multiple Importance Sampling (MIS). Combine samples from multiple strategies (e.g., BRDF sampling + light sampling) using the balance heuristic or power heuristic weights. For each sample from strategy i with PDF p_i, the MIS weight is w_i = p_i / Σ_j p_j. This nearly eliminates the catastrophic variance that occurs when a single strategy is a poor match for the integrand (e.g., a highly specular surface near a large area light).

Bidirectional Path Tracing (BDPT). Trace paths from both the camera and the light, connecting midpoints. Paths that are hard to find by camera-only tracing (caustics, inside a glass ball) are easy from the light side. BDPT combines all possible path lengths using MIS.

Metropolis Light Transport (MLT). Uses Metropolis-Hastings MCMC to sample paths proportional to their contribution. Once a bright path (e.g., a caustic) is found, mutations explore nearby paths — efficiently concentrating samples on high-contribution regions. Unbiased but requires careful implementation.

Denoising. Modern production renderers (e.g., OptiX AI Denoiser, Intel OIDN) use convolutional neural networks trained on clean/noisy image pairs to remove noise from as few as 4-8 samples per pixel. G-buffers (normals, depth, albedo) provide auxiliary inputs that help the network distinguish noise from signal. This is not unbiased but is immensely practical.

Low-discrepancy sequences. Replace pseudo-random numbers with quasi-random sequences (Halton, Sobol, scrambled Halton) that fill the sample space more uniformly. Quasi-Monte Carlo can reduce variance faster than 1/N for smooth integrands.

# MIS balance heuristic weight for n strategies
def mis_weight(pdf_a, pdf_b):
    """Weight for a sample from strategy A when combining A and B."""
    return pdf_a / (pdf_a + pdf_b)

# Usage in a path tracer with BRDF + light sampling:
# 1. Sample light: wi_light, pdf_light = sample_light(pt, scene)
# 2. Sample BRDF: wi_brdf,  pdf_brdf  = sample_brdf(ray_d, n, mat)
# 3. Evaluate both directions with MIS weights:
# L_light = brdf(wi_light)*Li(wi_light)*cos/pdf_light * mis_weight(pdf_light, brdf_pdf(wi_light))
# L_brdf  = brdf(wi_brdf) *Li(wi_brdf) *cos/pdf_brdf  * mis_weight(pdf_brdf,  light_pdf(wi_brdf))
# L_total = L_light + L_brdf

For decades, ray tracing was offline-only (hours per frame for film VFX). NVIDIA's Turing architecture (2018) introduced RT Cores — dedicated hardware for BVH traversal and ray-triangle intersection — enabling real-time ray tracing in games.

DXR / Vulkan Ray Tracing API. The shader pipeline adds new stages:

• Ray generation shader: runs once per pixel, calls TraceRay()
• Intersection shader: custom primitive intersection (optional; triangle is built-in)
• Any-hit shader: called for each candidate intersection (e.g., alpha test)
• Closest-hit shader: called once for the nearest hit (shading)
• Miss shader: called when no intersection is found (skybox, etc.)

The BVH (Acceleration Structure) is built by the driver. Traversal and leaf intersection happen in RT Core hardware. Shader invocations run on the standard CUDA/SM cores.

Hybrid rendering. Full path tracing at 60fps is still too expensive for complex scenes. Production games use hybrid approaches:

• Rasterize primary visibility (G-buffer)
• 1 ray per pixel for shadows or ambient occlusion
• 1-2 rays per pixel for reflections or diffuse GI
• AI denoise the result

This gives the visual quality of ray tracing at a fraction of the cost by concentrating rays only where rasterization falls short (reflections, shadows, GI).

Performance numbers. An RTX 4090 achieves ~100 gigarays/second for simple traversal. A fully path-traced 1080p frame at 1 spp takes ~2ms, so 4K/120fps with 64 spp requires enormous denoising and clever sampling.