Burst Photography

Burst photography takes several frames for one press of the shutter and computes a single picture from them. The frames are short and equally exposed; one of them is chosen as the reference, the others are aligned to it, and a merge combines them so that noise falls where the frames agree and the reference is kept where they do not. The result is an image with less noise and more bits per pixel than any frame, which a camera pipeline then finishes as if it were one exposure. The technique is also called burst denoising, multi-frame denoising or align-and-merge; Google’s HDR+ system, described in detail by Hasinoff et al. (2016), is the example followed here, and an independent reimplementation by Monod, Delon and Veit (2021) analyses it step by step.

Computational photography introduces burst capture alongside the field’s other methods, with Night Sight and a worked align-and-merge example based on phase correlation and pixel-wise rejection; this article goes into the HDR+ alignment and merge in detail, the choices they leave open, and burst super-resolution and learned burst denoising. Alignment as a general problem belongs to image registration.

Purpose

A small camera collects few photons per pixel, and a single exposure cannot fix that: a longer exposure blurs with hand shake and subject motion, and more gain amplifies noise. Hasinoff et al. (2016, Section 1) list the other ways to gather more light (a larger aperture, optical stabilization, exposure bracketing, flash) and the cost of each on a phone. A burst replaces one long exposure by several short ones. If $N$ frames of a static scene are aligned exactly and their noise is independent, their mean has $1/N$ of the noise variance of one frame, as computational photography derives; the work is in making the frames align, and in deciding what to do where they do not.

HDR+ uses the same idea against limited dynamic range. It deliberately underexposes every frame so that highlights do not clip, and merges the frames to recover the shadows, so that the authors treat high-dynamic-range imaging as denoising (Hasinoff et al., 2016, Section 3). The merged raw image has more bits than the input: 10-bit raw frames are merged to 14 bits in the authors’ corrected copy of the paper (Section 6), which an erratum of July 2018 corrects from the 12 bits stated in the original version. That deeper image is then tone mapped to bring the shadows up.

When to Use It

Burst merging helps whenever a single frame is noise-limited and the scene and camera hold still long enough, or move simply enough, for the frames to be aligned: low light, high-contrast scenes captured with underexposure, and hand-held shots where a long exposure would blur. Hasinoff et al. (2016, Section 3) describe their burst size as a trade-off: more frames in low light or in very high-contrast scenes, one or two in bright scenes, and in practice 2–8 frames chosen with their noise model.

It does not help where the frames cannot be related to the reference. Occlusions, non-rigid motion and changes of lighting defeat alignment (Section 5), and a robust merge then falls back to the reference frame, which leaves that region as noisy as one frame. If the scene moves so fast that features blur despite the short exposure, alignment may fail and leave excessive noise in the photograph (Hasinoff et al., 2016, Section 7). Exposure bracketing, which varies the exposure between frames to extend the range, is the subject of high dynamic range imaging; HDR+ avoids it because frames of different exposure are harder to align (Hasinoff et al., 2016, Section 1).

How It Works

Capture

HDR+ captures Bayer raw frames of constant exposure rather than the processed output of the camera’s image signal processor. Raw frames have more bits, are linear after black-level subtraction, and have a simple noise model; processing the raw data also avoids the ISP’s own tone mapping and spatial denoising, and the expensive demosaicking step runs once on the merged image rather than on every frame (Hasinoff et al., 2016, abstract and Section 2).

Exposure is chosen by example. A database of about 5,000 scenes captured with bracketing was hand-tuned for the best rendition; each new viewfinder frame is described by four sets of 64 quantiles of its brightness distribution (two spatial scales, the maximum and the average of the RGB channels), and the labels of the closest scenes are blended into a short exposure for capture and a longer synthetic exposure for tone mapping (Section 3 and supplement Section 5). The ratio between them, the dynamic-range compression, is capped at 8; in the authors’ verification set of 26,071 scenes, 4% would have needed more (Fig. 4). The exposure is then factored into time and gain on a fixed schedule: up to 8 ms at minimum gain, then gain up to 4×, then both together up to 100 ms and 96× gain (Section 3).

Noise model

On raw data the noise of each pixel is independent and its variance is affine in the signal $x$,

$$\sigma^2(x) = A x + B,$$

the photon (shot) noise contributing the term proportional to $x$ and read noise the constant; Hasinoff et al. (2016, Section 5) state that $A$ and $B$ depend only on the analog and digital gain, which their system controls. Monod et al. (2021, Section 4.1) write the same model as $\lambda_s x + \lambda_r$ and, for raw files without the gain settings, take the parameters from the DNG NoiseProfile tag or scale ISO 100 values by the ISO ratio. Image noise derives the model.

Reference frame

The merge does not average frames symmetrically: it denoises one reference frame with the help of the others. HDR+ takes the sharpest of the first three frames, measured by gradients in the green channel, a form of lucky imaging; limiting the choice to the first three keeps the perceived shutter lag short (Hasinoff et al., 2016, Section 4). Monod et al. (2021, Section 3.1) did not reimplement this step and use the reference recorded with each burst of the HDR+ dataset, or for other bursts a frame the user selects, by default the first.

Alignment

Alignment finds, for every tile of the reference, the matching tile in each alternate frame. HDR+ estimates only translations per tile and only to whole pixels of a downsampled image, which has a convenient consequence for Bayer data: each 2 × 2 RGGB block is averaged to one grey value, so the alignment works on a quarter-size grey image, and displacements are multiples of 2 raw pixels, so a red sample always lands on a red sample (Section 4). Undersampling of the colour planes is left to the merge, where aliasing is treated like any other misalignment.

The search runs coarse to fine on a four-level Gaussian image pyramid. In the supplement’s typical configuration (Section 2), the levels are downsampled by 2, 4 and 4 from the grey image; the coarsest level uses 8 × 8 tiles and the others 16 × 16; the three coarse levels search ±4 pixels minimizing the sum of squared differences ($L_2$), and the finest level searches ±1 pixel minimizing the sum of absolute differences ($L_1$). For a reference tile $T$ and a search area $I$ of the alternate frame, the squared-difference surface over all offsets is computed as

$$D_2 = \lVert T \rVert_2^2 + \operatorname{box}(I \circ I, n) - 2\, \mathcal{F}^{-1}\{ \mathcal{F}\{I\}^{*} \circ \mathcal{F}\{T\} \},$$

the tile’s energy plus an unnormalized $n \times n$ box filter of the squared search area minus twice the cross-correlation, computed with FFTs (Eq. 2, after Lewis’s fast normalized cross-correlation). At the coarse levels a quadratic $\tfrac12 \mathbf{u}^\top \mathbf{A} \mathbf{u} + \mathbf{b}^\top \mathbf{u} + c$ is fitted by weighted least squares to the $3 \times 3$ neighbourhood of the minimum, which reduces to six $3 \times 3$ filters, and its minimum $\boldsymbol{\mu} = -\mathbf{A}^{-1}\mathbf{b}$ gives a sub-pixel offset (Eqs. 3–4; the supplement, Section 4, derives the filters and forces $\mathbf{A}$ to be positive semi-definite). When a fine tile’s guess is inherited from the coarser level, the offsets of the three nearest coarse tiles are tried and the one with the smallest $L_1$ residual is kept, which helps where a coarse tile straddles a motion boundary (Section 4).

