Color Constancy

Color constancy, in computer vision, is the problem of estimating the colour of the light that illuminated a scene from an image of it, usually so that the light’s colour can be removed. A person in a room lit by coloured light unconsciously “discounts the illuminant”, in Helmholtz’s words as Barron (2015, Section 1) quotes them, and sees objects as if under white light; a camera that is to give the same impression has to infer the light and divide it out. Barron and Tsai (2017, Section 1) point to Edwin Land, who invented both Retinex theory and the Polaroid instant camera, as an example of the historical link between colour constancy and photography; Retinex as an image-enhancement filter is described in Frequency Domain Filtering. Barron (2015, Section 1) notes that the problem of inferring and correcting the illumination colour is commonly called “color constancy” or “white balance”; in cameras it runs as auto white balance (AWB). Barron and Tsai (2017, Section 1) draw a distinction between the two: colour constancy aims to recover the true scene behind the image, white balance to give the image a pleasant appearance, and for learned methods the difference comes down to the training data. The task is also called illuminant estimation, which is the half of the problem that is hard: once the light is known, removing it is a division (Digital Camera describes that step in the camera pipeline).

The problem is underconstrained. Barron (2015, Section 1) puts it as a question: a yellow pixel could be a white object under yellow light or a yellow object under white light. Every method therefore adds an assumption about the world, from “surfaces average to grey” to a model learned from thousands of labelled photographs, and each fails where its assumption fails. This article covers the image formation model behind the task, the classical statistical estimators and their unified form, gamut and probabilistic methods, the learned estimators that replaced them on benchmarks, the angular error by which all of them are measured, and a worked example on synthetic scenes where the assumption behind each estimator can be printed next to its error.

Problem Definition

Barron (2015, Section 2) states the model the field works with. Take a photometrically linear image $I$, black-level corrected, with no saturated pixels. Each pixel’s RGB value is the product of its “true” white-balanced value $W$ and an illumination $L$ shared by all pixels:

$$I = W \times L,$$

channel by channel. The goal is to estimate $L$ from $I$ and produce $W = I / L$. The per-channel product is the von Kries or diagonal model: correcting the illuminant means scaling each channel by its own gain. Barron (2015, Section 2) calls it an approximation, though an effective and widespread one, and a severe oversimplification of the full intrinsic-image problem, since it ignores shading, reflectance properties and spatially varying illumination. Global Image Descriptors uses the same diagonal model to build illuminant-invariant colour histograms, the other way of coping with an unknown light: describe the image in quantities that do not depend on it.

The magnitude of $L$ cannot be recovered. Finlayson (2013, Section 1) counts the unknowns: an image of $N$ pixels under one light has $3N + 3$ unknowns ($N$ surfaces and the light) and $3N$ measurements, and since a bright scene dimly lit cannot be told from a dim scene brightly lit, the unknowns reduce to $3N + 2$. Only the light’s chromaticity is sought. Barron (2015, Eqs. 2–6) makes this explicit by working in log-chrominance, $u = \log(g/r)$ and $v = \log(g/b)$, where the illuminant is two numbers $(L_u, L_v)$ and the image formation model becomes additive:

$$W_u = I_u - L_u, \qquad W_v = I_v - L_v .$$

A change of illuminant therefore translates every pixel’s log-chrominance by the same amount. After estimating $(L_u, L_v)$, Barron recovers an RGB illuminant by assuming $L$ has unit norm.

Input and Output

The input is a single image, or in a camera a stream of frames. Both Barron (2015, Section 2) and Barron and Tsai (2017, Section 2) assume a linear image with a black level of zero and no saturated pixels; Barron and Tsai add in a footnote that in practice saturated pixels are identified and removed from all later computation, in the same way that colour-checker pixels are ignored in benchmarks. Linear means raw or raw-like data: in a camera the estimate is made before the tone curve and gamma, from the values that Image Sensor and Digital Camera describe. An estimator applied to an sRGB JPEG first has to undo its transfer function (given in Color Spaces), and even then the camera’s own white balance and colour processing have already acted on it.

The output is an illuminant colour: an RGB triple up to scale, or a chromaticity. Barron and Tsai (2017, Section 1) list what a practical AWB needs besides accuracy: speed (they argue for 1–5 ms per frame, a small fraction of a phone’s compute budget), robustness to small low-bit-depth preview images (32 × 24 or 64 × 48 pixels, 8-bit), a confidence measure or a full posterior over illuminants, and temporal coherence so that colours in a video do not vary erratically. FFCC outputs a mean and a covariance for this reason.

Variants

