Introduction to 3D Rendering

The camera $z$ axis is $-\mathbf{f}$. The rotation $\mathbf{R}$ from camera space to world space has these axes as its columns (it is M_cam in the code):

A ray starts at $\mathbf{E}$, and its world-space direction is $\mathbf{R}$ times its camera-space direction. The camera-space direction comes from the pixel index, in four steps:

Ray Generation Code

DrJit ray generation code

A Camera class that generates one ray per pixel center. It uses Dr.Jit [14], the array library behind the Mitsuba renderer. Here fov is the vertical field of view.

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   # width / height: the reciprocal of a in the text
        self.H = H
        self.W = W
        self.f = 1.0   # distance to the image plane (D in the text, not the far plane)

    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))

        # Matrix3f takes rows, so transpose to get r, u, -f as columns (R in the text)
        M_cam = dr.transpose(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

        # Image-plane size at distance f: tan(fov / 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)

        # Rotate into world space
        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 

In the GLSL shader the pixel coordinate comes from gl_FragCoord. Dividing by u_resolution.y makes $y$ run from $-0.5$ to $0.5$ and scales $x$ by the aspect ratio, so with the image plane at distance $D = 1$ (f in the code) the vertical field of view is $2\tan^{-1}(0.5) \approx 53^\circ$. The shader shows each world-space direction as a color, mapped from $[-1, 1]$ to $[0, 1]$.

#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() {
    // y in [-0.5, 0.5]; x scaled by the aspect ratio
    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 * 0.5 + 0.5, 1.0);
}

A simpler camera for the remaining shaders

In the remaining shaders the camera is translated but not rotated, so it looks down $-z$, M_cam is the identity matrix and the ray direction can be built 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)); // already a world-space direction

Anti-Aliasing

One ray through the pixel center aliases hard edges into staircases (“jaggies”). A pixel should show the average color over its area. Averaging $N$ rays through random positions inside the pixel is a Monte Carlo estimate of that average.

Let $(p_x, p_y)$ be the pixel center on the image plane at distance $D$, as computed in the ray generation code, and let $s$ be the size of one pixel there (height / H in the code). Sample $i$ moves the point by offsets $\delta_{x,i}$ and $\delta_{y,i}$ drawn uniformly from $[-s/2, s/2]$:

Here $L(\mathbf{E}, \mathbf{d}_i)$ is the radiance (color) returned by the ray with origin $\mathbf{E}$ and direction $\mathbf{d}_i$. Uniform offsets estimate the plain average over the pixel square, a box filter.

Production renderers usually place these samples in a stratified or blue-noise pattern, or accumulate them over many frames. The shaders in this post use one ray through each pixel center.

Geometry Representations and Intersections

How a ray is intersected with the scene depends on how the geometry is represented. Three common representations are:

• Analytic primitives such as spheres, planes and cylinders are given by equations. Their intersection with a ray usually has a closed-form solution. They are cheap to intersect but can only model simple shapes.
• Triangle meshes approximate any surface with flat triangles. Rasterizers draw them, and ray tracers intersect them with a ray-triangle test.
• Signed distance functions (SDFs) give, for every point in space, the distance to the nearest surface: positive outside the object, zero on the surface and negative inside. A ray is intersected with an SDF iteratively, by sphere tracing [6].

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:

A point on the ray, with origin ${\color{#3273F6}\mathbf{o}}$ and unit direction ${\color{#FF4B4B}\mathbf{d}}$, is

Substituting the ray into the sphere equation gives

Let $\mathbf{q} = {\color{#3273F6}\mathbf{o}} - {\color{#109618}\mathbf{c}}$, the vector from the center to the ray origin. Writing the squared length as a dot product and expanding gives a quadratic equation in $t$:

This is $At^2 + Bt + C = 0$ with

• $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$

The roots are $t = \frac{-B \pm \sqrt{B^2 - 4AC}}{2A}$.

• If the discriminant $B^2 - 4AC$ is negative, the ray misses the sphere.
• If it is zero, the ray touches the sphere at one point.
• If it is positive, the ray crosses the sphere at two points. The smallest positive root is the first visible hit. If only the larger root is positive, the ray starts inside the sphere.

With $A = 1$ and $b = B/2 = \mathbf{q} \cdot {\color{#FF4B4B}\mathbf{d}}$, the roots simplify to $t = -b \pm \sqrt{h}$ with $h = b^2 - C$, a quarter of the discriminant. The shader below uses this form.

The hit point is $\mathbf{p} = {\color{#3273F6}\mathbf{o}} + t{\color{#FF4B4B}\mathbf{d}}$, and the unit normal there is $\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 the smallest positive t, or -1.0 on a miss. Also outputs the unit normal.
float intersectSphere(vec3 ro, vec3 rd, vec3 center, float radius, out vec3 normal) {
    vec3 q = ro - center;
    float b = dot(q, rd);                  // B / 2 (A = 1 because rd is normalized)
    float c = dot(q, q) - radius * radius; // C
    float h = b * b - c;                   // discriminant / 4
    if (h < 0.0) return -1.0;              // ray misses

    float t = -b - sqrt(h);                // nearer root
    if (t <= 0.0) t = -b + sqrt(h);        // ray starts inside the sphere
    if (t <= 0.0) return -1.0;             // sphere is behind the ray

    normal = (ro + t * rd - center) / radius;
    return t;
}

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)

Let a triangle have vertices ${\color{#FF9900}\mathbf{a}}$, ${\color{#0099C6}\mathbf{b}}$ and ${\color{#DD4477}\mathbf{c}}$. Any point $\mathbf{p}$ in its plane can be written with barycentric coordinates $(u, v)$:

With the edge vectors $\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 meets the plane of the triangle where

Moving the unknowns $t$, $u$, $v$ to the left gives three linear equations:

Let $\mathbf{s} = {\color{#3273F6}\mathbf{o}} - {\color{#FF9900}\mathbf{a}}$, the vector from $\mathbf{a}$ to the ray origin. Cramer’s rule solves the system with determinants. The determinant of the matrix with columns $\mathbf{x}$, $\mathbf{y}$, $\mathbf{z}$ is the scalar triple product $\det(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \mathbf{x} \cdot (\mathbf{y} \times \mathbf{z})$, the signed volume of the parallelepiped they span. The solution is

The Möller-Trumbore algorithm [5] evaluates these determinants with two cross products, using $\det(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \mathbf{x} \cdot (\mathbf{y} \times \mathbf{z}) = \mathbf{y} \cdot (\mathbf{z} \times \mathbf{x})$:

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

The ray hits the triangle (its interior or an edge) if $u \ge 0$, $v \ge 0$, $u + v \le 1$ and $t > 0$. The plane’s normal is never computed.

#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 g = cross(rd, e2);
    float det = dot(e1, g);

    // det near zero: the ray is parallel to the triangle's plane
    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, g) * invDet;
    if (u < 0.0 || u > 1.0) return false;

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

    // Calculate t
    t = dot(e2, k) * invDet;
    
    bary = vec3(1.0 - u - v, u, v);   // weights of a, b, c
    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)

A signed distance function (SDF) describes a surface implicitly: instead of storing vertices, it gives the distance from any point to the nearest surface.

Let $S \subset \mathbb{R}^3$ be a closed set (the solid object) with boundary surface $\partial S$. The signed distance function is the Euclidean distance from $\mathbf{p}$ to $\partial S$, with a sign that tells whether $\mathbf{p}$ is outside or inside:

The surface is the zero level set $\partial S = \{\mathbf{p} \in \mathbb{R}^3 \mid f(\mathbf{p}) = 0\}$.

If $f(\mathbf{p}) = d > 0$, the sphere of radius $d$ around $\mathbf{p}$ contains no surface. A ray can therefore advance by $d$ without passing through a surface, which is the idea behind sphere tracing.

The Eikonal Equation and Lipschitz Continuity

An exact signed distance function satisfies the eikonal equation

everywhere except where the gradient is undefined (for example on the medial axis, the points with more than one nearest surface point). Moving one unit directly away from the surface increases $f$ by exactly one.

Operations such as non-uniform scaling, twisting or blending usually break this property. Sphere tracing stays safe as long as $f$ is 1-Lipschitz:

Then $f$ never overestimates the distance to the surface. If $\|\nabla f\| > 1$ somewhere, $f$ can overestimate it, and the ray can step through the surface and miss geometry. If $\|\nabla f\|$ is much smaller than 1, the steps are safe but short, and the ray needs many more of them.

Constructive Solid Geometry (CSG)

Shapes are combined with Boolean operations. Quilez lists SDFs for many 2D [7] and 3D [8] primitives, and the operators below. If $a = f_1(\mathbf{p})$ and $b = f_2(\mathbf{p})$ are the signed distances to two shapes, the constructive solid geometry (CSG) operations are:

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;
}

• Union: the nearer shape decides, so the distance is the minimum.
• Intersection: a point is inside only if it is inside both shapes, so the larger value decides.
• Subtraction: negating $b$ swaps the inside and outside of the second shape, then the two are intersected.
• Smooth union: min has a crease where the two shapes meet. The smooth version blends the two distances with a quadratic wherever $|a - b| < 4k$, subtracting at most $k$ (at $a = b$). The result is no longer an exact distance near the blend, but it remains 1-Lipschitz.

Union is exact outside the shapes. Intersection and subtraction give a lower bound on the distance rather than the exact distance; a bound never overestimates the distance, which is enough for sphere tracing.

#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() {
    // Pixel coordinates, y from -1 to 1
    vec2 p = (2.0 * gl_FragCoord.xy - u_resolution.xy) / u_resolution.y;
    
    vec2 probe;
    
    // Without the mouse over the canvas, 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);
    
    // Outside blue, inside red (same colors as the figure)
    vec3 col = (d > 0.0) ? vec3(0.314, 0.627, 0.941) : vec3(0.922, 0.314, 0.235);

    // Darken toward the surface
    float decay = (d > 0.0) ? 3.0 : 4.0;
    col *= 1.0 - exp(-decay * abs(d));

    // Contour lines of constant distance
    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)));

    // Probe point and its empty circle
    float d_probe = sdCircle(probe, 0.5);
    
    vec3 yellow = vec3(1.0, 0.941, 0.0);

    // Circle of radius |f(probe)|: it touches the surface but does not cross it
    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

An implicit surface stores no normals, but the normal at a surface point $\mathbf{p}$ is the direction in which $f$ increases fastest, the gradient $\nabla f(\mathbf{p})$. For an exact SDF the eikonal equation makes it a unit vector:

Shaders approximate the gradient with finite differences. Central differences evaluate $f$ at a small offset $\varepsilon$ on either side along each axis (the common factor $1/2\varepsilon$ is removed by the normalization):

This takes six evaluations of $f$. Forward differences reuse $f(\mathbf{p})$ and take four, at the cost of a less accurate estimate:

The estimate is normalized because finite differences, smooth blends and distortions all make its length differ from 1.

The shader below renders a torus by sphere tracing and computes its normals with forward differences.

#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);
        
        // Close enough to the surface: 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;
        }
        // Too far: the ray missed
        if(t > max_d) break;
        
        // Advance by the distance to the nearest surface
        t += d;
    }

    fragColor = vec4(col, 1.0);
}

Shading Models

Intersection gives the hit point and the surface normal there. A shading model gives the color of that point.

Light reaching a surface is partly absorbed and partly scattered. The rendering equation [12] describes the light leaving a point:

The radiance $L_o$ leaving $\mathbf{x}$ in direction $\boldsymbol{\omega}_o$ is the emitted radiance $L_e$ plus the incoming radiance $L_i$ reflected toward $\boldsymbol{\omega}_o$, integrated over the hemisphere $\Omega$ above the surface. The bidirectional reflectance distribution function (BRDF) $f_s$ describes how the material reflects light, and the factor $\mathbf{n} \cdot \boldsymbol{\omega}_i$ is Lambert’s cosine law.

$L_i$ at one point is the $L_o$ of other points, so the equation is recursive and has no closed-form solution for real scenes. Real-time graphics therefore uses simpler empirical models that consider only the light arriving directly from the light sources. Whitted’s ray tracer [11] combined such a local model with traced shadow, reflection and refraction rays.

Gamma correction: lighting is computed in linear units, but displays expect values encoded with a power curve. The shaders below raise the final color to the power $1/2.2$, an approximation of the sRGB encoding.

Light Sources and Attenuation

Three common kinds of light source are:

• Directional lights are infinitely far away, like the sun. All rays are parallel and the intensity does not change with distance.
• Point lights emit equally in all directions from a point.
• Spot lights emit from a point into a cone.

For point and spot lights, the light arriving at a surface falls off with the square of the distance $d$. For artistic control, renderers often use constant ($k_c$), linear ($k_l$) and quadratic ($k_q$) terms instead:

The Phong Reflection Model

The Phong reflection model [9] is an empirical model with three terms: ambient, diffuse and specular.

1. Ambient Light

Light also arrives indirectly, after bouncing off other surfaces. Instead of computing it, the ambient term adds a constant amount, so that surfaces facing away from the light are not black.

$I_a$ is the ambient light intensity and $k_a$ the material’s ambient coefficient.

2. Lambertian Diffuse

A matte (Lambertian) surface, like chalk, reflects light equally in all directions, so its brightness does not depend on the viewer.

It does depend on the angle of the incoming light. A beam that hits the surface at an angle spreads over a larger area, so less light arrives per unit area. By Lambert’s cosine law the light received is proportional to $\cos\theta = \mathbf{n} \cdot {\color{#FF9800}\mathbf{l}}$, where $\mathbf{n}$ is the unit normal and ${\color{#FF9800}\mathbf{l}}$ the unit vector toward the light:

The clamp at zero removes light that arrives from behind the surface.

3. Specular Highlights

Smooth surfaces, like polished metal or plastic, also show specular highlights, bright spots where light is reflected mostly in one direction. Where they appear depends on the viewer.

A mirror reflects light along the reflection vector ${\color{#4CAF50}\mathbf{r}}$. Projecting ${\color{#FF9800}\mathbf{l}}$ onto $\mathbf{n}$ gives $(\mathbf{n} \cdot {\color{#FF9800}\mathbf{l}})\mathbf{n}$, and mirroring ${\color{#FF9800}\mathbf{l}}$ about the normal gives

A glossy surface reflects light into a lobe around ${\color{#4CAF50}\mathbf{r}}$. Phong models the falloff with the angle $\phi$ between ${\color{#4CAF50}\mathbf{r}}$ and the unit vector ${\color{#2196F3}\mathbf{v}}$ toward the viewer by a power of its cosine:

The exponent $\alpha$ is the shininess. Small values (around 5 to 10) give a broad highlight, large values (100 to 200) a small, sharp one. The specular term is usually set to zero where $\mathbf{n} \cdot {\color{#FF9800}\mathbf{l}} \le 0$, since the surface then faces away from the light; the widgets below do this.

Together:

The Phong model on a sphere. Drag to move the light, hover to read the terms at a pixel, and switch to Blinn-Phong.

Constructing $\mathbf{r} = 2(\mathbf{n} \cdot \mathbf{l})\mathbf{n} - \mathbf{l}$ step by step.

Below, the Phong model shades a test scene with two spheres on a checkerboard floor. Only local shading is computed, with no shadow or reflection rays, so the floor is lit even directly under the spheres.

#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

Blinn [10] replaced ${\color{#4CAF50}\mathbf{r}} \cdot {\color{#2196F3}\mathbf{v}}$ with the angle between the normal and the halfway vector, the unit vector midway between ${\color{#FF9800}\mathbf{l}}$ and ${\color{#2196F3}\mathbf{v}}$:

The highlight is brightest when ${\color{#9C27B0}\mathbf{h}} = \mathbf{n}$, which happens exactly when ${\color{#4CAF50}\mathbf{r}} = {\color{#2196F3}\mathbf{v}}$. The specular term becomes $\max(0, \mathbf{n} \cdot {\color{#9C27B0}\mathbf{h}})^\beta$. When $\mathbf{n}$, ${\color{#FF9800}\mathbf{l}}$ and ${\color{#2196F3}\mathbf{v}}$ lie in one plane, the angle between $\mathbf{n}$ and ${\color{#9C27B0}\mathbf{h}}$ is half of $\phi$, so $\beta \approx 4\alpha$ gives a highlight of about the same size.

When the light and the viewer are both infinitely far away, ${\color{#FF9800}\mathbf{l}}$ and ${\color{#2196F3}\mathbf{v}}$ are constant, so ${\color{#9C27B0}\mathbf{h}}$ is computed once instead of once per pixel. This made the model cheaper than Phong’s on early hardware.

Constructing $\mathbf{h}$ and comparing it with $\mathbf{n}$.

The same scene shaded with Blinn-Phong. Only the specular line changes: $\max(0, \mathbf{n} \cdot \mathbf{h})^\beta$ with $\beta = 4\alpha = 192$ replaces $\max(0, \mathbf{r} \cdot \mathbf{v})^\alpha$ with $\alpha = 48$. The highlights have about the same size; at grazing angles the Blinn-Phong highlight is stretched rather than round.

#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

A rasterizer also has to decide where on a triangle the model is evaluated:

1. Gouraud shading evaluates it at the three vertices and interpolates the resulting colors across the triangle. It is cheap, but a highlight that falls inside a large triangle, away from its vertices, is lost.
2. Phong shading interpolates the normal across the triangle and evaluates the model at every pixel. Highlights inside triangles are kept and curved surfaces look smooth. Lengyel calls this per-pixel evaluation, done with the Blinn-Phong model, Blinn-Phong shading [2].

The Phong reflection model (the formula above) and Phong shading (per-pixel evaluation with interpolated normals) are different things. The Phong reflection model can be evaluated with Gouraud shading.

Linear interpolation between two unit normals gives a vector shorter than 1, so the interpolated normal is normalized again in the fragment shader.

Limitations of Local Shading and Ambient Occlusion

The Phong and Blinn-Phong models look at one surface point at a time and know nothing about the rest of the scene. As a result:

• There are no cast shadows. A sphere on a floor casts no shadow, because the model never checks whether anything blocks the light.
• There is no contact darkening where two surfaces meet, such as under a sphere resting on a floor.
• The ambient term is constant. Every point gets the same $k_a I_a$, however enclosed it is.

Ambient occlusion (AO) addresses the last two without computing global illumination. At each point it measures how much of the hemisphere above the surface is blocked by nearby geometry, weighted by the cosine:

$V(\mathbf{p}, \boldsymbol{\omega})$ is $1$ if direction $\boldsymbol{\omega}$ is unblocked and $0$ otherwise, and the factor $1/\pi$ makes $A = 1$ for a fully open hemisphere. The shading is multiplied by $A$, which darkens corners and contact regions.

Screen-space ambient occlusion (SSAO) estimates $V$ from the depth buffer around each pixel, as a post-process. A ray tracer can compute it directly with short occlusion rays.

For a sphere of radius $R$ whose center is at distance $d$ from $\mathbf{p}$, in a direction at angle $\theta$ to the normal, the blocked fraction of $A$ is exactly $\cos\theta \, (R/d)^2$, as long as the whole sphere is above the surface’s tangent plane. The shader below subtracts this term for each sphere from $A = 1$ and multiplies all of the Phong terms by the result, which darkens the floor under the spheres.

#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;
}

// Occlusion by a sphere: cos(theta) * (R/d)^2 (exact while the sphere is above the tangent plane)
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

        // A = 1 minus the occlusion by each sphere
        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);
}

Real shadows need shadow rays, which come next.

Ray Tracing

The shading models above are local: they only use light that comes directly from the light sources. Whitted ray tracing [11] adds secondary rays. Where a primary ray hits a surface, new rays are traced toward the light, in the mirror direction and through transparent surfaces, recursively.

1. Shadow Rays (Visibility Testing)

A shadow ray is traced from $\mathbf{p}$ in the direction ${\color{#FF9800}\mathbf{l}}$ of the light.

• If it hits an object before reaching the light (${\color{#795548}t_{\text{hit}}} < {\color{#FF9800}t_{\text{light}}}$), the point is in shadow, and only the ambient term is used.
• Rounding errors can make the ray hit the surface it starts on (“shadow acne”), so its origin is moved slightly along the normal, to $\mathbf{p} + \epsilon\mathbf{n}$.

2. Reflection Rays

For a mirror-like surface the reflected direction is computed as in the Phong model, with the view vector in place of the light vector:

Here ${\color{#2196F3}\mathbf{v}} = -\mathbf{d}$ points from the surface back along the incoming ray. A new ray is traced from $\mathbf{p} + \epsilon\mathbf{n}$ in direction ${\color{#4CAF50}\mathbf{r}}$, and the color it returns is blended with the local shading according to the surface’s reflectivity.

3. Refraction Rays (Transmission)

Light entering a transparent material such as glass changes direction. By Snell’s law, the angle of incidence $\theta_L$ and the angle of transmission $\theta_T$ satisfy the following, where $\eta_L$ is the refractive index of the medium the light leaves and $\eta_T$ that of the medium it enters (about $1.0$ for air and $1.5$ for glass) [2]:

The transmitted direction ${\color{#9C27B0}\mathbf{t}}$ follows from splitting it into parts along and perpendicular to $\mathbf{n}$. With ${\color{#2196F3}\mathbf{v}} = -\mathbf{d}$ as above and $\cos\theta_L = \mathbf{n} \cdot {\color{#2196F3}\mathbf{v}}$,

where, by Snell’s law,

If the value under the square root is negative, there is no transmitted ray: all the light is reflected. This is total internal reflection, and it can only happen when light goes from a denser to a less dense medium ($\eta_L > \eta_T$).

The same scene with recursive ray tracing and an added glass sphere. Shadow rays now block the direct light, the chrome sphere reflects the scene, and rays refract through the glass sphere.

#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;
}

// Same as GLSL reflect(): r = d - 2 (d.n) n, equal to 2 (n.v) n - v with v = -d
vec3 myReflect(vec3 d, vec3 n) {
    return d - 2.0 * dot(d, n) * n;
}

// t = eta (-v) + (eta cosL - cosT) n with v = -d, eta = eta_L / eta_T, and n on the side of v.
// The caller checks for total internal reflection (disc < 0) first.
vec3 myRefract(vec3 d, vec3 n, float eta) {
    float cosL = dot(-d, n);
    float disc = 1.0 - eta * eta * (1.0 - cosL * cosL);   // cosT^2
    return eta * d + (eta * cosL - sqrt(disc)) * n;
}

// ior > 0.0 marks a dielectric (glass); ior == 0.0 means opaque
struct Hit { float t; vec3 n; vec3 col; float refl; float ior; };

Hit traceScene(vec3 ro, vec3 rd) {
    Hit best = Hit(1e20, vec3(0), vec3(0), 0.0, 0.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), 0.05, 0.0); }

    // 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, 0.0); }

    // Glass sphere (ior = 1.5, like crown glass)
    t = hitSphere(ro, rd, vec3(0.1, 0.7, -2.2), 0.6);
    if (t > 0.0 && t < best.t) { best = Hit(t, normalize(ro+rd*t - vec3(0.1,0.7,-2.2)), vec3(0.95,0.98,1.0), 0.0, 1.5); }

    // Floor plane
    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), 0.1, 0.0);
    }
    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 from p + eps n toward the light: in shadow, only the ambient term remains
    Hit sh = traceScene(p + hit.n * 0.01, l);
    float shadow = (sh.t < length(lightPos - p)) ? 0.0 : 1.0;

    float diff = max(dot(hit.n, l), 0.0) * shadow;
    float spec = pow(max(dot(myReflect(-l, hit.n), v), 0.0), 48.0) * shadow;
    return hit.col * 0.08 + 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;

    // 4 bounces: enough for glass entry + exit + one more interaction each
    for (int bounce = 0; bounce < 4; 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;

        if (hit.ior > 0.0) {
            // Glass: refract, or reflect on total internal reflection.
            // cosL < 0 means the ray is inside the sphere, leaving it.
            float cosL  = dot(hit.n, -rd);
            bool inside = cosL < 0.0;
            vec3  faceN = inside ? -hit.n :  hit.n;  // always points against rd
            float eta   = inside ?  hit.ior : 1.0 / hit.ior;
            float cosL2 = abs(cosL);
            float disc  = 1.0 - eta * eta * (1.0 - cosL2 * cosL2);

            if (disc < 0.0) {
                // Total internal reflection
                ro = p + faceN * 0.001;
                rd = myReflect(rd, faceN);
            } else {
                // Refract into / out of the medium
                // t = eta (-v) + (eta cosL - cosT) n, as in the text
                ro = p - faceN * 0.001;  // nudge into the new medium
                rd = myRefract(rd, faceN, eta);
            }
            atten *= 0.97; // a little absorption; no Fresnel reflection in this simple model
        } else {
            // Opaque: local shading, blended with a reflection ray by the reflectivity
            col  += atten * (1.0 - hit.refl) * shade(hit, p, rd);
            atten *= hit.refl;
            if (atten < 0.01) break;
            ro = p + hit.n * 0.001;
            rd = myReflect(rd, hit.n);
        }
    }

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

Path Tracing

Whitted ray tracing gives sharp shadows and mirror reflections, but no diffuse interreflection, such as red light from a red wall reaching a white floor. Path tracing solves the rendering equation with Monte Carlo integration.

At each hit, a path tracer continues the path in a direction sampled at random over the hemisphere. Each path is a noisy estimate of the pixel’s value, and the average over many paths per pixel converges to the solution of the rendering equation. Film renderers such as Arnold, RenderMan and Cycles are path tracers.

Monte Carlo integration, importance sampling and BRDFs are covered in the next post: Importance Sampling and Path Guiding for Path Tracing.

The widget below renders one scene four ways. Phong and Blinn-Phong use local shading only. Whitted ray tracing adds shadows and mirror reflections. Path tracing also includes diffuse interreflection (global illumination); its image gets less noisy as samples accumulate.

Sphere Tracing

Sphere tracing [6] is the form of ray marching that uses the SDF value as the step size. Instead of solving an equation for the hit, it evaluates $f$ at the current point on the ray, advances by that distance, and repeats:

The loop stops when $f$ falls below a small $\epsilon$ (a hit) or when $t$ exceeds a maximum distance (a miss):

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

For an exact SDF or any 1-Lipschitz bound, the steps never cross the surface. If the function may overestimate the distance, for example after a deformation, each step is multiplied by a factor below 1.

In the widget below, each circle is centered on a sample point and has the SDF value as its radius. The circles get smaller as the ray approaches a 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);

    // The ray direction follows the mouse, or animates without it
    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: one circle of radius f(p) per step (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);
}