Volume Rendering

Volume rendering forms an image from a field that gives every point a density and a colour, not from surfaces: each pixel integrates, along its camera ray, the light emitted at every point, attenuated by everything in front of it. It began as a way to display medical and scientific volumes without extracting surfaces, and it is now the image model of neural rendering methods such as neural radiance fields and, in discrete form, 3D Gaussian splatting.

Definition

Volume rendering computes the colour of a ray through a medium described by a density $\sigma(\mathbf{x}) \ge 0$, the rate at which it blocks light per unit length, and a colour $\mathbf{c}(\mathbf{x})$, the light it emits. The result sums the colour emitted at each depth, weighted by the density there and by the transmittance, the fraction of light surviving from that depth to the camera. This emission–absorption model ignores scattering.

Intuition

Picture the medium as a cloud of tiny opaque glowing particles, the model Max and Chen (2010, §2.1) start from. A ray either hits a particle and returns its colour, or passes through to the background. The density is the chance of a hit per unit length, the transmittance the chance of no hit so far, and the pixel the expected colour. Dense regions act like surfaces; thin ones, like smoke, let light through.

Formal Definition

The physical model

A pixel’s ray is $\mathbf{r}(t) = \mathbf{o} + t\,\mathbf{d}$: $\mathbf{o}$ is the camera centre $-R^\top \mathbf{t}$ and $\mathbf{d}$ the normalized $R^\top K^{-1}(u, v, 1)^\top$, for a pinhole camera with the world-to-camera convention $\mathbf{X}_c = R\,\mathbf{X}_w + \mathbf{t}$ of camera parameters. Write $\sigma(t)$ and $\mathbf{c}(t)$ for the field along it, between bounds $t_n$ and $t_f$.

Max and Chen (2010, §2.1) derive the transmittance from independence: not being blocked before $t$ and in $[t, t + dt]$ are independent events, so $T(t + dt) = T(t)\,(1 - \sigma(t)\,dt)$, which gives the ordinary differential equation

$$\frac{dT}{dt} = -\sigma(t)\,T(t), \qquad T(t_n) = 1, \qquad\Longrightarrow\qquad T(t) = \exp\!\Big(-\int_{t_n}^{t} \sigma(s)\,ds\Big).$$

For a homogeneous medium this reduces to the exponential attenuation of the Bouguer–Lambert law, which Beer extended to solutions (Max & Chen, §5.1). The probability that the ray stops in $[t, t + dt]$ is $T(t)\,\sigma(t)\,dt$, so the expected colour, including a background $\mathbf{c}_{\text{bg}}$ seen with probability $T(t_f)$, is

$$C = \int_{t_n}^{t_f} T(t)\,\sigma(t)\,\mathbf{c}(t)\,dt \;+\; T(t_f)\,\mathbf{c}_{\text{bg}}.$$

This is Max and Chen’s Equation (2). NeRF (Mildenhall et al., 2020, Eq. 1) states it without the background term and with the same particle reading of $\sigma$ and $T$.

Discretization

Split the ray at sample positions $t_n \le t_1 < \dots < t_N$, set $\delta_i = t_{i+1} - t_i$ with $t_{N+1} = t_f$, and treat $\sigma$ and $\mathbf{c}$ as constant, $\sigma_i$ and $\mathbf{c}_i$, on $[t_i, t_{i+1}]$, and zero before $t_1$. Each segment then integrates exactly (Max & Chen, 2010, §2.3), since the transmittance inside it is $T_i\,e^{-\sigma_i (t - t_i)}$:

$$\int_{t_i}^{t_{i+1}} T(t)\,\sigma_i\,\mathbf{c}_i\,dt = T_i\,\mathbf{c}_i \big(1 - e^{-\sigma_i \delta_i}\big).$$

Summing the segments gives the compositing equations:

