Immediate Radiosity

A real-time global illumination pipeline based on precomputed radiosity, implemented entirely on the GPU using Vulkan compute shaders.

Results

The following renders show the system running on four test scenes. Each row compares the same viewpoint rendered without global illumination, with diffuse GI only, and with the full pipeline (diffuse + specular indirect).

Sponza:

Box (Cornell-like scene with curved geometry):

Arena (cave/desert scene):

Abstract Park (multi-surface light bounce test):

Each triplet: no GI (left), diffuse GI only (center), full GI with specular indirect (right).

GI Component Decomposition

The forward shader composites the final image from diffuse and specular indirect contributions, each sampled independently from the radiance volume. Here is the decomposition for the Sponza scene:

Left: diffuse indirect (sampled from the radiance volume). Center: specular indirect (raytraced through the voxel volume). Right: full composite with PBR coefficients.

Demos

Key Performance Numbers

Introduction

Real-time global illumination remains one of the central unsolved problems in interactive rendering. Commercial engines approximate indirect light through a combination of voxel cone tracing, screen-space techniques, and radiance caches, but all of these approaches share a common cost: some form of ray traversal must occur every frame, either through hardware ray tracing, voxel cone marching, or probe resampling. That per-frame traversal cost scales with scene complexity and imposes a hard budget that competes directly with the broader rendering pipeline.

This project investigates a different formulation. Radiosity expresses global illumination as a linear system over surface patches: the radiance of each patch is the weighted sum of contributions from all mutually visible patches, where the weights (the form factors, or view factors) encode inter-patch geometry and visibility. Crucially, the form factors depend only on scene geometry, not on lighting. If the scene is static, they can be precomputed once and stored as a sparse matrix. The per-frame cost then reduces to a single sparse matrix-vector product (SpMV), a well-understood linear algebra primitive that maps efficiently onto GPU compute pipelines.

The central hypothesis is: precomputing visibility into a sparse radiosity matrix and reducing the runtime GI problem to an iterative SpMV can produce a competitive real-time indirect illumination result at lower per-frame cost than ray-traversal-based alternatives.

Theoretical Background

The Rendering Equation

The theoretical foundation of all global illumination algorithms is the rendering equation, formulated by Kajiya (1986):

$$L_o(\mathbf{x}, \omega_o) = L_e(\mathbf{x}, \omega_o) + \int_{\Omega} f_r(\mathbf{x}, \omega_i, \omega_o)\, L_i(\mathbf{x}, \omega_i)\, (\omega_i \cdot \mathbf{n})\, d\omega_i$$

where $L_o(\mathbf{x}, \omega_o)$ is the outgoing radiance at surface point $\mathbf{x}$ in direction $\omega_o$, $L_e$ is emitted radiance, $f_r$ is the bidirectional reflectance distribution function (BRDF), $L_i(\mathbf{x}, \omega_i)$ is the incident radiance arriving from direction $\omega_i$, and the integral is taken over the hemisphere $\Omega$ centered on the surface normal $\mathbf{n}$.

The fundamental difficulty stems from the recursive nature of this equation: $L_i(\mathbf{x}, \omega_i)$ is itself the outgoing radiance of whatever surface is visible from $\mathbf{x}$ along $\omega_i$, making the full solution a system of coupled integral equations over all surfaces in the scene. Direct illumination truncates this recursion at one bounce. Global illumination algorithms differ in how they approximate or solve the full recursive integral.

Radiosity Formulation

Radiosity simplifies the rendering equation by restricting to purely diffuse surfaces. Under this assumption, outgoing radiance is independent of viewing direction, and the continuous integral reduces to a finite sum over discrete patches:

$$B_i = E_i + \rho_i \sum_{j=1}^{n} F_{ij}\, B_j$$

where $B_i$ is the radiosity (total energy per unit area leaving) of patch $i$, $E_i$ is its emitted energy, $\rho_i$ its reflectivity, and $F_{ij}$ the form factor from patch $i$ to patch $j$. In matrix form:

$$\mathbf{B} = \mathbf{E} + \mathbf{M}\mathbf{B}, \quad \mathbf{M}_{ij} = \rho_i F_{ij}$$

The form factor $F_{ij}$ encodes the fraction of energy leaving patch $i$ that arrives at patch $j$, accounting for geometry, orientation, and mutual visibility. Since $\mathbf{M}$ depends only on scene geometry, it can be precomputed once for static scenes. The per-frame cost then reduces to iteratively solving the linear system, which in this implementation is done via Jacobi iteration dispatched as a GPU compute shader.

Prior Work and Context

Radiosity algorithms were first introduced by Goral et al. (1984) and refined by Cohen and Greenberg (1985) through the hemicube representation for efficient form factor computation, and further by Cohen et al. (1988) through progressive refinement. These contributions established the backbone for subsequent industry systems. However, the $O(n^2)$ complexity in patch count and the large memory requirements for the full form factor matrix made practical adoption limited. Stochastic ray-tracing lightmappers proved more popular due to their flexibility and better scalability.

