Monte Carlo denoising (feature-guided filtering)

Deep dive

How it works

Parts:

Feature buffers — albedo, normal, depth, motion vectors: noise-free first-hit data.
Demodulator — divides colour by albedo, multiplies back at the end.
History — last frame's irradiance and luminance moments, reprojected.
Variance estimator — per-pixel luminance variance from the moments.
A-trous filter — 5x5 B3-spline kernel [1/16, 1/4, 3/8, 1/4, 1/16], taps 2^i pixels apart in pass i.
Edge-stopping weights — product of normal, depth and luminance terms in [0, 1].

Render 1 to a few samples per pixel plus feature buffers; clamp extreme samples (fireflies).
Demodulate: irr = colour / albedo.
Reproject last frame's history, rejecting it where depth, normal or mesh id disagree; blend h = a*irr + (1-a)*h (a = 0.2 in SVGF).
Variance = second moment minus squared first moment of luminance; with short history use a spatial estimate.
For pass i = 0..k-1: out(p) = sum h(q) w(p,q) in(q) / sum h(q) w(p,q) over 25 taps at spacing 2^i; w = w_normal * w_depth * w_lum. SVGF scales w_lum by the variance; Dammertz et al. halve its sigma per pass.
Remodulate: colour = irr * albedo; tone-map afterwards.

noisy colour --> / albedo --> temporal accumulate --> variance --> a-trous x k --> * albedo --> out
 (1-4 spp)                        ^ reproject          (sigma_l)    steps 1,2,4,8,..
albedo, normal, depth, motion ----+----------------------+------------> edge-stopping weights w(p,q)

The code covers steps 2, 5, 6 (no history; sigma halves per pass instead of using variance).

Code

Synthetic scene: four surfaces, checker texture, light gradient, and a shadow the normals cannot see; 4 heavy-tailed samples per pixel.

use std::time::Instant;

type V3 = [f32; 3];
const NORMALS: [V3; 4] = [[0., 0., 1.], [0.6, 0., 0.8], [0., 0.6, 0.8], [-0.6, 0., 0.8]];
const IRR: [f32; 5] = [1.0, 0.25, 2.0, 0.5, 0.3]; // per surface; 4 = shadow on 2
const TINT: [V3; 4] = [[1., 0.6, 0.3], [0.3, 0.8, 0.4], [0.5, 0.5, 1.], [1., 1., 0.7]];

struct Frame { w: usize, h: usize, albedo: Vec<V3>, normal: Vec<V3>, truth: Vec<V3>, noisy: Vec<V3>, edge: Vec<bool> }

fn surface(x: usize, y: usize, w: usize, h: usize) -> usize {
    let (dx, dy) = (x as f32 - 0.3 * w as f32, y as f32 - 0.4 * h as f32);
    if dx * dx + dy * dy < (0.2 * h as f32).powi(2) { 3 } else if y * 4 > h * 3 { if x * 5 > w * 2 && x * 10 < w * 7 { 4 } else { 2 } } else if x * 2 > w { 1 } else { 0 }
}

fn make_frame(w: usize, h: usize, spp: u32) -> Frame {
    let mut rng = 0x9E3779B97F4A7C15u64;
    let mut uni = move || { rng ^= rng << 13; rng ^= rng >> 7; rng ^= rng << 17; ((rng >> 40) as f32 + 0.5) / 16777216.0 };
    let n = w * h;
    let mut f = Frame { w, h, albedo: vec![[0.; 3]; n], normal: vec![[0.; 3]; n], truth: vec![[0.; 3]; n], noisy: vec![[0.; 3]; n], edge: vec![false; n] };
    for y in 0..h { for x in 0..w {
        let (i, s) = (y * w + x, surface(x, y, w, h));
        let checker = if ((x / 8) + (y / 8)) % 2 == 0 { 0.9 } else { 0.25 }; // texture detail
        let e = IRR[s] * (0.6 + 0.4 * x as f32 / w as f32);                   // smooth lighting
        let m = if s == 4 { 2 } else { s }; // the shadow shares surface 2's normal and tint
        f.albedo[i] = TINT[m].map(|t| t * checker);
        f.normal[i] = NORMALS[m];
        let nb = [(x.saturating_sub(3), y), ((x + 3).min(w - 1), y), (x, y.saturating_sub(3)), (x, (y + 3).min(h - 1))];
        f.edge[i] = nb.iter().any(|&(a, b)| surface(a, b, w, h) != s); // within 3 px of a surface or shadow edge
        for c in 0..3 {
            f.truth[i][c] = f.albedo[i][c] * e;
            // heavy-tailed Monte Carlo samples: Pareto(alpha 3) scaled to mean 1, clamped at 8x the mean
            let m: f32 = (0..spp).map(|_| ((1.0 - uni()).powf(-1.0 / 3.0) * (2.0 / 3.0)).min(8.0)).sum::<f32>() / spp as f32;
            f.noisy[i][c] = f.albedo[i][c] * e * m;
        }
    }}
    f
}

