Tone Mapping

Tone mapping turns a high-dynamic-range image, whose values are proportional to scene radiance and may span more than four orders of magnitude, into an image that a display or a print can show, whose range is about two orders of magnitude. A global operator applies one curve to every pixel. A local operator lets the curve vary with each pixel’s surroundings, which keeps texture visible in both the shadows and the highlights, at the risk of halos: bright or dark fringes along strong edges.

High dynamic range imaging covers how a radiance map is captured and merged, summarizes tone mapping, and shows that one global curve cannot serve a dark room and a bright window at once. This article covers the operators themselves: Reinhard and colleagues’ photographic operator, in its global form and its local dodging-and-burning form, and Durand and Dorsey’s split of an image into a base layer and a detail layer. A global tone curve is a point operation, and the output still has to be encoded for the display with a transfer function, the subject of gamma correction.

Problem

Reinhard, Stark, Shirley and Ferwerda (2002, Section 1) put the problem in numbers: real-world light spans about ten orders of magnitude of absolute range, and “over four orders” of dynamic range within one scene, while prints and displays reproduce at best about two. They call the question of how to map measured or simulated scene luminances to display luminances, so that the image looks satisfactory, the tone reproduction problem.

The two papers that this article follows aim at a pleasing picture rather than a perceptual match. Durand and Dorsey (2002, Section 6) state that their method is not strictly a tone reproduction operator in the sense of Tumblin and Rushmeier (1993), because it does not try to imitate human vision; Reinhard et al. (2002, Section 1) turn to photographic practice because perceptual models of the time introduced artifacts such as ringing or visible clamping, and because appearance depends on more than contrast and brightness. The goal is to keep detail visible everywhere while the overall contrast shrinks by orders of magnitude.

Inputs and Outputs

Input: a radiance map, or the luminance of one, in floating point: values proportional to scene radiance, as produced by merging exposures or by a renderer. Absolute scale need not be known; the operators below normalize it.
Parameters: a key value or overall brightness, and operator-specific settings (a white point; a threshold and scales; spatial and range widths and a target contrast).
Output: display luminance in $[0, 1]$ (values above 1 are clipped), which is then encoded with a display transfer function and quantized, for example to 8 bits.

Both papers handle colour outside the operator. Reinhard et al. (2002, Section 4) compute luminance as $L = 0.27R + 0.67G + 0.06B$ ; Durand and Dorsey compress the intensity and recompose colour from the ratios of each channel to it (Section 6), and their project page gives the intensity as $(20R + 40G + B)/61$ and the output channel as that ratio times the new intensity, noting that any intensity formula may be used instead.

Intuition

Photographers faced the problem before computers. Reinhard et al. (2002, Section 2) build on Ansel Adams’s Zone System: the photographer measures a surface perceived as middle grey and places it on print zone V, about 18% reflectance; a scene that spans more zones than the print holds loses its extremes to black or white unless parts of the print are dodged (given less light during printing, so they print lighter) or burned (given more, so they print darker). A global tone curve is the exposure and development; dodging-and-burning is a local change of exposure.

A local operator can be seen as dividing each pixel by an estimate of its local adaptation level. Where that estimate is a smooth version of the image, a pixel keeps its contrast relative to its neighbours while large-scale differences, between a room and a window, are compressed. The difficulty is the estimate near a strong edge. If it averages across the edge, a dark pixel next to a bright window is divided by a level raised by the window and prints too dark, a dark halo; the bright side gets a bright halo. Reinhard’s operator avoids the problem by choosing, at each pixel, the largest neighbourhood that does not cross a large contrast; Durand and Dorsey avoid it with an edge-preserving filter, the bilateral filter.

Algorithm

Reinhard et al.: global photographic operator

Compute the log-average luminance $\bar L_w$ as an approximation of the scene’s key.
Scale the image so that $\bar L_w$ maps to the key value $a$ , by default 0.18.
Apply the compressive curve $L/(1+L)$ , or its version with a white point that lets the brightest values burn out.

Reinhard et al.: automatic dodging-and-burning

Scale as in steps 1–2 above.
At each of several scales $s$ , blur the scaled luminance with a centre Gaussian and a surround Gaussian 1.6 times as wide.
At each pixel, compute a normalized centre–surround difference at each scale, and pick the scale $s_m$ at which it falls below a threshold $\epsilon$ ; the paper’s rule is quoted under Mathematical Formulation.
Divide the pixel by $1 + V_1(s_m)$ , the centre blur at that scale, instead of $1 + L$ , and clip the result to 1.

Durand and Dorsey: bilateral base and detail

Compute the logarithm (base 10) of the input intensity.
Filter it with a bilateral filter to obtain the base layer; the detail layer is the log intensity minus the base.
Multiply the base by a compression factor chosen so that its range becomes a target contrast, by default 5 (a factor of 5 in luminance).
Add the detail back unchanged, shift so that the largest compressed base value maps to 1, exponentiate and recompose colour.

Mathematical Formulation

Global operator

Reinhard et al. (2002, Section 3.1) define the log-average luminance as

\bar L_w = \exp\!\left( \frac{1}{N} \sum_{x,y} \log\bigl(\delta + L_w(x,y)\bigr) \right),

