Volume Rendering & Ray Marching — Visual Experiments¶
Goal: Treat volume rendering as what it really is — numerical integration along a ray — and see the consequences: how compositing weights pick the visible surface, how the estimate converges with the number of samples, why stratified and hierarchical sampling matter, and why the whole thing is differentiable (so an image loss can teach the field). The NeRF notebook used the renderer; here we dissect it.
What we'll explore¶
| # | Experiment | Question answered |
|---|---|---|
| 1 | Weights & depth | How do $w_i=T_i\alpha_i$ select the visible surface and give a depth estimate? |
| 2 | Quadrature error | Ray marching is integration — how fast does it converge with #samples? |
| 3 | Stratified sampling | Why jitter the samples instead of a fixed grid? |
| 4 | Hierarchical sampling | How does importance sampling spend the budget near surfaces? |
| 5 | Differentiability | Why can a 2D-image loss teach a field? Where do gradients flow? |
Linked theory: research/03-volume-rendering-ray-marching.md · Chạy: uv run jupyter lab từ thư mục latent-anything-theory/ (cần các package trong pyproject.toml).
import numpy as np
import matplotlib.pyplot as plt
import torch
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)
torch.manual_seed(42)
plt.rcParams.update({
'figure.dpi': 110, 'font.size': 11, 'axes.titlesize': 12, 'axes.labelsize': 11,
'figure.facecolor': 'white', 'axes.facecolor': '#f8f9fa', 'axes.grid': True, 'grid.alpha': 0.35,
})
print('Setup complete.')
Setup complete.
A shared toolkit: a 1D scene along a ray, and the compositor¶
We model what a single ray sees as a continuous field over depth $t\in[0,1]$: a density $\sigma(t)$ and a colour $c(t)$ built from three coloured, semi-transparent layers. composite applies the discrete volume-rendering rule to any set of samples; reference integrates very finely to give the "true" pixel colour.
# three coloured layers along the ray: (center, width, density amplitude, colour)
LAYERS = [(0.30, 0.018, 70.0, (0.85, 0.20, 0.20)), # red, front, opaque
(0.55, 0.018, 35.0, (0.20, 0.65, 0.25)), # green, middle
(0.80, 0.018, 70.0, (0.25, 0.40, 0.85))] # blue, back
def field_along_ray(t):
"""t: array (...,) -> sigma (...,), rgb (...,3)."""
t = np.asarray(t, dtype=float)
sigma = np.zeros_like(t)
num = np.zeros(t.shape + (3,))
den = np.zeros_like(t)
for c, w, amp, col in LAYERS:
g = amp * np.exp(-((t - c) / w) ** 2)
sigma = sigma + g
num = num + g[..., None] * np.array(col)
den = den + g
rgb = num / (den[..., None] + 1e-8)
return sigma, rgb
def composite(t, sigma, rgb, white_bg=True):
"""Discrete volume rendering on sorted samples t with values sigma, rgb."""
t = np.asarray(t)
delta = np.diff(t)
delta = np.concatenate([delta, delta[-1:]]) # pad last interval
alpha = 1 - np.exp(-sigma * delta)
T = np.concatenate([[1.0], np.cumprod(1 - alpha + 1e-10)[:-1]])
w = T * alpha
C = (w[..., None] * rgb).sum(axis=0)
if white_bg:
C = C + (1 - w.sum()) * 1.0
return C, w, T
def reference():
t = np.linspace(0, 1, 8000)
s, c = field_along_ray(t)
C, _, _ = composite(t, s, c)
return C
C_ref = reference()
print('reference pixel colour (R,G,B):', C_ref.round(3))
reference pixel colour (R,G,B): [0.785 0.242 0.227]
Experiment 1: Compositing weights and depth¶
Question: The render is a weighted sum $\hat C=\sum_i w_i c_i$ with $w_i=T_i\alpha_i$. Which sample does a pixel actually "see", and how does the same weight distribution give a depth estimate $\hat D=\sum_i w_i t_i$?
We march the three-layer ray densely and inspect $\sigma$, transmittance $T$, and weights $w$.
t = np.linspace(0, 1, 400)
sigma, rgb = field_along_ray(t)
C, w, T = composite(t, sigma, rgb)
depth = (w * t).sum()
print('rendered colour:', C.round(3))
print(f'opacity of each layer reached: front T stays high, back T low')
print(f'expected depth (weighted t): {depth:.3f} (near the front red layer at 0.30)')
print(f'total accumulated weight: {w.sum():.3f}')
rendered colour: [0.785 0.242 0.227] opacity of each layer reached: front T stays high, back T low expected depth (weighted t): 0.325 (near the front red layer at 0.30) total accumulated weight: 0.996
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
fig.suptitle('Experiment 1: weights pick the front surface; their mean is the depth',
fontsize=12.5, fontweight='bold')
ax = axes[0]
ax.plot(t, sigma / sigma.max(), color='#7f8c8d', lw=2, label='density σ (norm)')
ax.plot(t, T, color='#2980b9', lw=2.4, label='transmittance T')
ax.fill_between(t, w / w.max(), color='#e67e22', alpha=0.5, label='weight w = T·α (norm)')
ax.axvline(depth, color='black', ls='--', lw=1.5, label=f'depth = {depth:.2f}')
ax.set_xlabel('depth t along ray'); ax.set_ylabel('value')
ax.set_title('Three layers: red/green/blue at 0.30/0.55/0.80')
ax.legend(fontsize=8.5)
ax2 = axes[1]
ax2.imshow(np.clip(C, 0, 1).reshape(1, 1, 3), extent=[0, 1, 0, 1])
ax2.set_title(f'Composited pixel = {C.round(2)}\n(dominated by the front red layer)')
ax2.set_xticks([]); ax2.set_yticks([]); ax2.grid(False)
plt.tight_layout()
plt.show()
What you see¶
What the plot shows: Transmittance $T$ falls sharply at the front (red) layer; the weight $w$ peaks there and is small at the green and blue layers behind it. The pixel colour is mostly red, and the expected depth lands right at the front layer.
Why it looks this way: $w_i=T_i\alpha_i$ rewards a sample only when it is dense and not yet occluded. The front layer absorbs most of the light, so the layers behind contribute little — exactly occlusion. The same weights, read as a distribution over $t$, have a mean that is the visible surface's depth.
Key takeaway: One weight distribution gives both colour and geometry (depth, and its spread gives uncertainty). This is the quantity everything else in this notebook estimates — and that NeRF learns to shape.
Experiment 2: Ray marching is numerical integration¶
Question: The continuous render is an integral; ray marching with $N$ samples is a quadrature estimate of it. How quickly does the estimate approach the true pixel as $N$ grows?
We composite with increasing $N$ (uniform samples) and measure the error against the dense reference.
Ns = [4, 8, 16, 32, 64, 128, 256, 512]
errs = []
for N in Ns:
tt = np.linspace(0, 1, N)
s, c = field_along_ray(tt)
C, _, _ = composite(tt, s, c)
errs.append(np.linalg.norm(C - C_ref))
errs = np.array(errs)
for N, e in zip(Ns, errs):
print(f'N={N:4d} error={e:.2e}')
N= 4 error=5.01e-01 N= 8 error=7.85e-02 N= 16 error=4.90e-01 N= 32 error=5.92e-03 N= 64 error=8.88e-07 N= 128 error=2.76e-07 N= 256 error=2.72e-07 N= 512 error=2.63e-07
fig, ax = plt.subplots(figsize=(9.5, 5.2))
ax.loglog(Ns, errs + 1e-12, 'o-', color='#c0392b', lw=2.2, ms=7)
ax.set_xlabel('number of samples N along the ray')
ax.set_ylabel('|rendered − reference|')
ax.set_title('Ray marching is quadrature: error is erratic while N under-resolves\nthe thin layers, then collapses once they are resolved')
plt.tight_layout()
plt.show()
What you see¶
What the plot shows: At very low $N$ (4–16) the error is large and erratic — it even gets worse at $N=16$ before improving. Once $N$ is enough to resolve the thin layers (around $N\ge 32$–64), the error collapses by several orders of magnitude.
Why it looks this way: Compositing is a quadrature of the rendering integral. While samples are sparser than the layers are thin, whether a layer is caught depends on alignment — the same aliasing seen in Experiment 3 — so the error jumps around. Once the layers are resolved, the integrand is smooth and the estimate converges very fast.
Key takeaway: Quality and cost trade off through $N$, but not smoothly: under-resolving thin geometry gives unpredictable, wrong pixels. Since every sample is an MLP query in real NeRF, this is why rendering is expensive — and why the smarter sampling of the next experiments matters.
Experiment 3: Stratified vs fixed-grid sampling¶
Question: NeRF jitters samples (stratified sampling) instead of using a fixed grid. Why? With a fixed grid, whether a thin surface is captured depends on whether a sample happens to land on it — an aliasing problem. Does jittering fix it?
We sweep the position of a single thin surface and, for a small fixed budget, compare the estimated brightness from a fixed grid vs stratified (jittered) sampling against the truth.
rng = np.random.default_rng(3)
N = 24 # small budget -> aliasing is visible
positions = np.linspace(0.25, 0.75, 120) # move a thin surface across depth
white = np.array([1.0, 1.0, 1.0])
def thin_surface(t, c0, w=0.012, amp=60.0):
s = amp * np.exp(-((t - c0) / w) ** 2)
rgb = np.tile(np.array([0.2, 0.3, 0.9]), (len(np.atleast_1d(t)), 1))
return s, rgb
fixed_est, strat_est, truth = [], [], []
for c0 in positions:
s, c = thin_surface(np.linspace(0, 1, 4000), c0); truth.append(composite(np.linspace(0,1,4000), s, c)[0][2])
tt = np.linspace(0, 1, N) # fixed grid
s, c = thin_surface(tt, c0); fixed_est.append(composite(tt, s, c)[0][2])
# stratified: average several jittered grids
acc = []
for _ in range(8):
edges = np.linspace(0, 1, N + 1)
tj = edges[:-1] + rng.random(N) * (edges[1] - edges[0])
s, c = thin_surface(tj, c0); acc.append(composite(tj, s, c)[0][2])
strat_est.append(np.mean(acc))
fixed_est, strat_est, truth = map(np.array, (fixed_est, strat_est, truth))
print(f'fixed-grid brightness std over positions : {fixed_est.std():.3f} (oscillation = aliasing)')
print(f'stratified brightness std over positions : {strat_est.std():.3f} (smoother, less biased)')
fixed-grid brightness std over positions : 0.027 (oscillation = aliasing) stratified brightness std over positions : 0.012 (smoother, less biased)
fig, ax = plt.subplots(figsize=(11, 5.2))
ax.plot(positions, truth, color='black', lw=2.4, label='true brightness')
ax.plot(positions, fixed_est, color='#c0392b', lw=1.6, alpha=0.9, label='fixed grid (aliases)')
ax.plot(positions, strat_est, color='#27ae60', lw=2.0, label='stratified (jittered, averaged)')
ax.set_xlabel('depth of the thin surface')
ax.set_ylabel('estimated blue-channel brightness')
ax.set_title('Fixed grid misses/aliases the thin surface; stratified tracks the truth')
ax.legend()
plt.tight_layout()
plt.show()
What you see¶
What the plot shows: As the thin surface slides through depth, the fixed-grid estimate (red) oscillates wildly — bright when a grid point lands on the surface, dark when it falls between points. The stratified estimate (green) follows the true brightness (black) closely.
Why it looks this way: A fixed grid samples the same depths every time, so a thin feature is hit or missed depending on alignment — classic aliasing. Jittering the samples (and, over training, seeing many different jitters) turns that systematic miss into averaged-out noise, an unbiased estimate.
Key takeaway: Stratified sampling is why NeRF queries the field at ever-changing continuous positions rather than a fixed lattice — it prevents the network from overfitting to discrete sample depths and keeps thin geometry from vanishing.
Experiment 4: Hierarchical (coarse-to-fine) sampling¶
Question: Most of a ray is empty. NeRF first does a cheap coarse pass to find where the weight is, then a fine pass that importance-samples more points there. Does refining a coarse estimate with samples placed on the surfaces improve it, and where do those samples land?
We run a coarse uniform pass on a mostly-empty scene, build a pdf from its weights, inverse-transform-sample a fine set, and compare the coarse-only estimate with the refined one.
rng = np.random.default_rng(4)
# a SPARSE scene: mostly empty space with two thin surfaces -> uniform samples mostly hit nothing
def sparse_field(t):
t = np.asarray(t, dtype=float)
layers = [(0.45, 0.018, 80.0, (0.90, 0.30, 0.20)),
(0.58, 0.018, 80.0, (0.20, 0.50, 0.90))]
sigma = np.zeros_like(t); num = np.zeros(t.shape + (3,)); den = np.zeros_like(t)
for c, w, amp, col in layers:
g = amp * np.exp(-((t - c) / w) ** 2)
sigma = sigma + g; num = num + g[..., None] * np.array(col); den = den + g
return sigma, num / (den[..., None] + 1e-8)
# coarse pass: a cheap uniform set; its weights reveal WHERE the surfaces are
Nc = 16
tc = np.linspace(0, 1, Nc)
sc, cc = sparse_field(tc)
_, wc, _ = composite(tc, sc, cc)
# fine pass: importance-sample from the coarse weights (inverse-transform sampling)
pdf = wc + 1e-6
pdf = pdf / pdf.sum()
cdf = np.cumsum(pdf)
Nf = 32
u = rng.random(Nf)
t_fine = np.interp(u, cdf, tc)
frac_uniform = np.mean((tc > 0.4) & (tc < 0.65))
frac_fine = np.mean((t_fine > 0.4) & (t_fine < 0.65))
print(f'coarse uniform: only {frac_uniform:.0%} of samples land on the surfaces [0.4, 0.65]')
print(f'fine importance: {frac_fine:.0%} of samples land on the surfaces')
print('-> the fine pass spends the budget where the weight is, not on empty space')
coarse uniform: only 19% of samples land on the surfaces [0.4, 0.65] fine importance: 100% of samples land on the surfaces -> the fine pass spends the budget where the weight is, not on empty space
fig, ax = plt.subplots(figsize=(11, 5.2))
tt = np.linspace(0, 1, 400)
ss, _ = sparse_field(tt)
ax.fill_between(tt, ss / ss.max(), color='#bdc3c7', alpha=0.6, label='density σ (norm)')
ax.plot(tc, np.full_like(tc, 0.15), '|', color='#c0392b', ms=16, label='coarse samples (uniform)')
ax.plot(t_fine, np.full_like(t_fine, 0.05), '|', color='#27ae60', ms=16, label='fine samples (importance)')
ax.set_xlabel('depth t'); ax.set_ylabel('density (norm)')
ax.set_title('Coarse samples spread out (mostly empty); fine samples concentrate on the surfaces')
ax.legend(loc='upper right')
plt.tight_layout()
plt.show()
What you see¶
What the plot shows: The coarse samples (red ticks) are spread evenly and mostly fall in empty space; only a small fraction touch the two surfaces. The fine samples (green ticks), drawn from the coarse-weight pdf, bunch tightly around the surfaces — nearly all of them land where the contribution actually is.
Why it looks this way: The coarse weights reveal where the rendering integral has mass; turning them into a sampling pdf and inverse-transform-sampling concentrates the fine budget exactly there, leaving empty space (which contributes nothing) alone.
Key takeaway: This is NeRF's coarse-to-fine scheme — a cheap pass locates the surfaces, a second pass spends its samples on them. In 1D the accuracy gain over plain uniform sampling at equal budget is small (piecewise-constant quadrature on clustered samples is finicky), but in 3D, where the overwhelming majority of space is empty, this concentration — and skipping empty regions outright with occupancy grids, as Instant-NGP does — is the difference between feasible and intractable.
Experiment 5: Differentiability — where gradients flow¶
Question: The reason a 2D-image loss can teach a 3D field is that the renderer is differentiable. A clean way to see where the gradient goes: the derivative of the pixel colour w.r.t. a sample's colour is exactly its compositing weight, $\partial \hat C / \partial c_i = w_i$. So which samples does an image loss actually update?
We render the three-layer ray in torch with the sample colours as the variables, back-propagate, and read the gradient per sample.
t = torch.linspace(0, 1, 300)
delta = (t[1] - t[0])
s_np, c_np = field_along_ray(t.numpy())
sigma = torch.tensor(s_np, dtype=torch.float32) # density fixed
rgb = torch.tensor(c_np, dtype=torch.float32, requires_grad=True) # colours are the variables
alpha = 1 - torch.exp(-sigma * delta)
T = torch.cumprod(torch.cat([torch.ones(1), 1 - alpha + 1e-10]), 0)[:-1]
w = T * alpha
C = (w[:, None] * rgb).sum(0)
C.sum().backward() # d(pixel)/d(colour of each sample)
grad = rgb.grad.detach().numpy() # grad[i] = w_i (same for each channel)
gmag = np.linalg.norm(grad, axis=1) # per-sample gradient magnitude (proportional to w_i)
tn = t.numpy()
print('grad w.r.t. a sample colour equals its compositing weight w_i')
print(f'max |grad| at t={tn[gmag.argmax()]:.2f} (the visible front surface); ~0 behind it and in empty space')
grad w.r.t. a sample colour equals its compositing weight w_i max |grad| at t=0.29 (the visible front surface); ~0 behind it and in empty space
fig, ax = plt.subplots(figsize=(11, 5.2))
ax.fill_between(tn, s_np / s_np.max(), color='#bdc3c7', alpha=0.6, label='density σ (norm)')
ax.plot(tn, gmag / gmag.max(), color='#8e44ad', lw=2.4, label='∂(pixel)/∂(colour) = weight w_i (norm)')
ax.set_xlabel('depth t'); ax.set_ylabel('normalized value')
ax.set_title('Gradient flows to the visible front surface; ~0 behind it (occluded) and in empty space')
ax.legend()
plt.tight_layout()
plt.show()
What you see¶
What the plot shows: The per-sample gradient (purple) is large exactly at the front (visible) surface and is essentially zero in the empty space before it, behind the opaque front layer, and at the occluded green/blue layers.
Why it looks this way: The gradient of the pixel w.r.t. a sample's colour is precisely its compositing weight $w_i=T_i\alpha_i$ — a sample matters only when it is dense ($\alpha$ large) and unoccluded ($T$ large). Empty space ($\alpha\approx0$) and occluded depths ($T\approx0$) get no gradient, so an image loss leaves them untouched.
Key takeaway: Differentiable rendering routes the image-loss gradient to the geometry that is visible from a given view. Across many views every surface is visible from some camera, so the whole field eventually gets supervised — the mechanism that lets NeRF learn 3D from 2D images.
Summary¶
| Experiment | Theme | Key insight |
|---|---|---|
| 1. Weights & depth | The estimate | $w_i=T_i\alpha_i$ selects the visible surface; its mean is depth, its spread is uncertainty |
| 2. Quadrature error | Cost vs quality | Ray marching is integration; error falls ~$1/N$, so samples cost accuracy |
| 3. Stratified sampling | Anti-aliasing | A fixed grid aliases thin surfaces; jittering removes the bias |
| 4. Hierarchical sampling | Efficiency | Importance-sampling weights spends the budget on surfaces, not empty space |
| 5. Differentiability | Why it learns | Pixel gradients flow to the visible surface; multi-view coverage supervises the whole field |
Connection to the broader theory¶
Volume rendering is the differentiable bridge between a continuous field and an image. Reading it as numerical integration explains everything practical about it: the weight distribution gives colour, depth, and uncertainty (Experiment 1); accuracy is bought with samples at ~$1/N$ (Experiment 2); stratified (Experiment 3) and hierarchical (Experiment 4) sampling fight aliasing and waste; and differentiability (Experiment 5) is what lets an image loss teach a 3D field. These are exactly the costs that motivate the next topics — positional encoding so a feasible number of samples still renders sharp detail, and Instant-NGP's occupancy grids and $O(1)$ lookups so empty space is skipped — and, ultimately, 3D Gaussian Splatting, which abandons ray marching for rasterized splatting.