Image losses and metrics

Every radiance-field paper ends with the same table: PSNR ↑, SSIM ↑, LPIPS ↓, averaged over held-out views. Every 3DGS trainer minimises the same loss, 0.8 · L1 + 0.2 · (1 − SSIM). These look like one-liners, and they are — but each has a convention (a log base, a data range, an averaging rule, a window, a padding mode, a colour encoding) that silently changes the number. A reconstruction engineer is expected to compute them exactly the way the reference code does, to know which differences are real, and to know why the field needed a learned metric on top. This lesson builds all of them in PyTorch, differentiably, in the (..., C, H, W) layout research code uses.

1. MSE and PSNR

The mean squared error between a render and the ground truth, and its logarithmic form, the peak signal-to-noise ratio:

MSE  = mean over C, H, W of (render − gt)²
PSNR = 10 · log10(max² / MSE)  =  20 · log10(max / RMSE)      [dB]

max is the peak of the data range: 1 for float images in [0, 1], 255 for 8-bit values. Decibels turn ratios into differences — dividing the MSE by 10 adds exactly 10 dB, halving the RMS error adds 6.02 dB — so "+0.5 dB" means the same relative improvement at 25 dB and at 35 dB. Rough anchors for novel view synthesis: below 20 dB is visibly wrong, 25–28 dB is typical for real-world captures (Mip-NeRF 360 scenes), 30–35 dB for synthetic scenes.

Three conventions to get right:

1. Per image, then average. 3DGS's psnr() reduces over everything but the batch axis (.view(img1.shape[0], -1).mean(1, keepdim=True)), and metrics.py averages the per-image values over the test set. The PSNR of the mean MSE is a different, smaller number, dominated by the worst image.
2. The data range matches the values. max_val = 255 with images in [0, 1] adds 48 dB of fiction.
3. Integers wrap. Images read back from PNGs are uint8; a - b in uint8 is modulo 256. Convert to float first.

import torch


def psnr(a, b, max_val=1.0):
    d = a.float() - b.float()                      # float first: uint8 would wrap
    return 10 * torch.log10(max_val ** 2 / (d * d).mean(dim=(-3, -2, -1)))


torch.manual_seed(0)
gt = torch.rand(2, 3, 16, 16)
noise = torch.tensor([0.1, 0.01]).view(2, 1, 1, 1)
render = gt + noise * torch.randn_like(gt)
p = psnr(render, gt)
print("per image:", [round(v, 2) for v in p.tolist()], "| mean:", round(p.mean().item(), 2))
mse = ((render - gt) ** 2).mean(dim=(-3, -2, -1))
print("PSNR of the mean MSE:", round((10 * torch.log10(1 / mse.mean())).item(), 2))

u8 = (gt * 255).round().to(torch.uint8)
off = (u8.int() + 3).clamp(0, 255).to(torch.uint8)   # 3 grey levels brighter
print("uint8 done right:", [round(v, 2) for v in psnr(off, u8, 255.0).tolist()])
wrapped = (u8 - off).float()                           # 0 - 3 = 253 in uint8
print("uint8 wrapped:   ", [round(v, 2) for v in (10 * torch.log10(255.0 ** 2 / (wrapped ** 2).mean(dim=(-3, -2, -1)))).tolist()])

2. SSIM: comparing local structure

PSNR treats every pixel independently. The structural similarity index (Wang, Bovik, Sheikh, Simoncelli 2004) compares neighbourhoods: around each pixel it computes Gaussian-weighted local means μ, variances σ² and the covariance σxy, and multiplies a luminance term and a contrast-structure term:

          2 μx μy + C1          2 σxy + C2
SSIM = ─────────────────── · ───────────────────        C1 = (0.01 L)², C2 = (0.03 L)²
        μx² + μy² + C1        σx² + σy² + C2