with $N$ pixels and a small $\delta$ that avoids the logarithm of zero. This is the corrected form: the authors’ project page gives an erratum stating that Equation 1 as printed, with $1/N$ outside the exponential, is wrong, and that the source code and all results in the paper are correct. The scaled luminance is

L(x,y) = \frac{a}{\bar L_w}\, L_w(x,y),

where the key value $a = 0.18$ for a normal-key scene; the authors vary it from 0.045 to 0.72 for low- and high-key scenes. The simple operator is

L_d = \frac{L}{1 + L},

which scales low luminances by about 1 and high luminances by about $1/L$ , and brings every value into $[0, 1)$ . Because a print need not show the brightest values, the authors extend it to

L_d = \frac{L \left( 1 + L / L_\text{white}^2 \right)}{1 + L},

where $L_\text{white}$ is the smallest luminance mapped to pure white; it reverts to $L/(1+L)$ as $L_\text{white} \to \infty$ , and by default $L_\text{white}$ is the largest luminance in the scene. In log terms the simple curve has slope $\mathrm{d}\log L_d / \mathrm{d}\log L = 1/(1+L)$ (our derivation): detail is kept where $L \ll 1$ and divided by about $L$ where $L \gg 1$ .

Dodging-and-burning

The local operator (Section 3.2) uses circularly symmetric Gaussian profiles

R_i(x, y, s) = \frac{1}{\pi (\alpha_i s)^2} \exp\!\left( -\frac{x^2 + y^2}{(\alpha_i s)^2} \right),

whose standard deviation is $\alpha_i s/\sqrt{2}$ , convolved with the scaled luminance to give a centre response $V_1$ ( $i = 1$ ) and a surround response $V_2$ ( $i = 2$ ). The normalized centre–surround function is

V(x, y, s) = \frac{V_1(x, y, s) - V_2(x, y, s)}{2^\phi a / s^2 + V_1(x, y, s)},

where dividing by $V_1$ makes the threshold independent of absolute luminance and the term $2^\phi a/s^2$ keeps $V$ from growing large when $V_1$ approaches zero (the paper’s text says “when $V$ approaches zero”; in our reading $V_1$ is meant). The paper selects the scale $s_m$ with $|V(x, y, s_m)| < \epsilon$ , and the local operator is

L_d(x, y) = \frac{L(x, y)}{1 + V_1\bigl(x, y, s_m(x, y)\bigr)}.

A dark pixel in a bright surround has $L < V_1$ and is darkened, and a pixel in a dark surround is compressed less: in both cases its contrast against its surround grows (Section 3.2). Because $V_1$ can be smaller than $L$ for very bright pixels, $L_d$ can exceed 1, and the authors clip it.

The text leaves two details open. It says “Starting at the lowest scale, we seek the first scale $s_m$ ” where $|V| < \epsilon$ , while the same section says the local area is in principle the largest one around the pixel without large contrast changes, and its Figure 9 shows that too small a scale loses detail. Read literally, in our reading, the first sentence would stop at the smallest scale wherever the image is smooth. The second detail is the set of scales: Section 4 gives “8 discrete scales increasing with a factor of 1.6 from 1 pixel wide to 43 pixels wide”, but eight scales from 1 by factors of 1.6 end at $1.6^7 \approx 26.8$ , and 43 is $1.6^8$ , a ninth scale. The example below runs both readings of each.

Base and detail

Durand and Dorsey (2002, Section 3.3) write the bilateral filter as

J_s = \frac{1}{k(s)} \sum_{p \in \Omega} f(p - s)\, g(I_p - I_s)\, I_p, \qquad k(s) = \sum_{p \in \Omega} f(p - s)\, g(I_p - I_s),

with a Gaussian spatial kernel $f$ and a Gaussian influence function $g$ in intensity; the bilateral filter article derives its behaviour. With $I$ the log intensity, the base is $B = J$ , the detail is $D = I - B$ , and the output log intensity is

O = c\,B + D - c \max B, \qquad c = \frac{\log_{10} 5}{\max B - \min B},

following the project page’s FAQ, which gives the target as log(5), the compression factor as the target divided by the range of the log base, the offset as the largest compressed base value, and states that all logarithms are in base 10. Writing the output as $O = I + (c - 1)B - c\max B$ (our rearrangement) shows where halos come from: with $c < 1$ , wherever the base is pulled up by a nearby bright region, the output falls, and wherever it is pulled down, the output rises. An edge-preserving base does not cross the edge, so the terms cancel.

Parameters