$$\alpha_i = 1 - \exp(-\sigma_i \delta_i), \qquad T_i = \exp\Big(-\sum_{j<i} \sigma_j \delta_j\Big) = \prod_{j<i} (1 - \alpha_j), \qquad \hat C = \sum_{i=1}^{N} T_i\,\alpha_i\,\mathbf{c}_i \;+\; T_{N+1}\,\mathbf{c}_{\text{bg}}.$$

$\alpha_i$ is the opacity of segment $i$, the probability that the ray stops in it once it has reached it. NeRF’s Equation (3) is this sum without the background term, taken, the paper says, from the quadrature in Max’s (1995) review; it notes the reduction to alpha compositing, and leaves $\delta_N$, which needs a $t_{N+1}$, undefined.

The sum is exact when $\sigma$ and $\mathbf{c}$ are constant on every segment, since each term is the exact integral; otherwise its error comes from approximating the field, not the exponential. Replacing $\alpha_i$ by $\sigma_i \delta_i$ gives a Riemann sum, accurate only while $\sigma_i \delta_i$ is small (Max & Chen, 2010, §2.2–2.3).

Weights, opacity and depth

The weights $w_i = T_i \alpha_i$ are the probabilities that the ray stops in segment $i$. Since $T_i - T_{i+1} = T_i \alpha_i$, the sum telescopes:

$$\sum_{i=1}^{N} w_i = 1 - T_{N+1},$$

the accumulated opacity of the ray; $T_{N+1}$ is the background’s weight. NeRF writes $\hat C = \sum_i w_i \mathbf{c}_i$ (Eq. 5). The same weights composite any per-sample quantity; the expected depth $\hat D = \sum_i w_i t_i / \sum_i w_i$, the mean stopping distance of rays that stop, is one way to read a depth map out of a radiance field. In a thick homogeneous slab the stopping distance beyond the front face is exponential with mean $1/\sigma$, so $\hat D$ lies $1/\sigma$ behind that surface (our derivation).

Front-to-back and back-to-front compositing

The sum can be evaluated in either order:

$$\text{front to back: } \; \mathbf{C} \leftarrow \mathbf{C} + T\,\alpha_i \mathbf{c}_i, \;\; T \leftarrow T\,(1 - \alpha_i); \qquad \text{back to front: } \; \mathbf{C} \leftarrow \alpha_i \mathbf{c}_i + (1 - \alpha_i)\,\mathbf{C}.$$

Back-to-front starts from $\mathbf{C} = \mathbf{c}_{\text{bg}}$ and needs no running transmittance; Levoy (1988) composited this way. Front-to-back starts from $\mathbf{C} = \mathbf{0}$, $T = 1$, adds $T\,\mathbf{c}_{\text{bg}}$ last, and allows early ray termination once $T$ is negligible (Max & Chen, 2010, §2.2); the remaining samples can then change the colour by at most $T$ times the largest colour (our bound). Both orders compute the same number.

The recurrence composites any depth-sorted list of semi-transparent primitives. 3D Gaussian splatting takes each $\alpha_i$ from a projected 2D Gaussian multiplied by a learned opacity and blends depth-sorted Gaussians front to back until the pixel saturates; Kerbl et al. (2023, Eqs. 1–3) show that this blending and NeRF’s sum share one image formation model (their product for $T_i$ misprints $\alpha_j$ as $\alpha_i$).

Gradients

NeRF notes that $\hat C$ is differentiable, which lets a field be fitted to photographs by gradient descent. Using $\partial \alpha_i / \partial \sigma_i = \delta_i (1 - \alpha_i)$ and $\partial T_k / \partial \sigma_i = -\delta_i T_k$ for $k > i$ (our derivation):

$$\frac{\partial \hat C}{\partial \mathbf{c}_i} = w_i, \qquad \frac{\partial \hat C}{\partial \sigma_i} = \delta_i \Big( T_{i+1}\,\mathbf{c}_i - \sum_{k>i} w_k \mathbf{c}_k - T_{N+1}\,\mathbf{c}_{\text{bg}} \Big).$$

Raising $\sigma_i$ adds sample $i$’s colour and dims everything behind it.

Properties

Sampling and step size