The authors report 24 ms per megapixel for alignment on their phone and displacements of up to 169 pixels (Section 4, Fig. 5); these figures are specific to their 2016 devices. Monod et al. (2021, Section 3.3) found that the sub-pixel step brought very few visual improvements in their implementation and do not apply it at the finest level, where sub-pixel offsets would mix Bayer colour planes. Alignment is allowed to be imperfect because the merge is built to tolerate it.

Merge

The merge works on tiles in the spatial frequency domain, separately for each Bayer colour plane: typically 16 × 16 tiles, or 32 × 32 in very dark scenes where low-frequency noise is objectionable (Section 5). Let $T_z(\boldsymbol{\omega})$ be the 2-D DFT of the tile of frame $z$, with $z = 0$ the reference and $N$ frames in all. Averaging every frequency, $\tilde T_0 = \frac1N \sum_z T_z$, is the temporal mean, which ghosts wherever alignment fails (Eq. 5, Fig. 6). HDR+ instead interpolates per frequency between each alternate and the reference:

$$\tilde T_0(\boldsymbol{\omega}) = \frac{1}{N} \sum_{z=0}^{N-1} \Big( T_z(\boldsymbol{\omega}) + A_z(\boldsymbol{\omega}) \big[ T_0(\boldsymbol{\omega}) - T_z(\boldsymbol{\omega}) \big] \Big), \qquad A_z(\boldsymbol{\omega}) = \frac{|D_z(\boldsymbol{\omega})|^2}{|D_z(\boldsymbol{\omega})|^2 + c\,\sigma^2},$$

with $D_z = T_0 - T_z$ (Eqs. 6–7). Where the difference between a frame and the reference is of the size of the noise, $A_z$ is small and the alternate contributes; where it is much larger, $A_z$ approaches 1 and the reference is used instead. Because the decision is made per frequency, a partly misaligned tile still contributes the frequencies that match. $A_z$ is a shrinkage operator of the Wiener filter’s form; the authors tried hard and soft thresholding and found this one the best balance between noise reduction and artifacts (Section 5).

Three details make the formula concrete:

• Noise per tile. The DFT cannot carry a signal-dependent covariance cheaply, so each tile is given one noise variance, the model evaluated at the root-mean-square of the tile’s samples; using the RMS rather than the mean filters high-contrast tiles more strongly (Section 5). Monod et al. (2021, Section 4.1) evaluate it on the reference tile.
• The constant $c$. It scales the spatial noise variance to the variance of $D_z$’s coefficients: a factor $n^2$ for the number of DFT samples, $1/4^2$ for the window and 2 because $D_z$ is a difference of two tiles, times a tuning factor that the authors fixed to 8 and that, they write, increases noise reduction at the expense of some robustness (Section 5). Monod et al. (2021, Section 4.2) write $c = k\tau$, call the justification of $k$ questionable, and keep it because $\tau$ can absorb it; as $\tau \to 0$ the merge returns the reference, and as $\tau \to \infty$ it returns the temporal mean. In their comparison with Google’s results they used $\tau = 75$, as they could not reproduce the paper’s denoising strength with $\tau = 8$ and their own noise parameters (Section 6.2.1).
• Window and overlap. Tiles overlap by half in each direction and are multiplied by the modified raised-cosine window $w(x) = \tfrac12 - \tfrac12 \cos\big(2\pi (x + \tfrac12)/n\big)$ for $0 \le x < n$; with denominator $n$ rather than $n - 1$, windows repeated at an offset of $n/2$ sum to one everywhere, and the half-sample shift avoids zeros (Section 5). Monod et al. (2021, Section 4.4) describe it as the blending window for the overlapped tiles, which also removes DFT artifacts at tile edges; the example below applies it once, before the DFT, and overlap-adds the merged tiles without a second window.

The paper explains why it merges pairwise rather than taking a DFT across frames as well, as earlier video denoisers did (Kokaram’s 3-D Wiener filtering, V-BM3D’s shrinkage). In a 3-D transform, one badly aligned frame spreads energy across all temporal frequencies, so the shrinkage rejects every alternate frame, not only the bad one; the pairwise filter judges each frame separately. The authors note that the pairwise filter gives up some noise reduction for well-aligned frames, in exchange for being cheaper and degrading more gracefully (Section 5, Fig. 7, a one-pixel toy example with $\sigma = 4$).

Because the temporal filter does no spatial filtering, HDR+ adds a spatial step in the 2-D DFT: the same shrinkage applied to the merged tile with the noise variance reduced to $\sigma^2/N$, as if all frames had been averaged perfectly, which the authors call conservative, and a piecewise linear “noise shaping” function of frequency that raises the assumed noise at high frequencies, tuned for subjective quality (Section 5). The function’s slopes are not published; Monod et al. (2021, Section 4.3) replace it with one proportional to $|\boldsymbol{\omega}|$.

Artifacts and finishing

Hasinoff et al. (2016, Section 5) report three kinds of artifact: noise left near strong high-contrast features, whose spectra are not sparse, so spatial shrinkage fails there (Fig. 8); mild ghosts, because the shrinkage never rejects a misaligned tile completely (Fig. 9); and occasional ringing, most often near poorly aligned clipped highlights. The merged raw image then passes through a software pipeline of thirteen steps, from black-level subtraction, lens shading, white balance and demosaicking through chroma denoising, colour correction, local tone mapping by a two-exposure variant of exposure fusion, dehazing, a global tone curve, chromatic aberration correction, sharpening, hue adjustments and dithering (Section 6). On a Qualcomm Snapdragon 810 the whole pipeline took 2.5 to 4 seconds per photograph depending on the number of frames (Section 7).

How to Apply It

The example implements the HDR+ merge and alignment on one image plane and compares the merge with the temporal mean and with 3-D DFT shrinkage, where the truth is known everywhere.

Data. The scene is scikit-image’s bundled camera image scaled to $[0, 1]$, taken as one colour plane of a raw image (there is no Bayer mosaic, so the 2 × 2 averaging step does not apply). A burst of $N = 8$ frames of $448 \times 448$ pixels is cut from it at eight fixed camera displacements of up to 10 pixels; frame 0 is the reference, and its noise-free crop is the ground truth. Two departures from a rigid scene were fixed when the example was planned:

• a misaligned block: in frame 5 only, a $128 \times 128$ region of buildings and horizon is displaced 6 pixels to the right, a local motion that the camera displacement does not describe;
• a moving object: a white disc of radius 12 pixels on the grass moves 4 pixels to the right per frame.

Each frame receives Poisson–Gaussian noise with variance $Ax + B$, $A = 0.01$ and $B = 0.0004$, for seeds 0–4. The error is the root-mean-square error (RMSE) against the truth in three regions: static (at least 32 pixels from the border and 16 pixels from the block and from every disc position), misaligned (the block) and ghost (the trail, the pixels covered by the disc in some alternate frame but not in the reference).