Geomerics Enlighten (later acquired by ARM and integrated into Unity and Unreal Engine) is the closest prior work. It also decouples expensive visibility computation from per-frame energy propagation, reducing runtime cost to a sparse matrix-vector product. However, Enlighten operates at a coarse patch resolution dictated by lightmap texel density and encodes only diffuse-to-diffuse (D→D) energy exchange. Specular indirect contributions are approximated separately through reflection captures rather than being derived from the transport matrix itself.

Voxel Cone Tracing (VCT), introduced by Crassin et al. (2011), represents the principal alternative for dynamic scenes. Commercial adaptations include NVIDIA VXGI (integrated into Unreal Engine 4) and Unreal Engine 5's Lumen, a hybrid software raymarcher against signed distance fields combined with hardware ray tracing. Unity's SDFGI and earlier Voxel GI baked probe systems represent a related family. The principal limitations of pure VCT are: voxel resolution imposes a lower bound on spatial frequency, causing light leaking through thin geometry; and the cone tracing step is bandwidth-bound. Sparse voxel representations and clipmap hierarchies partially mitigate the memory cost but do not eliminate the per-frame traversal cost.

This work differs from both families: unlike Enlighten, it operates at voxel resolution rather than patch resolution and targets real-time use without a separate offline bake step; unlike VCT, it eliminates per-frame ray traversal entirely, replacing it with a precomputed transport matrix.

Algorithm Overview

The pipeline is divided into two groups of passes:

Init passes (executed once at scene load, or whenever geometry changes):

Render passes (executed every frame, or whenever lighting changes):

The relight passes need only execute when the lighting configuration changes. If lighting remains static between frames, the radiance result can be cached, allowing the GPU to idle during steady-state rendering. This has a direct energy efficiency implication for applications with largely static lighting.

Initialization Pipeline

Voxel Image Mapping

The first step establishes the voxel resolution for the 3D image volumes. Given a chosen resolution parameter $N$, the scene's Axis-Aligned Bounding Box (AABB) is computed. For each axis $a \in \{x, y, z\}$, the extent is:

$$\Delta a = \max_i V_a[i] - \min_i V_a[i]$$

Rather than allocating a cubic $N \times N \times N$ grid (which wastes voxels on shorter axes), the grid is allocated proportionally:

$$N_a = \left\lfloor N \cdot \frac{\Delta a}{\max(\Delta x, \Delta y, \Delta z)} \right\rfloor$$

This guarantees every allocated voxel maps to a potentially occupied region of the scene. All three axes share the same isotropic voxel edge length:

$$\delta = \frac{\max(\Delta x, \Delta y, \Delta z)}{N}$$

The table below shows the resulting grid dimensions and voxel edge lengths for each test scene at three resolution levels.

The transformation from world coordinates to integer voxel coordinates is:

$$\mathbf{i} = \left\lfloor \frac{\mathbf{p} - \mathbf{o}}{\delta} \right\rfloor$$

where $\mathbf{o}$ is the scene AABB origin and $\mathbf{p}$ is the world-space position.

Voxelization

The voxelization pass converts continuous scene geometry into a discrete 3D volume through a single GPU draw call with a three-stage pipeline:

1. Vertex shader: transforms geometry to world space via the per-instance model matrix.
2. Geometry shader: selects the dominant projection axis per triangle to maximize rasterized footprint.
3. Fragment shader: writes albedo and encoded normal into the 3D image via imageStore.

Dominant-axis projection. A triangle with geometric normal $\mathbf{N} = (\mathbf{p}_1 - \mathbf{p}_0) \times (\mathbf{p}_2 - \mathbf{p}_0)$ is projected along the axis that maximizes its projected area:

$$k^* = \arg\max_{k \in \{x,y,z\}} |N_k|$$

Projecting along any other axis foreshortens the triangle, potentially reducing its rasterized footprint below one pixel and causing voxels to be missed.

Voxel write. For each rasterized fragment, the shader recovers the integer voxel index, samples the albedo from the material texture, compresses the surface normal into a single uint8 value $\eta$ (see Normal Encoding below), and writes:

$$\texttt{voxelAlbedo}[\mathbf{i}] \leftarrow (c_r, c_g, c_b, \eta)$$

Since imageStore is not guaranteed atomic for rgba8 images and multiple triangles may project to the same voxel, the write order is undefined. In practice this is acceptable because all triangles sharing a voxel are spatially adjacent and carry similar material properties.

Voxel Maps Visualization

The voxelization produces three volumes: the albedo map, the normal map, and (after SDF computation) the SDF map.

Box scene:

Sponza scene:

Albedo (left), normal direction encoded as color (center), and SDF distance field (right, brighter = farther from geometry).

