Volume rendering for NeRF

A neural radiance field (NeRF, Mildenhall et al. 2020) stores a scene as a function F(x, d) → (σ, c): at every 3-D point x, a density σ ≥ 0 (how much light is stopped per unit length) and, seen from direction d, an emitted colour c. There are no surfaces and no triangles. To make a pixel you shoot a ray through it, ask F at a few dozen points along the ray, and add up what reaches the camera. That last step — volume rendering — is all this lesson is about, and it is the same sum 3D Gaussian splatting evaluates per pixel, so it is worth knowing cold.

We'll derive the integral, turn it into the alpha-compositing sum the code uses, then cover where to put the samples (stratified and hierarchical sampling), the positional encoding that lets a small MLP represent sharp detail, and how the whole thing is trained.

1. The volume rendering integral

Take a ray r(t) = o + t·d from the camera centre o (see cameras and projection for how a pixel becomes a ray) between a near and a far bound. Model the medium as particles that absorb light and emit light of colour c — no scattering (Max 1995, the "emission–absorption" model).

A thin slab of thickness dt at depth t stops the fraction σ(t)·dt of the light crossing it. So the light surviving from the camera to depth t, the transmittance, obeys dT/dt = −σ(t)·T(t), T(t_n) = 1:

T(t) = exp(−∫_{t_n}^{t} σ(s) ds) (Beer–Lambert)

The slab at t emits σ(t)·c(t)·dt towards the camera (things that absorb more also emit more), and only T(t) of that arrives. Summing over the ray:

C(r) = ∫_{t_n}^{t_f} T(t) · σ(t) · c(t) dt (NeRF eq. 1)

There's a useful reading: T(t)·σ(t) is a probability density — the probability that the ray "terminates" (hits a particle) at depth t. C is then the expected colour at the termination point, and ∫ T σ t dt is the expected depth.

2. From the integral to a sum

The network is only evaluated at N sample depths t_0 < t_1 < … < t_{N−1}. Assume σ and c are constant on each interval [t_i, t_{i+1}), of length δ_i = t_{i+1} − t_i. Then every piece of the integral is exact:

• the light that gets through interval i is exp(−σ_i δ_i), so the interval stops the fraction α_i = 1 − exp(−σ_i δ_i);
• the light reaching interval i is T_i = Π_{j<i} (1 − α_j) = exp(−Σ_{j<i} σ_j δ_j);
• integrating T(t)·σ_i·c_i over the interval gives T_i · α_i · c_i.

C ≈ Σ_i T_i α_i c_i = Σ_i w_i c_i, w_i = T_i α_i (NeRF eq. 3)

This is alpha compositing (Porter & Duff 1984): each sample is a semi-transparent layer of opacity α_i, stacked front to back. Two details decide whether your renderer is right:

1. T_i is exclusive: the product runs over j < i. Sample i is not hidden by itself.
2. The last interval has no t_N. NeRF gives it δ = 1e10: it reaches to infinity, so a dense last sample is fully opaque. (Not INFINITY: 0 × ∞ is NaN when σ = 0.)

The weights have a budget: since T_{i+1} = T_i − w_i, Σ_{i<N} w_i = 1 − T_N. The sum of weights is the accumulated opacity — how much of the pixel the samples cover — and whatever is left, T_N, comes from behind them:

• colour = Σ w_i c_i + (1 − Σ w_i) · background (NeRF's Blender scenes use white: rgb + (1 − acc) in raw2outputs);
• depth = Σ w_i t_i, the expected termination depth. On a ray that misses everything it is ≈ 0, not "far": divide by the opacity if you want a depth map.

#include <cmath>
#include <print>
#include <vector>

int main() {
    std::vector<float> t{2.0f, 2.3f, 2.9f, 3.1f, 3.8f};
    std::vector<float> sigma{0.0f, 1.0f, 8.0f, 30.0f, 0.5f};

    float T = 1.0f, opacity = 0.0f, depth = 0.0f, optical = 0.0f;
    for (std::size_t i = 0; i < t.size(); ++i) {
        float delta = (i + 1 < t.size()) ? t[i + 1] - t[i] : 1e10f;
        float alpha = 1.0f - std::exp(-sigma[i] * delta);
        float w = T * alpha;                       // T_i: light that reached sample i
        std::println("i={} delta={:<6g} alpha={:.3f} T={:.3f} w={:.3f}", i, delta, alpha, T, w);
        opacity += w;
        depth += w * t[i];
        T *= 1.0f - alpha;                         // T_{i+1}
        optical += sigma[i] * delta;               // Σ σδ, for comparison
    }
    std::println("opacity={:.4f} depth={:.3f} T_left={:.4f}", opacity, depth, T);
    std::println("opacity + T_left = {:.4f}", opacity + T);
    // The product of (1 − α) equals exp(−Σσδ) over the finite intervals:
    std::println("T after 4 samples: {:.6f} vs exp(-sum) {:.6f}",
                 std::exp(-0.0f) * std::exp(-0.6f) * std::exp(-1.6f) * std::exp(-21.0f),
                 std::exp(-(0.0f + 0.6f + 1.6f + 21.0f)));
}

Most of the light stops at the samples at 2.3 and 2.9; the very dense sample at 3.1 takes whatever is left (about 0.11), and the last sample, even with δ = 1e10, sees practically nothing. The depth comes out ≈ 2.65: an expectation, so it can land between surfaces when a ray's weight is split between them — one reason NeRF depth maps look soft on semi-transparent geometry.

3. The over operator, two ways

Compositing can run in either direction:

• front to back: keep T; C += T·α_i·c_i; T *= (1 − α_i). Background last, weighted by T.
• back to front: start with C = background; for i from the back, C = α_i c_i + (1 − α_i) C — the over operator. No T needed: each layer covers the fraction α_i of whatever is behind it.

Both give the same colour (over is associative). Rasterisers that blend sorted transparent triangles use back to front (that's what glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA) computes). Ray marchers and 3DGS's tile rasteriser use front to back, because T tells you when to stop: once T < ε (3DGS uses 1e-4), nothing behind can change the pixel by more than ε. That's early ray termination, and it makes rendering a surface cost the samples in front of it, not the whole ray.

#include <cstddef>
#include <print>

int main() {
    const float alpha[] = {0.2f, 0.5f, 0.9f, 0.6f};
    const float color[] = {1.0f, 0.2f, 0.7f, 0.0f};   // one channel
    const float bg = 1.0f;
    const std::size_t n = 4;

    // front to back, stopping once T < eps
    float T = 1.0f, front = 0.0f;
    std::size_t used = 0;
    for (std::size_t i = 0; i < n && T >= 0.05f; ++i, ++used) {
        front += T * alpha[i] * color[i];
        T *= 1.0f - alpha[i];
    }
    front += T * bg;

    // back to front with the over operator
    float back = bg;
    for (std::size_t i = n; i-- > 0;) back = alpha[i] * color[i] + (1.0f - alpha[i]) * back;

    std::println("front-to-back (stopped after {} samples): {:.4f}", used, front);
    std::println("back-to-front over:                      {:.4f}", back);
}

With ε = 0.05 the loop stops after three samples. The skipped samples and the background can only change the result by at most T at the stop (0.04 here, times the largest colour difference), which is exactly the bound you choose with ε.

4. Where to put the samples

Stratified sampling. Evaluating F at the same fixed depths for every ray would let the network see only a discrete set of points and the scene would be learned at that resolution. Instead split [t_n, t_f] into N equal bins and draw one sample uniformly inside each:

t_i = t_n + (i + u_i) · (t_f − t_n) / N, u_i ~ U[0, 1)

Over many iterations every depth gets evaluated, so the MLP learns a continuous field, while each ray still covers the whole range evenly. At test time NeRF uses u_i = 0.5 (bin centres) for a deterministic image.

Hierarchical sampling. Most of a ray is empty space or hidden behind a surface; samples there are wasted. NeRF trains two networks. The coarse one is evaluated at N_c = 64 stratified samples; its weights w_i, normalised, are a piecewise-constant pdf over the ray — "where did the ray probably stop". Then N_f = 128 more samples are drawn from that pdf by inverse transform sampling, and the fine network is evaluated at all 64 + 128 samples, sorted by depth. Both renders are supervised.

Inverse transform sampling: the CDF of a piecewise-constant pdf is piecewise linear. For a uniform u ∈ [0, 1), find the bin whose CDF segment contains u (binary search) and interpolate linearly inside it. Where the pdf is high the CDF is steep, so many u values map into a short depth range.

Two numerical guards from the reference code:

• add 1e-5 to every weight before normalising. A ray through empty space has all-zero weights, and 0/0 = NaN would poison the batch. With the constant the pdf becomes uniform there;
• when a CDF segment is shorter than 1e-5, use 1 as the interpolation denominator instead of dividing by (almost) zero.

The bins are the midpoints between the coarse samples, and NeRF passes weights[1:-1] — the two end samples don't have a full bin.

#include <algorithm>
#include <cstddef>
#include <print>
#include <vector>

int main() {
    std::vector<float> bins{2, 3, 4, 5, 6};
    std::vector<float> w{0.1f, 0.1f, 0.6f, 0.2f};
    const std::size_t n = w.size();

    float total = 0.0f;
    for (float x : w) total += x + 1e-5f;
    std::vector<float> cdf(n + 1, 0.0f);
    for (std::size_t i = 0; i < n; ++i) cdf[i + 1] = cdf[i] + (w[i] + 1e-5f) / total;

    for (float u : {0.125f, 0.375f, 0.625f, 0.875f}) {
        std::size_t idx = static_cast<std::size_t>(std::upper_bound(cdf.begin(), cdf.end(), u) - cdf.begin());
        std::size_t below = idx > 0 ? idx - 1 : 0, above = std::min(n, idx);
        float denom = cdf[above] - cdf[below];
        if (denom < 1e-5f) denom = 1.0f;
        float t = bins[below] + (u - cdf[below]) / denom * (bins[above] - bins[below]);
        std::println("u = {:.3f} -> bin {} -> t = {:.3f}", u, below, t);
    }
}

5. Positional encoding

An MLP fed raw (x, y, z) produces blurry scenes: networks are biased towards low-frequency functions (spectral bias; Rahaman et al. 2019, explained via the neural tangent kernel in Tancik et al. 2020). NeRF maps each coordinate through sines and cosines of doubling frequencies first:

γ(p) = (sin(2⁰πp), cos(2⁰πp), …, sin(2^{L−1}πp), cos(2^{L−1}πp))

with L = 10 for positions and L = 4 for viewing directions. Nearby points with different fine detail now differ a lot in some input, and a smooth function of γ(p) can be a sharp function of p. The reference code differs from the paper in two ways you'll meet when reading it: it omits π (the scenes are normalised to a small box, so the lowest frequency still varies slowly over the scene), and it prepends the raw input, so a 3-D position becomes 3 + 3·2·10 = 63 values and a direction 3 + 3·2·4 = 27. Its layout is per frequency: all sines, then all cosines.

#include <cmath>
#include <print>
#include <vector>

std::vector<float> encode(const std::vector<float>& x, int L) {
    std::vector<float> out(x.begin(), x.end());
    for (int k = 0; k < L; ++k) {
        float f = std::ldexp(1.0f, k);           // 2^k, exact; never 1 << k
        for (float v : x) out.push_back(std::sin(f * v));
        for (float v : x) out.push_back(std::cos(f * v));
    }
    return out;
}

int main() {
    // Two points 0.01 apart: the raw coordinates barely differ, high frequencies differ a lot.
    auto a = encode({0.50f}, 10), b = encode({0.51f}, 10);
    for (int k = 0; k < 10; k += 3)
        std::println("k={} sin: {:+.3f} vs {:+.3f}", k, a[1 + 2 * k], b[1 + 2 * k]);
    std::println("size for xyz, L=10: {}", encode({0, 0, 0}, 10).size());
}

Precision note: in float, 2^k·x is exact (multiplying by a power of two only changes the exponent), but sin of a large argument loses relative precision; past k ≈ 15 the features are noise for float inputs. Instant-NGP replaces this encoding with a learned hash grid (next skills), and mip-NeRF with an integrated version that blurs frequencies finer than a pixel's footprint.

6. Training a NeRF

• Data: posed photos (camera intrinsics and extrinsics from COLMAP, or known for synthetic scenes). Each pixel is a ray with a known colour.
• Network: 8 fully connected layers of 256 with a skip connection; σ depends on x only (multi-view consistency), c on x and d (view-dependent effects such as specular highlights). σ goes through a ReLU, c a sigmoid.
• Loss: squared error between rendered and true colours, summed over the coarse and fine renders. Quality is reported as PSNR = −10·log₁₀(MSE) for colours in [0, 1] — 10 dB per 10× in MSE. NeRF reached ~31 dB on the synthetic Blender scenes.
• Optimiser: Adam, learning rate 5·10⁻⁴ decaying to 5·10⁻⁵, 100–300k iterations: one to two days on a 2020 GPU. That cost is what Instant-NGP (hash grids: minutes) and 3D Gaussian splatting (explicit primitives, rasterised: real time) attacked.

The gradient of the loss with respect to σ_i has two parts: raising σ_i makes sample i more opaque (more of c_i) but hides everything behind it (less of c_{j>i} and the background). Autodiff handles it, but you'll write it by hand in differentiable rendering.

7. Numerical pitfalls

The 1e-10 inside NeRF's cumprod keeps T from becoming exactly 0 after an opaque sample, which would kill the gradient flowing to samples behind it (and log-space implementations avoid it by summing σδ instead of multiplying).

In Python (NumPy / PyTorch)

dists  = z[..., 1:] - z[..., :-1]
dists  = cat([dists, 1e10 * ones_like(dists[..., :1])], -1) * norm(rays_d)
alpha  = 1 - exp(-relu(sigma) * dists)
T      = cumprod(cat([ones(N, 1), 1 - alpha + 1e-10], -1), -1)[:, :-1]   # exclusive!
w      = alpha * T
rgb    = (w[..., None] * c).sum(-2);  acc = w.sum(-1);  depth = (w * z).sum(-1)
rgb    = rgb + (1 - acc[..., None])                                       # white background

Note the * norm(rays_d): when ray directions aren't unit length, t is not a distance, and δ must be scaled to world units.

In an interview

• "Derive α from σ." Beer–Lambert on a constant-density slab: transmittance exp(−σδ), so α = 1 − exp(−σδ); T_i is the product of the slabs in front.
• "Why front to back?" Early termination, and it's how GPUs blend sorted splats in 3DGS (tile-based, stop at T < 1e-4).
• "Why the 1e10?" The last sample owns the rest of the ray; and 0·∞ = NaN.
• "How does hierarchical sampling work?" Coarse weights → piecewise-constant pdf → invert the piecewise-linear CDF with a binary search → fine samples where the surface is.
• "What does the positional encoding fix?" Spectral bias of MLPs: without it NeRF is blurry.