A fixed grid shows the field at the same points in every render. NeRF instead uses stratified sampling (Eq. 2), one uniform draw in each of $N$ equal bins, so that the network is evaluated at continuous positions over the course of optimization. A single stratified render is noisier than a uniform one; the benefit is that no region stays permanently between samples.

The step must be small compared with the thinnest structure. A layer thinner than $\delta$ is either missed or, when a sample lands in it, given the opacity of a whole step: aliasing along the ray. Jitter does not repair this; since $\hat C$ is nonlinear in the samples, the average of jittered renders stays biased until $\delta$ shrinks (our reasoning, shown in the examples).

Hierarchical sampling

Most samples land in empty space or behind opaque content. NeRF’s hierarchical volume sampling (§5.2) normalizes the weights of $N_c$ coarse samples into a piecewise-constant PDF, draws $N_f$ more by inverse transform sampling, and composites all $N_c + N_f$ ($N_c = 64$, $N_f = 128$, §5.3), treating them as a nonuniform discretization of the ray rather than as importance-sampling estimates.

Classical Direct Volume Rendering

In scientific and medical visualization, transfer functions map a scalar volume to colour and extinction, which are then integrated along rays (Max & Chen, 2010, §1). Levoy (1988; technical-report version) cast one ray per pixel through CT data and electron density maps, resampled colours and opacities at evenly spaced points by trilinear interpolation, and composited them back to front over an opaque background. Colours came from Phong shading with the data gradient as the normal; opacities from a separate classification, which for an isovalue surface falls off away from the chosen value at a rate inversely proportional to the gradient magnitude, to keep the transition layer’s thickness roughly constant.

Marching cubes instead commits to one level set and outputs a mesh. Levoy argued that such binary decisions create spurious surfaces or holes where features are small or poorly defined, while volume rendering can show weak or fuzzy surfaces and nested semi-transparent ones.

Densities from Signed Distance Functions

A density defines no surface (see 3D representations), so some methods render a signed distance function $f$, negative inside, through a density derived from it. VolSDF (Yariv et al., 2021, Eqs. 2–3) sets $\sigma = a\,\Psi_\beta(-f)$, with $\Psi_\beta$ the cumulative distribution function of a zero-mean Laplace distribution of scale $\beta$ and $a$ a learned scale (the paper’s $\alpha$). NeuS (Wang et al., 2021) states that taking a logistic density of $f$ directly as $\sigma$ puts the weight peak in front of the surface, and uses instead $\alpha_i = \max\big((\Phi_s(f(t_i)) - \Phi_s(f(t_{i+1})))/\Phi_s(f(t_i)), 0\big)$, with $\Phi_s$ the logistic CDF (Eq. 13); it states that its continuous weights peak on the surface where that is locally planar (Theorem 1, proved in the supplement).

Examples

The code renders media with known exact colour; composite returns $\hat C$, $w_i$ and $T_1, \dots, T_{N+1}$. The blocks run in sequence as one script (NumPy and SciPy, about 3 s).

The first block renders a ray on $[0, 4]$ through two slabs (red, $\sigma = 2$ on $[1, 1.5]$; blue, $\sigma = 4$ on $[2, 3]$; exact colour in closed form) and through a Gaussian blob whose colour changes along the ray (exact colour by adaptive quadrature), against white. Errors are the largest RGB error. “Uniform” samples the left end of each bin, “riemann” uses $\sigma_i \delta_i$ for $\alpha_i$, and “stratified” (NeRF’s Eq. 2) runs seeds 0–99; the last column is the error of their mean.

import numpy as np
from scipy.integrate import quad
from scipy.special import erf

T_N, T_F = 0.0, 4.0                       # near and far bounds of the ray segment
BG = np.array([1.0, 1.0, 1.0])            # white background

