Ray tracing volume densities

James T. KajiyaBrian P. Von Herzen

article1984SIGGRAPH1,505 citations

Presents the foundational ray tracing algorithms and radiative transfer solutions for rendering participating media such as clouds, fog, and flames defined within 3D volume grids.

Listen

Synthesizing realistic imagery of participating media and natural atmospheric phenomena—such as clouds, fog, flames, and dust—presents a major challenge in computer graphics. Earlier rendering techniques relied on restrictive plane-parallel models, low-reflectance assumptions, or fixed viewing and lighting geometries that fail when generating general, dynamic scenes.

The article demonstrates a generalized ray tracing framework capable of rendering arbitrary three-dimensional volume densities and introduces a mathematical approximation to handle multiple light scattering in highly reflective media like clouds.

The approach models volume densities on uniform 3D grids and splits the rendering workflow into two distinct stages: a precomputation phase that maps illumination through the volume, and an execution phase that integrates ray brightness and optical depth during ray tracing. To overcome the computational intractability of full radiative transfer, the authors develop a perturbation method expanded in spherical harmonics and truncate it to low orders. Additionally, they simulate cloud density evolution using a physical atmospheric convection model governed by partial differential equations.

The article establishes several key findings. First, single-scattering approximations cause severe visual defects, specifically artificial darkening on the shadowed sides of clouds, proving that realistic rendering of high-reflectance media requires accounting for multiple scattering. Second, separating the illumination precomputation from ray traversal allows unrestricted camera angles, internal viewing, and scene interactions such as cloud-induced shadows and reflections in other objects. Third, the spherical harmonic perturbation approach decouples directional scattering effectively, enabling numerical relaxation for high-reflectance scenarios. Finally, physically simulated convection grids (tested on 10x10x20 meshes) successfully generate dynamically realistic cloud formations across multi-minute lifespans, with render times ranging from 1 to 4 hours per frame on an IBM 4341.

These findings mean that computer-generated imagery can represent complex natural phenomena and particle systems without restrictive geometry, enhancing realism in visual simulations and animation. However, the computational cost remains high, requiring significant processing time per frame, and the current framework has a limitation where procedural objects do not cast reciprocal shadows back into the volume densities.

Teams implementing volumetric rendering should adopt precalculated illumination grids for ray tracing efficiency and apply low-order spherical harmonic approximations for bright scattering media. Future development should focus on optimizing grid computation, supporting bidirectional shadowing between solids and volumes, and evaluating higher-resolution physical simulations.

California Institute of Technology grid node numerical solution ray tracing ray tracing volume densities single scattering

Kajiya et al (1984).pdf

Paper: Particle systems—a technique for modeling a class of fuzzy objects, W. Reeves (1983). Reeves introduces particle systems to represent dynamic, non-rigid fuzzy phenomena such as fire and clouds, establishing the foundational representation that volume ray tracing models as continuous density fields.
Paper: Models of light reflection for computer synthesized pictures, J. Blinn (1977). Blinn establishes physics-based optical scattering and reflection models that underpin the radiative transfer and scattering approximations developed for participating media.
Paper: A reflectance model for computer graphics, Robert L. Cook et al. (1981). Cook and Torrance introduce a comprehensive microfacet-based physical reflection framework that provides the foundational optics principles for rendering light interaction in graphics.

Paper: The rendering equation, James T. Kajiya (1986). Kajiya formalizes the general rendering equation that generalizes early surface and volume scattering formulations into a unified framework solved via Monte Carlo path tracing.
Paper: Marching cubes: A high resolution 3D surface construction algorithm, William E. Lorensen et al. (1987). Lorensen and Cline introduce the marching cubes algorithm to extract explicit isosurfaces directly from the 3D volumetric density grids discussed in volume rendering.
Paper: Stable fluids, Jos Stam (1999). Stam develops an unconditionally stable numerical framework to simulate dynamic fluid and smoke advection through 3D density grids, substantially advancing dynamic volume modeling.
Paper: 3D Gaussian Splatting for Real-Time Radiance Field Rendering, Bernhard Kerbl et al. (2023). Kerbl et al. modernize volume and radiance field rendering by leveraging 3D Gaussian primitives to enable real-time novel view synthesis from volumetric representations.
Paper: TensoRF: Tensorial Radiance Fields, Anpei Chen et al. (2022). Chen et al. extend volumetric radiance and density fields through tensorial grid decompositions to dramatically accelerate volume reconstruction and rendering.