• Single global illuminant. The standard formulation and the one every benchmark below scores. FC4 states that while a scene can contain several illuminants, it addresses the traditional problem of a single global illumination colour (Hu et al., 2017, Section 3).
• Multiple illuminants. Indoor scenes with daylight through a window, or flash with ambient light, have different illuminant colours in different regions. No single diagonal correction can be right everywhere; the example below shows the bound this sets. Tone Mapping records that Durand and Dorsey (2002) listed white balance for scenes mixing indoor and outdoor light as an open problem.
• Video and viewfinder AWB. The estimate must be smoothed over time. Barron and Tsai (2017, Section 7) update a running Gaussian estimate with each frame’s posterior, a Kalman-filter-like step whose single parameter sets the expected variance of the illuminant over time (the Kalman filter is the general form).
• Camera-specific or cross-camera. Raw values depend on each sensor’s spectral sensitivities, so an estimator learned on one camera’s raw images does not transfer directly to another’s; Barron (2015, Section 6) learns one model per camera on the multi-camera dataset of Cheng et al. (2014).

Challenges

• Ambiguity of surface and light. Colour constancy is underconstrained by construction. FC4’s opening figure makes the point with patches: a uniform wall can be a yellow wall in white light or a white wall in yellow light, so a broad range of illuminants explains it, whereas a banana, which has an innate colour, does (Hu et al., 2017, Section 1).
• Scenes that break the prior. A statistic that assumes average grey fails when one colour dominates. Finlayson (2013, Section 1, Fig. 1) shows a scene with green leaves where dividing by the mean colour wrongly maps green toward grey.
• Clipping. Saturated pixels no longer obey $I = W \times L$. Finlayson (2013, Section 4) reports that the max-RGB estimate is known to perform less well with extensive clipping, and notes that a dataset in HDR format was compiled so that its data are not clipped.
• Nonlinear inputs. The estimators assume values proportional to light; a tone curve or gamma breaks the per-channel product structure. FC4, which fine-tunes networks pretrained on gamma-encoded ImageNet images, applies a gamma of $1/2.2$ to its linear inputs to make them resemble its pretraining data (Hu et al., 2017, Section 4.1), which shows that the choice is a design decision for learned methods too.
• Small datasets. Labelled data require a colour chart in each scene. Barron (2015, Section 6) considers his larger gain on the Gehler–Shi set likely due to its larger training folds (about 379 images against about 144), and his conclusion (Section 7) suggests that colour constancy may benefit from much larger datasets.

Approaches

Statistical estimators and their unified form

The oldest and simplest methods compute one statistic of the image and take it as the illuminant.

Gray world takes the mean of each channel. Finlayson (2013, Section 1) explains the name: if the average surface colour in the scene is grey, the mean RGB must have the illuminant’s colour, under the per-channel product model. Barron (2015, Section 1) and FC4 (Hu et al., 2017, Section 2) both describe it as assuming that reflectances average to grey; Barron cites Buchsbaum (1980) for it. Finlayson (2013, Section 1) adds the case for it: a strongly coloured light biases the first-order moments toward its own colour, so under a very yellow light the image mean is almost certainly yellow.

White patch, also called max-RGB, takes the maximum response in each channel (Finlayson, 2013, Section 4, calls it MaxRGB; Barron’s Table 1 cites Brainard and Wandell’s analysis of Retinex theory for it). It follows from the diagonal model that the estimate is exact when some surface in the scene reflects each channel as strongly as a white surface would, the case that check K2 in the example tests.

Shades of gray (Finlayson and Trezzi, 2004, as Barron 2015 cites it) replaces the mean by a Minkowski $p$-norm, and gray-edge (van de Weijer, Gevers and Gijsenij, 2007) applies the same statistic to image derivatives rather than values. Finlayson (2013, Eq. 6) summarizes the family, after van de Weijer et al., in one equation:

$$\left( \int \left| \frac{\partial^n \rho(\mathbf{x})}{\partial \mathbf{x}^n} \right|^p d\mathbf{x} \right)^{1/p} = k\, \hat{\rho}^{E}_{n,p,\sigma},$$

where $\rho(\mathbf{x})$ is the camera’s RGB response at location $\mathbf{x}$, first smoothed with a Gaussian of standard deviation $\sigma$, $n$ is the order of differentiation ($n = 0$ means none), and $k$ is the unknown scale of the light. In this notation:

• $n = 0$, $p = 1$ is gray world;
• $n = 0$, $p \to \infty$ is white patch, since the $p$-norm of nonnegative values tends to their maximum;
• $n = 0$, $1 < p < \infty$ is shades of gray;
• $n = 1$ and $n = 2$ are first- and second-order gray-edge, on the gradient magnitude and second derivatives.

Gray-edge reasons about differences between neighbouring surfaces rather than surfaces themselves. It follows from the diagonal model that each channel’s derivative is the illuminant’s gain for that channel times the derivative of the reflectance, so the illuminant still scales the statistic channel by channel. Finlayson (2013, Section 2) notes that arbitrary linear filters can replace the derivative operators. All members of the family share one property, which the example checks numerically: scaling a channel of the image by a gain scales the estimate’s channel by the same gain, so they are exact on any scene for which their statistic of the reflectance is grey.

Gamut mapping