Normal Encoding Strategies

Surface normals must be stored in the alpha channel of the R8G8B8A8_UNORM albedo volume. A zero alpha value serves as the voxel presence bit: a voxel is solid if and only if $\alpha \neq 0$. The remaining 255 non-zero values encode the surface normal direction. Three strategies were evaluated.

6-Direction (Face-Aligned)

The normal is quantized to the nearest axis-aligned face direction, yielding six possible values:

$$\hat{\mathbf{n}} \in \{(\pm1,0,0),\; (0,\pm1,0),\; (0,0,\pm1)\}$$

Encoding and decoding require no arithmetic beyond an index lookup. The quantization error is large for oblique normals.

26-Direction (Face, Edge, and Corner)

The normal is quantized to the nearest vertex of the $3\times3\times3$ integer lattice with the origin excluded, giving 26 candidate directions. Each is subsequently normalized:

$$\hat{\mathbf{n}} \in \{a,b,c\}^3 \setminus \{\mathbf{0}\},\quad a,b,c \in \{-1,0,1\}$$

This scheme provides better angular coverage than 6-direction at no additional decode cost, but still produces visible faceting on curved surfaces.

Octahedral Encoding

Octahedral encoding establishes a bijective mapping between $S^2$ and the unit square, distributing all 255 available levels uniformly over the sphere. The scheme projects a unit normal $\mathbf{n} = (n_x, n_y, n_z)$ onto the $\ell^1$ unit sphere:

$$\mathbf{p} = \frac{(n_x, n_y)}{|n_x| + |n_y| + |n_z|} \in [-1,1]^2$$

For the lower hemisphere ($n_z < 0$), the projection is folded into the unit square:

$$\mathbf{p} \leftarrow \bigl((1 - |p_y|)\,\text{sign}(p_x),\; (1 - |p_x|)\,\text{sign}(p_y)\bigr)$$

The two components are then quantized to 4-bit unsigned integers and packed into a single byte:

$$q_x = \text{round}\!\left(\frac{p_x+1}{2} \cdot 15\right), \quad q_y = \text{round}\!\left(\frac{p_y+1}{2} \cdot 15\right)$$

$$e = (q_x \ll 4) \mid q_y, \quad e \in [1, 255]$$

This makes full use of all 255 levels and distributes quantization error uniformly over $S^2$, at the cost of a small number of additional ALU instructions at decode time. Below is the GLSL implementation used in the project:

vec2 encodeOctahedral(vec3 n) {
    n = normalize(n);
    n /= (abs(n.x) + abs(n.y) + abs(n.z));
    vec2 oct = n.z >= 0.0 ? n.xy : (vec2(1.0) - abs(n.yx)) * sign(n.xy);
    return oct * 0.5 + 0.5;
}

float packNormal8(vec3 n) {
    vec2 oct = encodeOctahedral(n);
    uvec2 q = uvec2(clamp(oct, 0.0, 1.0) * 15.0 + 0.5);
    return float((q.x << 4) | q.y) / 255.0;
}

vec3 decodeOctahedral(vec2 oct) {
    oct = oct * 2.0 - 1.0;
    vec3 n = vec3(oct, 1.0 - abs(oct.x) - abs(oct.y));
    float t = max(-n.z, 0.0);
    n.xy += vec2(n.x >= 0.0 ? -t : t, n.y >= 0.0 ? -t : t);
    return normalize(n);
}

vec3 unpackNormal8(float pack) {
    uint bits = uint(pack * 255.0 + 0.5);
    vec2 q = vec2((bits >> 4) & 0xFu, bits & 0xFu) / 15.0;
    return decodeOctahedral(q);
}

Interactive Comparison

The following demo renders a sphere with normals encoded and decoded through each strategy. The color represents the decoded normal direction (RGB = XYZ mapped to [0,1]). The 6-direction mode shows extreme banding, 26-direction shows visible faceting, and octahedral produces a smooth gradient.

Experimental Comparison

The following renders show the effect on diffuse GI (Box scene, curved geometry):

Diffuse GI with 6-direction normals (left), 26-direction (center), and octahedral (right). Note banding on the curved sphere with coarser encodings.

Experimental comparison confirmed that octahedral encoding provides meaningfully better results for surfaces with high normal variation. For flat geometry the difference is minimal, but octahedral is recommended as the default since the decode overhead is negligible. There is no additional cost during matrix construction since both encoding and decoding are GPU-side operations.

Signed Distance Field Computation

Visibility queries during form factor computation require traversing the voxel volume. Raw voxel traversal (DDA) is correct but slow for rays crossing large empty regions. To accelerate this, a signed distance field is precomputed from the voxel volume, enabling rays to skip empty space in large steps during raymarching.

Jump Flood Algorithm