def composite(sigma, color, delta, bg=BG):
    """Front-to-back compositing. Returns colour, weights, and final transmittance."""
    alpha = 1.0 - np.exp(-sigma * delta)                              # alpha_i
    T = np.exp(-np.concatenate([[0.0], np.cumsum(sigma * delta)]))    # T_1..T_{N+1}
    w = T[:-1] * alpha                                                # w_i = T_i alpha_i
    return w @ color + T[-1] * bg, w, T

def samples(n, rng=None):
    """Uniform (left end of each bin) or stratified (one uniform draw per bin), as in NeRF Eq. 2."""
    edges = np.linspace(T_N, T_F, n + 1)
    t = edges[:-1] if rng is None else edges[:-1] + rng.uniform(0, 1, n) * (edges[1] - edges[0])
    delta = np.diff(np.append(t, T_F))                                # delta_i = t_{i+1} - t_i
    return t, delta

# Scene 1: two homogeneous slabs (red, then blue). Exact colour in closed form.
SLABS = [(1.0, 1.5, 2.0, [0.9, 0.1, 0.1]), (2.0, 3.0, 4.0, [0.1, 0.2, 0.9])]

def slab_field(t):
    s, c = np.zeros_like(t), np.zeros((len(t), 3))
    for a, b, sg, col in SLABS:
        m = (t >= a) & (t < b)
        s[m], c[m] = sg, col
    return s, c

def slab_exact():
    C, T = np.zeros(3), 1.0
    for a, b, sg, col in SLABS:
        C += T * (1 - np.exp(-sg * (b - a))) * np.array(col)
        T *= np.exp(-sg * (b - a))
    return C + T * BG

# Scene 2: a Gaussian blob of density whose colour changes along the ray.
A, MU, S = 3.0, 2.0, 0.3
blob_sigma = lambda t: A * np.exp(-0.5 * ((t - MU) / S) ** 2)
blob_tau = lambda t: A * S * np.sqrt(np.pi / 2) * (erf((t - MU) / (S * np.sqrt(2))) + erf(MU / (S * np.sqrt(2))))
blob_color = lambda t: np.stack([np.clip(0.5 + (t - MU), 0, 1), 0.3 + 0 * t, np.clip(0.5 - (t - MU), 0, 1)], -1)

def blob_field(t):
    return blob_sigma(t), blob_color(t)

def blob_exact():
    C = [quad(lambda t: np.exp(-blob_tau(t)) * blob_sigma(t) * blob_color(np.array([t]))[0, k],
              T_N, T_F, points=[MU], epsabs=1e-13, limit=200)[0] for k in range(3)]
    return np.array(C) + np.exp(-blob_tau(T_F)) * BG

def render(field, n, rng=None, riemann=False):
    t, d = samples(n, rng)
    s, c = field(t)
    if riemann:                           # Riemann sum: alpha_i replaced by sigma_i delta_i
        T = np.exp(-np.concatenate([[0.0], np.cumsum(s * d)]))
        C = (T[:-1] * s * d) @ c + T[-1] * BG
    else:
        C = composite(s, c, d)[0]
    return C

ex = {"slabs": slab_exact(), "blob": blob_exact()}
print("exact slabs", ex["slabs"].round(4), " exact blob", ex["blob"].round(4))
print("scene    N   uniform  riemann   stratified: median [min, max] over seeds 0-99, error of their mean")
for name, field in [("slabs", slab_field), ("blob", blob_field)]:
    for n in [10, 16, 50, 64, 250, 256]:
        runs = [render(field, n, np.random.default_rng(seed)) for seed in range(100)]
        st = [np.abs(C - ex[name]).max() for C in runs]
        print(f"{name:5s} {n:4d}  {np.abs(render(field, n) - ex[name]).max():.1e}"
              f"  {np.abs(render(field, n, riemann=True) - ex[name]).max():.1e}"
              f"   {np.median(st):.1e} [{min(st):.1e}, {max(st):.1e}]  {np.abs(np.mean(runs, 0) - ex[name]).max():.1e}")

Output:

exact slabs [0.6118 0.1422 0.395 ]  exact blob [0.4044 0.3733 0.7003]
scene    N   uniform  riemann   stratified: median [min, max] over seeds 0-99, error of their mean
slabs   10  6.8e-02  4.9e-01   8.2e-02 [4.1e-03, 5.1e-01]  2.3e-02
slabs   16  5.6e-17  2.1e-01   4.5e-02 [5.8e-04, 1.7e-01]  2.6e-03
slabs   50  1.3e-02  7.4e-02   1.8e-02 [3.1e-04, 7.3e-02]  5.8e-03
slabs   64  4.4e-16  4.6e-02   1.3e-02 [1.8e-04, 3.8e-02]  5.0e-04
slabs  250  2.5e-03  1.4e-02   4.0e-03 [8.5e-05, 1.2e-02]  3.0e-04
slabs  256  1.0e-15  1.1e-02   3.0e-03 [1.0e-04, 8.7e-03]  3.1e-04
blob    10  1.0e-02  2.3e-01   4.1e-02 [3.1e-03, 2.1e-01]  2.6e-02
blob    16  1.7e-03  1.4e-01   2.0e-02 [1.2e-03, 7.1e-02]  8.9e-03
blob    50  4.8e-04  4.4e-02   1.9e-03 [2.1e-04, 8.3e-03]  8.6e-04
blob    64  9.2e-05  3.4e-02   1.4e-03 [5.3e-05, 4.7e-03]  8.1e-04
blob   250  2.0e-05  8.6e-03   1.9e-04 [2.5e-05, 5.2e-04]  2.8e-05
blob   256  5.7e-06  8.4e-03   1.7e-04 [7.6e-06, 5.1e-04]  2.9e-05

For the slabs, the uniform sum is exact to rounding when $N$ is a multiple of 8, putting the boundaries on segment ends; otherwise the error falls about as $1/N$. The Riemann sum errs by 0.49 at $N = 10$, where $\sigma_i \delta_i$ reaches 1.6. Single stratified renders are less accurate than uniform ones (median 0.013 for the slabs at $N = 64$); the mean of 100 is closer.

The second block reads opacity, background weight and expected depth for a slab on $[2, 3.5]$ (exact mean $2 + 1/\sigma - L e^{-\sigma L}/(1 - e^{-\sigma L})$, $L = 1.5$), aliases a thin slab, and compares the compositing orders.

# Weights, opacity, background and expected depth for a slab whose front face (the "surface") is at t = 2.
print("sigma  opacity  T_bg      E[depth]  exact E[depth]")
for sg in [2.0, 5.0, 50.0, 500.0]:
    field = lambda t: (np.where((t >= 2.0) & (t < 3.5), sg, 0.0), np.full((len(t), 3), 0.5))
    t, d = samples(256)
    C, w, T = composite(*field(t), d)
    depth = w @ t / w.sum()                                   # expected depth of the ray's termination, given it hits
    L = 1.5
    exact = 2.0 + 1 / sg - L * np.exp(-sg * L) / (1 - np.exp(-sg * L))
    print(f"{sg:5.0f}  {w.sum():.4f}  {T[-1]:.2e}  {depth:.4f}    {exact:.4f}")

# A thin slab (width 0.05, optical depth 2) seen with 16 uniform samples: the result depends on the grid's phase.
thin = lambda t: (np.where((t >= 2.02) & (t < 2.07), 40.0, 0.0), np.full((len(t), 3), 0.5))
edges = np.linspace(T_N, T_F, 17)
op = []
for phase in np.arange(0, 1, 0.1):
    t = edges[:-1] + phase * 0.25
    op.append(composite(thin(t)[0], thin(t)[1], np.diff(np.append(t, T_F)))[1].sum())
strat = []
for seed in range(1000):
    t, d = samples(16, np.random.default_rng(seed))
    strat.append(composite(*thin(t), d)[1].sum())
print(f"thin slab: exact opacity {1 - np.exp(-2.0):.3f}; uniform, 10 grid phases: {np.round(op, 3)}")
print(f"stratified, seeds 0-999: mean {np.mean(strat):.3f}, fraction of renders that miss it {np.mean(np.array(strat) < 1e-9):.2f}")