Key value $a$ (Reinhard). Sets overall brightness; 0.18 by default, and in the paper 0.045 to 0.72 for dark or bright scenes.
White point $L_\text{white}$ (Reinhard, global). Default: the scene maximum, in scaled units. Smaller values burn out highlights deliberately.
Threshold $\epsilon$ and sharpening $\phi$ (Reinhard, local). The paper uses $\epsilon = 0.05$ and $\phi = 8$ . Decreasing $\epsilon$ forces larger scales, and increasing $\phi$ selects slightly larger scales at small $s$ ; both act as edge-enhancement controls (Section 3.2).
Scales (Reinhard, local). Centre-to-surround ratio 1.6, the next scale’s centre equal to the current surround; $\alpha_1 = 1/(2\sqrt{2}) \approx 0.35$ , chosen so that the smallest profile’s two standard deviations overlap one pixel, and $\alpha_2 = 1.6\,\alpha_1$ (Section 4). The authors computed the profile integrals with the error function, since the smallest profile is barely wider than a pixel.
Spatial width $\sigma_s$ (Durand and Dorsey). 2% of the image size; with the Gaussian or Tukey influence functions it had little influence on the result (Section 6).
Range width $\sigma_r$ (Durand and Dorsey). 0.4, which performed well in all their experiments; the authors speculate that it may relate to visual sensitivity or to the range of possible reflectances (Section 6). Since the project page states that all logarithms are base 10, in our reading 0.4 is in $\log_{10}$ units.
Base contrast (Durand and Dorsey). 5 worked well for all their examples; scenes with visible light sources may need a different value. With overall brightness, it is the only parameter the authors expose to the user.

Complexity

The global operator costs a few operations per pixel plus one pass for the log-average; Reinhard et al. (2002, Section 4, Table 1) report it as real-time on a 512 × 512 image on 2002 hardware. The local operator needs a centre and a surround Gaussian convolution per scale, and each scale’s surround is the next scale’s centre (Section 3.2); the authors computed them with the FFT and also tried Burt and Adelson’s multiscale spline approach, the construction behind the image pyramid, about 3.7 times faster than their FFT version, at the cost of blocky artifacts in smooth regions (Section 4, Fig. 14).

A brute-force bilateral filter costs $O(N r^2)$ for $N$ pixels and window radius $r$ . Durand and Dorsey (2002, Sections 5.1–5.2) made it fast with a piecewise-linear approximation in intensity, which turns the filter into a few convolutions that the FFT can compute, and by computing those convolutions on a downsampled image, reaching 1.4 s for a 700 × 480 image on a 700 MHz processor (Fig. 1). The later approximations are described under bilateral filter.

Implementation

The example builds a radiance map with known structure and compares eight operator settings on the same two measures: how much of a known texture survives in each region, and how large a halo appears beside a strong edge.

The scene is 320 × 384 pixels. A room is lit at level 1, a window at 1000 and a shadow at 0.01, and every surface carries a checkerboard reflectance texture of 4-pixel cells that alternate between 1.2 and 0.8. The radiance therefore spans a little over five orders of magnitude, and the texture has the same log contrast, $\log_{10}(1.2/0.8)$ , everywhere. Detail retention in a region is the output’s log contrast between the two kinds of cells divided by the input’s: 1 means the texture is kept, 0 that it is gone. The halo measure averages the log output over 64 rows that cross the window’s left edge, and reports the largest signed departure, within 64 columns of the edge, from the same profile far from the edge, on the room side and on the window side, in $\log_{10}$ units. A second table encodes each display image with the sRGB transfer function from gamma correction, after clipping to $[0, 1]$ , and reports the texture in 8-bit codes, the mean code of each region, and the share of window pixels above 1 before clipping.

The operators are Reinhard et al.‘s $L/(1+L)$ , the curve that the HDR imaging example also computes and which serves here as the baseline; their white-point curve with $L_\text{white}$ at the scene maximum; their dodging-and-burning operator under both readings of scale selection, rule A (the largest scale before $|V|$ first reaches $\epsilon$ ; the smallest if it already fails there) and rule B (the first scale at which $|V| < \epsilon$ ; the largest if none), each with 8 and with 9 scales; and Durand and Dorsey’s operator with a Gaussian base and with a brute-force bilateral base. The scene levels, the texture, the regions and these two conventions for when no scale qualifies were fixed when the example was planned; the operator parameters are the papers’ ( $a = 0.18$ , $\phi = 8$ , $\epsilon = 0.05$ , the $\alpha_i$ and the factor 1.6 from Reinhard et al.; $\sigma_s$ at 2% of the larger side, 7.68 pixels, $\sigma_r = 0.4$ and base contrast 5 from Durand and Dorsey; base-10 logarithms from their FAQ). The spatial window of both Durand bases is truncated at $3\sigma_s$ . Durand and Dorsey’s uncertainty correction (Section 5.3) is not used. Lines starting with C or L are checks fixed in advance; the last block compares the hand-written bilateral filter with scikit-image’s.

import os
for v in ["OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"]:
    os.environ[v] = "1"                       # one BLAS thread, set before NumPy loads
import numpy as np
from scipy.ndimage import gaussian_filter, convolve1d
from scipy.special import erf
from skimage.restoration import denoise_bilateral