The SDF is computed via the Jump Flood Algorithm (JFA), a standard $O(\log N)$ parallel algorithm that produces a Voronoi diagram (and therefore a distance field) on the GPU. JFA runs in a series of double-buffered passes over the volume; the number of passes is:

$$P = \lfloor \log_2(\max(N_x, N_y, N_z)) \rfloor + 1$$

In pass $k$ (counting from $k=0$), every voxel examines its neighbors at offset $2^{P-k-1}$ in each cardinal direction. If a neighbor holds a seed (solid voxel position) closer than the current best, the current voxel adopts that seed. After all $P$ passes, each voxel knows its nearest solid voxel and thus its unsigned distance to the nearest surface.

The key property is that JFA converges to the correct Voronoi diagram in logarithmic passes rather than the linear passes required by a naive flood fill. The resulting SDF is an approximation (not necessarily tight), but serves as a valid lower bound on the distance to the nearest occupied voxel along any ray, which is sufficient to preserve raymarching correctness at the cost of reduced skipping efficiency on strongly non-cubic AABBs.

Interactive JFA Demo

The following demo implements JFA on a 2D grid. Click to place seeds, then step through the algorithm to observe how the distance field propagates in logarithmic passes. Each cell is colored by its distance to the nearest seed.

Reference Implementation

The core evaluation logic of JFA in the 2D reference plotter:

evaluate(seed, neighbor, k) {
    if (neighbor[0] >= this.N || neighbor[0] < 0) return;
    if (neighbor[1] >= this.N || neighbor[1] < 0) return;

    var p = this.values[seed[0]][seed[1]];
    var q = this.values[neighbor[0]][neighbor[1]];

    if (p == null && q != null) {
        this.values[seed[0]][seed[1]] = this.values[neighbor[0]][neighbor[1]];
        this.values[seed[0]][seed[1]].distance = this.evaluateDistance(seed);
        return;
    }
    if (p != null && q != null) {
        var dq = this.evaluateDistance(neighbor);
        var dp = this.evaluateDistance(seed);
        if (dq < dp) {
            this.values[seed[0]][seed[1]] = this.values[neighbor[0]][neighbor[1]];
            this.values[seed[0]][seed[1]].distance = dq;
        }
    }
}

iterate() {
    for (var i = 0; i < this.N; i++) {
        for (var j = 0; j < this.N; j++) {
            var seed = [i, j];
            this.evaluate(seed, [i,     j + this.k], this.k);
            this.evaluate(seed, [i - this.k, j    ], this.k);
            this.evaluate(seed, [i + this.k, j    ], this.k);
            this.evaluate(seed, [i,     j - this.k], this.k);
        }
    }
    this.k = this.k / 2;
}

The final SDF volume is a single-channel 3D texture used exclusively to accelerate ray marching during form factor construction.

Solid Voxel List

Before constructing the form factor matrix, all solid voxels are compacted into a linear buffer called the SolidVoxelList. This serves two purposes:

1. Matrix construction iterates only over solid voxels, skipping empty space entirely.
2. The buffer establishes a stable mapping between entries of the radiance vector (indexed during Jacobi iteration) and voxel positions in the 3D radiance texture (indexed during scatter and forward shading).

The buffer layout (std430 alignment):

• Header (16 bytes): atomic counter for total solid voxel count.
• Entry $i$ (32 bytes each): voxel coordinates $(x, y, z)$ as uint, encoded normal as uint, albedo as vec3, padding.

A compute shader iterates over every voxel in the volume, tests the alpha-channel presence bit, and atomically appends solid voxels to the list. A companion ReverseMap (a 3D texture with R32_UINT format) is simultaneously populated, mapping each integer voxel position back to its index in the SolidVoxelList. This inverse map is used in the scatter pass to write radiance at the correct list position.

The solid voxel counts for the test scenes at $N=256$ are:

Form Factor Matrix Construction

The transport matrix $\mathbf{A} = \rho \cdot F$ is constructed stochastically on the GPU. For each solid voxel $i$, the shader fires $R$ cosine-weighted rays sampled from a jittered Hammersley sequence over the hemisphere defined by voxel $i$'s surface normal.

Hammersley Low-Discrepancy Sampling

The Hammersley sequence is a low-discrepancy point set constructed by pairing the Van der Corput radical-inverse function in base 2 with uniform stratification over $R$ samples. For hemisphere ray generation, each pair $(u_1, u_2) \in [0,1)^2$ is mapped to a cosine-weighted direction via:

$$\theta = \arccos\sqrt{1-u_1}, \quad \phi = 2\pi u_2$$

This provides superior hemisphere coverage compared to pseudo-random sampling, which is important since the matrix construction budget is fixed at $R$ rays per voxel.

Weight Computation and Storage

Each ray that hits a visible solid voxel $j$ (with incidence $\hat{\mathbf{d}} \cdot \hat{\mathbf{n}}_j < 0$) contributes an entry to row $i$ of the sparse matrix with weight:

$$w_{ij} = \frac{1}{R}$$

This is the unbiased Monte Carlo estimate of the geometric form factor $F_{ij}$ under cosine-weighted sampling. Visibility is evaluated by either DDA traversal on the integer voxel grid or SDF raymarching (see below).

The matrix is stored in a fixed-size CSR-like format: a flat buffer of (column-index, weight) pairs for each row, with a maximum of $R$ entries per row. The fixed-size layout avoids dynamic allocation on the GPU at the cost of wasted memory for voxels with fewer than $R$ visible neighbors.

DDA vs. SDF Raymarching

Two traversal strategies are available for the ray-voxel intersection during matrix construction:

DDA (Amanatides-Woo voxel traversal): Steps through the voxel grid one face at a time. Guaranteed to intersect every solid voxel the ray passes through. Reference-correct but traverses every empty voxel individually.

SDF raymarching: Uses the precomputed distance field to skip ahead by the minimum guaranteed distance to the nearest surface, approximately doubling traversal speed on sparse scenes. The SDF approximation from JFA is not tight, so the marcher may occasionally overshoot and miss a near-surface interaction. Experimentally, this produces slightly lower-intensity indirect shadows in edge cases, but no visible quality difference in typical scenes.

Matrix construction time across scenes. Blue: SDF raymarching. Red: DDA traversal. SDF is approximately 2× faster on average.

Relaxation Parameter

Multiple rays can land on the same destination voxel when source and destination are spatially close relative to voxel size. This produces duplicate column entries in a row, causing row weights to sum above 1.0, which violates energy conservation and prevents Jacobi convergence. Sorting rows and deduplicating would be expensive on the GPU. Instead, a relaxation parameter $\lambda \in (0, 1)$ scales the total propagation:

$$\mathbf{x}^{(k+1)} = \lambda\left(\mathbf{E} + \mathbf{A}\,\mathbf{x}^{(k)}\right)$$

Experimental results across all four test scenes showed optimal visual agreement with Cycles reference renders for $\lambda \in [0.7, 0.9]$. The parameter is scene-dependent and tunable at runtime.

Low-Probability High-Energy Paths

A known limitation of the stochastic sampling approach: rays that hit a very small but extremely energetic region may be missed if the row budget $R$ is exhausted before they are sampled. Because neighboring voxels may have different row entries, this can produce localized noise from matrix discontinuities. Increasing $R$ reduces this noise at linear cost to both memory and construction time.

Per-Frame Pipeline

Light Gathering

A compute kernel iterates over every solid voxel in the SolidVoxelList and gathers current direct lighting into an EmissionBuffer. The computation is identical to a standard forward shader for the supported light types:

• Spotlights / Sunlights: Shadow map visibility test per voxel, using the same PCF shadow map generated for the forward pass.
• Emissive geometry: Read directly from the emissive channel of the albedo volume, populated during voxelization.
• Omnilights: Cubemap shadow mapping (not implemented in this version due to time constraints).

This pass is the most extensible part of the algorithm. Any direct lighting source that can be evaluated at a 3D world position can be incorporated into the emission buffer without modifying any other pipeline stage.

Radiosity Solver (Jacobi SpMV)

The core solver performs one step of Jacobi power iteration per frame:

$$\mathbf{x}^{(k+1)} = \lambda\left(\mathbf{E} + \mathbf{A}\,\mathbf{x}^{(k)}\right)$$

where $\mathbf{x}$ is the per-voxel radiance vector (one vec4 per solid voxel), $\mathbf{E}$ is the current emission vector, and $\mathbf{A} = \rho \cdot F$ is the albedo-weighted form factor matrix. The relight converges to a steady state in several iterations, after which it can be cached until lighting changes.

Double buffering. Jacobi iteration reads the previous state while writing the next, requiring a ping-pong pair of radiance buffers. The GPU-equivalent of Gauss-Seidel (which updates in-place) would converge faster but requires reading and writing the same buffer with ordering guarantees, which is not efficiently achievable in a compute shader.

Workgroup SpMV. Each workgroup is responsible for computing a single dot product (one row of the matrix-vector product). The workgroup size is fixed to 64 threads, matching the RDNA3 wavefront width.

For a row of size $R$, each thread $t$ is responsible for entries at indices $\{t,\; t + 64,\; t + 128,\; \ldots\}$ within the row. Thread $t$ reads its assigned entries, accumulates partial results into shared memory, and then participates in a barrier-synchronized tree reduction:

1. Each thread accumulates its partial sum into sharedMem[localID].
2. A sequence of barrier-synchronized reduction steps halves the active thread count each time.
3. Thread 0 writes the final scalar result to the output radiance buffer at the correct solid-voxel index.

This pattern is efficient because all threads in a workgroup read contiguous ranges of the row, maintaining good cache locality on the SSBO.