fn lum(c: V3) -> f32 { 0.2126 * c[0] + 0.7152 * c[1] + 0.0722 * c[2] }

/// One pass of the 5x5 B3-spline a-trous filter with taps `step` pixels apart.
fn atrous_pass(src: &[V3], dst: &mut [V3], normal: &[V3], w: usize, h: usize, step: i32, sigma_c: f32, edge_aware: bool) {
    const K: [f32; 5] = [1. / 16., 1. / 4., 3. / 8., 1. / 4., 1. / 16.];
    for y in 0..h as i32 { for x in 0..w as i32 {
        let p = y as usize * w + x as usize;
        let (mut acc, mut wsum, lp) = ([0f32; 3], 0f32, lum(src[p]));
        for j in -2..=2i32 { for i in -2..=2i32 {
            let (qx, qy) = ((x + i * step).clamp(0, w as i32 - 1), (y + j * step).clamp(0, h as i32 - 1));
            let q = qy as usize * w + qx as usize;
            let mut wt = K[(i + 2) as usize] * K[(j + 2) as usize];
            if edge_aware {
                let nd = normal[p].iter().zip(&normal[q]).map(|(a, b)| a * b).sum::<f32>().max(0.0);
                wt *= nd.powi(64) * (-(lp - lum(src[q])).abs() / sigma_c).exp(); // normal and colour edge-stopping
            }
            for c in 0..3 { acc[c] += wt * src[q][c]; }
            wsum += wt;
        }}
        dst[p] = acc.map(|a| a / wsum);
    }}
}

/// Demodulate albedo, run `iters` passes (step 1,2,4,...; sigma_c halves each pass), remodulate.
fn denoise(f: &Frame, iters: u32, edge_aware: bool, demodulate: bool) -> Vec<V3> {
    let alb = |i: usize, c: usize| if demodulate { f.albedo[i][c].max(1e-3) } else { 1.0 };
    let mut a: Vec<V3> = (0..f.w * f.h).map(|i| [0, 1, 2].map(|c| f.noisy[i][c] / alb(i, c))).collect();
    let mut b = a.clone();
    let mut sigma = 1.0;
    for k in 0..iters {
        atrous_pass(&a, &mut b, &f.normal, f.w, f.h, 1 << k, sigma, edge_aware);
        std::mem::swap(&mut a, &mut b);
        sigma *= 0.5;
    }
    (0..f.w * f.h).map(|i| [0, 1, 2].map(|c| a[i][c] * alb(i, c))).collect()
}

fn rmse(f: &Frame, img: &[V3], edges_only: bool) -> f32 {
    let (mut s, mut n) = (0f32, 0f32);
    for i in 0..f.w * f.h { if !edges_only || f.edge[i] {
        for c in 0..3 { s += (img[i][c] - f.truth[i][c]).powi(2); n += 1.0; } } }
    (s / n).sqrt()
}

fn main() {
    let f = make_frame(256, 256, 4);
    let (before, before_e) = (rmse(&f, &f.noisy, false), rmse(&f, &f.noisy, true));
    let ours = denoise(&f, 5, true, true);
    let blur = denoise(&f, 5, false, true);      // same footprint, no edge-stopping
    let raw = denoise(&f, 5, true, false);       // edge-aware but no albedo demodulation
    println!("RMSE all/edge: noisy {before:.3}/{before_e:.3}  a-trous {:.3}/{:.3}  blur {:.3}/{:.3}  no-demod {:.3}/{:.3}",
        rmse(&f, &ours, false), rmse(&f, &ours, true), rmse(&f, &blur, false), rmse(&f, &blur, true), rmse(&f, &raw, false), rmse(&f, &raw, true));
    assert!(rmse(&f, &ours, false) < before / 8.0);                // solid reduction overall
    assert!(rmse(&f, &ours, true) < before_e / 4.0);               // edges stay sharp
    assert!(rmse(&f, &blur, true) > 10.0 * rmse(&f, &ours, true));   // same footprint, no edge-stopping smears them
    assert!(rmse(&f, &raw, false) > 10.0 * rmse(&f, &ours, false));  // texture detail needs demodulation
    let big = make_frame(1920, 1080, 4);
    let t = Instant::now();
    let out = denoise(&big, 5, true, true);
    println!("1080p, 5 passes, 1 thread: {:?} (rmse {:.3} vs noisy {:.3})", t.elapsed(), rmse(&big, &out, false), rmse(&big, &big.noisy, false));
}

Measured on this Mac (Apple M4 Pro, release build, one thread, 256x256 frame). RMSE over the frame: noisy 0.092, filtered 0.007, no edge-stopping 0.091, no demodulation 0.153. Near edges: 0.096, 0.012, 0.185, 0.164. 1080p, 5 passes: 739 to 885 ms over nine runs.

Complexity and constants

Time O(k * 25 * N); a dense filter of the same reach (125x125 at k = 5) needs 15,625 taps.
Cost drivers: taps, passes, and the exp and power per tap.