Abstract

This paper presents new algorithms to trace objects represented by densities within a volume grid, e.g. clouds, fog, flames, dust, particle systems. We develop the light scattering equations, discuss previous methods of solution, and present a new approximate solution to the full three-dimensional radiative scattering problem suitable for use in computer graphics. Additionally we review dynamical models for clouds used to make an animated movie.

Topics

grid node numerical solution ray tracing scattering problems 3d graphics

Keywords

ray tracing volume densities single scattering optical densities spherical harmonics Legendre polynomials plane parallel atmospheres albedo approximation particle systems image synthesis

§1 Introduction
§2 The scattering equation
§3 Solving the scattering equation
3.1 Blinn's Low Albedo approximation
3.2 A ray tracing algorithm for the low albedo case
§4 High Albedo approximation.
4.1 A Perturbation solution, conservative systems
4.2 The Scattering equation expressed in Spherical Harmonics
4.3 Matrix elements for the position
4.4 Matrix elements for the phase integral
§5 Generating density models
5.1 A Cloud Model for Generating Density Functions
§6 Computer Results
§7 Summary
§8 References

Knowls

Knowl 1 — Two-Pass Ray Tracing Algorithm for Volume Densities
algorithm
The two-pass volume rendering algorithm computes the single-scattering radiance of volume densities on a 3D grid, allowing arbitrary viewing angles, internal/external light sources, and intersections with solid procedural geometry.
Input: Density grid $\rho(x,y,z)$ , extinction coefficient $\tau$ , phase function $p(\cos\Theta)$ , light sources $\{L_i\}$ , camera ray parameters, solid scene database
Output: Pixel radiance $B$
// Pass 1: Precompute light illumination array (performed at most once per frame)
for each light source $i$ do
    for each grid element $(x,y,z)$ in volume extent do
        Compute ray path $\Gamma_{x,y,z}(t) = (x(t), y(t), z(t))$ from light $i$ to $(x,y,z)$
         $I_i(x,y,z) \leftarrow \exp\left(-\tau \int_{\Gamma_{x,y,z}} \rho(\gamma)\,d\gamma\right)$
    end for
end for
// Pass 2: Ray trace eye rays
for each viewing ray $R(t) = (x(t), y(t), z(t))$ from the eye do
    Intersect $R(t)$ with volume bounding box to get entry distance $d_1$ and exit distance $d_2$
    Intersect $R(t)$ with solid scene geometry to get nearest solid intersection distance $d_{\text{global}}$
     $\lambda_1 \leftarrow \max(0, d_1)$
     $\lambda_2 \leftarrow \min(d_{\text{global}}, d_2)$
    if $\lambda_1 < \lambda_2$ then
        Evaluate the brightness integral over $[\lambda_1, \lambda_2]$ using Romberg numerical integration:
         $B \leftarrow \int_{\lambda_1}^{\lambda_2} \exp\left(-\tau \int_{\lambda_1}^t \rho(x(\mu), y(\mu), z(\mu))\,d\mu\right) \cdot \left[\sum_i I_i(x(t), y(t), z(t)) p(\cos\Theta_i)\right] \cdot \rho(x(t), y(t), z(t))\,dt$
    else
         $B \leftarrow 0$
    end if
    Composite $B$ with light reflected from solid surface at $d_{\text{global}}$