Gamut mapping, for which Barron (2015, Section 1) cites Forsyth (1990) among others, uses calibration rather than a single statistic. As Finlayson (2013, Section 2) describes it, it learns the gamut of camera RGBs that surfaces can produce under a known reference light, and looks for an illuminant that, once divided out, puts all of the image’s colours inside that gamut. The constraint usually leaves a set of feasible lights rather than one, and a single answer is selected by a statistical argument; gamut mapping can also work on colour derivatives. In Barron’s Table 1 the pixel-based and intersection-based gamut variants have medians of 2.33° and 2.39°, well below those of gray world and white patch (6.28° and 5.68°), and means of 4.20° against their 6.36° and 7.55°.

Probabilistic and corrected-moment methods

Adding probabilities to the gamut idea gives Bayesian colour constancy: with likelihoods of RGBs under each light and a prior over lights, the illuminant is estimated by maximum likelihood or a posterior (Finlayson, 2013, Section 2). The reprocessed Color Checker dataset that benchmarks still use came from Gehler et al.’s (2008) Bayesian Color Constancy Revisited (not read; its role as the source of the images is from Barron, 2015, Section 6, and Hu et al., 2017, Section 4.1).

Finlayson (2013) argued that the simple moments are useful but biased, and corrected them with a learned matrix. His corrected gray world is $\hat{\rho}^E \approx C \underline{\mu}$, with $\underline{\mu}$ the mean RGB and $C$ a $3 \times 3$ matrix (Eq. 5); the general form maps a vector of $M$ moments, including second- and higher-order moments with cross-channel terms and moments of colour-edge images, through an $M \times 3$ matrix. Because the light’s intensity is unknown, he fits $C$ to map the orientation of the moments onto the orientation of the true illuminant, which turns a least-squares problem into an alternating least-squares procedure (Section 3). In Barron’s Table 1 the corrected-moment variants were the strongest prior method: 2.8° mean and 2.0° median for the 19-edge-moment version, as reported in Finlayson’s paper, with Barron’s own reruns of the authors’ code slightly different (2.96° to 3.12° mean).

Learned estimators

Barron (2015) reframed colour constancy as detection. Because a change of illuminant translates the log-chrominance histogram, scoring every possible illuminant with a linear function of the histogram is a convolution, and the estimate is the location of the maximum response. Convolutional Color Constancy (CCC) builds a luminance-weighted log-chroma histogram (Eq. 7), convolves it with a learned filter and takes the argmax (Eq. 9). The filter is trained discriminatively with a softmax over histogram bins and a loss equal to the angle between each candidate and the true illuminant (Eqs. 10–11). Barron contrasts this with the generative framing of earlier statistical and learned methods, which model what white-balanced images look like and then use that model to correct new ones: CCC instead learns to tell well-balanced from badly balanced versions of the same image. His extended model adds histograms of three augmented images (a sharpened, rectified image, a soft max filter and a local standard deviation), each computed channel by channel and preserving scalar multiplication (Eq. 14), so that the translation property still holds (Section 5). In his experiments, discriminative training improved on generative training by 30–40% and the extended features on the simple ones by 10–20% (Section 6).

Fast Fourier Color Constancy (FFCC) (Barron and Tsai, 2017) makes CCC fast. It wraps the log-chroma histogram onto a small torus (64 × 64 bins of width 1/32), performs the convolution with FFTs, and fits a bivariate von Mises distribution to the softmax output to get a mean and covariance (Sections 3–4). Wrapping makes the estimate ambiguous up to the torus period, a problem the authors call illuminant aliasing; they resolve it either with the image’s mean chroma (“gray world de-aliasing”) or by placing the histogram so that all plausible illuminants fall inside one period (“gray light de-aliasing”). A learned per-illuminant gain and bias map lets the model prefer some illuminants over others, which a torus-shift-invariant convolution cannot do on its own (Section 5). The abstract claims 13–20% lower error than the previous state of the art and a 250–3000× speed-up; on a 2016 Google Pixel XL, a Halide implementation on the camera’s thumbnails ran at 1.44 ms per image (Section 8). Computational Photography notes the learned white balance that Google’s Night Sight adds for night scenes.

FC4 (Hu, Wang and Lin, 2017) is a fully convolutional network built on AlexNet or SqueezeNet layers pretrained on ImageNet, which outputs a local illuminant estimate and a confidence for every region of the image and pools them by a confidence-weighted sum (Eqs. 2–4). The authors show that the pooling’s gradients reach each local estimate in proportion to its confidence (Eq. 5), so the network learns which regions to trust, such as faces and objects with an innate colour, without extra supervision. Its loss is the angular error (Eq. 1). The abstract claims better accuracy than the previous state of the art at 120× greater efficiency.

These methods are built on convolutional neural networks, softmax outputs and, for FC4, transfer learning; Barron’s framing borrows directly from object detection.

Example: Four Estimators on Synthetic Linear Mondrians

