Introduction to 3D Rendering

We will use this $M_{camera}$ matrix to transform these rays into the world space. For generating the rays we will use the pixel index.

Ray Generation Code

DrJit ray generation code

The following Python code shows how we can implement a camera class in order to create rays:

The code below uses drjit, the same array library that powers the Mitsuba renderer.

import math
import drjit as dr

# can switch between 'cuda' or 'llvm'
from drjit.cuda import Float, Int, Array3f, Matrix3f

class Ray:
    def __init__(self, origin, direction):
        self.o = origin
        self.d = direction

class Camera:
    def __init__(self, pos, target, fov, aspect=None, H=800, W=800):
        self.pos = Array3f(*pos)
        self.target = Array3f(*target)
        self.world_up = Array3f(0, 1, 0)
        self.fov = fov
        self.aspect = aspect if aspect is not None else W / H
        self.H = H
        self.W = W
        self.f = 1.0

    def lookAt(self):
        forward = dr.normalize(self.target - self.pos)
        right = dr.normalize(dr.cross(forward, self.world_up))
        up = dr.normalize(dr.cross(right, forward))

        M_cam = Matrix3f(right, up, -forward)
        return M_cam


    def rays(self):
        
        # Generate 1D pixel indices
        num_pixels = self.H * self.W
        idx = dr.arange(Int, num_pixels)

        # Generate the indices of pixels as (x, y)
        x_ind = idx % self.W
        y_ind = idx // self.W

        # Scale to [0, 1] range with +0.5 to reach center 
        # or change 0.5 to randomness in [0, 1] to jitter within pixel
        x_center = (Float(x_ind) + 0.5) / self.W
        y_center = (Float(y_ind) + 0.5) / self.H

        xs = x_center * 2.0 - 1.0
        ys = y_center * -2.0 + 1.0

        # Get the height and width in physical units
        # tan(beta / 2) = height /2 * f
        height = 2.0 * self.f * math.tan(math.radians(self.fov) / 2.0)
        width = height * self.aspect

        # Scale to [-width/2, width/2] x [-height/2, height/2]
        px = xs * (width / 2.0)
        py = ys * (height / 2.0)

        # Create the Camera-to-World transformation matrix.
        M_cam = self.lookAt()

        # Pack the camera-space ray directions into an Array3f
        dirs_cam = Array3f(px, py, -self.f)

        # 5. Matrix multiplication (Matrix @ Vector)
        rd = M_cam @ dirs_cam
        rd = dr.normalize(rd) # normalize the ray direction

        ro = self.pos 

        return Ray(ro, rd) # Convert to Ray class and return 

For the glsl shader, the pixel coordinate itself comes from gl_FragCoord. Dividing by u_resolution.y keeps the field of view stable while correcting the aspect ratio:

(We can also see the color of pixels in Fig. is same as the image formed)

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

// Builds the Camera-to-World transformation matrix
mat3 lookAt(vec3 pos, vec3 target, vec3 world_up) {
    vec3 forward = normalize(target - pos);
    vec3 right   = normalize(cross(forward, world_up));
    vec3 up      = cross(right, forward);
    
    // GLSL matrices are column-major. 
    // We negate 'forward' because the camera looks down -Z in camera space.
    return mat3(right, up, -forward);
}

void main() {
    // Pixel coordinates mapped to[-1, 1] range (roughly)
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    
    // Camera Setup
    // Place the camera at the origin, looking straight down the -Z axis
    vec3 ro = vec3(0.0, 0.0, 0.0);         // Ray Origin (Eye)
    vec3 target = vec3(0.0, 0.0, -1.0);    // Look at -Z
    vec3 world_up = vec3(0.0, 1.0, 0.0);
    
    mat3 M_cam = lookAt(ro, target, world_up);
    
    // Focal length (controls Field of View)
    float f = 1.0; 
    
    // Camera-space ray direction
    vec3 rd_cam = normalize(vec3(uv, -f)); 
    
    // Transform ray direction to World Space
    vec3 rd = M_cam * rd_cam;

    // Visualize direction: [-1,1] -> [0,1].
    fragColor = vec4(rd, 1.0);
}

The Simplified Camera Shortcut

For the rest of the GLSL examples in this post, we will use a common “shortcut”. Because the camera stays at the origin vec3(0,0,0) (or is shifted only along the Z-axis without rotation) and always looks straight down the $-Z$ axis, the M_cam matrix is simply the identity matrix.

This means we can skip the matrix multiplication entirely and construct our ray direction directly from the pixel coordinates:

vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
vec3 ro = vec3(0.0, 0.0, 3.0);       // Shifted back slightly, but not rotated
vec3 rd = normalize(vec3(uv, -1.0)); // Direct to world-space!

Anti-Aliasing

Sampling only the pixel center is simple, but it aliases hard edges. The jaggies (also called staircase artifacts) are clearly visible in the polygon render. Ideally, we want the average colors of objects within the pixel area. A common fix is to take several samples at slightly perturbed positions inside the pixel and average the results (A Monte Carlo estimate of the weighted area integration):

In a production renderer those offsets are usually stratified, blue-noise, or accumulated over many frames. For this introductory shader, one center ray is enough to keep the moving parts visible.

Geometry Representations and Intersections

Once we have rays shooting into our scene, the immediate next question is: how is the scene represented? Based on how we define our geometry, the intersection logic—whether we compute exact intersections or iteratively sample along the ray—will fundamentally change.

In computer graphics, geometry is typically represented in three main ways:

• Analytic primitives such as spheres, planes, and cylinders are defined by exact mathematical equations. Their intersections usually have closed-form algebraic solutions. While they are lightning-fast to intersect, they are not very flexible for modeling complex real-world objects.
• Triangles (Explicit Geometry) are the fundamental building blocks of 3D graphics. They are explicit surface primitives that can approximate any shape. They are widely used in rasterization because they project cleanly to screen space, but ray tracers also intersect them directly using highly optimized ray-triangle tests.
• Signed Distance Functions (Implicit Geometry) are continuous mathematical fields. Instead of storing a mesh of explicit vertices, an SDF returns the shortest distance from any point in space to the nearest surface. It returns a positive value outside the object, zero exactly on the surface, and a negative value inside. Intersecting with an SDF is done through an iterative process called ray marching, as there is no explicit surface to algebraically intersect with.

Let’s dive into the mathematics and code for each of these representations.

1. Ray-Sphere Intersection