Efficiency

Cost: one 1080p frame, five passes: 259 million taps in about 0.75 to 0.9 s, 3 ns per tap (measured on this Mac, M4 Pro, one scalar thread; division mine). Bound: compute (a power and an exp per tap); each buffer is about 25 MB.

Levers:

Run on a GPU — 10.2 ms per 1080p frame including temporal and variance steps (Schied et al. 2017, Titan X Pascal) — real time.
Fewer weight buffers — 42.6 ms with three, 9.6 ms with none, 1080p, 5 passes (Dammertz et al. 2010, GTX 285) — when edges need no colour or position.
Fewer passes — time is linear in passes (Dammertz et al. 2010) — higher sample counts.

Naive → standard → state of the art:

More samples · four times the samples, half the noise.
Edge-avoiding a-trous, normal, position and colour weights · 42.6 ms at 1080p (Dammertz et al. 2010, GTX 285).
SVGF · 10.2 ms at 1080p from 1 sample per pixel; learned denoisers for final frames (Bako et al. 2017; Open Image Denoise).

Implementation notes

Filter linear HDR irradiance, before tone mapping; floor albedo.
Normal weight: a high power of the dot product (code: 64; SVGF: 128). SVGF's depth weight uses the screen-space depth gradient, so slanted planes are not cut.

Pitfalls

Fireflies: an extreme pixel differs from all neighbours, so a colour weight leaves it alone. Clamp samples.
Temporal lag: a small a is stable but lags; fast geometry near a light leaves detached shadows (SVGF).
Luminance-only variance over-blurs coloured reflections (SVGF).
Demodulation fits diffuse surfaces only (Dammertz et al.).
Bias: the result is a blur, not a sample mean.

Variants and alternatives

Regression on feature bufferslocal fits (Bitterli et al. 2016); higher-order fits overfit noise (Bako et al. 2017).

Kernel-predicting networka CNN outputs a weight kernel per pixel (Bako et al. 2017).

Recurrent autoencodera network with memory over frames, for interactive sequences (Chaitanya et al. 2017).

Joint bilateral upsamplingguided weights used to upsample (Kopf et al. 2007).

Adaptive samplingsend more samples where noise is high (survey: Zwicker et al. 2015).

In production

Blender Cycles offers Open Image Denoise and OptiX denoisers (Blender manual).
Open Image Denoise: learned filters from 1 sample per pixel to converged, optional albedo and normal inputs (README).
NVIDIA Real-Time Denoisers: RELAX (a-trous based) and ReBLUR, in 15+ AAA titles; RELAX 3.25 ms at 1440p, RTX 4080 (README).
OptiX 8 ships an AI denoiser based on the recurrent autoencoder (NVIDIA).

Implementations

name	lang	repo	install	verdict	licence	health	as_of	notes
Intel Open Image Denoise	C++	RenderKit/oidn	`brew install open-image-denoise`	—	—	—	—	learned filters, CPU and GPU
NVIDIA Real-Time Denoisers (NRD)	C++	NVIDIA-RTX/NRD	`git clone https://github.com/NVIDIA-RTX/NRD`	—	—	—	—	RELAX (a-trous), ReBLUR, SIGMA
Falcor SVGF pass	C++	NVIDIAGameWorks/Falcor	`git clone https://github.com/NVIDIAGameWorks/Falcor`	—	—	—	—	`Source/RenderPasses/SVGFPass`
Cycles	C++	blender/cycles	`git clone https://github.com/blender/cycles`	—	—	—	—	Open Image Denoise and OptiX integration

Sources

Dammertz, Sewtz, Hanika, Lensch, Edge-Avoiding A-Trous Wavelet Transform for Fast Global Illumination Filtering, HPG 2010. paper · pdf
Schied et al., Spatiotemporal Variance-Guided Filtering: Real-Time Reconstruction for Path-Traced Global Illumination, HPG 2017. paper · page
Bako et al., Kernel-Predicting Convolutional Networks for Denoising Monte Carlo Renderings, ACM TOG 36(4), 2017. paper
Chaitanya et al., Interactive Reconstruction of Monte Carlo Image Sequences using a Recurrent Denoising Autoencoder, ACM TOG 36(4), 2017. paper
Bitterli et al., Nonlinearly Weighted First-order Regression for Denoising Monte Carlo Renderings, Computer Graphics Forum 2016. paper
Zwicker et al., Recent Advances in Adaptive Sampling and Reconstruction for Monte Carlo Rendering, Computer Graphics Forum 2015. paper
Kopf, Cohen, Lischinski, Uyttendaele, Joint Bilateral Upsampling, ACM TOG 26(3), 2007. paper
Hasinoff et al., Burst Photography for High Dynamic Range and Low-Light Imaging on Mobile Cameras, ACM TOG 35(6), 2016. paper
Open Image Denoise · OptiX denoiser · Cycles sampling manual