The example compares gray world, white patch, shades of gray ($p = 6$) and first-order gray-edge ($n = 1$, $\sigma = 1$ pixel, $p = 6$) on synthetic “Mondrian” scenes: 128 × 128 images made of a background and 60 rectangles of uniform reflectance. All parameters were fixed before any run and not tuned.

Why the data are linear. The radiance is generated as $E = R \times L$, reflectance times illuminant per channel, so the diagonal model holds exactly. The sensor adds Gaussian noise of standard deviation 0.005 and clips to $[0, 1]$, and nothing else: no tone curve, no gamma. Values are therefore proportional to light except for noise and clipping, and the code prints a check that doubling the light doubles the noise-free radiance. A value of 1 is saturation. Every estimator ignores pixels with any saturated channel, as Barron and Tsai (2017, Section 2) describe; gray-edge also ignores pixels within three 4-connected steps of a saturated one (three binary erosions), because its Gaussian derivative mixes values from a neighbourhood; SciPy’s kernel for $\sigma = 1$ reaches 4 pixels, so the margin removes most but not all of that neighbourhood.

Two reflectance models. The ranking of these estimators depends on how surface colours are distributed, which the sources leave open, so the results are printed under two models fixed in advance:

• M1, independent uniform: each channel of each surface is drawn independently from $U(0.05, 0.9)$.
• M2, smooth spectra: each surface has a reflectance spectrum on 400–700 nm, the sigmoid of a random smooth function (a constant plus three cosine and three sine harmonics), seen through Gaussian red, green and blue sensitivities centred at 610, 540 and 460 nm (standard deviation 30 nm) and normalized so that a perfect white maps to $(1, 1, 1)$. The light still acts as a per-channel gain, so only the statistics of surface colours differ from M1.

The illuminant is drawn per scene from $U(0.3, 1)^3$. Each model is run on 100 scenes under five conditions:

• base;
• dominant colour: the left half of the image becomes one green surface $(0.10, 0.60, 0.15)$;
• clipped highlights: four discs of radius 4 pixels with radiance $4L$, which saturates every channel, scored with the saturation mask and, as a separate row, without it;
• two illuminants: the left half lit by $L$, the right half by a second light $L_2$, scored as the worse of the two errors against $L$ and $L_2$.

Measure. The angular error in degrees between the estimated and true illuminant, as Barron (2015, Section 6) and FC4 (Hu et al., 2017, Eq. 1) define it, computed in the numerically stable form $\operatorname{atan2}(\lVert a \times b \rVert, a \cdot b)$. Each cell prints the median and the [minimum, maximum] over the 100 scenes.

Checks. The code first prints exact checks: gray world on a noise-free scene whose mean reflectance is grey (K1); white patch on a noise-free scene containing an unclipped white surface (K2); for every method, that the estimate of $R \times L$ divided by the estimate of $R$ is parallel to $L$, the per-channel equivariance described above (K3); and that shades of gray with $p = 1$ equals gray world (K3b). For two illuminants there is an exact bound. The angular error obeys the triangle inequality on the sphere, so for any estimate $e$, $\angle(e, L) + \angle(e, L_2) \ge \angle(L, L_2)$, and the worse of the two errors is at least half the angle between the lights. K4 prints the smallest margin over that bound across all scenes and methods. Data checks print the angle of each model’s mean reflectance and per-channel maximum reflectance from grey (the gray-world and white-patch assumptions), that the base scenes contain no saturated pixels, and that every highlight pixel is saturated.

The plan fixed these hypotheses before the run: (H1) in the base condition white patch has the lowest median error and gray world the highest, under both models; (H2) from base to dominant colour, gray world’s median rises more than any other method’s; (H3) with clipped highlights, white patch and shades of gray have higher medians without the mask than with it; (H4) every method’s median worst-half error with two illuminants exceeds its own base median. A disclosure: the data-property lines (the two reflectance gaps, the saturated-pixel counts and the angle between the two lights) and the K4 margin were printed by a planning script on one to three scenes before these hypotheses were written, and the hypotheses took the two reflectance gaps into account; the models were not changed after that. A 20-scene pilot of the same script preceded the run; the editor approved it unchanged and declined a third reflectance model that the pilot report offered as an option, and nothing was changed after the pilot or after the run. The script runs in about 3 seconds (NumPy 2.3.5, SciPy 1.16.3), and two runs gave identical output.

import os, sys
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage

H = W = 128
SIGMA_NOISE = 0.005          # additive Gaussian sensor noise, linear units
DOMINANT = np.array([0.10, 0.60, 0.15])
HIGHLIGHT_GAIN = 4.0         # specular highlight = 4 x illuminant, clips every channel

# ---- reflectance model M2: smooth random spectra seen by Gaussian RGB sensitivities ----
lam = np.arange(400, 701, 10.0)
t = (lam - 400) / 300
sens = np.stack([np.exp(-0.5 * ((lam - c) / 30.0) ** 2) for c in (610, 540, 460)])
sens /= sens.sum(1, keepdims=True)          # a perfect white (s = 1) maps to (1, 1, 1)