# Front-to-back with early ray termination versus back-to-front: red slab in front of the opaque slab (sigma 50).
t, d = samples(256)
s = np.where((t >= 1.0) & (t < 1.5), 2.0, 0.0) + np.where((t >= 2.0) & (t < 3.5), 50.0, 0.0)
c = np.where((t < 1.5)[:, None], [0.9, 0.1, 0.1], [0.5, 0.5, 0.5])
a = 1 - np.exp(-s * d)
C_btf = BG.copy()
for i in reversed(range(len(t))):                            # back to front: C <- alpha c + (1 - alpha) C
    C_btf = a[i] * c[i] + (1 - a[i]) * C_btf
for eps in [0.0, 1e-4, 1e-2]:
    C_ftb, T, used = np.zeros(3), 1.0, 0
    for i in range(len(t)):                                  # front to back: C <- C + T alpha c,  T <- T (1 - alpha)
        if T <= eps:
            break
        C_ftb += T * a[i] * c[i]
        T *= 1 - a[i]
        used += 1
    C_ftb += T * BG
    print(f"stop at T<={eps:g}: {used:3d}/256 samples, |front-to-back - back-to-front| = {np.abs(C_ftb - C_btf).max():.1e}")

Output:

sigma  opacity  T_bg      E[depth]  exact E[depth]
    2  0.9502  4.98e-02  2.4136    2.4214
    5  0.9994  5.53e-04  2.1915    2.1992
   50  1.0000  2.68e-33  2.0132    2.0200
  500  1.0000  0.00e+00  2.0000    2.0020
thin slab: exact opacity 0.865; uniform, 10 grid phases: [0. 1. 1. 0. 0. 0. 0. 0. 0. 0.]
stratified, seeds 0-999: mean 0.215, fraction of renders that miss it 0.79
stop at T<=0: 256/256 samples, |front-to-back - back-to-front| = 8.9e-16
stop at T<=0.0001: 139/256 samples, |front-to-back - back-to-front| = 3.4e-05
stop at T<=0.01: 133/256 samples, |front-to-back - back-to-front| = 3.7e-03

Opacity and background weight sum to 1. The expected depth approaches $2 + 1/\sigma$ as $\sigma L$ grows (2.020 at $\sigma = 50$); the discrete values sit up to $\delta/2 = 0.008$ low because each weight sits at the start of its segment. The thin slab (opacity 0.865) is missed by 8 of the 10 uniform grids and rendered opaque by the other 2; 79% of stratified renders miss it, and their mean opacity is 0.215. The two orders agree to $10^{-15}$. Stopping at $T \le 10^{-4}$ skips 117 samples deep inside and behind the opaque slab, with an error of $3.4 \times 10^{-5}$; the empty space in front still costs its samples.

The third block checks the closed-form gradients on a stratified 64-sample render of the blob; behind is $\sum_{k>i} w_k \mathbf{c}_k + T_{N+1} \mathbf{c}_{\text{bg}}$.

# Closed-form gradients of C_hat with respect to sigma_i and c_i, checked by central differences.
t, d = samples(64, np.random.default_rng(0))
s, c = blob_field(t)
C, w, T = composite(s, c, d)
behind = np.cumsum((w[:, None] * c)[::-1], 0)[::-1] - w[:, None] * c + T[-1] * BG   # colour arriving from behind sample i
dC_dsigma = d[:, None] * (T[1:, None] * c - behind)        # (N, 3)
dC_dc = w                                                    # dC/dc_i = w_i for each channel
h, fd_s, fd_c = 1e-6, np.zeros((64, 3)), np.zeros(64)
for i in range(64):
    e = np.zeros(64); e[i] = h
    fd_s[i] = (composite(s + e, c, d)[0] - composite(s - e, c, d)[0]) / (2 * h)
    cp, cm = c.copy(), c.copy(); cp[i, 0] += h; cm[i, 0] -= h
    fd_c[i] = (composite(s, cp, d)[0][0] - composite(s, cm, d)[0][0]) / (2 * h)