end for
The inner line integral $\int_{\lambda_1}^t \rho\,d\mu$ computes the optical depth between the viewer and the scattering element, while $I_i(x(t), y(t), z(t))$ provides the precomputed incident light attenuation from source $i$ to that point. The phase angle $\Theta_i$ is the angle between the incoming light direction and the ray direction towards the eye.
Knowl 2 — Radiative Transfer Scattering Equation for Volume Densities
equation
The radiative transfer equation governing the intensity field $I(x, s)$ at position $x \in \mathbb{R}^3$ along direction unit vector $s$ (where $\|s\| = 1$ ) in a participating medium with local density $\rho(x)$ , optical depth per unit density $\kappa$ , and no intrinsic thermal emission is:

$-\frac{1}{\kappa \rho(x)} s \cdot \nabla_x I(x, s) - I(x, s) + \frac{1}{4\pi} \int_{\|\tilde{s}\|=1} p(s, \tilde{s}) I(x, \tilde{s}) \, d\tilde{s} = 0$

where:
- $s \cdot \nabla_x I(x, s) = \frac{dI}{ds}$ is the directional derivative of intensity along the direction $s$ .
- $\kappa \rho(x)$ is the local extinction coefficient.
- $p(s, \tilde{s})$ is the scattering phase function defining the proportion of light scattered from incident direction $\tilde{s}$ into direction $s$ .
- For isotropic media, $p(s, \tilde{s})$ depends only on the phase angle $\Theta$ between $s$ and $\tilde{s}$ (where $\cos\Theta = s \cdot \tilde{s}$ ), with perfectly diffuse scattering given by $p(\cos\Theta) = \omega_0$ and Rayleigh scattering given by $p(\cos\Theta) = \frac{3}{4}\omega_0(1 + \cos^2\Theta)$ , where $\omega_0$ is the single-scattering albedo.
Knowl 3 — Perturbation Expansion for High-Albedo Volume Radiative Transfer
model/method
For high-albedo scattering media (such as clouds where single-scattering albedo $\omega_0 \approx 1$ ), standard Neumann series expansions converge slowly. The solution is obtained via a perturbation expansion in the absorption parameter $\beta = 1 - \omega_0$ .

Normalizing the phase function by $p(\Theta) = \omega_0 \bar{p}(\Theta) = (1 - \beta)\bar{p}(\Theta)$ , the radiative transfer equation is expressed in operator form:

$L I + (1 - \beta) M I = 0$

where the linear operators are defined as:

$L I = -\frac{1}{\kappa \rho(x)} s \cdot \nabla_x I(x, s) - I(x, s)$

$M I = \frac{1}{4\pi} \int_{\|\tilde{s}\|=1} \bar{p}(s, \tilde{s}) I(x, \tilde{s}) \, d\tilde{s}$

Expanding the radiation intensity into a power series in $\beta$ ,

$I = \sum_{k=0}^\infty \beta^k I_k$

and matching equal powers of $\beta$ yields a recursive sequence of forced conservative ( $\omega_0 = 1$ ) scattering equations:

$L I_k + M I_k = -M I_{k-1}$

where $I_{-1} = 0$ .
Knowl 4 — Spherical Harmonics Representation and P-Wave Truncation for Radiative Transfer
model/method
The angular dependence of the intensity field $I(x, s)$ is expanded in normalized spherical harmonics $Y_{lm}(s) = P_{lm}(\cos\theta) e^{im\phi}$ :

$I(x, s) = \sum_{l=0}^\infty \sum_{m=-l}^l I^{lm}(x) Y_{lm}(s)$

where $P_{lm}$ are the associated Legendre polynomials of degree $l$ and order $m$ . Projecting the radiative transfer equation onto each spherical harmonic basis function using the inner product $\langle X | O | Y \rangle = \int_0^{2\pi} \int_0^\pi X^*(s) O Y(s) \sin\theta \, d\theta \, d\phi$ produces a coupled system of first-order partial differential equations for the spatial intensity coefficients $I^{lm}(x)$ :