Methods. All merges use 16 × 16 tiles with half overlap and the window above, applied once before the DFT; the noise variance of each tile is the model evaluated at the RMS of the reference tile. The pairwise merge follows Eqs. 6–7 with $c = n^2 \cdot \tfrac{1}{4^2} \cdot 2 \cdot \tau$. The 3-D DFT shrinkage applies the same form, $|X|^2 / (|X|^2 + c_3 \sigma^2)$, to every coefficient of the DFT across tiles and frames and keeps frame 0; its constant $c_3 = n^2 \cdot \tfrac{1}{4^2} \cdot N \cdot \tau$, with $N$ in place of the factor 2 for an $N$-point temporal DFT, is our own analogue, not a setting taken from Hasinoff et al.’s Fig. 7. The tuning factor takes three values fixed in advance: 1, the paper’s 8, and Monod et al.’s 75. HDR+’s spatial denoising step is left out of every method, because its noise-shaping function is not published. The temporal mean and the reference frame pass through the same tiles.

Alignment. Every merge runs twice: with the true camera displacements for every tile (“global”, an oracle that is wrong by design in the misaligned block and at the disc), and with offsets from our implementation of HDR+’s alignment (“HDR+ tiles”): the four-level pyramid with steps 2, 4 and 4, tiles of 8 and 16, the FFT form of Eq. 2 with ±4 search at the three coarse levels, the supplement’s sub-pixel fit there, the three-candidate upsampling, and a ±1 $L_1$ search at the finest level, whose tiles are the merge tiles.

Checks and predictions. The plan fixed seven checks before the first run: the overlapped windows sum to one (C1); on a noise-free burst without the block or the moving disc, the pairwise merge, the mean and the 3-D shrinkage with zero noise variance reproduce the truth (C2); $\tau = 0$ returns the reference and a huge $\tau$ the mean (C3); Eq. 2 by FFT equals the brute-force sum (C4); the alignment on noise-free frames is exact on at least 92% of static tiles and 95% of the displaced block’s tiles in frame 5 (C5); the generated noise has variance $Ax + B$ within 5% (C6); and the mean’s static RMSE is within 1% of $\sqrt{\overline{(Ax + B)}/N}$ (C7). It also fixed eight predictions, H1–H8, stated with the results below; each comparison is printed as its range of paired differences over the five seeds and holds only if its sign is the same for every seed.

import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[_v] = "1"                      # one BLAS thread, set before NumPy loads
import numpy as np
from numpy.lib.stride_tricks import sliding_window_view as windows
from scipy import ndimage
from skimage import data

# ---------------- set-up (fixed at planning) ----------------
N, n, HALF = 8, 16, 8              # frames, merge tile size, half overlap
A_NOISE, B_NOISE = 0.01, 0.0004    # noise variance A*x + B (signal x in [0, 1])
H = W = 448                        # reference frame size; scene is the 512x512 camera image
OFF = 32                           # reference frame = scene[32:480, 32:480]
SHIFTS = [(0, 0), (3, -7), (-9, 4), (6, 10), (-2, -5), (8, -3), (-6, -9), (10, 6)]  # camera (dy, dx)
MIS_FRAME, MIS_EXTRA = 5, 6        # in frame 5, scene rows 128:256, cols 320:448 move 6 px right
DISC_R, DISC_V, DISC_STEP = 12, 1.0, 4   # disc radius, value, px per frame (rightwards)
TAUS = (1, 8, 75)                  # tuning factors: 1, paper's 8, Monod et al.'s 75
CAM = data.camera().astype(float) / 255.0


def scene(k, moving=True, misalign=True):
    s = CAM.copy()
    if misalign and k == MIS_FRAME:
        s[128:256, 320:448] = CAM[128:256, 320 - MIS_EXTRA:448 - MIS_EXTRA]
    yy, xx = np.mgrid[0:512, 0:512]
    cx = 176 + DISC_STEP * (k if moving else 0)
    s[(yy - 432) ** 2 + (xx - cx) ** 2 <= DISC_R ** 2] = DISC_V
    return s


def clean_frames(moving=True, misalign=True):
    return np.array([scene(k, moving, misalign)[OFF + dy:OFF + dy + H, OFF + dx:OFF + dx + W]
                     for k, (dy, dx) in enumerate(SHIFTS)])


def add_noise(clean, rng):
    return rng.poisson(clean / A_NOISE) * A_NOISE + rng.normal(0, np.sqrt(B_NOISE), clean.shape)


# ---------------- regions (reference coordinates) ----------------
def disc_mask(k):
    yy, xx = np.mgrid[0:H, 0:W]
    return (yy - 400) ** 2 + (xx - (144 + DISC_STEP * k)) ** 2 <= DISC_R ** 2


VALID = np.zeros((H, W), bool); VALID[32:416, 32:416] = True
MIS = np.zeros((H, W), bool); MIS[96:224, 288:416] = True
DISCS = np.array([disc_mask(k) for k in range(N)])
TRAIL = DISCS[1:].any(0) & ~DISCS[0]
GROW = lambda m: ndimage.binary_dilation(m, np.ones((2 * n + 1, 2 * n + 1), bool))
STATIC = VALID & ~GROW(MIS) & ~GROW(DISCS.any(0))
REGIONS = {"static": STATIC, "misaligned": MIS & VALID, "ghost": TRAIL & VALID}


def rmse(est, truth, m):
    return np.sqrt(((est - truth)[m] ** 2).mean())


# ---------------- tiles, window, overlap-add ----------------
w1 = 0.5 - 0.5 * np.cos(2 * np.pi * (np.arange(n) + 0.5) / n)   # Hasinoff et al., Section 5
WIN = np.outer(w1, w1)
P = H + 2 * HALF                                    # padded reference size, 464
ORG = np.arange(0, P - n + 1, HALF)                 # tile origins, 57 per axis
APAD = 40                                           # alternates padded by 40 (= HALF + 32)


def pad_ref(f):
    return np.pad(f, HALF, mode="reflect")


def pad_alt(f):
    return np.pad(f, APAD, mode="reflect")


def fetch(img, oy, ox, size, base):
    """Tiles of img at origins (oy, ox) (arrays of the same shape) shifted by base; edge-clamped."""
    r = np.clip(oy[..., None] + base + np.arange(size), 0, img.shape[0] - 1)
    c = np.clip(ox[..., None] + base + np.arange(size), 0, img.shape[1] - 1)
    return img[r[..., :, None], c[..., None, :]]


def stack_tiles(frames, offsets):
    """(N, T, T, n, n) tiles: frame 0 from the padded reference, others at the given offsets."""
    oy, ox = np.meshgrid(ORG, ORG, indexing="ij")
    out = [fetch(pad_ref(frames[0]), oy, ox, n, 0)]
    for k in range(1, N):
        out.append(fetch(pad_alt(frames[k]), oy + offsets[k, ..., 0], ox + offsets[k, ..., 1],
                         n, APAD - HALF))
    return np.array(out)


def overlap_add(tiles):
    out = np.zeros((P, P))
    for i, ty in enumerate(ORG):
        for j, tx in enumerate(ORG):
            out[ty:ty + n, tx:tx + n] += tiles[i, j]
    return out[HALF:HALF + H, HALF:HALF + W]