Light Scattering

After the Jacobi iteration, the updated radiance lives in the flat SSBO RadianceBuffer. For efficient sampling in the forward shader, it must be written into a 3D texture (RadianceTexture, RGBA16F) indexed by voxel position. A compute shader dispatches one invocation per solid voxel and performs this scatter using the voxel coordinates stored in the SolidVoxelList.

The RadianceTexture is cleared at the start of each scatter pass to remove stale radiance from the previous iteration.

Shadow hard-edge correction. The RadianceBuffer contains embedded direct light gathered in the EmissionBuffer. If written directly to the RadianceTexture and sampled in the forward shader, direct light would be applied twice: once from the voxel volume sample and once from the forward shader's own direct lighting pass. To correct this, the scatter pass subtracts the emission before writing:

$$\mathbf{L}_{\text{vox}}(\mathbf{x}) = \mathbf{R}(\mathbf{x}) - \lambda\, \mathbf{E}(\mathbf{x})$$

where $\mathbf{L}_{\text{vox}}$ is the value written to the texture, $\mathbf{R}$ is the accumulated radiosity result, and $\mathbf{E}$ is the direct emission. A side effect of skipping this subtraction is that direct light values remain embedded in the RadianceTexture, allowing them to appear in specular reflections computed by the raytracer in the forward pass. This creates a design trade-off between shadow boundary accuracy and reflection quality.

Pre-multiplied alpha sampling. When the forward shader samples the RadianceTexture using a fragment's world position, the sample may land on a neighboring empty voxel (black in the texture). Linear filtering with pre-multiplied alpha corrects this: the alpha channel stores a presence value (1.0 for solid, 0.0 for empty), and the sampled color is divided by the interpolated alpha:

$$\mathbf{c}_{\text{GI}} = \frac{\mathbf{L}_{RGB}}{\alpha}$$

An additional blur kernel over the 3D sample further smooths radiance and suppresses noise from low-probability high-energy paths. Since radiosity radiance has inherently low spatial frequency, this introduces minimal error while meaningfully reducing visible artifacts.

Specular Extensions

Specular Light Tracing

Since the voxelization stage produces all data necessary for DDA ray traversal, a voxel-space ray tracer in the forward shader can evaluate specular indirect paths at shading time.

For each fragment, rays are stochastically generated within the specular cone defined by the fragment's surface normal, jittered by an amount proportional to the PBR roughness parameter (wider cone = rougher surface). Each ray is traced using DDA through the voxel volume until it hits a solid voxel or escapes the volume. The radiance at the hit voxel is read from the RadianceTexture. Results from all rays are averaged and weighted by PBR coefficients (Cook-Torrance BRDF with GGX NDF) to produce the specular indirect contribution.

The final GI value per fragment is the sum of:

• Diffuse indirect: sampled directly from the RadianceTexture at the fragment's world position.
• Specular indirect: averaged over stochastically traced rays in the specular cone.

A practical note: this approach can produce significant specular aliasing for surfaces with high-frequency normal variation. Using geometric normals instead of interpolated normals for the specular cone direction reduces this, at the cost of some shading smoothness.

SD* — Specular Writeback

As a novel contribution, the specular value computed during forward shading can be written back into the emission volume via imageStore, allowing specular contributions to re-enter the radiosity solver in subsequent iterations. This creates a specular-to-diffuse (SD*) transport path without requiring per-frame ray traversal in the GI solver itself.

The primary challenge is memory contention: many fragment shader threads may attempt to write to the same voxel simultaneously. A proximity mask limits writes to at most one representative fragment per voxel cell.

Let $\mathbf{p}_w$ be the fragment world position, $\mathbf{c}$ the camera position, and $\hat{n}$ the surface normal.

Step 1 — Camera distance weighting. A sub-linear distance factor:

$$d = \|\mathbf{c} - \mathbf{p}_w\|^{\alpha}, \quad \alpha < 1$$

softens the falloff with camera distance.

Step 2 — Voxel-space projection. The fragment is projected into voxel space $\mathbf{v} = \text{world2volume}(\mathbf{p}_w)$, and its fractional cell position $\mathbf{g} = \text{fract}(\mathbf{v}) \in [0,1)^3$ gives its position within the voxel cell.

Step 3 — Normal-weighted centering. Each axis of $\mathbf{g}$ is blended toward the voxel center (0.5) proportionally to the surface normal alignment with that axis:

$$g_i' = \text{mix}(g_i,\ 0.5,\ |\hat{n}_i|), \quad i \in \{x, y, z\}$$

A face perpendicular to axis $i$ ($|\hat{n}_i| \approx 1$) collapses all its fragments to the same voxel-center coordinate along that axis.

Step 4 — Proximity score. Distance from the blended grid point to the voxel center:

$$s = \|\mathbf{g}' - 0.5\|$$