$\sum_{l, m} \frac{1}{\kappa \rho(x)} \nabla I^{lm}(x) \cdot \langle Y_{l'm'}(s) | s | Y_{lm}(s) \rangle - I^{l'm'}(x) + \sum_{l, m} I^{lm}(x) \langle Y_{l'm'}(s) | p | Y_{lm}(s) \rangle = 0$

For computer graphics rendering, the spherical harmonic series is truncated at degree $l = 1$ ("p-wave" truncation). This retains four spatial fields ( $l=0, m=0$ and $l=1, m=-1,0,1$ ), reducing the high-dimensional angular transport problem to a small system of coupled 3D PDEs solvable across the volume grid by numerical relaxation.
Knowl 5 — Matrix Elements for Direction and Phase Operators in Spherical Harmonics
theoretical result
Projection of the radiative transport operators onto the spherical harmonic basis $Y_{lm}$ yields exact coupling matrix elements:
1. Direction Operator $s$ : Using the complex variable $u = x + iy = \sin\theta e^{i\phi}$ and $z = \cos\theta$ :
$\langle Y_{l'm'}(s) | u | Y_{lm}(s) \rangle = (k_0 \delta_{l', l+1} - k_1 \delta_{l', l-1}) \delta_{m', m+1} \cdot 2\pi$

where $k_0 = \frac{(l - m + 1)(l - m + 2)}{2l + 1}$ and $k_1 = \frac{(l + m - 1)(l + m)}{2l + 1}$ . For the $z$ -component:

$\langle Y_{l'm'}(s) | z | Y_{lm}(s) \rangle = (k_2 \delta_{l', l+1} - k_3 \delta_{l', l-1}) \delta_{m' m} \cdot 2\pi$

where $k_2 = \frac{l + m}{2l + 1}$ and $k_3 = \frac{l - m + 1}{2l + 1}$ . Consequently, the spatial gradient operator couples only adjacent spherical harmonic degrees $l \pm 1$ .
1. Phase Function Operator $p$ : Expanding an isotropic phase function $p(\cos\Theta)$ in Legendre polynomials $p(\cos\Theta) = \sum_{k=0}^\infty w_k P_k(\cos\Theta)$ and applying Laplace's addition theorem yields:
$\langle Y_{l'm'}(s) | p | Y_{lm}(s) \rangle = \frac{4\pi}{2l+1} w_l \delta_{l l'} \delta_{m m'}$

Hence, the phase function scattering operator is strictly diagonal in the spherical harmonic representation, with no scattering coupling between different spherical harmonic modes.
Knowl 6 — Dynamical Convection Model for Cumulus Cloud Density Generation
model/method
A physical simulation of 3D cumulus convection produces dynamic cloud volume optical densities using nine coupled differential and algebraic equations on a uniform grid:
1. Momentum equations (with advection, friction, and buoyancy): $\frac{\partial u}{\partial t} = -V \cdot \nabla u - F_x$ $\frac{\partial v}{\partial t} = -V \cdot \nabla v - F_y$ $\frac{\partial w}{\partial t} = -V \cdot \nabla w - F_z + \theta$ where $V = (u, v, w)$ is wind velocity, $F = (F_x, F_y, F_z) = \frac{1}{t_f} V$ is the friction vector with friction timescale $t_f$ , and $\theta$ is potential temperature providing upward buoyancy in the vertical velocity $w$ .
2. Potential temperature conservation (with advection and condensation heating): $\frac{\partial \theta}{\partial t} = -V \cdot \nabla \theta + \frac{L_{vw}}{c_p} \frac{\partial q_l}{\partial t} + Q$ where $L_{vw}$ is the latent heat of vaporization of water, $c_p$ is the specific heat of air at constant pressure, $q_l$ is the liquid water mixing ratio, and $Q$ represents external heat sources (e.g. ground solar heating).
3. Continuity (incompressible air flow): $\nabla \cdot V = 0$
4. Moisture and cloud condensation: $\frac{\partial q}{\partial t} = -V \cdot \nabla q$ $q_s(z) = A \exp(-\alpha z)$ $q_l = \max(q - q_s, 0)$ where $q$ is the total water mixing ratio, $q_s(z)$ is the saturation mixing ratio decaying exponentially with altitude $z$ (fixed by boundary conditions $q_s = 0.02$ at the base and $q_s = 0.002$ at the top), and excess condensed liquid water $q_l$ is directly assigned as the local optical volume density $\rho(x,y,z)$ .
Knowl 7 — Limitations in Single-Scattering Darkening and Asymmetric Shadowing
limitation
The volume density ray tracing framework has two key limitations:
1. Low-Albedo Shadowing Artifacts: Single-scattering approximations break down in high-albedo media like clouds (where $\omega_0 \approx 1$ ). When clouds are viewed from the side such that self-shadowing occurs, the shadowed regions appear unnaturally dark because second- and higher-order internal scattering, which physically illuminates these regions, is omitted.
2. Asymmetric Procedural Shadowing: Although volume density clouds cast accurate shadows on other procedural surfaces and on themselves, other procedural solid objects in the scene do not cast shadows into or upon the volume densities.
Knowl 8 — Empirical Computational Performance of Cloud Simulation and Volume Ray Tracing
empirical result
The rendering and dynamical modeling methods yielded the following performance benchmarks:
- Fractal Volume Densities: Rendering clouds synthesized from 3D FFT $1/f$ noise at $512 \times 512$ image resolution on an IBM 4341 processor required 1 to 4 hours of CPU time for grid dimensions ranging from $16 \times 16 \times 16$ to $128 \times 128 \times 16$ . Rendering a fractal cloud intersecting a fractal mountain database at $256 \times 256$ resolution required 6 hours of CPU time on the same machine.
- Hydrodynamic Cloud Simulation: Integrating the 9-equation cumulus model on a $10 \times 10 \times 20$ grid using a forward-differencing scheme on a VAX 11/780 took approximately 10 CPU seconds per time step, corresponding to approximately 1 second of physical cloud evolution. Rendering the resulting density grids at $512 \times 512$ resolution on an IBM 4341 required 2 CPU hours per frame.

Coverage note — Procedural density generation via 3D FFT 1/f noise (Voss method) and particle system rasterization (Reeves method) were omitted as standalone knowls because they represent existing techniques cited and applied as inputs rather than primary theoretical contributions of this paper.

References

1.Anselone, P.M., and Gibbs, A.G., 1974: Convergence of the discrete ordinates method for the transport equation, Constructive and Computational methods for differential and integral equations, Springer Verlag Lecture notes in math 430.
2.Appel, A., 1968: Some techniques for shading machine renderings of solids, 1968 SJCC, 37-45.
3.Blinn, J.F., 1982: Light reflection functions for simulation of clouds and dusty surfaces. Proc. SIGGRAPH82. In Comput. Gr. 16,3, 21-29.
4.Chandrasekhar, S., 1950: Radiative Transfer, Oxford University Press.
5.Clark, T.L., 1979: Numerical Simulations with a threedimensional cloud model: lateral boundary condition experiments and multicellular severe storm simulations. J. of the Atmospheric Sciences, 36, 2191.
6.Courant, R. and Hilbert, D., 1953: Methods of Mathematical Physics v.1, Interscience, New York.
7.Dahlquist, G., and Bjork, A., 1974: Numerical Methods, Prentice Hall, New York.
8.Goldstein, E. and Nagle, R. 1971: 3D visual simulation, Simulation 16, 25-31.
9.Kajiya, J.T., 1983: Ray tracing procedurally defined objects, SIGGRAPH83, Comput. Gr. 17,3, 91-102.
10.Kajiya, J.T., 1982: Ray tracing parametric patches, SIGGRAPH82, Comput. Gr. 16,3, 245-254.
11.Keller, H.B., 1960a: Approximate solutions of transport problems, SIAM J. Appl. Math. 8, 43-73.
12.Keller, H.B., 1960b: On the pointwise convergence of the discrete ordinates method, SIAM J. Appl. Math. 8, 560-567.
13.Max, N., 1983: Panel on the simulation of natural phenomena, Proc. SIGGRAPH83, In Comput. Gr. 17,3, 137-139.
14.Schlesinger, R.E., 1975: A three-dimensional numerical model of an isolated deep convective cloud: Preliminary results. J. of the Atmospheric Sciences, 32, 934-957.
15.Schlesinger, R.E., 1978: A three-dimensional numerical model of an isolated thunderstorm, part I: comparative experiments for variable ambient wind shear. J. of the Atmospheric Sciences, 35, 690-713.
16.Schlesinger, R.E., 1980: A three-dimensional numerical model of an isolated thunderstorm, part II: dynamics of updraft splitting and mesovortex couplet evolution. J of the Atmospheric Sciences, 37, 395.
17.Simpson, J., Van Helvoirt, G., McCumber, M., 1982: Three-dimensional simulations of cumulus congestus clouds on GATE day 261. J. of the Atmospheric Sciences, 39, 126.
18.Reeves, W.T., 1983: Particle systems—a technique for modeling a class of fuzzy objects, ACM Trans. on Graphics, 2,2.
19.Voss, R., 1983: Fourier synthesis of gaussian fractals: 1/f noises, landscapes, and flakes, Tutorial on State of the Art Image Synthesis v.10, SIGGRAPH83.
20.Wallace, J. M., and Hobbs, P. V., 1977: Atmospheric Science, Academic Press, pp.359-407.
21.Whitted, T., 1980: An improved illumination model for shaded display, Comm. ACM 23, 343-349.

Citation

MLA

Kajiya, J. T., and B. P. Von Herzen. “Ray Tracing Volume Densities”. Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques, 1984, pp. 165–74, https://doi.org/10.1145/800031.808594.

APA

Kajiya, J. T., & Von Herzen, B. P. (1984). Ray tracing volume densities. Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques, 165–174. https://doi.org/10.1145/800031.808594

Chicago

Kajiya, J. T., and B. P. Von Herzen. 1984. “Ray Tracing Volume Densities”. Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques, 165–74. https://doi.org/10.1145/800031.808594.

Harvard

Kajiya, J.T. and Von Herzen, B.P. (1984) “Ray tracing volume densities”, Proceedings of the 11th annual conference on Computer graphics and interactive techniques. ACM, pp. 165–174. Available at: https://doi.org/10.1145/800031.808594.

Vancouver

1. Kajiya JT, Von Herzen BP (1984) Ray tracing volume densities. In: Proceedings of the 11th annual conference on Computer graphics and interactive techniques. ACM, pp 165–174

BibTeX

@inproceedings{Kajiya_1984, series={SIGGRAPH ’84}, title={Ray tracing volume densities}, url={http://dx.doi.org/10.1145/800031.808594}, DOI={10.1145/800031.808594}, booktitle={Proceedings of the 11th annual conference on Computer graphics and interactive techniques}, publisher={ACM}, author={Kajiya, James T. and Von Herzen, Brian P}, year={1984}, month=Jan, pages={165–174}, collection={SIGGRAPH ’84} }

Metadata:Crossref

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF

Ray tracing volume densities

Abstract

Topics

Keywords

Table of Contents

Knowls

Knowl 1 — Two-Pass Ray Tracing Algorithm for Volume Densities

Knowl 2 — Radiative Transfer Scattering Equation for Volume Densities

Knowl 3 — Perturbation Expansion for High-Albedo Volume Radiative Transfer

Knowl 4 — Spherical Harmonics Representation and P-Wave Truncation for Radiative Transfer

Knowl 5 — Matrix Elements for Direction and Phase Operators in Spherical Harmonics

Knowl 6 — Dynamical Convection Model for Cumulus Cloud Density Generation

Knowl 7 — Limitations in Single-Scattering Darkening and Asymmetric Shadowing

Knowl 8 — Empirical Computational Performance of Cloud Simulation and Volume Ray Tracing

References

Citation

BibTeX

Access the Paper