def refl_m1(rng):                           # M1: independent uniform channels
    return rng.uniform(0.05, 0.9, 3)

def refl_m2(rng):                           # M2: sigmoid of a random smooth function
    g = rng.normal(-1.0, 1.0) + sum(rng.normal(0, 1.5) * np.cos(np.pi * m * t)
                                    + rng.normal(0, 1.5) * np.sin(np.pi * m * t) for m in (1, 2, 3))
    s = 1 / (1 + np.exp(-g))
    return sens @ s

MODELS = {"M1 independent uniform": refl_m1, "M2 smooth spectra": refl_m2}

def mondrian(rng, refl, n=60):
    R = np.empty((H, W, 3)); R[:] = refl(rng)
    for _ in range(n):
        y, x = rng.integers(0, H, 2); h, w = rng.integers(8, 40, 2)
        R[y:y + h, x:x + w] = refl(rng)
    return R

def illuminant(rng):
    return rng.uniform(0.3, 1.0, 3)

def capture(E, rng):
    """Linear sensor: radiance E plus Gaussian noise, clipped to [0, 1]; 1 = saturation."""
    return np.clip(E + rng.normal(0, SIGMA_NOISE, E.shape), 0, 1)

# ---- estimators (Finlayson 2013, Eq. 6 family) ----
def minkowski(V, p):
    return np.mean(V ** p, 0) ** (1 / p)

def gray_world(I, m):  return I[m].mean(0)
def white_patch(I, m): return I[m].max(0)
def shades_of_gray(I, m, p=6): return minkowski(np.clip(I[m], 0, None), p)
def gray_edge(I, m, p=6, s=1.0):
    G = np.stack([np.hypot(ndimage.gaussian_filter(I[..., c], s, order=(0, 1)),
                           ndimage.gaussian_filter(I[..., c], s, order=(1, 0))) for c in range(3)], -1)
    me = ndimage.binary_erosion(m, iterations=3, border_value=1)   # drop edges within 3 px of a masked pixel
    return minkowski(G[me], p)

METHODS = {"gray world": gray_world, "white patch": white_patch,
           "shades of gray p=6": shades_of_gray, "gray-edge n=1 p=6 s=1": gray_edge}

def ang(a, b):
    """Angular error in degrees, numerically stable form."""
    return np.degrees(np.arctan2(np.linalg.norm(np.cross(a, b)), np.dot(a, b)))

def valid_mask(I):
    return ~(I >= 1.0).any(-1)               # mask pixels with any saturated channel

def scene(model_i, i, refl):
    rng = np.random.default_rng([model_i, i])
    R = mondrian(rng, refl)
    L, L2 = illuminant(rng), illuminant(rng)
    hy, hx = rng.integers(8, H - 8, (2, 4))  # four highlight centres
    return R, L, L2, hy, hx