Normalized by camera distance and clamped:

$$\hat{s} = \text{clamp}\!\left(1 - \frac{k \cdot s}{d},\ 0,\ 1\right)$$

Step 5 — Mask condition. A fragment issues the imageStore only if $\hat{s} \geq \tau$ for threshold $\tau \in (0,1)$. This heuristically selects at most one representative fragment per voxel cell under normal viewing conditions.

Experimentally, SD* adds perceptual realism through IBL-like indirect specular contributions but introduces significant temporal flickering and carries a performance cost several times higher than regular shading. A more robust implementation would write specular data to a color attachment and process writes in a dedicated compute pass, avoiding fragment shader contention entirely. This remains as future work.

Memory Analysis

Let $V = N_x N_y N_z$ be the total voxel count, $s$ the number of solid voxels, $M_s = 2^{\lceil\log_2 s\rceil}$ the GPU solid-list capacity rounded up to the next power of two (for lock-free atomic allocation), and $R$ the per-row form-factor budget.

Volume buffers (seven 3D textures over the full grid):

Three textures at 8 B/voxel (R16G16B16A16_SFLOAT): radiance ping, radiance pong, JFA work volume. Three at 4 B/voxel (R8G8B8A8_UNORM or R32_UINT): albedo, JFA seed scratch, ReverseMap:

$$\mathcal{M}_{\text{vol}}(V) = 36\,V \quad \text{[bytes]}$$

Matrix GI buffers (allocated against $M_s$):

The solid list stores 8 uint32 fields per entry (32 B) plus a 16 B header; each form-factor entry stores a column index and quantized weight (2 uint32 = 8 B); row counters use 4 B; and three vec4 arrays of 16 B each hold the Jacobi solution ping-pong pair and the emission cache:

$$\mathcal{M}_{\text{mat}}(M_s, R) = M_s(84 + 8R) + 16 \quad \text{[bytes]}$$

The dominant term is $8R \cdot M_s$ (the matrixHits buffer), accounting for approximately 97% of $\mathcal{M}_{\text{mat}}$ at the default parameters ($R = 384$).

Total:

$$\mathcal{M}_{\text{total}}(V, M_s, R) = 36V + M_s(84 + 8R) + 16 \quad \text{[bytes]}$$

Both terms scale as $O(N^3)$, but $\mathcal{M}_{\text{mat}}$ carries the coefficient $8R$, making it the dominant cost in practice. For the Sponza scene at $N = 256$ ($V = 4{,}300{,}544$, $s = 441{,}511$, $M_s = 2^{20}$, $R = 384$):

$$\mathcal{M}_{\text{total}} = 3\,287.72 \text{ MiB}$$

The matrixHits buffer alone accounts for 93% of total GPU memory consumption at these parameters. Matrix compression (e.g., block SVD, hierarchization) is the most impactful direction for reducing this footprint.

Experimentation

Test Platform

Default experimental parameters:

Voxel Resolution

$N = 128$ (left), $N = 256$ (center), $N = 384$ (right). Computation time scales cubically with $N$. Lower resolutions produce smoother results with fewer artifacts from high-frequency normals, at the cost of spatial detail and light-leak through thin geometry.

At $N = 128$, the Sponza voxel edge length is approximately 29 cm — coarser than the scene's thinnest geometry, which allows some light leaking. At $N = 256$ (14.5 cm), the main visible artifacts are reduced. At $N = 384$, diminishing returns are observed in quality while construction time grows significantly.

Relaxation Parameter

$\lambda = 0.6$ (left), $\lambda = 0.8$ (center), $\lambda = 1.0$ (right). The optimal range $\lambda \in [0.7, 0.9]$ matches Cycles reference renders most closely. Values approaching 1.0 may fail to converge if row weights sum above 1 due to duplicate entries.

DDA vs. SDF Matrix Construction

Across all four test scenes, SDF raymarching reduces matrix construction time by approximately 2× compared to DDA, with no perceptible difference in final GI quality for typical lighting configurations. In edge cases with tightly-spaced geometry, SDF can miss near-surface interactions due to the approximate JFA distance field, producing slightly lower-intensity indirect shadows.

Normal Encoding

The 6- and 26-direction schemes require no arithmetic at decode time and trivially preserve the presence bit. The octahedral scheme makes full use of all 255 levels and distributes quantization error uniformly over $S^2$ at the cost of a small number of additional ALU instructions. For surfaces with significant normal variation (curved geometry), octahedral encoding is visibly superior. For flat planar geometry, the difference is negligible. Octahedral encoding is recommended as the default.

Relight Performance

One Jacobi iteration (the entire relight step) runs in 2–4 ms on the AMD RX 7700 XT. The bottleneck is memory bandwidth rather than compute: each iteration reads the full matrixHits buffer (3 GiB at default parameters) to perform one SpMV pass. Performance scales directly with memory bandwidth. The AMD 7700 XT is particularly strong in this regard; performance on bandwidth-limited mobile or integrated GPUs would be substantially worse.