print(f"max |dC/dsigma| {np.abs(dC_dsigma).max():.3f}, max error vs finite differences {np.abs(dC_dsigma - fd_s).max():.1e}")
print(f"max |dC/dc|     {np.abs(dC_dc).max():.3f}, max error vs finite differences {np.abs(dC_dc - fd_c).max():.1e}")

Output:

max |dC/dsigma| 0.047, max error vs finite differences 2.3e-10
max |dC/dc|     0.095, max error vs finite differences 1.8e-11

Both agree with central differences to about $10^{-10}$, the differences’ own accuracy.

The fourth block renders a 64 × 64 image of a dense orange ball ($\sigma = 50$) and a thin blue one ($\sigma = 3$) partly in front of it; the exact image follows from the ray–sphere chords.

# A 64 x 64 image of two balls of constant density seen by a pinhole camera, X_c = R X_w + t.
BALLS = [(np.array([0.0, 0.0, 0.0]), 0.6, 50.0, np.array([0.9, 0.5, 0.1])),     # dense: nearly a surface
         (np.array([-0.45, 0.1, -1.0]), 0.35, 3.0, np.array([0.1, 0.3, 0.9]))]  # thin fog in front of it
cam = np.array([0.0, -0.8, -3.0])                     # camera centre in the world frame
z = -cam / np.linalg.norm(cam)                        # look at the origin
x = np.cross([0.0, -1.0, 0.0], z); x /= np.linalg.norm(x)
R = np.stack([x, np.cross(z, x), z])                  # rows: camera axes in world coordinates
tvec = -R @ cam
W, f = 64, 64 / (2 * np.tan(np.radians(20)))          # 40 degree field of view
K = np.array([[f, 0, (W - 1) / 2], [0, f, (W - 1) / 2], [0, 0, 1]])
u, v = np.meshgrid(np.arange(W), np.arange(W))
dirs = np.stack([u, v, np.ones_like(u)], -1).reshape(-1, 3) @ np.linalg.inv(K).T @ R   # R^T K^-1 x per pixel
dirs /= np.linalg.norm(dirs, axis=1, keepdims=True)
origin = -R.T @ tvec                                  # equals cam
NEAR, FAR = 1.0, 5.0

def chord(ctr, r):
    """Entry and exit distances of every ray through a ball (nan where it misses)."""
    oc = origin - ctr
    b = dirs @ oc
    disc = b ** 2 - (oc @ oc - r ** 2)
    sq = np.sqrt(np.where(disc > 0, disc, np.nan))
    return -b - sq, -b + sq

def exact_image():
    segs = [(*chord(c0, r), sg, col) for c0, r, sg, col in BALLS]
    order = np.argsort(np.stack([np.nan_to_num(s[0], nan=np.inf) for s in segs], 1), 1)
    C, T = np.zeros((len(dirs), 3)), np.ones(len(dirs))
    for k in range(len(BALLS)):                       # composite the chords in depth order, per ray
        for j, (t0, t1, sg, col) in enumerate(segs):
            m = (order[:, k] == j) & ~np.isnan(t0)
            a = 1 - np.exp(-sg * (t1[m] - t0[m]))
            C[m] += T[m, None] * a[:, None] * col
            T[m] *= 1 - a
    return C + T[:, None] * BG

def render_image(n, rng=None):
    edges = np.linspace(NEAR, FAR, n + 1)
    jitter = np.zeros((len(dirs), n)) if rng is None else rng.uniform(0, 1, (len(dirs), n))
    t = edges[:-1] + jitter * (edges[1] - edges[0])  # (rays, n)
    delta = np.diff(np.concatenate([t, np.full((len(dirs), 1), FAR)], 1), axis=1)
    p = origin + t[..., None] * dirs[:, None, :]
    sigma, color = np.zeros(t.shape), np.zeros(t.shape + (3,))
    for c0, r, sg, col in BALLS:
        inside = np.linalg.norm(p - c0, axis=-1) < r
        sigma[inside], color[inside] = sg, col
    alpha = 1 - np.exp(-sigma * delta)
    T = np.exp(-np.concatenate([np.zeros((len(dirs), 1)), np.cumsum(sigma * delta, 1)], 1))
    w = T[:, :-1] * alpha
    C = np.einsum("rn,rnk->rk", w, color) + T[:, -1:] * BG
    return C, w, t, T