Each factor is 1 when the two neighbourhoods agree and falls towards 0 as they differ; the small constants keep flat, dark patches (μ ≈ 0, σ ≈ 0) from dividing by zero. L is the data range, so with images in [0, 1], C1 = 1e-4 and C2 = 9e-4. The local statistics come from five blurs with the same window: μx = blur(x), μy = blur(y), and the second moments blur(x²), blur(y²), blur(xy), from which the means' products are subtracted.

The window is the one from the paper: an 11 × 11 Gaussian with σ = 1.5, normalised to sum to 1 — the outer product of a 1-D Gaussian with itself.

In PyTorch the blur is F.conv2d with groups = C: a weight of shape (C, 1, 11, 11) is C independent one-channel filters, so the red channel is blurred only with red. (Without groups, conv2d wants a (C_out, C_in, k, k) weight and sums across input channels.) This is 3DGS's utils/loss_utils.py, almost line for line:

import torch
import torch.nn.functional as F


def gaussian(size, sigma, dtype=torch.float32):
    x = torch.arange(size, dtype=dtype) - size // 2
    g = torch.exp(-x ** 2 / (2 * sigma ** 2))
    return g / g.sum()


def ssim(x, y, size=11, sigma=1.5):                  # (N, C, H, W)
    C = x.shape[-3]
    g = gaussian(size, sigma, x.dtype)               # 3DGS: window.type_as(img1)
    win = torch.outer(g, g).expand(C, 1, size, size).contiguous()
    blur = lambda t: F.conv2d(t, win, padding=size // 2, groups=C)
    mx, my = blur(x), blur(y)
    vx = blur(x * x) - mx * mx
    vy = blur(y * y) - my * my
    cxy = blur(x * y) - mx * my
    c1, c2 = 0.01 ** 2, 0.03 ** 2
    m = (2 * mx * my + c1) * (2 * cxy + c2) / ((mx * mx + my * my + c1) * (vx + vy + c2))
    return m.mean()


torch.manual_seed(0)
gt = torch.rand(1, 3, 32, 32)
print("window sum:", round(gaussian(11, 1.5).sum().item(), 6), "| centre:", round(gaussian(11, 1.5)[5].item(), 4))
print("ssim(gt, gt)        =", round(ssim(gt, gt).item(), 6))
print("ssim(gt, gt + 0.05) =", round(ssim(gt, gt + 0.05).item(), 4), " (a uniform shift barely matters)")
print("ssim(gt, noisy)     =", round(ssim(gt, gt + 0.1 * torch.randn_like(gt)).item(), 4))
print("ssim(gt, 1 - gt)    =", round(ssim(gt, 1 - gt).item(), 4), " (inverted structure: negative)")

Implementation choices change the number, and papers mix them more than they should:

• Padding. 3DGS pads with zeros (padding = window_size // 2), so the output map is H × W and border pixels see a neighbourhood darkened by the padding. pytorch-msssim uses a valid convolution ((H − 10) × (W − 10) map); scikit-image crops the borders.
• Window. scikit-image's structural_similarity defaults to a 7 × 7 uniform window; it matches Wang et al. only with gaussian_weights=True, sigma=1.5, use_sample_covariance=False, data_range=1.
• Averaging. 3DGS's size_average=True averages over the batch; size_average=False gives one value per image.

When reproducing a paper's table, use its evaluation script.

3. The 3DGS training loss

L = (1 − λ) · L1 + λ · (1 − SSIM),      λ = 0.2  (lambda_dssim)

L1 rather than L2 for the per-pixel term: its gradient doesn't grow with the error, so a few badly wrong pixels (a floater, a specular highlight the model can't represent yet) don't dominate the update — even though the metric you report is PSNR, an L2 quantity. The D-SSIM term (1 − SSIM) adds a gradient that cares about local contrast and edges, which keeps textures sharp. SSIM is built from convolutions and arithmetic, so autograd differentiates it for free; the price is five blurs per step, which is why recent trainers swap in a fused CUDA SSIM with a hand-written backward (fused-ssim).

import torch
import torch.nn.functional as F


def ssim(x, y, size=11, sigma=1.5):
    C = x.shape[-3]
    r = torch.arange(size, dtype=x.dtype) - size // 2
    g = torch.exp(-r ** 2 / (2 * sigma ** 2))
    g = g / g.sum()
    win = torch.outer(g, g).expand(C, 1, size, size).contiguous()
    blur = lambda t: F.conv2d(t, win, padding=size // 2, groups=C)
    mx, my = blur(x), blur(y)
    vx, vy, cxy = blur(x * x) - mx ** 2, blur(y * y) - my ** 2, blur(x * y) - mx * my
    c1, c2 = 0.01 ** 2, 0.03 ** 2
    return ((2 * mx * my + c1) * (2 * cxy + c2) / ((mx ** 2 + my ** 2 + c1) * (vx + vy + c2))).mean()


torch.manual_seed(0)
gt = torch.rand(3, 24, 24)                          # unbatched (C, H, W), like train.py
render = torch.full_like(gt, 0.5).requires_grad_(True)
opt = torch.optim.Adam([render], lr=0.05)
for step in range(201):
    opt.zero_grad(set_to_none=True)
    l1 = (render - gt).abs().mean()
    s = ssim(render, gt)
    loss = 0.8 * l1 + 0.2 * (1 - s)
    loss.backward()
    opt.step()
    if step % 50 == 0:
        print(f"step {step:3}: L1 {l1.item():.4f}  SSIM {s.item():.4f}  loss {loss.item():.4f}")

(conv2d accepts an unbatched (C, H, W) input, which is how 3DGS calls it on a single rendered image.)

4. sRGB and linear light

Pixel values in a PNG or JPEG are sRGB-encoded: roughly linear^(1/2.2), with a short linear segment near black. The encoding spends more code values on dark tones, where vision is more sensitive. The exact transfer function (IEC 61966-2-1):

decode: linear = s / 12.92                      if s ≤ 0.04045
                 ((s + 0.055) / 1.055)^2.4      otherwise
encode: s      = 12.92 · l                      if l ≤ 0.0031308
                 1.055 · l^(1/2.4) − 0.055      otherwise

The thresholds belong to different variables: 0.0031308 = 0.04045 / 12.92.

Where it matters for radiance fields:

• Standard training doesn't convert. NeRF, 3DGS and nerfstudio train directly on the dataset's sRGB values; the model learns display-referred colour, and published metrics are computed on those sRGB values in [0, 1].
• Physics needs linear light. Light adds linearly, so averaging pixels (downsampling, blur, a pixel's footprint over an edge) or changing exposure should happen in linear space. RawNeRF trains on linear raw data with a loss weighted to mimic a tone-mapped comparison; HDR and relighting work stays linear until display.
• Blending encoded values darkens. The average of sRGB black and white is 0.5, which displays at 21 % of white's light; the physical average encodes to 0.735.

In PyTorch the piecewise curve is a torch.where, with one gradient trap: torch.where differentiates both branches, and the power branch has an infinite slope at 0. Its contribution is multiplied by zero where it isn't selected — and 0 · inf = NaN:

import torch


def linear_to_srgb(l):
    return torch.where(l <= 0.0031308, 12.92 * l, 1.055 * l ** (1 / 2.4) - 0.055)


def srgb_to_linear(s):
    return torch.where(s <= 0.04045, s / 12.92, ((s + 0.055) / 1.055) ** 2.4)


s = torch.linspace(0, 1, 11, dtype=torch.float64)
print("round trip error:", (linear_to_srgb(srgb_to_linear(s)) - s).abs().max().item())
print("mid-grey 0.5 decodes to", round(srgb_to_linear(torch.tensor(0.5)).item(), 4))

l = torch.tensor([0.0, 0.5], requires_grad=True)
linear_to_srgb(l).sum().backward()
print("gradient at 0 without care:", l.grad.tolist())

l = torch.tensor([0.0, 0.5], requires_grad=True)
safe = torch.where(l <= 0.0031308, 12.92 * l, 1.055 * l.clamp_min(0.0031308) ** (1 / 2.4) - 0.055)
safe.sum().backward()
print("with the power branch clamped:", [round(v, 4) for v in l.grad.tolist()])

5. Why LPIPS exists

PSNR and SSIM are both aligned comparisons: pixel (i, j) of the render against pixel (i, j) of the ground truth. Two things people barely notice cost them dearly — a one-pixel misregistration of sharp detail — and one thing people hate costs them little: blur, because a blurred prediction is the minimum-squared-error answer when detail is uncertain.

import torch
import torch.nn.functional as F


def ssim(x, y, size=11, sigma=1.5):
    C = x.shape[-3]
    r = torch.arange(size, dtype=x.dtype) - size // 2
    g = torch.exp(-r ** 2 / (2 * sigma ** 2))
    g = g / g.sum()
    win = torch.outer(g, g).expand(C, 1, size, size).contiguous()
    blur = lambda t: F.conv2d(t, win, padding=size // 2, groups=C)
    mx, my = blur(x), blur(y)
    vx, vy, cxy = blur(x * x) - mx ** 2, blur(y * y) - my ** 2, blur(x * y) - mx * my
    c1, c2 = 0.01 ** 2, 0.03 ** 2
    return ((2 * mx * my + c1) * (2 * cxy + c2) / ((mx ** 2 + my ** 2 + c1) * (vx + vy + c2))).mean()


def psnr(x, y):
    return 10 * torch.log10(1 / ((x - y) ** 2).mean())


torch.manual_seed(0)
blocks = (torch.rand(1, 1, 8, 8) > 0.5).float()
gt = blocks.repeat_interleave(4, -1).repeat_interleave(4, -2) * 0.8 + 0.1   # sharp texture
blurred = F.avg_pool2d(F.pad(gt, (2, 2, 2, 2), mode="replicate"), 5, stride=1)
shifted = torch.roll(gt, 1, dims=-1)                                        # one pixel to the right
for name, img in [("blurred", blurred), ("shifted", shifted)]:
    print(f"{name:8} PSNR {psnr(img, gt).item():5.2f} dB   SSIM {ssim(img, gt).item():.3f}")

The blurred image wins on PSNR although it has lost all its detail; the shifted one, perfect to a viewer, scores no better. LPIPS (Zhang et al. 2018, "The Unreasonable Effectiveness of Deep Features as a Perceptual Metric") compares images in the feature space of a pretrained network instead:

1. run both images (scaled to [−1, 1]) through AlexNet or VGG and take the activations of five layers;
2. unit-normalise each feature vector across channels, take the squared difference, weight channels with learned linear weights w_l;
3. average spatially and sum over layers.

The weights were fitted to human two-alternative judgements ("which of these two distortions is closer to the reference?"), and the resulting distance agrees with people far more often than PSNR or SSIM. It isn't available in this sandbox (it needs pretrained weights), but you'll see it in every evaluation script: torchmetrics' LearnedPerceptualImagePatchSimilarity (AlexNet by default) in nerfstudio, lpips(…, net_type="vgg") in 3DGS's metrics.py — and AlexNet and VGG LPIPS values are not comparable with each other, one more reason a table needs its evaluation code.

In an interview

Expect "define PSNR and say what +1 dB means" (a 21 % lower MSE), "how is SSIM computed?" (local Gaussian-weighted means, variances and covariance; luminance × contrast-structure; the constants), "why does 3DGS train with L1

• D-SSIM rather than MSE?" (robustness of L1, structure from SSIM, λ = 0.2), and "what's wrong with PSNR as a metric?" (aligned, rewards blur — hence LPIPS). Practical follow-ups: why your SSIM differs from a paper's (window, padding, data range, averaging), what goes wrong with uint8 images, and when you would work in linear rather than sRGB values. Be ready to write F.conv2d(img, window, padding=5, groups=C) from memory and explain groups.