def noise_var(stack):
    rms = np.sqrt((stack[0] ** 2).mean(axis=(-1, -2)))    # RMS of the reference tile
    return (A_NOISE * rms + B_NOISE)[..., None, None]


def safe_ratio(num, den):
    return np.divide(num, den, out=np.ones_like(num), where=den > 0)


def merge(stack, method, tau=None, var=None):
    """method: 'reference', 'mean', 'pairwise' (Eqs. 6-7), '3d' (3-D DFT shrinkage)."""
    FT = np.fft.fft2(stack * WIN)
    s2 = noise_var(stack) if var is None else var
    if method == "reference":
        M = FT[0]
    elif method == "mean":
        M = FT.mean(0)
    elif method == "pairwise":
        c = n ** 2 / 4 ** 2 * 2 * tau
        D = FT[0] - FT
        Az = safe_ratio(np.abs(D) ** 2, np.abs(D) ** 2 + c * s2)
        M = (FT + Az * D).mean(0)                 # (1 - Az) Tz + Az T0, averaged over z
    elif method == "3d":
        c3 = n ** 2 / 4 ** 2 * N * tau           # our analogue: factor N for an N-point temporal DFT
        X = np.fft.fft(FT, axis=0)
        G = safe_ratio(np.abs(X) ** 2, np.abs(X) ** 2 + c3 * s2)
        M = np.fft.ifft(G * X, axis=0)[0]
    return overlap_add(np.fft.ifft2(M).real)


# ---------------- HDR+ tile alignment (Section 4 and supplement Section 2) ----------------
FACT = [1, 2, 8, 32]                 # cumulative downsampling of levels 0..3 (steps 2, 4, 4)
TSIZE = [16, 16, 16, 8]
RAD = [1, 4, 4, 4]
NORM = [1, 2, 2, 2]


def pyramid(img):
    lv = [img]
    for l in (1, 2, 3):
        f = FACT[l] // FACT[l - 1]
        lv.append(ndimage.gaussian_filter(lv[-1], f / 2, mode="reflect")[::f, ::f])
    return lv


def d2_fft(T, I, r):
    """Eq. 2: ||T||^2 + box(I*I, n) - 2 * cross-correlation, for all (2r+1)^2 offsets."""
    m = T.shape[-1] + 2 * r
    Tp = np.zeros(T.shape[:-2] + (m, m)); Tp[..., :T.shape[-2], :T.shape[-1]] = T
    corr = np.fft.ifft2(np.conj(np.fft.fft2(Tp)) * np.fft.fft2(I)).real[..., :2 * r + 1, :2 * r + 1]
    box = windows(I ** 2, T.shape[-2:], axis=(-2, -1)).sum(axis=(-1, -2))
    return (T ** 2).sum(axis=(-1, -2))[..., None, None] + box - 2 * corr


def dp_brute(T, I, p):
    diff = windows(I, T.shape[-2:], axis=(-2, -1)) - T[..., None, None, :, :]
    return (np.abs(diff) ** p).sum(axis=(-1, -2))


def grid(level_size, l):
    return np.arange(0, level_size - TSIZE[l] + 1, TSIZE[l] // 2)


# Sub-pixel fit (supplement Section 4, Eqs. 13-33): weighted least squares of a quadratic
# to the 3x3 window of D2 around its minimum; u is the column (x) and v the row (y) offset.
_uv = [(u, v) for v in (-1, 0, 1) for u in (-1, 0, 1)]           # row-major vec order
_X = np.array([[0.5 * u * u, u * v, 0.5 * v * v, u, v, 1.0] for u, v in _uv])   # Eq. 16
_Wd = np.diag(np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]], float).ravel())     # Eqs. 15, 17
FSUB = np.linalg.solve(_X.T @ _Wd @ _X, _X.T @ _Wd)                           # Eq. 21


def subpixel(D, k, r):
    """mu = -A^{-1} b (Eq. 30) per tile, as (dy, dx); zero where it cannot be used."""
    iy, ix = k // (2 * r + 1), k % (2 * r + 1)
    inner = (iy > 0) & (iy < 2 * r) & (ix > 0) & (ix < 2 * r)
    iy, ix = np.clip(iy, 1, 2 * r - 1), np.clip(ix, 1, 2 * r - 1)
    a, b = np.indices(D.shape[:2])
    sub = np.stack([D[a, b, iy + dy, ix + dx] for dy in (-1, 0, 1) for dx in (-1, 0, 1)], -1)
    A11, A12, A22, b1, b2, _ = np.moveaxis(sub @ FSUB.T, -1, 0)
    A11, A22 = np.maximum(A11, 0), np.maximum(A22, 0)                         # Eqs. 26-27
    A12 = np.where(A11 * A22 - A12 ** 2 < 0, 0, A12)                          # Eq. 28
    det = A11 * A22 - A12 ** 2
    ok = inner & (det > 0)
    sdet = np.where(ok, det, 1.0)
    mu_u = -(A22 * b1 - A12 * b2) / sdet                                       # Eq. 32
    mu_v = -(A11 * b2 - A12 * b1) / sdet
    ok &= np.hypot(mu_u, mu_v) <= 1                                            # Monod et al.: ||mu|| <= 1
    return np.stack([np.where(ok, mu_v, 0), np.where(ok, mu_u, 0)], -1)