def conditions(R, L, L2, hy, hx):
    yy, xx = np.mgrid[:H, :W]
    Rd = R.copy(); Rd[:, :W // 2] = DOMINANT
    E_hl = R * L
    hl = np.zeros((H, W), bool)
    for y, x in zip(hy, hx): hl |= (yy - y) ** 2 + (xx - x) ** 2 <= 16
    E_hl[hl] = HIGHLIGHT_GAIN * L
    E_two = R * L; E_two[:, W // 2:] = R[:, W // 2:] * L2
    return {"base": R * L, "dominant colour": Rd * L, "clipped highlights": E_hl,
            "two illuminants": E_two}, hl

def run(n_scenes):
    out = {}
    print("checks")
    # K1 gray world exact on a grey-mean scene (noise-free)
    Rg = np.full((32, 32, 3), 0.4); Rg[:16] = [0.2, 0.4, 0.6]; Rg[16:] = [0.6, 0.4, 0.2]
    Lc = np.array([0.9, 0.7, 0.4]); m32 = np.ones((32, 32), bool)
    print("K1 gray world, grey-mean scene, error (deg): %.1e" % ang(gray_world(Rg * Lc, m32), Lc))
    # K2 white patch exact with a white patch present (noise-free, unclipped)
    rng = np.random.default_rng(99)
    Rw = mondrian(rng, refl_m1); Rw[10:26, 10:26] = 0.95
    m = np.ones((H, W), bool)
    print("K2 white patch, white patch present, error (deg): %.1e" % ang(white_patch(Rw * Lc, m), Lc))
    # K3 equivariance under the diagonal model: est(R*L) / est(R) is parallel to L
    e3 = max(ang(f(Rw * Lc, m) / f(Rw, m), Lc) for f in METHODS.values())
    print("K3 max over methods of angle(est(R*L)/est(R), L) (deg): %.1e" % e3)
    print("K3b shades of gray p=1 vs gray world (deg): %.1e" % ang(shades_of_gray(Rw * Lc, m, p=1), gray_world(Rw * Lc, m)))
    # linearity of the synthetic data (noise-free, unclipped)
    print("linearity max|E(2L) - 2E(L)|: %.1e" % np.abs(Rw * (2 * Lc) - 2 * (Rw * Lc)).max())
    for mi, (mname, refl) in enumerate(MODELS.items()):
        errs = {}; props = {"gw gap": [], "wp gap": [], "sat base": [], "sat clip": [], "hl px": [], "sat minus hl": [],
                            "bound margin": [], "L1-L2 angle": []}
        for i in range(n_scenes):
            R, L, L2, hy, hx = scene(mi, i, refl)
            props["gw gap"].append(ang(R.reshape(-1, 3).mean(0), np.ones(3)))
            props["wp gap"].append(ang(R.reshape(-1, 3).max(0), np.ones(3)))
            Es, hl = conditions(R, L, L2, hy, hx)
            nrng = np.random.default_rng([mi, i, 7])
            for cname, E in Es.items():
                I = capture(E, nrng)
                mk = valid_mask(I)
                variants = [(cname, mk)]
                if cname == "clipped highlights":
                    variants.append(("clipped highlights, no mask", np.ones((H, W), bool)))
                    props["sat clip"].append((~mk).sum()); props["hl px"].append(hl.sum())
                    props["sat minus hl"].append((~mk[hl]).sum() - hl.sum())
                if cname == "base": props["sat base"].append((~mk).sum())
                for vname, mm in variants:
                    for meth, f in METHODS.items():
                        e = f(I, mm)
                        if cname == "two illuminants":
                            e1, e2 = ang(e, L), ang(e, L2)
                            val = max(e1, e2)
                            props["bound margin"].append(val - ang(L, L2) / 2)
                        else:
                            val = ang(e, L)
                        errs.setdefault(vname, {}).setdefault(meth, []).append(val)
            props["L1-L2 angle"].append(ang(L, L2))
        print()
        print("model", mname, "| scenes:", n_scenes)
        print("data: median angle(mean reflectance, grey) %.2f deg; median angle(channel max reflectance, grey) %.2f deg"
              % (np.median(props["gw gap"]), np.median(props["wp gap"])))
        print("data: saturated pixels, base: max %d; clipped: min %d max %d (highlight discs: min %d max %d px)"
              % (max(props["sat base"]), min(props["sat clip"]), max(props["sat clip"]),
                 min(props["hl px"]), max(props["hl px"])))
        print("data: min over scenes of (saturated highlight pixels - highlight pixels): %d" % min(props["sat minus hl"]))
        print("data: median angle between the two illuminants %.2f deg" % np.median(props["L1-L2 angle"]))
        print("K4 min over scenes and methods of worst-half error - angle(L1,L2)/2 (deg): %.2e"
              % min(props["bound margin"]))
        print("angular error (deg): median [min, max] over scenes; two illuminants: worse of the two halves")
        print("%-30s" % "condition" + "".join("%-24s" % k for k in METHODS))
        for vname, d in errs.items():
            print("%-30s" % vname + "".join("%-24s" % ("%.2f [%.2f, %.2f]" % (np.median(v), min(v), max(v)))
                                            for v in d.values()))
        out[mname] = errs
    return out

if __name__ == "__main__":
    run(int(sys.argv[1]) if len(sys.argv) > 1 else 100)

Output:

checks
K1 gray world, grey-mean scene, error (deg): 6.1e-13
K2 white patch, white patch present, error (deg): 0.0e+00
K3 max over methods of angle(est(R*L)/est(R), L) (deg): 3.2e-12
K3b shades of gray p=1 vs gray world (deg): 0.0e+00
linearity max|E(2L) - 2E(L)|: 0.0e+00

model M1 independent uniform | scenes: 100
data: median angle(mean reflectance, grey) 6.63 deg; median angle(channel max reflectance, grey) 0.65 deg
data: saturated pixels, base: max 0; clipped: min 178 max 196 (highlight discs: min 178 max 196 px)
data: min over scenes of (saturated highlight pixels - highlight pixels): 0
data: median angle between the two illuminants 17.50 deg
K4 min over scenes and methods of worst-half error - angle(L1,L2)/2 (deg): 1.60e-01
angular error (deg): median [min, max] over scenes; two illuminants: worse of the two halves
condition                     gray world              white patch             shades of gray p=6      gray-edge n=1 p=6 s=1   
base                          6.26 [1.41, 12.65]      0.68 [0.02, 2.23]       2.93 [0.23, 6.96]       3.35 [0.37, 10.90]      
dominant colour               15.91 [8.45, 24.00]     0.87 [0.12, 4.72]       3.38 [0.05, 8.77]       3.17 [0.42, 11.16]      
clipped highlights            6.16 [1.43, 12.71]      0.66 [0.04, 2.18]       2.92 [0.14, 6.99]       3.30 [0.37, 10.76]      
clipped highlights, no mask   5.96 [0.68, 12.20]      13.75 [1.13, 25.26]     10.27 [0.63, 20.56]     15.49 [1.33, 29.96]     
two illuminants               14.21 [2.09, 28.46]     13.84 [3.43, 26.74]     13.05 [3.75, 29.06]     13.95 [4.05, 30.27]     

model M2 smooth spectra | scenes: 100
data: median angle(mean reflectance, grey) 8.38 deg; median angle(channel max reflectance, grey) 1.32 deg
data: saturated pixels, base: max 0; clipped: min 163 max 196 (highlight discs: min 163 max 196 px)
data: min over scenes of (saturated highlight pixels - highlight pixels): 0
data: median angle between the two illuminants 17.19 deg
K4 min over scenes and methods of worst-half error - angle(L1,L2)/2 (deg): 6.61e-02
angular error (deg): median [min, max] over scenes; two illuminants: worse of the two halves
condition                     gray world              white patch             shades of gray p=6      gray-edge n=1 p=6 s=1   
base                          7.42 [0.38, 22.02]      1.14 [0.18, 6.07]       3.40 [0.50, 10.81]      3.13 [0.40, 10.90]      
dominant colour               17.81 [11.38, 25.99]    1.67 [0.09, 7.16]       3.52 [0.40, 10.28]      3.78 [0.24, 12.16]      
clipped highlights            7.33 [0.17, 21.92]      1.07 [0.16, 6.15]       3.35 [0.54, 10.75]      3.29 [0.35, 10.86]      
clipped highlights, no mask   7.23 [0.16, 20.96]      14.03 [1.13, 27.93]     8.87 [0.89, 22.95]      12.13 [0.68, 29.97]     
two illuminants               15.04 [3.23, 31.43]     13.20 [3.63, 27.81]     12.68 [3.50, 28.80]     13.49 [3.30, 29.60]     

Checks. K1, K2, K3, K3b and the linearity check held, at $3.2 \times 10^{-12}$ degrees or below against a tolerance of $10^{-9}$. The two-illuminant bound K4 held in every scene and method, with smallest margins of 0.16° (M1) and 0.066° (M2). No base scene had a saturated pixel and every highlight pixel was saturated, in both models.

H1 held under both models, and it answers the question the two models were set up to ask: changing the reflectance model did not change the winner. White patch had the lowest base median (0.68° in M1, 1.14° in M2) and gray world the highest (6.26°, 7.42°). The data lines show why. In both models the per-channel maximum reflectance lies within about a degree of grey (0.65°, 1.32°), while the mean reflectance lies several degrees from it (6.63°, 8.38°). Since the estimators are exact when their statistic of the reflectance is grey (K1–K3), gray world’s and white patch’s errors follow the gaps of their own statistics, the two gaps the code prints. In our reading, the small channel-maximum gap comes from drawing 61 surfaces independently per scene, which leaves some surface near the top of each channel’s range in both models; it is a property of many-surface synthetic scenes rather than evidence about real photographs. The models do change the order of the two middle methods, a result that was not hypothesized: shades of gray beat gray-edge under M1 (2.93° against 3.35°) and lost to it under M2 (3.40° against 3.13°).

H2 held under both models. Making half the image green raised gray world’s median by 9.65° (M1) and 10.39° (M2), while the other methods moved by less than 0.7° either way; gray-edge’s median even fell slightly in M1. A green half pulls the mean toward green by construction, and a single uniform region adds edges only along its border; in our reading, white patch barely moves because the green surface’s largest value (0.60) stays below the channel maxima that the other surfaces reach.

H3 held under both models. Without the mask, white patch’s median rose from 0.66° to 13.75° (M1) and from 1.07° to 14.03° (M2); shades of gray’s rose from 2.92° to 10.27° and from 3.35° to 8.87°. A fully saturated pixel reads $(1, 1, 1)$ whatever the light, so it reports a grey light to any statistic that weights bright pixels heavily, and white patch’s estimate becomes exactly grey. Two further results were not hypothesized. Gray-edge was also hurt strongly without the mask (15.49° in M1, the largest of all, and 12.13° in M2): in our reading, the step from a disc at $(1, 1, 1)$ to its surround has size $1 - R_c L_c$ in channel $c$, which does not scale with the light as the diagonal model requires. Gray world’s median was slightly lower without the mask (5.96° against 6.16°, 7.23° against 7.33°): with the mask the clipped discs are removed, and without it they add grey pixels to the mean; in our reading, adding grey pixels can pull a mean that was off-grey slightly toward the truth or away from it depending on the scene, and the medians show only a small net effect. With the mask, every method’s median stayed within 0.2° of its base value.

H4 held for all four methods under both models: every worst-half median with two illuminants (12.68° to 15.04°) exceeded the method’s base median. Here the methods are nearly level; in our reading this is because the bound K4 forces each of them to be wrong by at least half the angle between the lights (median angles 17.50° and 17.19°) on one side or the other, whatever statistic it computes. A single global estimate cannot fix a scene lit by two lights.

The comparison is between known assumptions on synthetic data, not between methods on photographs. Real benchmark results differ in kind: on the Gehler–Shi set (below), white patch has the highest mean error of all the methods in Barron’s Table 1.

Datasets and Benchmarks

Benchmark scenes contain a Macbeth Color Checker from which the ground-truth illuminant is obtained; FC4 masks the chart by setting its region to zero for both training and testing (Hu et al., 2017, Section 4.1).

• Color Checker Dataset, reprocessed by Shi and Funt (often called Gehler–Shi): 568 raw images from Gehler et al. (2008), in the version reprocessed by Shi and Funt (Barron, 2015, Section 6; Hu et al., 2017, Section 4.1). Barron (2015, Section 6) describes it as from a single camera, while FC4 (Hu et al., 2017, Section 4.1) and FFCC (Barron and Tsai, 2017, Fig. 4) state that it contains images from two cameras.
• NUS 8-Camera Dataset of Cheng et al. (2014): 1,736 raw images from 8 cameras, about 210 per camera, with the same scenes taken by each camera (Barron, 2015, Section 6; Hu et al., 2017, Section 4.1). Results are computed per camera and summarized by the geometric mean across cameras.

Barron (2015, Section 6) and FC4 (Hu et al., 2017, Section 4.1) evaluate by three-fold cross-validation, and Barron and Tsai (2017, Fig. 5) show it for their training on the Gehler–Shi set. Results reported on the Gehler–Shi set in Barron’s (2015) Table 1, read from the page image (angular error in degrees, mean / median):

Method	Mean	Median
White-Patch	7.55	5.68
Gray-World	6.36	6.28
1st-order Gray-Edge	5.33	4.52
Shades-of-Gray	4.93	4.01
Bayesian	4.82	3.46
Pixel-based Gamut	4.20	2.33
Corrected-Moment (19 Edge), as reported by Finlayson	2.8	2.0
CCC (disc+ext)	1.95	1.22

Barron notes that several baseline numbers were taken from earlier papers that used different error statistics, which leaves some entries in the table missing. Within this table, white patch has the highest mean but a lower median than gray world, a sign of a heavy tail of large failures; the gap between the classical statistics and the learned and corrected methods is larger than the gaps among the statistics. This is a different data set and protocol from the synthetic example above, so the two are not compared number for number.

Evaluation Metrics

The standard measure is the angular error: the angle in degrees between the estimated and true illuminant RGB vectors,

$$\varepsilon = \frac{180}{\pi} \arccos \frac{\hat{\mathbf{L}} \cdot \mathbf{L}}{\lVert \hat{\mathbf{L}} \rVert \, \lVert \mathbf{L} \rVert},$$

which ignores the unrecoverable scale of the light (Hu et al., 2017, Eq. 1; Barron, 2015, Eq. 11 and Section 6). Because a few scenes can produce very large errors, papers report a set of statistics over the dataset: mean, median, tri-mean, the means of the best and worst 25%, and on the Gehler–Shi set the 95th percentile; on the NUS set FC4 also reports the geometric mean of the others (Hu et al., 2017, Section 4.1; Barron, 2015, Section 6). The arccos form loses precision for very small angles; the example uses an equivalent arctangent form. Other error measures exist; Barron (2015, Table 1) notes that missing entries in his table come from earlier work that used different error measures.

Applications

• Cameras. Automatic white balance estimates the illuminant and applies the gains to the raw data before the colour correction matrix and tone curve (Digital Camera, Image Formation). FFCC was designed for this setting, running in real time on a phone’s viewfinder thumbnails (Barron and Tsai, 2017).
• Recognition with colour. FC4’s authors motivate colour correction by its benefit to visual recognition, where colour distinguishes objects (Hu et al., 2017, Section 1). The alternative is to make features invariant to the illuminant, as Global Image Descriptors and the chromaticity spaces in Color Spaces describe.
• Consistency across frames and cameras. Because AWB is estimated per image, the gains applied to the same surface can change from frame to frame, which breaks the brightness constancy assumption that optical flow and direct methods rely on; on inspection lines white balance is fixed rather than estimated, as Industrial Visual Inspection describes.

Open Problems

• Multiple and spatially varying illuminants. The benchmarks and most methods assume one global light; the example shows that any global estimate then has an unavoidable error of at least half the angle between the lights. Estimating a per-region or per-pixel illuminant brings back the full intrinsic-image problem that the diagonal model sets aside (Barron, 2015, Section 2).
• Data. Labelled raw datasets are small and camera-specific. Barron (2015, Section 7) suggests that colour constancy may benefit from much larger datasets, as object detection and recognition did.
• Truth or preference. Barron and Tsai (2017, Section 1) note that colour constancy and white balance have come to mean different things, recovering the true scene against producing a pleasing image, and that for learned methods the choice is made by the training data. What the “right” illuminant is for a photograph, as opposed to a measurement, is not settled by the angular error against a chart.
• Processed images. Many images available to vision systems are already white-balanced and tone-mapped JPEGs, for which the linear model above no longer holds; in our reading, correcting them after the fact is a different problem from estimating the light on raw data.