# --- 1. Radiance map: three illumination levels times a checker texture -----
H, W = 320, 384
yy, xx = np.mgrid[0:H, 0:W]
illum = np.ones((H, W))                                    # room
illum[(xx >= 192) & (yy >= 48) & (yy < 208)] = 1000.0      # window
illum[(xx >= 16) & (xx < 144) & (yy >= 240)] = 0.01        # shadow
m = 0.2                                                    # reflectance 1 +/- m
hi = ((xx // 4 + yy // 4) % 2) == 0                        # 4-px checker cells
Lw = illum * np.where(hi, 1 + m, 1 - m)
regions = {"room": (slice(16, 80), slice(16, 112)),
           "window": (slice(96, 160), slice(256, 384)),
           "shadow": (slice(272, 320), slice(48, 112))}
band = slice(96, 160)                                      # rows crossing the window edge
d_in = np.log10((1 + m) / (1 - m))
print(f"radiance {Lw.min():.4f} to {Lw.max():.1f}: {np.log10(Lw.max() / Lw.min()):.3f} decades")
for k, R in regions.items():
    print(f"  {k:6s} illumination {np.unique(illum[R])}, cells {hi[R].sum()} bright / {(~hi[R]).sum()} dark")

# --- 2. Measures ---------------------------------------------------------------
def retention(lo, R):          # output log contrast of the texture / input log contrast
    return (lo[R][hi[R]].mean() - lo[R][~hi[R]].mean()) / d_in

def halo(lo):                  # signed largest deviation of the band profile near the edge
    p = lo[band].mean(0)
    dr = p[128:192] - p[0:32].mean()          # room side, against the far room
    dw = p[192:256] - p[320:384].mean()       # window side, against the far window
    return dr[np.argmax(np.abs(dr))], dw[np.argmax(np.abs(dw))]

def srgb(x):
    x = np.clip(x, 0, 1)
    return np.where(x <= 0.0031308, 12.92 * x, 1.055 * x ** (1 / 2.4) - 0.055)

def codes(disp, R):            # 8-bit detail and mean code after clipping and sRGB
    c = np.round(255 * srgb(disp))
    return c[R][hi[R]].mean() - c[R][~hi[R]].mean(), c[R].mean()

lin = np.log10(Lw)
print("C0 input: |1 - retention| max "
      f"{max(abs(1 - retention(lin, R)) for R in regions.values()):.2e}, halo {halo(lin)[0]:.2e} {halo(lin)[1]:.2e}")

# --- 3. Reinhard et al. (2002): global curves, Eqs. 1-4 -----------------------
a, delta = 0.18, 1e-6
Lbar = np.exp(np.mean(np.log(delta + Lw)))   # Eq. 1 with the authors' erratum (1/N inside exp)
L = a / Lbar * Lw                            # Eq. 2
Lwhite = L.max()                             # the paper's default white point
eq3 = lambda x: x / (1 + x)
eq4 = lambda x, Lwh=Lwhite: x * (1 + x / Lwh**2) / (1 + x)
print(f"log-average {Lbar:.4f}, scaled levels room {a/Lbar:.4f} window {1000*a/Lbar:.2f} "
      f"shadow {0.01*a/Lbar:.6f}, L_white {Lwhite:.2f}")
res = {}
for name, f in [("Eq. 3  L/(1+L)", eq3), ("Eq. 4  L_white = max", eq4)]:
    res[name] = (np.log10(f(L)), f(L))
    cf = max(abs(retention(np.log10(f(L)), R)
                 - (np.log10(f(a / Lbar * illum[R][0, 0] * (1 + m)))
                    - np.log10(f(a / Lbar * illum[R][0, 0] * (1 - m)))) / d_in) for R in regions.values())
    print(f"C5 {name}: |image - closed form| retention max {cf:.2e}")
print(f"C4 Eq. 4 at L_white {eq4(Lwhite):.15f}; L_white=1e12 vs Eq. 3 max |diff| "
      f"{np.abs(eq4(L, 1e12) - eq3(L)).max():.2e}")

# --- 4. Reinhard et al. (2002): dodging-and-burning, Eqs. 5-9 ------------------
phi, eps = 8.0, 0.05
a1 = 1 / (2 * np.sqrt(2)); a2 = 1.6 * a1

def taps(sig):                 # Eq. 5 profile integrated over each pixel (erf), 1-D factor
    rad = int(np.ceil(4 * sig)) + 1
    k = np.arange(-rad, rad + 1)
    return 0.5 * (erf((k + 0.5) / (sig * np.sqrt(2))) - erf((k - 0.5) / (sig * np.sqrt(2))))

def blur(img, sig):
    t = taps(sig); t = t / t.sum()
    return convolve1d(convolve1d(img, t, axis=0, mode="reflect"), t, axis=1, mode="reflect")

def dodge_burn(L, nscales, rule):
    V1s, Vs = [], []
    for i in range(nscales):
        s = 1.6 ** i
        V1 = blur(L, a1 * s / np.sqrt(2)); V2 = blur(L, a2 * s / np.sqrt(2))   # Eq. 6
        V1s.append(V1); Vs.append((V1 - V2) / (2**phi * a / s**2 + V1))       # Eq. 7
    ok = np.abs(np.stack(Vs)) < eps                                           # Eq. 8
    if rule == "A":   # largest scale before |V| first reaches eps
        first_bad = np.where((~ok).any(0), (~ok).argmax(0), nscales)
        sm = np.maximum(first_bad - 1, 0)
    else:             # "B": the first scale at which |V| < eps
        sm = np.where(ok.any(0), ok.argmax(0), nscales - 1)
    V1 = np.take_along_axis(np.stack(V1s), sm[None], 0)[0]
    return L / (1 + V1)                                                       # Eq. 9

print("C7 tap sums before normalising, min over s=1, 1.6^8 and both alphas: "
      f"{min(taps(al * s / np.sqrt(2)).sum() for s in (1, 1.6**8) for al in (a1, a2)):.6f}")
const = np.full((H, W), 3.7)
for ns in (8, 9):
    for rule in "AB":
        name = f"local rule {rule}, {ns} scales"
        dev = np.abs(dodge_burn(const, ns, rule) - eq3(const)).max()
        print(f"C6 {name}: constant image vs Eq. 3 max |diff| {dev:.2e}")
        Ld = dodge_burn(L, ns, rule)
        res[name] = (np.log10(Ld), Ld)

# --- 5. Durand and Dorsey (2002): base-detail split on log10 luminance ---------
ss = 0.02 * max(H, W); r = int(np.ceil(3 * ss)); sr = 0.4

def bilateral(img, ss, sr, r):  # brute force over a (2r+1)^2 window
    P = np.pad(img, r, mode="symmetric")
    num = np.zeros_like(img); den = np.zeros_like(img)
    for dy in range(-r, r + 1):
        for dx in range(-r, r + 1):
            q = P[r + dy:r + dy + img.shape[0], r + dx:r + dx + img.shape[1]]
            w = np.exp(-(dy * dy + dx * dx) / (2 * ss * ss) - (q - img) ** 2 / (2 * sr * sr))
            num += w * q; den += w
    return num / den

lI = np.log10(Lw)
B_gauss = gaussian_filter(lI, ss, mode="reflect", truncate=r / ss)
B_bilat = bilateral(lI, ss, sr, r)
print(f"C3 bilateral(sigma_r=1e6) vs Gaussian base max |diff| "
      f"{np.abs(bilateral(lI, ss, 1e6, r) - B_gauss).max():.2e}")
for name, B in [("Durand, Gaussian base", B_gauss), ("Durand, bilateral base", B_bilat)]:
    c = np.log10(5) / (B.max() - B.min())          # base contrast 5 (log10 5)
    out = c * B + (lI - B) - c * B.max()           # largest base value -> 1
    res[name] = (out, 10 ** out)
    print(f"{name}: compression factor c = {c:.4f}")

# --- 6. Results -------------------------------------------------------------------
print("\nC1 halo of the global curves:",
      " ".join(f"{h:.2e}" for n in list(res)[:2] for h in halo(res[n][0])))
print("\n                         detail retention (log)       halo (log10)")
print("                          room  window  shadow    room side  window side")
for name, (lo, disp) in res.items():
    hr, hw = halo(lo)
    print(f"{name:24s} " + " ".join(f"{retention(lo, R):6.4f}" for R in regions.values())
          + f"    {hr:+.4f}     {hw:+.4f}")
print("\n                         8-bit detail (codes)      mean code            window")
print("                          room  window shadow     room window shadow   share > 1")
for name, (lo, disp) in res.items():
    cs = [codes(disp, R) for R in regions.values()]
    print(f"{name:24s} " + " ".join(f"{d:6.2f}" for d, _ in cs) + "   "
          + " ".join(f"{mc:6.1f}" for _, mc in cs) + f"     {(disp[regions['window']] > 1).mean():.3f}")

# --- 7. Library check: scikit-image 0.25.2 denoise_bilateral ---------------------
off = lI.min()                                    # its colour table assumes values >= 0
B_sk = denoise_bilateral(lI - off, win_size=2 * r + 1, sigma_color=sr, sigma_spatial=ss,
                         bins=100000, mode="symmetric") + off
inner = (slice(r, H - r), slice(r, W - r))
print(f"\nL1 skimage vs hand bilateral, interior max |diff| {np.abs(B_sk - B_bilat)[inner].max():.2e} (log10)")
win = 5
imp = np.ones((21, 21)); imp[10, 10] = 2.0
resp = denoise_bilateral(imp, win_size=win, sigma_color=1e3, sigma_spatial=1.0, mode="constant", cval=1.0) - 1.0
wts = np.array([[resp[10 - dy, 10 - dx] for dx in range(-2, 3)] for dy in range(-2, 3)])
g = np.arange(-win // 2, win // 2 + 1)                      # the grid of _compute_spatial_lut
table = np.exp(-0.5 * np.hypot(*np.meshgrid(g, g, indexing="ij")) ** 2).ravel()
as_read = table[:win * win].reshape(win, win)               # read as range_lut[kr * win_size + kc]
k1 = np.exp(-0.5 * np.arange(-2, 3) ** 2); centred = np.outer(k1, k1)
print(f"L2 impulse weights: table has {table.size} entries for a {win}x{win} window; "
      f"vs table read with stride {win}: {np.abs(wts - as_read / as_read.sum()).max():.2e}; "
      f"vs centred Gaussian: {np.abs(wts - centred / centred.sum()).max():.4f}")
print("   weight at offset (0, 0):", f"{wts[2, 2]:.4f}", "centred Gaussian:", f"{centred[2, 2] / centred.sum():.4f}")

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; identical on two runs; about 4 s):

radiance 0.0080 to 1200.0: 5.176 decades
  room   illumination [1.], cells 3072 bright / 3072 dark
  window illumination [1000.], cells 4096 bright / 4096 dark
  shadow illumination [0.01], cells 1536 bright / 1536 dark
C0 input: |1 - retention| max 1.55e-15, halo 0.00e+00 0.00e+00
log-average 3.7538, scaled levels room 0.0480 window 47.95 shadow 0.000480, L_white 57.54
C5 Eq. 3  L/(1+L): |image - closed form| retention max 7.55e-15
C5 Eq. 4  L_white = max: |image - closed form| retention max 1.22e-15
C4 Eq. 4 at L_white 1.000000000000000; L_white=1e12 vs Eq. 3 max |diff| 0.00e+00
C7 tap sums before normalising, min over s=1, 1.6^8 and both alphas: 0.999959
C6 local rule A, 8 scales: constant image vs Eq. 3 max |diff| 0.00e+00
C6 local rule B, 8 scales: constant image vs Eq. 3 max |diff| 0.00e+00
C6 local rule A, 9 scales: constant image vs Eq. 3 max |diff| 0.00e+00
C6 local rule B, 9 scales: constant image vs Eq. 3 max |diff| 0.00e+00
C3 bilateral(sigma_r=1e6) vs Gaussian base max |diff| 3.60e-12
Durand, Gaussian base: compression factor c = 0.1398
Durand, bilateral base: compression factor c = 0.1393

C1 halo of the global curves: 0.00e+00 0.00e+00 0.00e+00 -4.34e-19

                         detail retention (log)       halo (log10)
                          room  window  shadow    room side  window side
Eq. 3  L/(1+L)           0.9549 0.0210 0.9995    +0.0000     +0.0000
Eq. 4  L_white = max     0.9549 0.0351 0.9995    +0.0000     -0.0000
local rule A, 8 scales   1.0000 0.7554 1.0000    -0.3076     +0.0154
local rule B, 8 scales   0.9569 0.0657 0.9995    -1.3281     +0.0898
local rule A, 9 scales   1.0000 0.7554 1.0000    -0.3076     +0.0154
local rule B, 9 scales   0.9569 0.0657 0.9995    -1.3281     +0.0898
Durand, Gaussian base    1.0000 1.0000 1.0000    -1.2232     +1.2232
Durand, bilateral base   0.9583 0.9583 0.9583    -0.0000     +0.0000

                         8-bit detail (codes)      mean code            window
                          room  window shadow     room window shadow   share > 1
Eq. 3  L/(1+L)            12.00   1.00   1.00     60.0  252.5    1.5     0.000
Eq. 4  L_white = max      12.00   2.00   1.00     60.0  254.0    1.5     0.000
local rule A, 8 scales    12.00  19.90   1.00     60.0  244.8    1.5     0.373
local rule B, 8 scales    12.00   2.98   1.00     60.0  252.5    1.5     0.000
local rule A, 9 scales    12.00  19.90   1.00     60.0  244.8    1.5     0.373
local rule B, 9 scales    12.00   2.98   1.00     60.0  252.5    1.5     0.000
Durand, Gaussian base     31.00  22.00  24.00    166.5  244.0  124.0     0.500
Durand, bilateral base    29.00  21.00  22.00    166.5  244.5  124.0     0.500

L1 skimage vs hand bilateral, interior max |diff| 7.07e-03 (log10)
L2 impulse weights: table has 36 entries for a 5x5 window; vs table read with stride 5: 8.81e-08; vs centred Gaussian: 0.2150
   weight at offset (0, 0): 0.0015 centred Gaussian: 0.1621

The checks held. The texture measure gives exactly 1 and the halo measure exactly 0 on the input itself (C0), and the global curves have no halo, to rounding (C1): a point operation maps every copy of a radiance level to the same output, so both sides of the edge look as they do far from it. The scene has the planned levels, range and equal numbers of bright and dark cells in every region (C2). The retention that the image gives for the global curves matches the closed form from the two radiance levels of each region (C5); the white-point curve sends $L_\text{white}$ to 1 and becomes $L/(1+L)$ when $L_\text{white}$ is huge (C4); the integrated Gaussian taps lose less than $10^{-4}$ of their mass to truncation (C7); the dodging-and-burning operator gives exactly $L/(1+L)$ on a constant image under all four variants, as its normalized kernels require (C6); and the hand-written bilateral filter with $\sigma_r = 10^6$ reproduces SciPy’s Gaussian base (C3).

Global curves. The log-average puts the room at a scaled level of 0.048, the window at 47.95 and the shadow at 0.00048. The slope $1/(1+L)$ predicts the retention: about $1/1.048 \approx 0.95$ in the room, $1/48.95 \approx 0.02$ in the window and almost exactly 1 in the shadow, close to the printed values (the finite ±20% texture shifts them slightly). The planned outcome, below 0.05 in the window, above 0.99 in the shadow and between 0.90 and 1.00 in the room, held for both curves. The white point gives the window a little more contrast, but in 8-bit codes both curves leave the window’s texture 1 or 2 codes apart and the shadow at a mean code of 1.5, where its texture is 1 code; the planned limits (at most 2 codes in the window, at most 1 in the shadow, mean below 5) held. In log terms the shadow keeps its detail; on the display it is black. That is the case for local operators, which the HDR imaging example reaches with stripes in a room and a window.

Dodging-and-burning. Under rule A, the window’s retention rises from 0.0351 to 0.7554, more than 0.5 and about 21 times the white-point curve’s, and the room and shadow keep 1.0000; in codes the window texture spans 19.90 codes. These planned outcomes held under both scale sets. Rule A’s window retention stays below 1, and in our reading that is because the selected scale varies from pixel to pixel, so that at some pixels $V_1$ still follows the texture. Under rule B the window keeps only 0.0657, below the planned 0.1 (held), close to the global curve: in our reading the smallest scale already passes the threshold in smooth textured areas, so $V_1 \approx L$ and the operator reduces to $L/(1+L)$ , the failure Reinhard et al. describe for too small a scale. The 8- and 9-scale rows are identical as printed, so on this scene the disagreement about the number of scales does not matter (planned to differ by less than 0.01 in every retention and halo value, held). The shadow is a different matter under every variant: its mean code stays 1.5 and its texture 1 code (planned at most 1 code and a mean below 5, held). Dividing by $1 + V_1$ changes nothing where $V_1 \ll 1$ , so the operator brightens no dark region; it compresses bright surroundings.

The dodging-and-burning operator does produce a halo on the room side of the window: $-0.3076$ under rule A, smaller in size than the Gaussian base’s (held), and $-1.3281$ under rule B, larger than it (held); the window side stays within $0.05$ under rule A ( $+0.0154$ , held). A construction argument accounts for rule A’s value, in our reading at the column next to the edge. There the smallest centre profile, with standard deviation $\alpha_1/\sqrt{2} = 0.25$ pixels, puts $(1 - \operatorname{erf}(0.5/(0.25\sqrt{2})))/2 \approx 2.3\%$ of its weight on the window, which is 1000 times brighter, so $V_1 \approx 0.048 \times (1 + 0.023 \times 999) \approx 1.14$ ; the surround, with standard deviation 0.4 pixels, puts about 10.6% there, so $V_2 \approx 0.048 \times (1 + 0.106 \times 999) \approx 5.1$ , and at $s = 1$ , $V \approx (1.14 - 5.1)/(2^8 \times 0.18 + 1.14) \approx -0.08$ . Even the smallest scale fails the threshold, so rule A uses it, and the output falls by $\log_{10}(2.14/1.048) \approx 0.31$ , as printed. Rule B’s larger halo comes, in our reading, from a larger profile being selected next to the edge.

Base and detail. Both Durand variants keep the texture everywhere, as planned: retention between 0.95 and 1.05 in all three regions held for the Gaussian base (1.0000) and for the bilateral base (0.9583). A Gaussian of 7.68 pixels removes a 4-pixel checkerboard almost completely from the base, so the texture goes into the detail layer untouched. The bilateral base keeps a little of it, because its range weights favour neighbours on cells of the same colour (our reading); that texture is then compressed with the base, and the bilateral variant retains less than the Gaussian one, as planned (held). In codes, both show 21 to 31 codes of texture in every region and lift the shadow to a mean code of 124 (planned at least 15 codes and a mean above 100, held). The window’s largest base value is mapped to 1, so the bright cells there exceed 1 and half the window clips (share 0.500, about half as planned, held); Reinhard’s rule A clips part of the window too.

The difference between the bases is the halo. With a Gaussian base the room side drops by $1.2232$ in $\log_{10}$ , a factor of about 17, and the window side rises by the same amount (planned beyond $\mp 0.3$ , held). The rearranged formula predicts it: in the column next to the edge, a Gaussian of 7.68 pixels puts about 47% of its weight on the window, which is 3 decades brighter, so the base rises by about $3 \times 0.47 = 1.42$ , and the output changes by $(c - 1) \times 1.42 = (0.1398 - 1) \times 1.42 \approx -1.22$ . The bilateral base gives no measurable halo at four decimals (planned within 0.02, held): a neighbour three decades away has a range weight of $e^{-(3 - 0.18)^2 / (2 \times 0.4^2)}$ , about $10^{-11}$ even for the closest texture values (our arithmetic), so the base does not cross the edge. Durand and Dorsey (2002, Section 1) use an edge-preserving filter for the base in order to avoid such halo artifacts, and the survey by Paris, Kornprobst, Tumblin and Durand (2009, Fig. 3.4) shows the same contrast on a photograph.

scikit-image’s bilateral filter. Check L1, agreement between scikit-image 0.25.2’s denoise_bilateral and the hand-written filter within $10^{-3}$ in the interior, failed, as planned: the two differ by $7.07 \times 10^{-3}$ . In our reading, the source explains it. The Python helper _compute_spatial_lut builds the spatial weights on np.arange(-win_size // 2, win_size // 2 + 1), which for an odd window of 5 runs from −3 to 2: six values per axis, 36 entries, not centred. The Cython loop in _denoise_cy.pyx reads the table as range_lut[kr * win_size + kc], with kr and kc from 0 to win_size - 1. The impulse response measured in L2 matches that table read row by row with a stride of 5 to $8.8 \times 10^{-8}$ (held, criterion $10^{-6}$ ), and differs from a centred Gaussian by 0.2150 (held, criterion more than 0.1); the weight at the centre is 0.0015 instead of 0.1621. That is the library’s behaviour at version 0.25.2, as printed; the example’s Durand base uses the hand-written filter.

Every design value above was fixed when the example was planned, and nothing changed after the first run.

Properties and Behavior

Global curves cannot create halos. A point operation sends equal inputs to equal outputs, so a uniform region looks the same next to an edge and far from it (C1 above). The price is that detail survives only where the curve is steep; with $L/(1+L)$ that is the dark end in log terms, and on a display the dark end may still be black.
Halos come from a base that crosses edges. Durand and Dorsey (2002, Section 2) call haloing the major pitfall of local methods and note that it becomes more critical with high dynamic range: in 8-bit images edge contrast is limited to about two orders of magnitude, which limits the strength of halos. Chiu et al.’s low-pass gain produced pronounced halos, and multiscale methods reduce but do not always remove them, in Durand and Dorsey’s review. The example’s Gaussian base shows a halo that grows with the size of the edge’s step in log radiance.
Edge-preserving bases leave small halos of their own. Durand and Dorsey (2002, Section 5.3) observed small halos at antialiased edges and around flare, where few neighbours fall inside the range kernel and the base is a statistical estimate from little data; they use the bilateral normalization factor $k$ as a measure of uncertainty and blend the output toward a slightly blurred version where it is small. Amplifying the detail layer can also produce gradient reversals next to edges, as the bilateral filter article describes.
Scale selection is the core of dodging-and-burning. Too small a scale reduces the local operator to the global one, and too large a scale forms dark rings around bright areas (Reinhard et al., 2002, Section 3.2, Fig. 9). The example shows both the reduction (rule B) and a dark fringe beside the edge even under rule A.
The log domain. Paris et al. (2009, Section 3.2.7) explain the choice of logarithms for contrast management: the range kernel then acts on multiplicative contrast, the same everywhere in the range. The cost is that zero maps to minus infinity and noise in dark regions is amplified, and many users add a small constant on the order of the noise level before taking the logarithm.
Measurement against picture. Every operator here is nonlinear and local operators vary across the image, so tone-mapped values are no longer proportional to radiance. HDR imaging discusses why a vision system should usually keep the linear radiance map and tone-map only the copy people look at.

Limitations

No objective target. Reinhard et al. (2002, Section 4) compared operators by a “knock-out race” of progressively harder images, for lack of well-tested quantitative methods; Durand’s project page calls its own comparison of the three SIGGRAPH 2002 operators informal. The measures in the example quantify detail and halos on a synthetic scene; they do not say which picture looks better.
Parameters are image-dependent. Reinhard et al. set the key per image; Durand and Dorsey note that scenes with visible light sources may need another base contrast; the compression factor and offset of Durand and Dorsey’s method are computed from each image, so they change from frame to frame in a video unless they are smoothed (our reasoning).
Clipping of the brightest region. Both local operators in the example clip part of the brightest region: Reinhard et al. by design, to let light sources burn out, and Durand and Dorsey because the largest base value maps to white while detail rides on top of it.
Colour. Operating on luminance and recomposing colour by ratios keeps hue but can oversaturate compressed regions (our reasoning); Durand and Dorsey (2002, Section 7) list white balance for scenes mixing indoor and outdoor light as open.
Not a measurement. A tone-mapped image should not feed methods that assume brightness proportional to light, such as photometric stereo or brightness constancy across frames, as digital camera explains for camera tone curves.

Variants

Other global curves. The S-shaped functions of Schlick and of Tumblin et al., in Durand and Dorsey’s review (Section 2), keep some highlight and shadow detail but reduce contrast severely there. Ward Larson et al.’s histogram adjustment, also in their review, discards the empty parts of the histogram, a relative of histogram equalization on the image histogram; its global nature shows when the histogram is uniform.
Multiscale and perceptual local operators. Durand and Dorsey review Jobson et al.’s multiscale retinex, Pattanaik et al.’s psychophysically derived filter banks, and Tumblin and Turk’s low-curvature image simplifier, an anisotropic-diffusion decomposition that extracts fine detail but is slow and needs per-image coefficients. Reinhard et al. compare their operator with homomorphic filtering, which suppresses low frequencies and enhances high ones.
Other base filters. Durand’s project page notes that the bilateral base can serve as the adaptation level of most global operators, and that the photographic local operator, which blurs but stops at edges, could use the bilateral filter instead. Edge-preserving alternatives to the bilateral filter, such as the guided filter, are described under bilateral filter.
Gradient-domain compression. Fattal et al. (2002) published a gradient-domain method at the same SIGGRAPH as the two operators above, as Durand’s project page notes. Durand’s informal impression there is that it preserves detail best of the three, “maybe sometimes too good”, and that the photographic operator best gives a photographic appearance.
Exposure fusion. Blending the bracketed exposures directly, without a radiance map, produces a display image in one step; HDR imaging describes it, together with the Laplacian-pyramid blending it relies on.
Inverse tone mapping. The reverse problem, expanding one low-dynamic-range image into a plausible radiance map, is covered under single-image HDR reconstruction in HDR imaging.
In cameras. Cameras tone-map inside their image signal processors. The digital camera article models the step as a global curve and notes that real cameras may use edge-aware local tone mapping; computational photography describes HDR+, a burst pipeline that tone-maps its merged frames with a local operator.