def align(ref, alt):
    """Offsets (T, T, 2) per level-0 tile, from coarse to fine."""
    R, Al = pyramid(pad_ref(ref)), pyramid(pad_alt(alt))
    off = None
    for l in (3, 2, 1, 0):
        g = grid(R[l].shape[0], l)
        oy, ox = np.meshgrid(g, g, indexing="ij")
        base = (APAD - HALF) // FACT[l]
        T = fetch(R[l], oy, ox, TSIZE[l], 0)
        if off is None:
            cand = np.zeros(oy.shape + (2,), int)
        else:   # three candidates from the nearest coarse tiles, chosen by L1 residual
            F = FACT[l + 1] // FACT[l]
            gc = grid(R[l + 1].shape[0], l + 1)
            cc = gc + TSIZE[l + 1] / 2
            x = (g + TSIZE[l] / 2) / F
            i0 = np.abs(x[:, None] - cc[None, :]).argmin(1)
            i1 = np.clip(i0 + np.sign(x - cc[i0]).astype(int), 0, len(cc) - 1)
            cands = [off[i0[:, None], i0[None, :]], off[i1[:, None], i0[None, :]],
                     off[i0[:, None], i1[None, :]]]
            cands = [np.rint(c * F).astype(int) for c in cands]
            errs = [np.abs(fetch(Al[l], oy + c[..., 0], ox + c[..., 1], TSIZE[l], base) - T).sum((-1, -2))
                    for c in cands]
            cand = np.choose(np.argmin(errs, 0)[..., None], cands)
        r = RAD[l]
        I = fetch(Al[l], oy + cand[..., 0] - r, ox + cand[..., 1] - r, TSIZE[l] + 2 * r, base)
        D = d2_fft(T, I, r) if NORM[l] == 2 else dp_brute(T, I, 1)
        k = D.reshape(D.shape[:2] + (-1,)).argmin(-1)
        off = cand + np.stack([k // (2 * r + 1) - r, k % (2 * r + 1) - r], -1)
        if l > 0:                 # sub-pixel refinement at the L2 levels only
            off = off + subpixel(D, k, r)
    return off


def true_offsets():
    """Global (camera) offsets for every tile and frame, and the local truth in R_mis."""
    T = len(ORG)
    glob = np.zeros((N, T, T, 2), int)
    for k, (dy, dx) in enumerate(SHIFTS):
        glob[k] = (-dy, -dx)
    local = glob.copy()
    oy, ox = np.meshgrid(ORG - HALF, ORG - HALF, indexing="ij")    # tile origins in ref coords
    inside = (oy >= 96) & (oy + n <= 224) & (ox >= 288) & (ox + n <= 410)
    local[MIS_FRAME][inside] += (0, MIS_EXTRA)
    return glob, local, inside


def tile_sets():
    """Static level-0 tiles (no overlap with R_mis or any disc, inside the 32-px margin)."""
    T = len(ORG)
    ok = np.zeros((T, T), bool)
    bad = MIS | DISCS.any(0)
    for i, ty in enumerate(ORG - HALF):
        for j, tx in enumerate(ORG - HALF):
            if ty >= 32 and tx >= 32 and ty + n <= 416 and tx + n <= 416:
                ok[i, j] = not bad[ty:ty + n, tx:tx + n].any()
    return ok


METHODS = [("reference", None), ("mean", None)] + [("pairwise", t) for t in TAUS] + \
          [("3d", t) for t in TAUS]


def label(m, t):
    return m if t is None else f"{m} tau={t}"


# ================================ the example ================================
import scipy, skimage
SEEDS = range(5)
print(f"NumPy {np.__version__}, SciPy {scipy.__version__}, scikit-image {skimage.__version__}")
print(f"{N} frames of camera ({H}x{W}), noise variance {A_NOISE}*x + {B_NOISE}; tiles {n}x{n}, half overlap")
truth = clean_frames()[0]
glob, local, inside = true_offsets()
static_tiles = tile_sets()

# ---- checks ----
ones = overlap_add(np.broadcast_to(WIN, (len(ORG), len(ORG), n, n)))
print(f"C1 overlapped windows sum to one: max |sum - 1| = {np.abs(ones - 1)[VALID].max():.1e}")
print(f"   window sum of squares / n^2 = {(WIN ** 2).sum() / n ** 2:.6f} (paper's window factor 1/4^2 = {1 / 16:.6f})")
st = stack_tiles(clean_frames(moving=False, misalign=False), glob)
zero = np.zeros_like(noise_var(st))
c2 = [np.abs(merge(st, "pairwise", 8) - truth)[VALID].max(), np.abs(merge(st, "mean") - truth)[VALID].max(),
      np.abs(merge(st, "3d", 8, var=zero) - truth)[VALID].max()]
print(f"C2 noise-free static burst, max |merge - truth|: pairwise tau=8 {c2[0]:.1e}, mean {c2[1]:.1e}, "
      f"3-D with sigma^2 = 0 {c2[2]:.1e}")
cfull = clean_frames()
nf = np.array([np.zeros_like(glob[0])] + [align(cfull[0], cfull[k]) for k in range(1, N)])
nf_static = np.mean([(nf[k] == glob[k]).all(-1)[static_tiles].mean() for k in range(1, N)])
nf_mis = (nf[MIS_FRAME] == local[MIS_FRAME]).all(-1)[inside].mean()
print(f"C5 alignment on noise-free frames: exact on static tiles {nf_static:.4f}, "
      f"frame {MIS_FRAME} inside R_mis {nf_mis:.4f}")

res = {}          # (alignment, method label, region) -> list over seeds
align_stats = []
c3a = c3b = c4 = 0.0
c6, c7 = [], []
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    fr = add_noise(cfull, rng)
    nz = fr - cfull
    slope, icpt = np.polyfit(cfull[:, VALID].ravel(), (nz[:, VALID] ** 2).ravel(), 1)
    c6 += [slope / A_NOISE, icpt / B_NOISE]
    est = np.array([np.zeros_like(glob[0])] + [align(fr[0], fr[k]) for k in range(1, N)])
    ex = np.mean([(est[k] == glob[k]).all(-1)[static_tiles].mean() for k in range(1, N)])
    w1px = np.mean([(np.abs(est[k] - glob[k]).max(-1) <= 1)[static_tiles].mean() for k in range(1, N)])
    exm = (est[MIS_FRAME] == local[MIS_FRAME]).all(-1)[inside].mean()
    align_stats.append((ex, w1px, exm))
    if seed == 0:
        R1, A1 = pyramid(pad_ref(fr[0])), pyramid(pad_alt(fr[1]))
        g = grid(R1[1].shape[0], 1); oy, ox = np.meshgrid(g, g, indexing="ij")
        Tt, It = fetch(R1[1], oy, ox, 16, 0), fetch(A1[1], oy - 4, ox - 4, 24, 16)
        b = dp_brute(Tt, It, 2); c4 = np.abs(d2_fft(Tt, It, 4) - b).max() / np.abs(b).max()
    for aname, offs in (("global", glob), ("HDR+ tiles", est)):
        st = stack_tiles(fr, offs)
        if aname == "global":
            c3a = max(c3a, np.abs(merge(st, "pairwise", 0) - fr[0])[VALID].max())
            c3b = max(c3b, np.abs(merge(st, "pairwise", 1e12) - merge(st, "mean"))[VALID].max())
            pred = np.sqrt(((A_NOISE * truth + B_NOISE) / N)[STATIC].mean())
        for m, t in METHODS:
            out = merge(st, m, t)
            for rname, mask in REGIONS.items():
                res.setdefault((aname, label(m, t), rname), []).append(rmse(out, truth, mask))
            if aname == "global" and m == "mean":
                c7.append(rmse(out, truth, STATIC) / pred)
print(f"C3 limits (global alignment): tau = 0 gives the reference frame, max diff {c3a:.1e}; "
      f"tau = 1e12 gives the mean, max diff {c3b:.1e}")
print(f"C4 Eq. 2 by FFT against brute force (level 1, seed 0): max relative diff {c4:.1e}")
print(f"C6 noise regression over seeds: slope/A {min(c6[0::2]):.4f}-{max(c6[0::2]):.4f}, "
      f"intercept/B {min(c6[1::2]):.4f}-{max(c6[1::2]):.4f}")
print(f"C7 mean RMSE (static) / predicted sigma/sqrt(N): {min(c7):.4f}-{max(c7):.4f}")

a = np.array(align_stats)
print("\nHDR+ tile alignment on the noisy bursts, seeds 0-4 (min-max):")
print(f"  static tiles exact {a[:, 0].min():.4f}-{a[:, 0].max():.4f}, within 1 px "
      f"{a[:, 1].min():.4f}-{a[:, 1].max():.4f}; frame {MIS_FRAME} inside R_mis exact "
      f"{a[:, 2].min():.4f}-{a[:, 2].max():.4f}")

print("\nRMSE, mean over seeds 0-4 [min, max]")
print(f"{'alignment':11s} {'method':16s} {'static':>24s} {'misaligned':>24s} {'ghost':>24s}")
for aname in ("global", "HDR+ tiles"):
    for m, t in METHODS:
        cells = []
        for rname in REGIONS:
            v = np.array(res[(aname, label(m, t), rname)])
            cells.append(f"{v.mean():.5f} [{v.min():.5f}, {v.max():.5f}]")
        print(f"{aname:11s} {label(m, t):16s} " + " ".join(f"{c:>24s}" for c in cells))

print("\nPaired differences over seeds 0-4 (min, max)")
def diff(a1, l1, a2, l2, region):
    d = np.array(res[(a1, l1, region)]) - np.array(res[(a2, l2, region)])
    return f"{d.min():+.5f}, {d.max():+.5f}"
for aname in ("global", "HDR+ tiles"):
    P = lambda t: f"pairwise tau={t}"
    Q = lambda t: f"3d tau={t}"
    print(f"[{aname}]")
    print(f"  H1 static: tau1-tau8 {diff(aname, P(1), aname, P(8), 'static')}; tau8-tau75 "
          f"{diff(aname, P(8), aname, P(75), 'static')}; tau75-mean {diff(aname, P(75), aname, 'mean', 'static')}")
    print("  H2 static, 3-D minus pairwise: " + "; ".join(
        f"tau={t} {diff(aname, Q(t), aname, P(t), 'static')}" for t in TAUS))
    print("  H3 misaligned, pairwise minus 3-D: " + "; ".join(
        f"tau={t} {diff(aname, P(t), aname, Q(t), 'misaligned')}" for t in TAUS))
    print(f"  H4 misaligned, pairwise tau=8 minus mean: {diff(aname, P(8), aname, 'mean', 'misaligned')}")
    print("  H5 ghost, pairwise minus mean: " + "; ".join(
        f"tau={t} {diff(aname, P(t), aname, 'mean', 'ghost')}" for t in TAUS)
          + f"; tau1-tau8 {diff(aname, P(1), aname, P(8), 'ghost')}; tau8-tau75 {diff(aname, P(8), aname, P(75), 'ghost')}")
print("[HDR+ tiles minus global]")
print(f"  H7 misaligned: mean {diff('HDR+ tiles', 'mean', 'global', 'mean', 'misaligned')}; pairwise tau=8 "
      f"{diff('HDR+ tiles', 'pairwise tau=8', 'global', 'pairwise tau=8', 'misaligned')}")
print(f"  H8 static: mean {diff('HDR+ tiles', 'mean', 'global', 'mean', 'static')}")

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; about 16–18 seconds; two runs printed identical output):

NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2
8 frames of camera (448x448), noise variance 0.01*x + 0.0004; tiles 16x16, half overlap
C1 overlapped windows sum to one: max |sum - 1| = 2.2e-16
   window sum of squares / n^2 = 0.140625 (paper's window factor 1/4^2 = 0.062500)
C2 noise-free static burst, max |merge - truth|: pairwise tau=8 7.8e-16, mean 7.8e-16, 3-D with sigma^2 = 0 6.7e-16
C5 alignment on noise-free frames: exact on static tiles 0.9283, frame 5 inside R_mis 0.9619
C3 limits (global alignment): tau = 0 gives the reference frame, max diff 8.9e-16; tau = 1e12 gives the mean, max diff 5.2e-10
C4 Eq. 2 by FFT against brute force (level 1, seed 0): max relative diff 4.3e-15
C6 noise regression over seeds: slope/A 0.9961-1.0003, intercept/B 0.9825-1.0272
C7 mean RMSE (static) / predicted sigma/sqrt(N): 0.9955-1.0032

HDR+ tile alignment on the noisy bursts, seeds 0-4 (min-max):
  static tiles exact 0.5179-0.5204, within 1 px 0.7619-0.7680; frame 5 inside R_mis exact 0.4238-0.5095

RMSE, mean over seeds 0-4 [min, max]
alignment   method                             static               misaligned                    ghost
global      reference        0.06715 [0.06694, 0.06740] 0.08451 [0.08383, 0.08522] 0.07978 [0.07614, 0.08347]
global      mean             0.02373 [0.02362, 0.02380] 0.03363 [0.03351, 0.03377] 0.18439 [0.18323, 0.18541]
global      pairwise tau=1   0.05262 [0.05247, 0.05285] 0.06669 [0.06614, 0.06737] 0.06714 [0.06358, 0.07052]
global      pairwise tau=8   0.03139 [0.03129, 0.03156] 0.04014 [0.03988, 0.04051] 0.05099 [0.04851, 0.05379]
global      pairwise tau=75  0.02400 [0.02389, 0.02407] 0.03135 [0.03130, 0.03142] 0.05964 [0.05783, 0.06214]
global      3d tau=1         0.05145 [0.05131, 0.05163] 0.06546 [0.06494, 0.06605] 0.06814 [0.06499, 0.07186]
global      3d tau=8         0.02718 [0.02711, 0.02724] 0.03461 [0.03423, 0.03494] 0.05374 [0.05144, 0.05620]
global      3d tau=75        0.02544 [0.02540, 0.02549] 0.03075 [0.03060, 0.03091] 0.07064 [0.06977, 0.07190]
HDR+ tiles  reference        0.06715 [0.06694, 0.06740] 0.08451 [0.08383, 0.08522] 0.07978 [0.07614, 0.08347]
HDR+ tiles  mean             0.02397 [0.02390, 0.02404] 0.02888 [0.02871, 0.02902] 0.05585 [0.05402, 0.05777]
HDR+ tiles  pairwise tau=1   0.05305 [0.05288, 0.05329] 0.06682 [0.06632, 0.06749] 0.06713 [0.06403, 0.07021]
HDR+ tiles  pairwise tau=8   0.03249 [0.03240, 0.03267] 0.04061 [0.04023, 0.04098] 0.05145 [0.05038, 0.05348]
HDR+ tiles  pairwise tau=75  0.02455 [0.02449, 0.02463] 0.02987 [0.02968, 0.03003] 0.05404 [0.05303, 0.05550]
HDR+ tiles  3d tau=1         0.05179 [0.05164, 0.05202] 0.06482 [0.06427, 0.06541] 0.06679 [0.06331, 0.07001]
HDR+ tiles  3d tau=8         0.02804 [0.02794, 0.02815] 0.03284 [0.03253, 0.03321] 0.05163 [0.05051, 0.05290]
HDR+ tiles  3d tau=75        0.02617 [0.02615, 0.02621] 0.02736 [0.02710, 0.02759] 0.06089 [0.06049, 0.06134]

Paired differences over seeds 0-4 (min, max)
[global]
  H1 static: tau1-tau8 +0.02110, +0.02131; tau8-tau75 +0.00733, +0.00749; tau75-mean +0.00026, +0.00028
  H2 static, 3-D minus pairwise: tau=1 -0.00122, -0.00114; tau=8 -0.00432, -0.00408; tau=75 +0.00133, +0.00160
  H3 misaligned, pairwise minus 3-D: tau=1 +0.00109, +0.00139; tau=8 +0.00530, +0.00573; tau=75 +0.00039, +0.00082
  H4 misaligned, pairwise tau=8 minus mean: +0.00619, +0.00696
  H5 ghost, pairwise minus mean: tau=1 -0.12163, -0.11271; tau=8 -0.13670, -0.12943; tau=75 -0.12686, -0.12237; tau1-tau8 +0.01507, +0.01727; tau8-tau75 -0.01051, -0.00706
[HDR+ tiles]
  H1 static: tau1-tau8 +0.02045, +0.02063; tau8-tau75 +0.00786, +0.00803; tau75-mean +0.00057, +0.00060
  H2 static, 3-D minus pairwise: tau=1 -0.00132, -0.00124; tau=8 -0.00451, -0.00436; tau=75 +0.00152, +0.00172
  H3 misaligned, pairwise minus 3-D: tau=1 +0.00185, +0.00207; tau=8 +0.00757, +0.00790; tau=75 +0.00227, +0.00293
  H4 misaligned, pairwise tau=8 minus mean: +0.01152, +0.01199
  H5 ghost, pairwise minus mean: tau=1 +0.00743, +0.01370; tau=8 -0.00622, -0.00296; tau=75 -0.00228, -0.00098; tau1-tau8 +0.01365, +0.01673; tau8-tau75 -0.00425, -0.00153
[HDR+ tiles minus global]
  H7 misaligned: mean -0.00506, -0.00457; pairwise tau=8 +0.00027, +0.00059
  H8 static: mean +0.00016, +0.00028

Checks. All seven held. The windows sum to one to rounding error, and on a noise-free static burst every merge returns the truth, since all differences $D_z$ vanish. The $\tau$ limits behave as Monod et al. describe. The alignment line C5 shows that our implementation is not exact even without noise; in our reading, the misses lie in the two situations that Lafenetre et al. (2023, Section 2.1) name as limits of local registration, near edges (the aperture problem) and in flat areas, where the residual surface is a valley or nearly flat and the coarse-to-fine search can settle on a wrong offset. The window line shows that $\sum w^2 / n^2$ for the 2-D raised cosine is $0.140625$, not the paper’s $1/4^2$: along each axis $\sum_x w(x)^2 = 3n/8$, so the 2-D value is $(3/8)^2 = 9/64$, whereas $1/4^2$ is the square of the 2-D window’s mean value, $1/4$. Since a DFT coefficient of windowed white noise of variance $\sigma^2$ has variance $\sigma^2 \sum w^2$, in our reading the paper’s $c$ assumes a noise variance in $D_z$ that is $9/4$ times smaller than the window gives, which a larger tuning factor offsets, consistent with Monod et al.’s remark that $\tau$ absorbs $k$.

Static region. With the oracle alignment, everything in the static region differs only by noise, and the mean is the best of the methods, at the noise level C7 predicts. The pairwise merge approaches it as $\tau$ grows, and at the paper’s $\tau = 8$ it leaves roughly a third more error than the mean: that is the noise reduction the robust merge gives up where alignment succeeds. The 3-D shrinkage beats the pairwise merge at $\tau = 1$ and 8, as the paper states, but not at 75: by construction, its gain goes to zero for every coefficient as $c_3$ grows, signal included, while the pairwise merge tends to the mean.

Misaligned block. One frame in eight is displaced by 6 pixels inside the block, so the mean takes one eighth of a shifted image there; its error rises only moderately above the static level. The pairwise merge at $\tau = 8$ does worse than the mean, and the 3-D shrinkage does better than the pairwise merge at every $\tau$. Hasinoff et al.’s Fig. 7 argument does not appear in this set-up. Two points from the construction bear on it. First, our 3-D constant $c_3$ is our own analogue, so H3 tests that shrinkage and not the 3-D setting of the paper’s figure. Second, at the same $\tau$ the constants differ by the factor $N/2 = 4$, so equal $\tau$ is not equal shrinkage strength for the two methods.

Ghost trail. The mean with oracle alignment leaves a strong ghost, since seven of eight frames put a white disc or its trail where the reference shows grass. With that alignment, every pairwise merge cuts the ghost error to well under half of the mean’s. The pairwise ghost error is not smallest at the smallest $\tau$: the trail is small, and a merge that keeps close to the reference keeps its noise, whose level the reference row shows. Between $\tau = 8$ and $\tau = 75$ the ghost grows, as the merge leans more on the alternates.

HDR+ alignment. On noisy frames only about half the static tiles get exactly the true offset, and about three quarters are within one pixel. The static error of the mean nevertheless rises only slightly over the oracle. In our reading, a one-pixel misplacement barely changes a tile with little texture, so such misses add little error. In the misaligned block, tile alignment helps the mean, because tiles inside the block can follow frame 5’s local displacement; it does not help the pairwise merge at $\tau = 8$, which in our reading already leaned on the reference there. The ghost of the mean falls sharply with tile alignment; in our reading, tiles near the disc partly follow it.

The predictions, each read from the paired differences:

• H1 (static: pairwise error falls with $\tau$ and stays above the mean): held under both alignments.
• H2 (static: 3-D below pairwise at the same $\tau$): held at $\tau = 1$ and 8 and failed at $\tau = 75$, under both alignments.
• H3 (misaligned block: pairwise below our 3-D shrinkage at the same $\tau$): failed at every $\tau$ under both alignments.
• H4 (misaligned block: pairwise at $\tau = 8$ below the mean): failed under both alignments.
• H5 (ghost: pairwise below the mean at every $\tau$, and the pairwise ghost rising with $\tau$): with global alignment, pairwise below the mean held at every $\tau$; with HDR+ tiles it failed at $\tau = 1$ and held at 8 and 75. The rise from $\tau = 1$ to 8 failed under both alignments; the rise from 8 to 75 held under both.
• H6 (noisy alignment less often exact than noise-free): held.
• H7 (misaligned block: lower error with HDR+ tiles than with global alignment): held for the mean and failed for the pairwise merge at $\tau = 8$.
• H8 (static: the mean’s error higher with HDR+ tiles than with the oracle): held.

All settings, regions, seeds and values of $\tau$ were fixed before the first run, and nothing was changed after it.

The example is an illustration, not a benchmark: the frames share one plane, the noise model is known exactly, and the misalignment and motion are of one size each. It shows what Hasinoff et al. and Monod et al. describe, that $\tau$ trades noise reduction against robustness, and also, as the reversal of H2 between $\tau = 8$ and $\tau = 75$ shows, that the ranking of two merges at equal $\tau$ depends on how their constants are scaled, which the sources leave open for the 3-D filter.

Variants

• Frequency-domain and 3-D stacks. Video denoisers that filter 3-D stacks of matched patches, from Kokaram’s 3-D DFT Wiener filtering to V-BM3D, are the methods HDR+ departs from (Hasinoff et al., 2016, Section 5); BM3D describes the single-image collaborative filter on which V-BM3D builds.
• Fourier burst accumulation. Delbracio and Sapiro’s method, as Hasinoff et al. (2016, Section 5) describe it, combines frequency content across a burst more aggressively to remove motion blur from long exposures, in the extreme taking the maximum of each frequency over the burst. HDR+ instead keeps the reference frame’s blur. Hasinoff et al. cite the video-deblurring paper of that year in IEEE Transactions on Computational Imaging and Monod et al. the camera-shake paper in IEEE Transactions on Image Processing, two works by the same authors in 2015. Monod et al. (2021, Section 4.2) add that the HDR+ filter neither introduces blur from other frames nor removes it.
• Burst super-resolution. When frames are offset by fractions of a pixel, a merge onto a finer grid recovers detail beyond one frame’s sampling. Wronski et al. (2019) merge raw colour-filter-array frames directly onto an RGB grid, replacing demosaicking; computational photography describes their use of hand tremor, and Lafenetre, Facciolo and Eboli (2023) reimplement the method. As Lafenetre et al. describe it, alignment is block matching on a four-level pyramid, as in HDR+, refined to sub-pixel accuracy by three Lucas–Kanade-type iterations; the merge is a memory-efficient variant of normalized convolution, a kernel regression whose anisotropic Gaussian kernels follow local structure tensors, weighted by a robustness mask against misalignment; and the zoom factor is limited to about 2. On a synthetic set, their implementation reached about 0.1 pixel of registration error after three iterations when the frames were first low-pass filtered rather than averaged in 2 × 2 Bayer blocks (Section 2.1, Fig. 1). Super-resolution covers single- and multi-image methods in general.
• Learned burst denoising. Mildenhall et al. (2018) trained a kernel prediction network that outputs, for every pixel, a stack of $N$ filters of $5 \times 5$ taps applied to the frames and averaged, so that the network aligns, weights and denoises at once (Section 4, Eq. 2). It was trained on synthetic bursts made from internet photographs: gamma undone, misalignments of up to ±2 pixels, occasional large shifts to simulate alignment failure, and shot and read noise drawn from the range measured on a phone (Section 3.2). An annealed loss that first asks each filtered frame to match the truth keeps the network from ignoring the alternates (Section 4.2). On their linear synthetic test set at the lowest gain, Table 1 reports PSNR 36.47 dB for the network, 34.60 for VBM4D, 31.96 for HDR+, 28.70 for the reference frame and 24.70 for the plain burst average; at the highest gain, which lies outside the noise levels seen in training, 27.97, 26.52, 20.05, 15.76 and 20.00. The authors state that they ran HDR+ with its spatial denoising disabled “by setting c = 0 in Eq. 7”; taken literally, $c = 0$ makes $A_z = 1$ and returns the reference frame (check C3 above), yet their Table 1 scores HDR+ above the reference frame at every gain, so in our reading the remark refers to the spatial step’s constant, which the paper does not spell out.
• Very low light. Night Sight extends the HDR+ merge to longer frames chosen by measured motion; computational photography describes it.

Best Practices

• Work on linear raw data with a calibrated noise model. The merge’s decisions are only as good as $\sigma^2$; Hasinoff et al. (2016, Section 2) give linearity and an accurate noise model as reasons to start from raw, and Monod et al. (2021, Section 4.1) found metadata-based or ISO-scaled parameters more consistent than blind noise estimation.
• Keep the exposure constant across the burst. Equal exposures align more reliably than brackets (Hasinoff et al., 2016, Section 1).
• Treat the tuning factor as a tuning factor. $c$ mixes a scaling of noise variances with a free factor, and Monod et al. needed $\tau = 75$ where the paper used 8; the example shows the static noise and the ghost moving in opposite directions as $\tau$ changes.
• Overlap tiles and use a window that sums to one. Without them, tiles merged independently show seams and DFT edge artifacts (Monod et al., 2021, Section 4.4, Figs. 6 and 8).
• Judge a merge by region. A single RMSE over the image hides the trade: the example’s static, misaligned and ghost regions rank the methods differently. Monod et al. (2021, abstract) note that HDR+ does not necessarily aim at maximizing standard denoising metrics, but at natural, visually pleasing images.

Pitfalls

• Averaging without robustness. The temporal mean is optimal only where frames agree; with motion it ghosts, as the example’s oracle-aligned mean shows on the moving disc.
• Expecting a robust merge to remove noise everywhere. Where alternates are rejected, the output stays close to the reference and its noise; HDR+ accepts that some regions look noisier than others (Hasinoff et al., 2016, Section 1).
• Reading $c$ as a calibrated constant. The window factor $1/4^2$ does not equal the window’s $\sum w^2 / n^2$, as the example prints, and Monod et al. call its justification questionable.
• Sub-pixel offsets on Bayer data. A sub-pixel shift of a downsampled grey image mixes colour planes in the raw mosaic; Monod et al. (2021, Section 3.3) observed significant colour shifts and skip sub-pixel refinement at the finest level.
• Comparing merges at “the same” parameter. Constants scale differently in different filters (2 against $N$ in the example), so equal $\tau$ is not equal strength; in the example the static ranking of the 3-D and pairwise merges reverses between $\tau = 8$ and $\tau = 75$.
• Trusting timing figures across devices. The 4 seconds per 12-megapixel photograph and 24 ms per megapixel for alignment are Hasinoff et al.’s measurements on 2016 phones.

Evidence and Impact

Hasinoff et al. (2016, Section 8) describe HDR+ as deployed on the Nexus 6, 5X and 6P, with bursts of up to 8 frames on 12–13-megapixel sensors (Section 7). Their evidence is mostly qualitative: side-by-side comparisons with a single-exposure pipeline on the same phone (Fig. 10), examples of failure (Fig. 11), and, in the supplement (Section 6), full-resolution results for 30 bursts processed by their merge and by JPEG-based burst fusion methods (two variants of Liu et al.’s 2014 method, CV-BM3D with a hand-picked noise level, and Photoshop’s Merge to HDR Pro) through a common dcraw finishing pipeline. The authors write that a quantitative comparison with competing commercial systems is hard because their algorithms are proprietary, and released an archive of raw bursts with their merged and final results.

Monod et al. (2021, Section 6.1) describe that archive, the HDR+ dataset, as 3,640 bursts (28,461 images, 765 GB) with a curated subset of 153 bursts and results from 2016 and from a 2017 version with “algorithm refinements and updated tuning”. Comparing their reimplementation with Google’s 2017 merged images on the subset bursts whose sizes matched, they report an average PSNR of 55.2 dB between the two merges, against 48.6 dB between Google’s merge and the noisy reference, which they read as evidence that their alignment and temporal denoising behave similarly to Google’s (Section 6.2.1, Fig. 12); they also found that the 2016 merged files probably include spatial denoising. The GitHub copy of their open-source implementation is under the AGPL-3.0 licence as of November 2026.

Learned methods were compared against HDR+ on synthetic data with known truth: Mildenhall et al.’s kernel prediction network scored above HDR+ and the other baselines at every gain in their Table 1. Burst super-resolution of the kind Wronski et al. introduced was, by Lafenetre et al.’s account (2023, Section 1), used for the digital zoom of the Google Pixel 3, and the code accompanying their IPOL article is in the first author’s GitHub repository, under the MIT licence as of November 2026.