Comparison Against Reference

The diffuse indirect component was compared against Cycles (path-traced ground truth) and EEVEE (real-time) for the same lighting configuration. The main visible difference against Cycles is shadow boundary quality (PCF shadow maps vs. soft ray-traced shadows), not the GI computation itself. The indirect color bleeding and overall energy distribution match Cycles closely for all four test scenes.

Limitations

GPU memory. The view factor matrix dominates memory consumption. At $N = 256$ with $R = 384$, the Sponza scene requires over 3 GiB of VRAM for the GI subsystem alone. Scaling to $N = 512$ or larger scenes rapidly exceeds available VRAM on consumer hardware.

Static geometry only. The form factor matrix encodes visibility between all patch pairs at construction time. Any change in scene geometry (moving an object, opening a door) invalidates the matrix and requires reconstruction. Partial matrix updates for small changes are theoretically possible but not implemented.

Memory bandwidth bottleneck. The per-frame relight reads the entire matrix per iteration. At 3 GiB and 2–4 ms, this saturates the memory bus at ~750 GB/s, which is close to the theoretical peak bandwidth of the RX 7700 XT. Performance on bandwidth-limited hardware degrades proportionally.

Specular writeback (SD).* The proximity mask strategy reduces contention but does not eliminate it. The result suffers from temporal flickering due to frame-to-frame variation in which fragments are selected. The correct solution (writing to a color attachment and processing in compute) was not implemented within the scope of this project.

Specular aliasing. The voxel-space DDA ray tracer in the forward shader produces aliasing for high-frequency normal variation. Using geometric normals for the ray direction mitigates this at the cost of shading accuracy.

Conclusions and Future Work

The experimental results support the central hypothesis: a precomputed sparse radiosity matrix enables competitive real-time global illumination with a per-frame cost of 2–4 ms on current hardware, without per-frame ray traversal. The diffuse indirect component closely matches path-traced references. The specular extension, while approximate, adds perceptual realism at low additional cost.

The primary constraint is GPU memory. At practical resolutions, the transport matrix requires several gigabytes of VRAM, which limits scalability to scenes that fit within the memory budget of high-end consumer hardware.

Future directions:

• Matrix compression: randomized block SVD, hierarchical radiosity, or grouped-patch compression to reduce the dominant $O(R \cdot M_s)$ memory term while preserving transport accuracy.
• Variable-resolution patches: sparse voxel representations with non-uniform in-world patch sizes, allowing high-detail regions near the camera and coarser patches farther away.
• Dynamic partial updates: selective matrix recomputation for moving geometry, reusing unchanged rows and reconstructing only affected regions.
• Async compute overlap: the relight pass is memory-bound while the forward pass is compute-bound; scheduling them in parallel on separate GPU engines could improve overall frame time.
• SD via compute*: replacing the fragment-shader writeback with a color attachment and a dedicated compute pass would eliminate contention and enable robust specular-to-diffuse transport.

The most practical near-term application is as an optional offline-to-realtime GI backend for renderers such as EEVEE, bridging the quality gap with Cycles at a fraction of the render time (Cycles takes minutes per frame; Immediate Radiosity runs in under 4 ms after a one-time setup cost).

Novel Contributions

• SD* (Specular Writeback): a new transport path that writes specular contributions from the forward shader back into the diffuse radiance cache, enabling indirect specular propagation through subsequent radiosity iterations without per-frame ray traversal in the solver.
• Stochastic SDF raymarching for form factor computation: using the JFA-computed approximate SDF to accelerate matrix construction, achieving approximately 2× speedup over DDA with negligible quality loss.
• Comprehensive parameter study across four scenes: voxel resolution ($N \in \{128, 256, 384\}$), matrix row size ($R \in \{128, 256, 384\}$), normal encoding strategy (6-direction, 26-direction, octahedral), traversal mode (DDA vs. SDF), relaxation parameter ($\lambda \in [0.6, 1.0]$), sampling mode (pre-multiplied alpha, blur kernel).
• Full GPU-only implementation: all passes including voxelization, JFA, matrix construction, Jacobi iteration, and scatter are compute shaders with no CPU-side iteration. The only CPU involvement is dispatch call submission.

Technologies

• C++ / Vulkan / GLSL
• GPU compute shaders (SSBO-based SpMV, Jacobi power iteration)
• Signed distance fields via Jump Flood Algorithm
• Voxelization pipeline with geometry-shader dominant-axis projection
• Physically-based rendering (Cook-Torrance microfacet BRDF, GGX NDF, Disney parameterization)
• Low-discrepancy sampling (Hammersley / Van der Corput sequence)
• Sparse matrix storage (fixed-row CSR-like format)
• Pre-multiplied alpha volume sampling