Tone mapping (HDR to display range)

Deep dive

How it works

In: linear scene-referred RGB floats (any positive range), an exposure (manual, or key value a) and a curve · Out: display-encoded 8-bit (or PQ/HLG) RGB · State: none for stills; one smoothed average luminance for video · Cost: O(1) ALU per pixel plus one O(N) average-luminance reduction; about 10 ns/px scalar here (measured below)

1. Linearise the input (render output, or camera raw after white balance and a matrix to RGB primaries). All later steps assume linear light (colour conversion).
2. Exposure. Log-average Lw = exp(mean(ln(δ + L))), small δ avoids ln 0. Scale L' = (a / Lw) · L, key value a = 0.18 by default, 0.18 to 0.36 for high-key scenes (Reinhard et al. 2002). For video smooth Lw over time, or set exposure by hand in stops.
3. Curve Ld = f(L'), the display luminance:
   • Reinhard: Ld = L' / (1 + L'). Extended form Ld = L'(1 + L'/Lwhite²) / (1 + L'), where Lwhite is the smallest luminance mapped to pure white (default: the scene maximum).
   • Hable (Uncharted 2): f(x) = (x(Ax + CB) + DE) / (x(Ax + B) + DF) − E/F, with A..F = 0.15, 0.50, 0.10, 0.20, 0.02, 0.30 and white W = 11.2 as published, divided by f(W) so W maps to 1; he applies exposure bias 2 inside.
   • ACES: the full RRT + ODT is heavy. Narkowicz's fit saturate(x(2.51x + 0.03) / (x(2.43x + 0.59) + 0.14)) has max error 0.0138, pre-exposed (1 maps to about 0.8; multiply input by 0.6 for the original curve); it over-saturates brights.
4. Colour. Luminance-based: rgb · Ld/L' keeps ratios (hue) but a saturated bright pixel can exceed 1. Per-channel: apply f to R, G, B; always in range, but highlights desaturate and hues shift.
5. Gamut. Clip, or compress out-of-range values before encoding.
6. Encode. SDR sRGB OETF (IEC 61966-2-1). HDR: PQ ((c1 + c2 Y^m1)/(1 + c3 Y^m1))^m2 with Y = L/10000, m1 = 2610/16384, m2 = 2523/32, c1 = 107/128, c2 = 2413/128, c3 = 2392/128, or HLG (BT.2100).
7. Dither, then quantise. Add a threshold in (0,1) codes (ordered Bayer or noise) to the encoded value, then floor: banding becomes fine noise and sub-code darks survive on average.

Code

use std::time::Instant;

type Rgb = [f32; 3];
const EPS: f32 = 1e-4;

fn luminance(c: Rgb) -> f32 { 0.2126 * c[0] + 0.7152 * c[1] + 0.0722 * c[2] } // Rec.709 / sRGB primaries

/// Reinhard's log-average ("key") luminance: exp(mean(ln(eps + L))).
fn log_average(img: &[Rgb]) -> f32 {
    let s: f64 = img.iter().map(|&c| (EPS + luminance(c)).ln() as f64).sum();
    (s / img.len() as f64).exp() as f32
}

/// Reinhard extended: L (1 + L/Lw^2) / (1 + L); maps L = lw to exactly 1.
fn reinhard_ext(l: f32, lw: f32) -> f32 { l * (1.0 + l / (lw * lw)) / (1.0 + l) }

/// Hable "Uncharted 2" curve (A..F, W from his blog), exposure bias 2, scaled so W maps to 1.
fn hable(x: f32) -> f32 {
    let f = |x: f32| ((x * (0.15 * x + 0.10 * 0.50) + 0.20 * 0.02) / (x * (0.15 * x + 0.50) + 0.20 * 0.30)) - 0.02 / 0.30;
    f(2.0 * x) / f(11.2)
}

fn srgb_oetf(l: f32) -> f32 {
    if l <= 0.0031308 { 12.92 * l } else { 1.055 * l.powf(1.0 / 2.4) - 0.055 }
}

const BAYER4: [[u8; 4]; 4] = [[0, 8, 2, 10], [12, 4, 14, 6], [3, 11, 1, 9], [15, 7, 13, 5]];

/// Tone map on luminance, keep the RGB ratios: out = rgb * (curve(L) / L). Returns linear display RGB (clamped to 0..1).
fn tonemap(img: &[Rgb], exposure: f32, curve: &dyn Fn(f32) -> f32) -> Vec<Rgb> {
    img.iter().map(|&c| {
        let c = c.map(|v| v * exposure);
        let l = luminance(c);
        let k = if l > 0.0 { curve(l) / l } else { 0.0 };
        c.map(|v| (v * k).clamp(0.0, 1.0))
    }).collect()
}

/// sRGB-encode and quantise to 8 bit with a 4x4 ordered dither (threshold in (0,1) codes).
fn encode8(lin: &[Rgb], w: usize) -> Vec<u8> {
    let mut out = Vec::with_capacity(lin.len() * 3);
    for (i, c) in lin.iter().enumerate() {
        let t = (BAYER4[(i / w) % 4][(i % w) % 4] as f32 + 0.5) / 16.0;
        for &v in c { out.push((srgb_oetf(v) * 255.0 + t).floor().clamp(0.0, 255.0) as u8); }
    }
    out
}

/// Synthetic scene: luminance 1e-3..1e4 across x (log scale), orange tint, plus bright spots.
fn scene(w: usize, h: usize) -> Vec<Rgb> {
    let tint = [1.0f32, 0.5, 0.25];
    let n = luminance(tint);
    let mut img = Vec::with_capacity(w * h);
    for y in 0..h {
        for x in 0..w {
            let l = 1e-3f32 * 1e7f32.powf(x as f32 / (w - 1) as f32);
            let spot = (y % 64 == 32 && x % 97 == 0) as u32 as f32 * 5e4; // row 0 has no spots
            img.push(tint.map(|t| (l + spot) * t / n));
        }
    }
    img
}

fn main() {
    let (w, h) = (4096, 2048);
    let img = scene(w, h);
    let exposure = 0.18 / log_average(&img); // key value a = 0.18
    let lmax = img.iter().map(|&c| luminance(c)).fold(0.0, f32::max) * exposure;
    let reinhard = |l: f32| reinhard_ext(l, lmax);
    assert!((reinhard(lmax) - 1.0).abs() < 1e-5 && (hable(11.2 / 2.0) - 1.0).abs() < 1e-5);

    for (name, curve) in [("reinhard-ext", &reinhard as &dyn Fn(f32) -> f32), ("hable", &hable)] {
        let t = Instant::now();
        let lin = tonemap(&img, exposure, curve);
        let px = encode8(&lin, w);
        let mps = (w * h) as f64 / t.elapsed().as_secs_f64() / 1e6;
        // monotonic in input luminance along the (single-tint) ramp
        let row = &lin[..w];
        assert!(row.windows(2).all(|p| luminance(p[1]) >= luminance(p[0]) - 1e-6));
        // dark detail survives: dithering lifts sub-code values into nonzero codes
        assert!(px[..3 * w / 10].iter().any(|&v| v > 0));
        // colour ratio kept wherever no channel clipped
        for (a, b) in row.iter().zip(&img[..w]) {
            if a[0] < 1.0 && a[1] > 1e-6 { assert!(((a[0] / a[1]) / (b[0] / b[1]) - 1.0).abs() < 1e-3); }
        }
        println!("{name}: {mps:.0} Mpixel/s (1 thread, {w}x{h}, tonemap + encode + dither)");
    }
}

Asserts: white point maps to 1, monotone luminance, dark detail survives dithering, RGB ratio kept where nothing clips. The orange tint (red is 1.70 times luminance) clips red above Ld of about 0.59. Measured on this Mac (Apple M4 Pro), one thread, release build, 4096x2048: Reinhard-extended 88 to 95 and Hable 98 to 108 Mpixel/s over several runs, including allocation.

Complexity and constants

• Per pixel: luminance (3 FMA), one curve (a division, a few FMAs), 3 multiplies, 3 sRGB encodes (powf unless tabulated), 3 dither adds.
• Average luminance: one ln per pixel; 417 Mpixel/s here (scalar).
• Memory: f32 RGB is 12 B/px in, 3 B/px out; 8.4 Mpixel is 100 MB in. The exposure needs the whole frame: an extra pass or last frame's value.
• Local operators add one blur per scale (convolution, pyramid).

Efficiency

• Cost: measured on this Mac (Apple M4 Pro), scalar Rust, one thread: tone map + sRGB encode + dither about 9 to 11 ns/px (88 to 108 Mpixel/s); encode with dither alone 118 to 143 Mpixel/s; log-average 417 Mpixel/s. Three powf per pixel dominate: compute-bound scalar; with a table or on a GPU, memory-bandwidth bound.
• Levers:
  • OETF as a table — 10x less encode time (6.4 vs 0.61 ns/px, measured on colour conversion) — 8-bit SDR output
  • fuse exposure, curve and encode into the last pass or tile — skips a 12 B/px round trip — always
  • 3D LUT for a heavy view transform — one fetch instead of the transform; Bevy ships LUTs for AgX — GPU
  • estimate Lw from a subsample or small mip — fewer ln calls (not measured) — exposure need not be exact
  • filter a downsampled base layer — Durand and Dorsey report about two orders of magnitude from piecewise-linear approximation plus subsampling — bilateral local operator
• Naive → standard → state of the art:
  • scalar powf per channel, whole-frame passes · about 10 ns/px here
  • one fused shader, closed-form curve, hardware sRGB targets or LUT · bandwidth-bound (not measured)
  • fitted view transform as 3D LUT with gamut compression, PQ/HLG output · one LUT fetch per pixel

Implementation notes

• Keep the working buffer float (f16 on GPUs); 8-bit linear bands in the darks.
• Luminance weights follow the working primaries (Rec.709 here); drop NaN and Inf from the log-average.
• PQ is absolute, so there is no exposure: map SDR white to a chosen nit level first.
• Check the curve is monotone and ends at 1, as the asserts do.

Pitfalls

• Curve applied to sRGB-encoded values, or twice (engine plus OS or TV HDR mapping).
• Per-channel mapping of saturated lights: orange fire goes yellow-white.
• A few extreme pixels (sun, fireflies from Monte Carlo noise) drag the average and darken the frame: clamp or use a percentile.
• Local operators: halos at strong edges when the base/detail split or scale fails; Reinhard et al. choose per pixel the largest neighbourhood without strong contrast.

Variants and alternatives

In production

• Real-time rendering: Bevy ships Reinhard, ACES-fitted, AgX and others; Narkowicz moved a game renderer from Reinhard to a filmic curve.
• Film and VFX: OpenColorIO with ACES transforms. Photography: OpenCV Reinhard, Drago, Mantiuk and Mertens fusion.
• Medical: DICOM Window Center and Width. Astronomy: Astropy AsinhStretch. Video players: libplacebo.

Implementations

Sources

• Reinhard, Stark, Shirley, Ferwerda, Photographic Tone Reproduction for Digital Images, ACM TOG 21(3), SIGGRAPH 2002. DOI · tech report
• Durand, Dorsey, Fast Bilateral Filtering for the Display of High-Dynamic-Range Images, ACM TOG 21(3), SIGGRAPH 2002. DOI · PDF
• Mertens, Kautz, Van Reeth, Exposure Fusion, Pacific Graphics 2007. DOI · PDF
• Hable, Filmic Tonemapping Operators, Filmic Worlds blog. post
• Narkowicz, ACES Filmic Tone Mapping Curve, 2016. post
• ITU-R Recommendation BT.2100 (PQ and HLG for HDR television). rec; SMPTE ST 2084:2014 (PQ). DOI
• Lupton et al., Preparing Red-Green-Blue Images from CCD Data, PASP 116, 2004. DOI
• DICOM PS3.3 C.11.2 (Window Center, Window Width). standard · Astropy, Image stretching and normalization. docs
• Sobotka, AgX. repo · Khronos, PBR Neutral. repo