A sphere is one of the simplest analytic primitives to intersect. It is defined by its center point ${\color{#109618}\mathbf{c}}$ and its radius ${\color{#990099}R}$. Any point $\mathbf{p}$ on the surface of the sphere satisfies the equation:

We want to find if our ray ever reaches a point $\mathbf{p}$ that satisfies this equation. Recall our ray equation, defined by an origin ${\color{#3273F6}\mathbf{o}}$ and a normalized direction ${\color{#FF4B4B}\mathbf{d}}$:

To find the intersection, we simply substitute the ray equation into the sphere equation:

To make the math cleaner, let $\mathbf{q} = {\color{#3273F6}\mathbf{o}} - {\color{#109618}\mathbf{c}}$ (the vector from the sphere’s center to the ray’s origin). Expanding the squared magnitude (which is the dot product of the vector with itself) gives us a classic quadratic equation in terms of $t$:

This is in the standard quadratic form $At^2 + Bt + C = 0$, where:

• $A = {\color{#FF4B4B}\mathbf{d}} \cdot {\color{#FF4B4B}\mathbf{d}}$ (If our ray direction is normalized, $A = 1$)
• $B = 2(\mathbf{q} \cdot {\color{#FF4B4B}\mathbf{d}})$
• $C = \mathbf{q} \cdot \mathbf{q} - {\color{#990099}R}^2$

Using the quadratic formula $t = \frac{-B \pm \sqrt{B^2 - 4AC}}{2A}$, we can find the roots.

• If the discriminant ($B^2 - 4AC$) is negative, the ray misses the sphere.
• If it is zero, the ray grazes the sphere (one hit).
• If it is positive, the ray goes through the sphere (two hits). The smallest positive root $t$ is our first visible hit.

Once we find the hit distance $t$, the hit point is $\mathbf{p} = {\color{#3273F6}\mathbf{o}} + t{\color{#FF4B4B}\mathbf{d}}$, and the surface normal is simply the normalized vector from the center to the hit point: $\mathbf{n} = (\mathbf{p} - {\color{#109618}\mathbf{c}}) / {\color{#990099}R}$.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

// Returns t if it hits, -1.0 otherwise. Also outputs the normal.
float intersectSphere(vec3 ro, vec3 rd, vec3 center, float radius, out vec3 normal) {
    vec3 q = ro - center;
    float b = dot(q, rd);
    float c = dot(q, q) - radius * radius;
    
    // Discriminant (A = 1 because rd is normalized)
    float h = b * b - c;
    if (h < 0.0) return -1.0; // Ray misses
    
    // Smallest positive root
    float t = -b - sqrt(h);
    if (t > 0.0) {
        vec3 p = ro + t * rd;
        normal = normalize(p - center);
        return t;
    }
    return -1.0;
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 0.0, 3.0);
    vec3 rd = normalize(vec3(uv, -1.0));

    vec3 center = vec3(0.0, 0.0, 0.0);
    float radius = 1.0;
    
    vec3 normal;
    float t = intersectSphere(ro, rd, center, radius, normal);
    
    vec3 col = vec3(0.05); // Background color
    if (t > 0.0) {
        // Map normal from [-1, 1] to [0, 1] for visualization
        col = normal * 0.5 + 0.5;
    }

    fragColor = vec4(col, 1.0);
}

2. Ray-Triangle Intersection (Möller-Trumbore)

Triangles are the backbone of 3D modeling. Let a triangle be defined by its three vertices ${\color{#FF9900}\mathbf{a}}$, ${\color{#0099C6}\mathbf{b}}$, and ${\color{#DD4477}\mathbf{c}}$. Any point $\mathbf{p}$ lying on the plane of this triangle can be written using barycentric coordinates $(u, v)$:

Let’s define the two edge vectors of the triangle originating from ${\color{#FF9900}\mathbf{a}}$ as $\mathbf{e}_1 = {\color{#0099C6}\mathbf{b}} - {\color{#FF9900}\mathbf{a}}$ and $\mathbf{e}_2 = {\color{#DD4477}\mathbf{c}} - {\color{#FF9900}\mathbf{a}}$. The ray hits the triangle when the ray equation equals the triangle plane equation:

Rearranging the terms to isolate our unknowns ($t, u, v$) on the left gives us a system of three linear equations:

Let $\mathbf{s} = {\color{#3273F6}\mathbf{o}} - {\color{#FF9900}\mathbf{a}}$ (the vector from vertex $\mathbf{a}$ to the ray origin). We can solve this $3\times3$ system elegantly using Cramer’s Rule, which expresses the solution in terms of determinants.

Geometrically, the determinant of a $3\times3$ matrix composed of three vectors $\det(\mathbf{u}, \mathbf{v}, \mathbf{w})$ is the scalar triple product, written as $\mathbf{u} \cdot (\mathbf{v} \times \mathbf{w})$. It represents the volume of the parallelepiped formed by the three vectors.

Using Cramer’s Rule, the solution is:

This translates to the famous Möller-Trumbore algorithm. To make it computationally efficient, we define cross products that we can reuse:

1. Let $\mathbf{p} = {\color{#FF4B4B}\mathbf{d}} \times \mathbf{e}_2$. The denominator becomes $\mathbf{e}_1 \cdot \mathbf{p}$. If this is $0$, the ray is parallel to the triangle.
2. The $u$ coordinate becomes $u = (\mathbf{s} \cdot \mathbf{p}) / (\mathbf{e}_1 \cdot \mathbf{p})$.
3. Let $\mathbf{q} = \mathbf{s} \times \mathbf{e}_1$. The $v$ coordinate becomes $v = ({\color{#FF4B4B}\mathbf{d}} \cdot \mathbf{q}) / (\mathbf{e}_1 \cdot \mathbf{p})$.
4. The hit distance becomes $t = (\mathbf{e}_2 \cdot \mathbf{q}) / (\mathbf{e}_1 \cdot \mathbf{p})$.

The intersection is strictly inside the triangle if $u \ge 0$, $v \ge 0$, $u + v \le 1$, and $t > 0$. Notice how this brilliantly evaluates the hit without ever explicitly computing the plane’s normal!

#version 300 es
precision highp float;

uniform vec2 u_resolution;
uniform float u_time;
out vec4 fragColor;

// Möller-Trumbore Ray-Triangle Intersection
bool intersectTriangle(vec3 ro, vec3 rd, vec3 a, vec3 b, vec3 c, out float t, out vec3 bary) {
    vec3 e1 = b - a;
    vec3 e2 = c - a;
    vec3 pvec = cross(rd, e2);
    float det = dot(e1, pvec);

    // If det is near zero, ray lies in the plane of the triangle
    if (abs(det) < 0.0001) return false;

    float invDet = 1.0 / det;
    vec3 svec = ro - a;
    
    // Calculate u and test bounds
    float u = dot(svec, pvec) * invDet;
    if (u < 0.0 || u > 1.0) return false;

    vec3 qvec = cross(svec, e1);
    
    // Calculate v and test bounds
    float v = dot(rd, qvec) * invDet;
    if (v < 0.0 || u + v > 1.0) return false;

    // Calculate t
    t = dot(e2, qvec) * invDet;
    
    bary = vec3(1.0 - u - v, u, v);
    return t > 0.0;
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 0.0, 3.0);
    vec3 rd = normalize(vec3(uv, -1.0));

    float wobble = 0.25 * sin(u_time * 2.0);
    vec3 a = vec3(-1.0, -0.7, 0.0);
    vec3 b = vec3( 1.0, -0.6, 0.0);
    vec3 c = vec3( wobble, 0.9, 0.0);

    float t;
    vec3 bary;
    vec3 col = vec3(0.05); // Background
    
    if (intersectTriangle(ro, rd, a, b, c, t, bary)) {
        col = bary; // Visualize barycentric coordinates directly as RGB
    }

    fragColor = vec4(col, 1.0);
}

3. Signed Distance Functions (SDFs)

In contrast to explicit geometry representations like triangle meshes, a Signed Distance Function (SDF) describes geometry implicitly through a continuous mathematical field. Rather than storing discrete vertices, an SDF evaluates the shortest distance from any given point in space to the nearest surface.

Formally, let $\Omega \subset \mathbb{R}^3$ be a closed set representing a solid object, with its boundary surface denoted by $\partial \Omega$. The signed distance function $f(\mathbf{p})$ at any point $\mathbf{p} \in \mathbb{R}^3$ is defined as the minimum metric distance from $\mathbf{p}$ to the boundary $\partial \Omega$, signed by whether $\mathbf{p}$ is in the interior or exterior of the object:

The sign convention dictates that $f(\mathbf{p})$ is positive outside the object, negative inside, and exactly zero on the surface $\partial \Omega$. Consequently, the surface itself is defined implicitly as the zero-level set of the function: $\partial \Omega = \{\mathbf{p} \in \mathbb{R}^3 \mid f(\mathbf{p}) = 0\}$.

Because $f(\mathbf{p})$ guarantees the exact distance to the nearest surface, it provides a strictly bound spatial radius. If $f(\mathbf{p}) = d$, a sphere of radius $d$ centered at $\mathbf{p}$ is guaranteed to be entirely empty. This property allows ray marching algorithms to safely advance along a ray by distance $d$ without prematurely overshooting the surface.

The Eikonal Equation and Lipschitz Continuity

For a function to be a valid Euclidean signed distance field, it must satisfy the Eikonal equation:

This nonlinear partial differential equation states that the magnitude of the gradient of the function must equal exactly 1 everywhere (except at singularities, such as the medial axis, where the gradient is undefined). Geometrically, this implies that the field’s rate of change perfectly matches physical space: moving one unit of distance directly away from the surface increases the SDF value by exactly one.

In practice, when constructing complex 3D scenes, SDFs are often subjected to spatial deformations—such as scaling, twisting, or blending—that violate the strict Eikonal property. For rendering algorithms like sphere tracing to remain mathematically safe, the function must at least maintain a Lipschitz bound of $1$:

If the gradient magnitude exceeds 1 ($\|\nabla f\| > 1$), the function overestimates the true distance. A sphere tracer would march too far and step through the surface, causing visual artifacts or missing geometry entirely. Conversely, if the gradient magnitude is substantially less than 1, the algorithm remains mathematically safe but becomes highly inefficient, forcing the ray to take excessively small, conservative steps. Maintaining the Eikonal property, or a strict Lipschitz bound, is a fundamental constraint when designing operations for SDFs.

Constructive Solid Geometry (CSG)

Complex scenes can be constructed from primitive SDFs using Boolean operations. If $a = f_1(\mathbf{p})$ and $b = f_2(\mathbf{p})$ represent the signed distances to two individual shapes at a given point, standard Constructive Solid Geometry (CSG) operations can be formulated as combinations of min and max functions:

float opUnion(float a, float b) {
    return min(a, b);
}

float opIntersection(float a, float b) {
    return max(a, b);
}

float opSubtraction(float a, float b) {
    return max(a, -b);
}

float opSmoothUnion(float a, float b, float k) {
    k *= 4.0;
    float h = max(k - abs(a - b), 0.0);
    return min(a, b) - h * h * 0.25 / k;
}

These operators directly manipulate the distance fields based on their spatial definitions:

• Union (opUnion): The combined surface is determined by the closest object, yielding the minimum distance.
• Intersection (opIntersection): A point is only inside the intersection if it is inside both shapes, meaning the limiting boundary is the one furthest away (the maximum).
• Subtraction (opSubtraction): Subtracting the second shape from the first involves inverting the second field (turning inside to outside by negating $b$) and then intersecting them.
• Smooth Union (opSmoothUnion): While standard CSG operators introduce non-differentiable $C^0$ seams at the intersection boundaries, smooth operators replace this abrupt transition with a polynomial blending curve of width $k$. While smooth blending slightly violates the Eikonal equation near the transition, it reliably maintains the Lipschitz bound.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
uniform float u_time;
uniform vec2 u_mouse;

out vec4 fragColor;

// Signed distance to a 2D circle
float sdCircle(in vec2 p, in float r) {
    return length(p) - r;
}

void main() {
    // Normalize pixel coordinates (from -1 to 1 on the y-axis)
    vec2 p = (2.0 * gl_FragCoord.xy - u_resolution.xy) / u_resolution.y;
    
    vec2 probe;
    
    // If mouse is at (0,0)[meaning not hovered], animate the probe
    if (u_mouse.x <= 0.0 && u_mouse.y <= 0.0) {
        probe = vec2(cos(u_time * 0.8) * 0.8, sin(u_time * 0.5) * 0.6);
    } else {
        // Map mouse pixel coordinates to normalized screen space
        probe = (2.0 * u_mouse - u_resolution.xy) / u_resolution.y;
    }

    // Evaluate the SDF for the current pixel
    float d = sdCircle(p, 0.5);
    
    // Base Colors (Matched exactly to LaTeX RGB values)
    // Outside (Blue): RGB(80, 160, 240) -> vec3(0.314, 0.627, 0.941)
    // Inside (Red): RGB(235, 80, 60) -> vec3(0.922, 0.314, 0.235)
    vec3 col = (d > 0.0) ? vec3(0.314, 0.627, 0.941) : vec3(0.922, 0.314, 0.235);
    
    // Add dark gradients near the surface (decay factor matches LaTeX)
    float decay = (d > 0.0) ? 3.0 : 4.0;
    col *= 1.0 - exp(-decay * abs(d));
    
    // Add repeating isoline waves (Math matched to LaTeX)
    col *= 0.85 + 0.15 * cos(150.0 * d);
    
    // Draw the exact surface (d = 0) as a white outline
    col = mix(col, vec3(1.0), 1.0 - smoothstep(0.0, 0.01, abs(d)));

    // Draw the SDF probe demonstrating the safe-stepping ring
    float d_probe = sdCircle(probe, 0.5);
    
    // Yellow Probe: RGB(255, 240, 0) -> vec3(1.0, 0.941, 0.0)
    vec3 yellow = vec3(1.0, 0.941, 0.0);
    
    // Draw the yellow ring (radius equals the exact distance to the surface)
    col = mix(col, yellow, 1.0 - smoothstep(0.0, 0.005, abs(length(p - probe) - abs(d_probe)) - 0.0025));
    
    // Draw the yellow center point
    col = mix(col, yellow, 1.0 - smoothstep(0.0, 0.005, length(p - probe) - 0.015));

    fragColor = vec4(col, 1.0);
}

Calculating Normals for SDFs

Unlike explicit meshes, implicit surfaces do not inherently store normal vectors. However, the surface normal at any point $\mathbf{p}$ on the boundary $\partial \Omega$ can be derived directly from the distance field. The normal vector $\mathbf{n}$ points in the direction of steepest ascent, which corresponds to the gradient of the SDF, denoted as $\nabla f(\mathbf{p})$.

By definition of the Eikonal equation ($\|\nabla f(\mathbf{p})\| = 1$), the gradient of a mathematically perfect SDF is already a unit vector:

In a shader implementation, this continuous derivative is typically approximated using finite differences. The Central Difference method evaluates the SDF at infinitesimally small offsets $\varepsilon$ along each Cartesian axis:

This formulation requires six evaluations of the distance function. For real-time applications where rendering performance is critical, a Forward Difference method is often preferred. This approach reuses the evaluation at the center point $\mathbf{p}$, reducing the total number of evaluations to four:

(Note: Although a pure SDF yields a unit-length gradient, numerical precision limits, smooth CSG operations, and domain distortions introduce slight deviations from the perfect Eikonal property. Consequently, the resulting normal vector must be explicitly normalized prior to shading calculations).

The following shader demonstrates these concepts by rendering a Torus using sphere tracing, with surface normals calculated dynamically via the forward difference method.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
uniform float u_time;
out vec4 fragColor;

// SDF for a Torus
float sdTorus(vec3 p, vec2 t) {
    vec2 q = vec2(length(p.xz) - t.x, p.y);
    return length(q) - t.y;
}

// Scene SDF (Combines all objects)
float map(vec3 p) {
    // Rotate the point slightly over time for animation
    float s = sin(u_time);
    float c = cos(u_time);
    mat3 rotX = mat3(1, 0, 0, 0, c, -s, 0, s, c);
    mat3 rotY = mat3(c, 0, s, 0, 1, 0, -s, 0, c);
    
    vec3 pt = rotY * rotX * p;
    return sdTorus(pt, vec2(1.0, 0.4));
}

// Forward difference normal approximation
vec3 calcNormal(vec3 p) {
    const float eps = 0.001;
    const vec2 h = vec2(eps, 0);
    
    // Evaluate the SDF at the point and at small offsets
    float base = map(p);
    return normalize(vec3(
        map(p + h.xyy) - base,
        map(p + h.yxy) - base,
        map(p + h.yyx) - base
    ));
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 0.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));

    // Sphere Tracing loop
    float t = 0.0;
    float max_d = 20.0;
    
    vec3 col = vec3(0.05); // Background
    
    for(int i = 0; i < 100; i++) {
        vec3 p = ro + rd * t;
        float d = map(p);
        
        // If we are close enough to the surface, it is a hit
        if(abs(d) < 0.001) {
            vec3 n = calcNormal(p);
            // Map normal to [0,1] range for visualization
            col = n * 0.5 + 0.5;
            break;
        }
        // If we marched too far, the ray missed the object
        if(t > max_d) break;
        
        // Step forward by the safe distance bound
        t += d;
    }

    fragColor = vec4(col, 1.0);
}

Shading Models

Intersection algorithms tell us where a ray hits and which way the surface faces. Shading models answer the next question: what color should that point be?

When light strikes a surface, part of its energy is absorbed and part is scattered. The physics of this interaction is captured by the Rendering Equation:

The outgoing radiance $L_o$ toward the viewer equals the emitted light $L_e$ plus the integral of all incoming light $L_i$ scattered by the surface. We integrate over $\Omega$ (the upper hemisphere defined by the surface normal) because opaque surfaces only receive light from their exterior side. The function $f_s$ is the Bidirectional Reflectance Distribution Function (BRDF), which encodes how a material distributes incoming light, and the dot product $(\mathbf{n} \cdot \boldsymbol{\omega}_i)$ accounts for foreshortening via Lambert’s cosine law.

Evaluating this integral is computationally immense — it involves infinite recursive scattering between all objects in the scene. Real-time graphics have therefore historically relied on simpler, empirical models that approximate the result locally.

Gamma Correction: Before we look at the code, it is important to remember that physical light calculations are linear, but our computer monitors display light non-linearly. To fix this, we apply a $1/2.2$ power function to our final color at the very end of our shaders—a process known as gamma correction.

Light Sources and Attenuation

Before applying any reflection model, we first need to know how much light energy actually reaches the surface.

• Directional (Infinite) Lights: Placed infinitely far away (like the sun). All rays are parallel, and intensity does not diminish over distance.
• Point Lights: Radiate equally in all directions from a specific point.
• Spot Lights: Radiate from a point, but restricted to a specific cone/direction.

For point and spot lights, energy disperses as it travels through space. Physically, light intensity falls off according to the inverse-square law ($1/d^2$). However, for artistic control in rendering pipelines, this distance attenuation is generally modeled using constant ($k_c$), linear ($k_l$), and quadratic ($k_q$) terms:

The Phong Reflection Model

With the incoming light intensity established, we can now compute how the surface responds. The most famous empirical model is the Phong reflection model, which decomposes local illumination into three physically motivated components: Ambient, Diffuse, and Specular.

1. Ambient Light

In the real world, light bounces off walls, floors, and other objects, filling the environment with indirect illumination. Instead of simulating these complex global bounces, the ambient term provides a constant base level of light so that objects in shadow are not completely pitch black.

Here, $I_a$ is the intensity of the ambient light, and $k_a$ is the material’s ambient reflection coefficient.

2. Lambertian Diffuse

A perfectly matte (diffuse) surface, like chalk or unpolished wood, scatters light equally in all directions. Because the scattered light is uniform, the appearance of a diffuse surface does not depend on where the camera is positioned.

However, it does depend heavily on the surface orientation relative to the light. If you hold a flashlight directly above a surface, the beam is concentrated. If you tilt the surface, the same beam spreads over a larger area, decreasing the intensity per unit area. This is governed by Lambert’s cosine law, which we calculate using the dot product between the normalized surface normal $\mathbf{n}$ and the normalized light direction ${\color{#FF9800}\mathbf{l}}$:

We clamp the dot product at zero because a negative value means the light is hitting the back of the surface, which should contribute no illumination.

3. Specular Highlights

Smooth surfaces, like polished metal or wet plastic, are not perfect diffusers. They exhibit specular highlights—bright spots where the light reflects strongly in a specific direction. Unlike diffuse reflection, specular reflection is highly dependent on the viewer’s position.

A perfect mirror reflects light exactly along the reflection vector ${\color{#4CAF50}\mathbf{r}}$. Using basic vector projection, we can compute ${\color{#4CAF50}\mathbf{r}}$ by taking the incoming light vector ${\color{#FF9800}\mathbf{l}}$, projecting it onto the normal $\mathbf{n}$, and reflecting it across the normal:

For surfaces that are shiny but not perfect mirrors, the reflected light scatters in a tight lobe around ${\color{#4CAF50}\mathbf{r}}$. As the angle $\phi$ between the reflection vector ${\color{#4CAF50}\mathbf{r}}$ and the view direction ${\color{#2196F3}\mathbf{v}}$ increases, the intensity drops. Phong modeled this drop-off using a cosine power function:

The exponent $\alpha$ is the shininess coefficient. A low value (e.g., $5 - 10$) produces a broad, soft highlight like plastic, while a high value (e.g., $100 - 200$) produces a sharp, tight highlight typical of polished metals.

Combining all three gives the complete Phong reflection model:

Interactive Phong reflection model. Drag on the sphere to reposition the light; toggle between Phong and Blinn-Phong specular models.

Below is the Phong reflection model applied to a test scene with two spheres on a checkerboard floor. Only local shading is computed here — no shadow rays, no reflections. Notice how the floor receives uniform light everywhere, even directly beneath the spheres where a shadow should be. This is the key limitation of local-only models.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

float hitSphere(vec3 ro, vec3 rd, vec3 c, float r) {
    vec3 q = ro - c;
    float b = dot(q, rd);
    float disc = b * b - dot(q, q) + r * r;
    if (disc < 0.0) return -1.0;
    float t = -b - sqrt(disc);
    return t > 0.001 ? t : -1.0;
}

float hitPlane(vec3 ro, vec3 rd, float y) {
    float t = (y - ro.y) / rd.y;
    return t > 0.001 ? t : -1.0;
}

struct Hit { float t; vec3 n; vec3 col; };

Hit traceScene(vec3 ro, vec3 rd) {
    Hit best = Hit(1e20, vec3(0), vec3(0));
    float t;

    // Red sphere
    t = hitSphere(ro, rd, vec3(-0.8, 0.0, -1.0), 1.0);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(-0.8, 0, -1)), vec3(0.85, 0.15, 0.1));

    // Silver sphere
    t = hitSphere(ro, rd, vec3(1.0, -0.3, 0.0), 0.7);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(1, -0.3, 0)), vec3(0.9));

    // Floor
    t = hitPlane(ro, rd, -1.0);
    if (t > 0.0 && t < best.t) {
        vec3 p = ro + rd * t;
        float ck = mod(floor(p.x) + floor(p.z), 2.0);
        best = Hit(t, vec3(0, 1, 0), mix(vec3(0.3), vec3(0.7), ck));
    }
    return best;
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 1.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));
    vec3 lightPos = vec3(3.0, 4.0, 2.0);

    Hit hit = traceScene(ro, rd);
    vec3 col = vec3(0.03, 0.04, 0.06);

    if (hit.t < 1e19) {
        vec3 p = ro + rd * hit.t;
        vec3 l = normalize(lightPos - p); // Light vector
        vec3 v = normalize(ro - p);       // View vector
        vec3 r = reflect(-l, hit.n);      // Reflection vector

        // Phong: Ambient + Diffuse + dot(r, v)^alpha
        vec3 ambient  = hit.col * 0.08;
        vec3 diffuse  = hit.col * max(dot(hit.n, l), 0.0);
        float spec    = pow(max(dot(r, v), 0.0), 48.0);
        col = ambient + diffuse + vec3(0.7) * spec;
    }

    col = pow(col, vec3(1.0 / 2.2));
    fragColor = vec4(col, 1.0);
}

The Blinn-Phong Modification

Calculating the exact reflection vector $\mathbf{r}$ at every pixel requires several vector operations. Jim Blinn proposed a cheaper approximation that produces nearly identical results: the halfway vector.

Instead of finding the angle between the reflection ${\color{#4CAF50}\mathbf{r}}$ and the viewer ${\color{#2196F3}\mathbf{v}}$, we calculate the vector that sits exactly halfway between the light ${\color{#FF9800}\mathbf{l}}$ and the viewer ${\color{#2196F3}\mathbf{v}}$:

If the halfway vector ${\color{#9C27B0}\mathbf{h}}$ perfectly aligns with the surface normal $\mathbf{n}$, the highlight is at its maximum. The specular term becomes $\max(0, \mathbf{n} \cdot {\color{#9C27B0}\mathbf{h}})^\beta$. (Note that to match the visual size of a Phong highlight, the Blinn-Phong exponent $\beta$ must be roughly $4$ times larger than the Phong exponent $\alpha$).

This was a massive optimization for early rendering pipelines: if the light and the camera are infinitely far away, ${\color{#FF9800}\mathbf{l}}$ and ${\color{#2196F3}\mathbf{v}}$ are constant, meaning ${\color{#9C27B0}\mathbf{h}}$ only needs to be computed once per scene, not once per pixel.

The same scene as above, now shaded with Blinn-Phong. The only code change is in the specular calculation — dot(n, h)^β replaces dot(r, v)^α. The highlights are slightly rounder and wider.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

float hitSphere(vec3 ro, vec3 rd, vec3 c, float r) {
    vec3 q = ro - c;
    float b = dot(q, rd);
    float disc = b * b - dot(q, q) + r * r;
    if (disc < 0.0) return -1.0;
    float t = -b - sqrt(disc);
    return t > 0.001 ? t : -1.0;
}

float hitPlane(vec3 ro, vec3 rd, float y) {
    float t = (y - ro.y) / rd.y;
    return t > 0.001 ? t : -1.0;
}

struct Hit { float t; vec3 n; vec3 col; };

Hit traceScene(vec3 ro, vec3 rd) {
    Hit best = Hit(1e20, vec3(0), vec3(0));
    float t;

    // Red sphere
    t = hitSphere(ro, rd, vec3(-0.8, 0.0, -1.0), 1.0);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(-0.8, 0, -1)), vec3(0.85, 0.15, 0.1));

    // Silver sphere
    t = hitSphere(ro, rd, vec3(1.0, -0.3, 0.0), 0.7);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(1, -0.3, 0)), vec3(0.9));

    // Floor
    t = hitPlane(ro, rd, -1.0);
    if (t > 0.0 && t < best.t) {
        vec3 p = ro + rd * t;
        float ck = mod(floor(p.x) + floor(p.z), 2.0);
        best = Hit(t, vec3(0, 1, 0), mix(vec3(0.3), vec3(0.7), ck));
    }
    return best;
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 1.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));
    vec3 lightPos = vec3(3.0, 4.0, 2.0);

    Hit hit = traceScene(ro, rd);
    vec3 col = vec3(0.03, 0.04, 0.06);

    if (hit.t < 1e19) {
        vec3 p = ro + rd * hit.t;
        vec3 l = normalize(lightPos - p); // Light vector
        vec3 v = normalize(ro - p);       // View vector
        vec3 h = normalize(l + v);        // Halfway vector

        // Blinn-Phong: Ambient + Diffuse + dot(n, h)^beta
        vec3 ambient  = hit.col * 0.08;
        vec3 diffuse  = hit.col * max(dot(hit.n, l), 0.0);
        float spec    = pow(max(dot(hit.n, h), 0.0), 192.0);
        col = ambient + diffuse + vec3(0.7) * spec;
    }

    col = pow(col, vec3(1.0 / 2.2));
    fragColor = vec4(col, 1.0);
}

Shading Paradigms: Gouraud vs. Phong

The math above tells us how to calculate color, but the pipeline must decide where to calculate it. In rasterization, there are two primary ways to apply these equations across a triangle mesh:

1. Gouraud Shading (Vertex Shading): The lighting equations are evaluated only at the three vertices of a triangle. The resulting colors are then linearly interpolated across the pixels of the triangle. This is very fast, but it suffers from severe artifacts: if a sharp specular highlight falls in the exact center of a large triangle, the vertices will all be dark, and the highlight will completely disappear!
2. Phong Shading (Pixel/Fragment Shading): Instead of interpolating colors, we linearly interpolate the normal vector $\mathbf{n}$ across the triangle’s surface. The full lighting equation is then evaluated at every single pixel. This solves the highlight problem and produces beautifully smooth curved surfaces.

Note the terminology: The Phong reflection model (the math equation for ambient/diffuse/specular light) is fundamentally different from Phong shading (the technique of interpolating normals per-pixel). You can technically evaluate the Phong reflection model using Gouraud shading, though it will look poor!

Furthermore, when linearly interpolating unit normal vectors across a triangle in a rasterizer, the resulting vectors will slightly shrink in length. Therefore, they must be re-normalized dynamically in the fragment shader before applying the reflection model.

Limitations: Self-Shadowing and Ambient Occlusion

The Phong and Blinn-Phong models evaluate lighting at each surface point independently — they have no knowledge of the rest of the scene. This means:

• No self-shadowing. A sphere sitting on a floor will not cast a shadow, because the model never asks whether anything blocks the light.
• No contact darkening. Where two surfaces meet (e.g., the base of a sphere on a floor), real-world light is partially occluded by nearby geometry. Local models miss this entirely.
• The ambient term is a flat constant. Every surface point receives the same $k_a I_a$, regardless of how enclosed or exposed it is.

Ambient Occlusion (AO) is a popular technique that addresses the third problem without requiring full global illumination. The idea is simple: at each surface point, estimate what fraction of the surrounding hemisphere is blocked by nearby geometry. Points in corners, crevices, and contact regions score low (more occluded), while exposed surfaces score high.

Here $V(\mathbf{p}, \boldsymbol{\omega})$ is a binary visibility function — $1$ if the hemisphere direction $\boldsymbol{\omega}$ is unblocked, $0$ if nearby geometry occludes it. The result $A(\mathbf{p})$ is then multiplied into the ambient (and sometimes diffuse) term, darkening the tucked-away areas that a flat ambient constant would otherwise over-illuminate.

In screen-space implementations like SSAO (Screen-Space Ambient Occlusion), the visibility check is approximated by sampling the depth buffer around each pixel — a fast, purely post-processing approach that requires no scene knowledge at all. In ray-traced pipelines, short-range occlusion rays can compute AO directly.

Here is our test scene rendered with local Phong shading plus an Ambient Occlusion term computed analytically for each sphere. Notice how contact shadows naturally appear underneath the spheres where they approach the floor plane.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

float hitSphere(vec3 ro, vec3 rd, vec3 c, float r) {
    vec3 q = ro - c;
    float b = dot(q, rd);
    float disc = b * b - dot(q, q) + r * r;
    if (disc < 0.0) return -1.0;
    float t = -b - sqrt(disc);
    return t > 0.001 ? t : -1.0;
}

float hitPlane(vec3 ro, vec3 rd, float y) {
    float t = (y - ro.y) / rd.y;
    return t > 0.001 ? t : -1.0;
}

struct Hit { float t; vec3 n; vec3 col; };

Hit traceScene(vec3 ro, vec3 rd) {
    Hit best = Hit(1e20, vec3(0), vec3(0));
    float t;

    // Red sphere
    t = hitSphere(ro, rd, vec3(-0.8, 0.0, -1.0), 1.0);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(-0.8, 0, -1)), vec3(0.85, 0.15, 0.1));

    // Silver sphere
    t = hitSphere(ro, rd, vec3(1.0, -0.3, 0.0), 0.7);
    if (t > 0.0 && t < best.t)
        best = Hit(t, normalize(ro + rd*t - vec3(1, -0.3, 0)), vec3(0.9));

    // Floor
    t = hitPlane(ro, rd, -1.0);
    if (t > 0.0 && t < best.t) {
        vec3 p = ro + rd * t;
        float ck = mod(floor(p.x) + floor(p.z), 2.0);
        best = Hit(t, vec3(0, 1, 0), mix(vec3(0.3), vec3(0.7), ck));
    }
    return best;
}

// Analytic Ambient Occlusion of a sphere
float sphereAO(vec3 p, vec3 n, vec3 c, float r) {
    vec3 v = c - p;
    float d = length(v);
    return max(0.0, dot(n, v / d)) * (r * r) / (d * d);
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 1.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));
    vec3 lightPos = vec3(3.0, 4.0, 2.0);

    Hit hit = traceScene(ro, rd);
    vec3 col = vec3(0.03, 0.04, 0.06);

    if (hit.t < 1e19) {
        vec3 p = ro + rd * hit.t;
        vec3 l = normalize(lightPos - p); // Light vector
        vec3 v = normalize(ro - p);       // View vector
        vec3 r = reflect(-l, hit.n);      // Reflection vector

        // Compute analytical AO from both spheres
        float ao = 1.0;
        ao -= sphereAO(p, hit.n, vec3(-0.8, 0.0, -1.0), 1.0);
        ao -= sphereAO(p, hit.n, vec3(1.0, -0.3, 0.0), 0.7);
        ao = clamp(ao, 0.0, 1.0);

        // Combine with local Phong lighting
        vec3 ambient  = hit.col * 0.08 * ao;
        vec3 diffuse  = hit.col * max(dot(hit.n, l), 0.0) * ao;
        float spec    = pow(max(dot(r, v), 0.0), 48.0) * ao;
        col = ambient + diffuse + vec3(0.7) * spec;
    }

    col = pow(col, vec3(1.0 / 2.2));
    fragColor = vec4(col, 1.0);
}

We will see full shadow and reflection rays in the next section.

Ray Tracing

The shading models above are strictly local—they only consider the light hitting a surface directly from a source. Whitted ray tracing extends this idea into a global context by tracing visibility recursively. Once a primary ray hits a surface, the rendering algorithm spawns new secondary rays to evaluate the global environment.

Instead of just intersection math, Ray Tracing is defined by this recursive routing logic:

1. Shadow Rays (Visibility Testing)

To determine if a point $\mathbf{p}$ is in shadow, we don’t guess. We cast a shadow ray from $\mathbf{p}$ directly toward the light source ${\color{#FF9800}\mathbf{l}}$.

• If the ray hits an object before it reaches the light (${\color{#795548}t_{\text{hit}}} < {\color{#FF9800}t_{\text{light}}}$), the point is occluded. We skip the diffuse and specular calculations and only apply ambient light.
• To prevent a surface from shadowing itself due to floating-point inaccuracies (Shadow Acne), we offset the ray origin slightly along the normal: $\mathbf{p}_{\text{origin}} = \mathbf{p} + \epsilon\mathbf{n}$.

2. Reflection Rays

For shiny surfaces like mirrors or water, we calculate the reflection vector ${\color{#4CAF50}\mathbf{r}}$ exactly as we did in the Phong model:

(Where ${\color{#2196F3}\mathbf{v}}$ points from the surface back to the camera). We then cast a new ray from $\mathbf{p} + \epsilon\mathbf{n}$ in the direction of ${\color{#4CAF50}\mathbf{r}}$. Whatever color that ray eventually hits is added to the current pixel’s color.

3. Refraction Rays (Transmission)

For transparent materials like glass, light bends as it passes through the surface boundary. According to Snell’s Law, the angle of incidence $\theta_L$ and the angle of transmission $\theta_T$ are related by the indices of refraction ($\eta$):

The exact transmission vector ${\color{#9C27B0}\mathbf{t}}$ can be derived from the normal $\mathbf{n}$ and the incoming view vector ${\color{#2196F3}\mathbf{v}}$. If the light enters a medium with a lower index of refraction at a steep enough angle, the math under the square root becomes negative. This physical phenomenon is Total Internal Reflection—no light is refracted; it is all reflected inside the object.

Here is the same scene rendered with full recursive ray tracing. Shadow rays now block direct light where occluded, and reflection rays on the chrome sphere pick up the environment. Compare with the local-only shaders above.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

float hitSphere(vec3 ro, vec3 rd, vec3 c, float r) {
    vec3 q = ro - c;
    float b = dot(q, rd);
    float disc = b * b - dot(q, q) + r * r;
    if (disc < 0.0) return -1.0;
    float t = -b - sqrt(disc);
    if (t < 0.001) t = -b + sqrt(disc);
    return t > 0.001 ? t : -1.0;
}

float hitPlane(vec3 ro, vec3 rd, float y) {
    float t = (y - ro.y) / rd.y;
    return t > 0.001 ? t : -1.0;
}

struct Hit { float t; vec3 n; vec3 col; float refl; };

Hit traceScene(vec3 ro, vec3 rd) {
    Hit best = Hit(1e20, vec3(0), vec3(0), 0.0);
    vec3 n; float t;

    // Red sphere
    t = hitSphere(ro, rd, vec3(-0.8, 0.0, -1.0), 1.0);
    if (t > 0.0 && t < best.t) { best = Hit(t, normalize(ro+rd*t - vec3(-0.8,0,-1)), vec3(0.85,0.15,0.1), 0.05); }

    // Chrome sphere
    t = hitSphere(ro, rd, vec3(1.0, -0.3, 0.0), 0.7);
    if (t > 0.0 && t < best.t) { best = Hit(t, normalize(ro+rd*t - vec3(1.0,-0.3,0.0)), vec3(0.9), 0.8); }

    // Floor plane
    t = hitPlane(ro, rd, -1.0);
    if (t > 0.0 && t < best.t) {
        vec3 p = ro + rd * t;
        float check = step(0.5, fract(p.x * 0.5)) + step(0.5, fract(p.z * 0.5));
        check = mod(check, 2.0);
        best = Hit(t, vec3(0,1,0), mix(vec3(0.3), vec3(0.7), check), 0.1);
    }
    return best;
}

vec3 shade(Hit hit, vec3 p, vec3 rd) {
    vec3 lightPos = vec3(3.0, 4.0, 2.0);
    vec3 l = normalize(lightPos - p);
    vec3 v = -rd;

    // Shadow ray
    Hit sh = traceScene(p + hit.n * 0.01, l);
    float shadow = (sh.t < length(lightPos - p)) ? 0.15 : 1.0;

    float diff = max(dot(hit.n, l), 0.0) * shadow;
    float spec = pow(max(dot(reflect(-l, hit.n), v), 0.0), 64.0) * shadow;
    return vec3(0.06) + hit.col * diff + vec3(0.7) * spec;
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 1.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));

    vec3 col = vec3(0.0);
    float atten = 1.0;

    // 2 bounces: primary + one reflection
    for (int bounce = 0; bounce < 2; bounce++) {
        Hit hit = traceScene(ro, rd);
        if (hit.t > 1e19) {
            col += atten * vec3(0.03, 0.04, 0.06);
            break;
        }
        vec3 p = ro + rd * hit.t;
        col += atten * (1.0 - hit.refl) * shade(hit, p, rd);

        // Reflection ray for next bounce
        atten *= hit.refl;
        if (atten < 0.01) break;
        ro = p + hit.n * 0.01;
        rd = reflect(rd, hit.n);
    }

    col = pow(col, vec3(1.0 / 2.2));
    fragColor = vec4(col, 1.0);
}

Path Tracing

Whitted ray tracing handles mirrors and hard shadows well, but it cannot capture the soft, diffuse inter-reflections that make real scenes look natural — colored light bleeding from a red wall onto a white floor, for example. Path tracing solves this by actually evaluating the rendering equation via Monte Carlo integration.

Instead of spawning deterministic shadow and reflection rays, the path tracer shoots the ray in a random direction sampled from the hemisphere at each bounce. Over hundreds or thousands of samples per pixel, the noisy estimates converge to the true solution of the rendering equation. This is the same technique used by production renderers in film (Arnold, RenderMan, Cycles).

The key ideas — Monte Carlo integration, importance sampling, and the mathematics of BRDFs — are covered in depth in the next post in this series: A comprehensive look into the Importance Sampling and Path Guiding for Path Tracing

The interactive below shows the same scene under different rendering paradigms. You can compare the local shading models (Phong and Blinn-Phong) — which calculate light per surface and ignore inter-object visibility — with global illumination models (Ray Tracing and Path Tracing), which trace ray paths recursively to simulate shadows, glossy reflections, and diffuse light transport.

Sphere Tracing

As described in the SDF section above, Sphere Tracing (or Ray Marching) is the algorithm that turns the distance returned by an SDF into safe steps along a ray. Instead of solving an analytic equation, we evaluate the field at the current point, advance by that distance, and repeat.

At the level of one march step, the update rule is:

That equation is the heart of the method: the SDF tells us how far we can move along the ray without crossing a surface, so each iteration makes a safe, adaptive jump.

The core loop is simple: march until the distance becomes smaller than an epsilon, or until the ray travels farther than the scene limit. The shader below is just one concrete implementation of that loop for a small scene.

In pseudocode, the algorithm is easy to think about:

SphereTrace(ray origin ro, ray direction rd):
    t = 0
    for i in 1..maxSteps:
        p = ro + t * rd
        d = sceneSDF(p)
        if d < epsilon:
            return hit at p
        t = t + d
        if t > maxDistance:
            break
    return miss

The key idea is that the current SDF value is a guaranteed empty radius around the sample point, so the next evaluation can jump forward by exactly that amount. If the field is exact, this is extremely efficient. If the field is only a bound or has been distorted, the same loop still works, but you usually introduce a safety factor so the steps stay conservative.

If you want to visualize what is happening, it helps to think of the ray as a sequence of expanding circles centered on the current sample point. Each march step uses the SDF value as the radius of the next safe jump, so the circles shrink as the ray gets close to the surface.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
uniform float u_time;
uniform vec2 u_mouse;

out vec4 fragColor;

float sdCircle(vec2 p, float r) {
    return length(p) - r;
}

float sceneSDF(vec2 p) {
    float c1 = sdCircle(p - vec2(-0.45, 0.10), 0.35);
    float c2 = sdCircle(p - vec2(0.42, -0.24), 0.22);
    float fl  = p.y + 0.72;
    return min(min(c1, c2), fl);
}

void main() {
    vec2 p = (2.0 * gl_FragCoord.xy - u_resolution.xy) / u_resolution.x;

    // Fixed ray origin at the top center
    vec2 ro = vec2(0.0, 0.8);

    // Only ray direction animates / follows mouse
    vec2 target = (u_mouse.x <= 0.0 && u_mouse.y <= 0.0)
        ? vec2(0.7 * cos(u_time * 0.55), 0.15 + 0.3 * sin(u_time * 0.8))
        : (2.0 * u_mouse - u_resolution.xy) / u_resolution.x;

    vec2 rd = normalize(target - ro);

    // --- SDF field visualization ---
    float d = sceneSDF(p);
    vec3 col = (d > 0.0) ? vec3(0.314, 0.627, 0.941) : vec3(0.922, 0.314, 0.235);
    float decay = (d > 0.0) ? 3.0 : 4.0;
    col *= 1.0 - exp(-decay * abs(d));
    col *= 0.85 + 0.15 * cos(150.0 * d);
    col = mix(col, vec3(1.0), 1.0 - smoothstep(0.0, 0.01, abs(d)));

    // --- Ray line ---
    float rayLine = 1.0 - smoothstep(0.0, 0.012,
        length(p - (ro + clamp(dot(p - ro, rd), 0.0, 4.0) * rd)));
    col = mix(col, vec3(1.0), rayLine * 0.18);

    // --- Sphere tracing circles (20 steps) ---
    vec2 marchP = ro;
    float t = 0.0;
    for (int i = 0; i < 20; i++) {
        float sd  = sceneSDF(marchP);
        float asd = abs(sd);

        float ring     = abs(length(p - marchP) - asd);
        float ringGlow = 1.0 - smoothstep(0.0, 0.008, ring - 0.003);
        float frac     = float(i) / 19.0;
        vec3  ringCol  = mix(vec3(0.3, 0.75, 1.0), vec3(1.0, 0.85, 0.2), frac);
        col = mix(col, ringCol, ringGlow * 0.9);

        float dot1 = 1.0 - smoothstep(0.0, 0.025, length(p - marchP) - 0.012);
        col = mix(col, vec3(1.0), dot1 * 0.9);

        if (sd < 0.001 || t > 4.0) break;
        t     += asd;
        marchP = ro + t * rd;
    }

    // --- Ray origin dot ---
    float roGlow = 1.0 - smoothstep(0.0, 0.03, length(p - ro) - 0.018);
    col = mix(col, vec3(1.0, 0.9, 0.3), roGlow);

    fragColor = vec4(col, 1.0);
}

Here is the same two-sphere scene, shown through that marching loop.

#version 300 es
precision highp float;

uniform vec2 u_resolution;
out vec4 fragColor;

// Scene defined via Signed Distance Functions
float map(vec3 p) {
    float s1 = length(p - vec3(-0.8, 0.0, -1.0)) - 1.0;  // Red sphere
    float s2 = length(p - vec3(1.0, -0.3, 0.0)) - 0.7;   // Silver sphere
    float fl = p.y + 1.0;                                  // Floor
    return min(min(s1, s2), fl);
}

// Gradient-based normal estimation
vec3 getNormal(vec3 p) {
    vec2 e = vec2(0.001, 0.0);
    return normalize(vec3(
        map(p + e.xyy) - map(p - e.xyy),
        map(p + e.yxy) - map(p - e.yxy),
        map(p + e.yyx) - map(p - e.yyx)
    ));
}

// Shadow ray via sphere tracing
float calcShadow(vec3 ro, vec3 rd) {
    float t = 0.02;
    for (int i = 0; i < 30; i++) {
        float h = map(ro + rd * t);
        if (h < 0.001) return 0.15;
        t += h;
    }
    return 1.0;
}

// Identify which object was hit for coloring
vec3 getColor(vec3 p) {
    float s1 = length(p - vec3(-0.8, 0.0, -1.0)) - 1.0;
    float s2 = length(p - vec3(1.0, -0.3, 0.0)) - 0.7;
    float fl = p.y + 1.0;
    if (s1 < s2 && s1 < fl) return vec3(0.85, 0.15, 0.1);
    if (s2 < fl) return vec3(0.9);
    float ck = mod(floor(p.x) + floor(p.z), 2.0);
    return mix(vec3(0.3), vec3(0.7), ck);
}

void main() {
    vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution.xy) / u_resolution.y;
    vec3 ro = vec3(0.0, 1.0, 4.0);
    vec3 rd = normalize(vec3(uv, -1.0));
    vec3 lightPos = vec3(3.0, 4.0, 2.0);

    vec3 col = vec3(0.03, 0.04, 0.06);

    // Sphere Tracing
    float t = 0.0;
    bool hit = false;
    for (int i = 0; i < 80; i++) {
        float d = map(ro + rd * t);
        if (d < 0.001) { hit = true; break; }
        if (t > 20.0) break;
        t += d;
    }

    if (hit) {
        vec3 p = ro + rd * t;
        vec3 n = getNormal(p);
        vec3 c = getColor(p);
        vec3 l = normalize(lightPos - p);
        vec3 v = normalize(ro - p);
        vec3 r = reflect(-l, n);

        float shadow = calcShadow(p + n * 0.01, l);
        float diff = max(dot(n, l), 0.0) * shadow;
        float spec = pow(max(dot(r, v), 0.0), 48.0) * shadow;
        col = c * 0.08 + c * diff + vec3(0.7) * spec;
    }

    col = pow(col, vec3(1.0 / 2.2));
    fragColor = vec4(col, 1.0);
}