psnr = lambda a, b: 10 * np.log10(1.0 / np.mean((a - b) ** 2))
gt = exact_image()
print(f"image {W}x{W}: {1 - (gt == BG).all(1).mean():.0%} of pixels see a ball, "
      f"{(~np.isnan(chord(*BALLS[0][:2])[0])).mean():.0%} the dense one")
hit = ~np.isnan(chord(*BALLS[0][:2])[0])               # rays that reach the dense ball
print("  N   PSNR uniform  PSNR stratified, seeds 0-9 [min, max]   samples used on those rays if stopping at T<=1e-4")
for n in [32, 64, 128, 256]:
    Cu, w, t, T = render_image(n)
    ps = [psnr(render_image(n, np.random.default_rng(s))[0], gt) for s in range(10)]
    used = (T[hit, :-1] > 1e-4).sum(1).mean() / n
    print(f"{n:4d}   {psnr(Cu, gt):5.1f} dB      {min(ps):5.1f} - {max(ps):5.1f} dB      {used:.0%}")

# Depth from the weights against the true entry point of the dense ball (rays that miss the fog ball).
t_in, _ = chord(*BALLS[0][:2])
fog_in, _ = chord(*BALLS[1][:2])
Cu, w, t, T = render_image(256)
m = ~np.isnan(t_in) & np.isnan(fog_in)
depth = (w * t).sum(1)[m] / w.sum(1)[m]
print(f"expected depth - true depth over {m.sum()} rays: median {np.median(depth - t_in[m]):+.4f}, "
      f"max |.| {np.abs(depth - t_in[m]).max():.4f}  (1/sigma = {1 / BALLS[0][2]:.3f}, delta = {4 / 256:.4f})")

Output:

image 64x64: 36% of pixels see a ball, 23% the dense one
  N   PSNR uniform  PSNR stratified, seeds 0-9 [min, max]   samples used on those rays if stopping at T<=1e-4
  32    40.1 dB       32.0 -  34.8 dB      50%
  64    42.9 dB       39.9 -  41.7 dB      48%
 128    50.0 dB       46.6 -  47.8 dB      47%
 256    56.7 dB       52.6 -  53.5 dB      48%
expected depth - true depth over 814 rays: median +0.0206, max |.| 0.0287  (1/sigma = 0.020, delta = 0.0156)

PSNR rises by about 6–7 dB per doubling of $N$ beyond 64, and uniform sampling beats single stratified renders by 1–8 dB. Early termination evaluates about half of the samples on rays that reach the dense ball. On 814 rays that see only the dense ball, the expected depth is a median 0.0206 behind the true surface, near $1/\sigma = 0.02$; the second block’s low bias is offset because the first sample inside the ball lies up to $\delta$ past its entry (our reasoning).

Common Misconceptions

• “$\alpha_i = \sigma_i \delta_i$.” That is a first-order approximation; $1 - e^{-\sigma_i \delta_i}$ stays below 1 and is exact for piecewise-constant media.
• “The weights sum to one.” They sum to the opacity $1 - T_{N+1}$; the rest belongs to the background.
• “The expected depth is the surface.” It is the mean stopping distance: behind the front of a homogeneous medium, and between layers when a ray meets several.

Where It Is Used

• Scientific and medical visualization: scanned and simulated volumes, through transfer functions.
• Neural radiance fields: density and colour fitted to photographs through the sum’s gradients.
• Neural surface reconstruction: NeuS and VolSDF, for 3D reconstruction.
• 3D Gaussian splatting: front-to-back compositing of sorted Gaussians.