Image Signal Processing Pipeline

An image signal processing pipeline is the sequence of operations that turns the raw values of a camera’s image sensor into a finished picture: black-level subtraction, lens shading correction, white balance, demosaicing, denoising, colour correction, tone mapping and gamma, sharpening and compression. The hardware that runs it is the image signal processor (ISP), and the digital camera article describes each of its stages in turn. This article treats the pipeline as a designed system. Its stages do not commute, so the order in which they run, and how clipping, noise and colour pass from one stage to the next, decide what the output looks like. The same chain can run in fixed hardware or in software, and it can be inverted approximately to turn processed photographs back into raw-like training data.

The pipeline is also called the camera pipeline, the imaging pipeline or the raw processing pipeline; “ISP” names both the chip and, loosely, the chain it runs.

Definition

An image signal processing pipeline is an ordered composition of operations, each of which maps an image to an image, that converts a mosaic of linear sensor values into a displayable, usually compressed, colour image in a standard colour space. A stage acts on one pixel at a time (offsets, gains, colour matrices, tone curves) or on a neighbourhood (demosaicing, denoising, sharpening). It may use parameters estimated from the image or from the camera’s state, such as the white-balance gains or the exposure.

Dewil (2022, Section 1) summarizes what pipelines share and where they differ: stages and their order differ between camera makers and software, but all share essential elements such as demosaicking, denoising, colour correction and tone mapping.

Intuition

Each stage is designed on an assumption about its input. A demosaicer that borrows detail across channels assumes the channels look alike. A denoiser tuned for one noise level assumes the noise has that level. A white-balance gain assumes that a pixel’s value is proportional to the light that reached it. An earlier stage can break the assumption of a later one. The red and blue white-balance gains in the Darmstadt Noise Dataset’s metadata all exceed 1 (Brooks et al., 2019, Section 3.4), so on a neutral surface the raw red and blue channels are dimmer copies of green, and a cross-channel demosaicer run before white balance borrows mismatched detail. After the tone curve, noise is no longer of one level. A clipped pixel is no longer proportional to its light, so a gain applied to it invents a colour. In our reading, designing a pipeline is largely a matter of putting each stage where its assumption still holds, at a cost the hardware can afford.

Formal Definition

Write the pipeline as a composition of stage maps applied to the raw mosaic $\mathbf{r}$:

$$\mathbf{y} = (S_n \circ \cdots \circ S_2 \circ S_1)(\mathbf{r}).$$

Swapping two adjacent stages leaves the output unchanged only if they commute, $S_i \circ S_j = S_j \circ S_i$. A few cases are decided by their form; we state them as reasoning, not as statements from a source:

• Per-channel gains and per-channel linear filters commute. White balance multiplies every sample of channel $c$ by a constant $g_c$, and a demosaicer that interpolates each channel only from its own samples is linear in that channel’s samples, so scaling the samples scales the result. Bilinear demosaicing is of this kind.
• Per-channel gains and cross-channel filters do not. A demosaicer whose estimate of one channel adds a weighted Laplacian of another channel’s samples, as Malvar–He–Cutler’s does (Getreuer, 2011, Section 2), mixes samples that white balance would have scaled by different constants.
• Pointwise non-linear curves and spatial filters do not. A tone curve $T$ and a filter $h$ satisfy $T(h * x) \ne h * T(x)$ unless $T$ is affine.
• Clipping and gains do not. Clipping at the white level and then multiplying by $g > 1$ gives $g \min(x, 1)$, which exceeds 1 for any clipped pixel, while multiplying first and clipping gives at most 1.

The example below measures the first, second and fourth cases and a consequence of the third.

Stage Orders in the Sources

Sources agree on the main stages but not on their order. The table lists each source’s order as the source gives it. None of them is the order of every camera, and this article picks none as correct.

Source	Order of the main stages
Karaimer and Brown (2016, Fig. 1), as Digital Camera reports it	pre-processing (black-level offset, normalization, bad-pixel mask); white balance; demosaicing; post-processing (noise reduction, sharpening); colour transform; colour rendering (tone mapping, sRGB gamma, colour manipulation); display or compression
Delbracio et al. (2021, Fig. 6), as Digital Camera reports it	demosaicing before white balance; noise reduction after tone mapping
HDR+ finishing pipeline (Hasinoff et al., 2016, Section 6)	black level; lens shading; white balance; demosaicking; chroma denoising; colour correction; dynamic range compression; dehazing; global tone curve and sRGB gamma; chromatic aberration correction; sharpening; hue-specific adjustments; dithering
Brooks et al. (2019, Section 4.2), processing applied before their loss	white balance; bilinear demosaicing; colour correction to sRGB; gamma compression

Two disagreements matter most. The first is white balance against demosaicing. Karaimer and Brown, HDR+ and Brooks et al.’s processing put white balance first, and Delbracio et al. put demosaicing first. Brooks et al.’s own sections differ on this: their Section 3, which describes a conventional pipeline from sensor measurement to final image, takes demosaicing (Section 3.2) before white balance (Section 3.4), and their Section 4.2 applies white balance first. For their bilinear demosaicer the difference vanishes, as the first case above shows.

The second is denoising, before or after the non-linear photo-finishing steps; Digital Camera reports where Karaimer and Brown and Delbracio et al. place it. HDR+ does both in its own way. Its merge is a temporal denoiser on raw data (burst photography explains it), and its finishing pipeline adds chroma denoising after demosaicking.

A third difference is completeness. Brooks et al.’s processing has no tone curve, although their unprocessing inverts one (Section 4.1). They describe their simplified processing as the one the Darmstadt Noise Dataset benchmark uses, and as a good approximation of general pipelines (Section 4.2).

Properties

White balance and demosaicing

Color constancy estimates the illuminant, and white balance divides it out by per-channel gains, which are large (the Darmstadt ranges are given under inverse pipelines below). Whether the demosaicer sees the channels before or after these gains makes no difference for per-channel methods. It does make a difference for methods that borrow detail from green, as Malvar–He–Cutler and the other cross-channel methods that the demosaicing article compares do. In our reading, this is one argument for the white-balance-first order: the methods’ inter-channel assumption refers to the scene’s colours, which white balance restores, not to the illuminant’s tint. HDR+ scales the four RGGB channels before demosaicking, with gains supplied by the hardware ISP (Section 6, steps 3–4). The example below measures the effect for Malvar–He–Cutler.

Clipping, highlights and white balance

The sensor clips at its white level separately in each channel. A neutral surface brighter than the white level is recorded as the same clipped value in every channel that saturated. White balance then multiplies those values by gains that differ between channels. Unless the pipeline clips again after white balance, or handles highlights in some other way, a clipped neutral highlight takes the colour of the gains themselves. With green as the reference channel (gain 1) and red and blue gains above 1, green is the lowest channel, and the highlight renders magenta. Brooks et al. (2019, Section 3.4) make the reverse point for their inverse pipeline. Their inverse gains are almost always below one, so plain division leaves a synthetic raw image with almost no clipped pixels, which they consider a problem because handling saturated intensities correctly is critical when denoising. Their highlight-preserving inverse gain is described under inverse pipelines below.

Noise through the pipeline

Raw noise is close to independent between pixels and has a variance affine in the signal (image noise). Brooks et al. (2019, Section 3.1) model it as a heteroscedastic Gaussian with variance $\lambda_{\text{read}} + \lambda_{\text{shot}}\,x$ and contrast it with the noise of processed images, which may have very complex characteristics because of non-linearities and correlation across pixel values. Digital camera describes how demosaicing, the colour matrix and the curves correlate, amplify and reshape it. A denoiser placed late in the chain faces noise that is spatially correlated and whose level depends on brightness and colour. A denoiser placed early faces the simple raw model, but on a mosaic, before the detail the demosaicer will reconstruct. Image denoising notes the consequence for methods: raw denoising has a simple noise model, and sRGB denoising has a complex one.

Bit depth

Precision changes along the chain. In HDR+, 10-bit raw frames are merged to 14 bits (burst photography gives the correction of July 2018), and the errata state that intermediate results in the finishing pipeline are 12 bits. The final step adds blue noise from a precomputed table before reducing from 12 to 8 bits, to avoid quantization artifacts (Section 6, step 13). Hasinoff et al. (Section 2) put the advantage of 10-bit raw over the 8-bit output of mobile ISPs at less than 2 bits, because raw is linear while the 8-bit output already carries a gamma curve.

Stages that depend on the scene

White balance, exposure and local tone mapping are chosen per image. Hasinoff et al. (2016, Section 3) describe a mobile ISP as continuously adjusting exposure time, gain, focus and white balance as the user aims the camera. A pipeline is therefore not one fixed function of its input; Digital Camera describes what this breaks for methods that compare intensities across images.

Hardware and Software Pipelines

A hardware ISP is dedicated hardware that runs the whole pipeline in real time (digital camera gives its speed and how manufacturers tune it); in the phones HDR+ ran on, it is part of the system on a chip (Hasinoff et al., 2016, Section 2). A software pipeline runs the same kind of stages on a general-purpose processor. It is slower, but it is free to change the stages and their order.

HDR+ uses both in one camera (Hasinoff et al., 2016, Section 2, Fig. 3):

• Viewfinder. A real-time pipeline produces a low-resolution stream on the hardware ISP of the phone’s system on a chip. In the authors’ implementation the stream is 1.6 megapixels at 15–30 frames per second. They use the ISP because it is power efficient.
• Photograph. After the shutter press, a burst of raw frames is captured at constant exposure. A non-real-time pipeline, written in the Halide language and running on the application processor, aligns and merges the frames and then finishes the merged image. The result is a single 12-megapixel 8-bit photograph, computed in about 4 seconds.

The two pipelines’ images look different; the authors write that their viewfinder is not WYSIWYG. The software pipeline depends on the camera interface exposing raw frames and per-frame exposure and gain (Android’s Camera2 API in their system). It still takes some parameters from the hardware ISP: the lens shading correction image, the white-balance gains and the colour correction matrix (Section 6, steps 2, 3 and 6).

Hasinoff et al. (2016, Section 2) give three reasons for starting from raw data rather than from the ISP’s output:

• More bits. Raw frames have more bits per pixel than the ISP’s output.
• Linearity. Raw values are linear, which lets the noise be modelled accurately, makes alignment and merging more reliable and makes auto-exposure easier.
• Portability. Merging the ISP’s output would mean modelling and reversing processing that is proprietary and scene-dependent.

The cost is that the whole photographic pipeline, including lens shading and chromatic aberration correction and demosaicking, must be implemented in software. Merging on raw data has one saving: demosaicking runs once, on the merged image, rather than on every frame.

The software pipeline is still constrained by the device. To make it fast enough, the authors avoided algorithms with non-local communication or data dependencies that prevent parallelization, kept memory use low and used fixed-point arithmetic where possible (Section 7). For the same reasons they ruled out three kinds of algorithm:

• globally iterative ones, such as FlexISP;
• ones with large or dynamic spatial support, such as BM3D;
• expensive tone mapping, such as local Laplacian filters.

Up to 104 megapixels of raw data can go into one photograph (Section 7). Hardware and software also remain coupled in failure: the authors’ most serious failure mode was that at very low light the ISP’s autofocus and white-balance estimation begin to fail, which can leave a photograph out of focus or with a colour cast even when merging works.

The HDR+ Finishing Pipeline

Hasinoff et al. (2016, Section 1) note that academic papers in this area describe only algorithms, not complete systems, and that the algorithms in commercial systems are proprietary and not easily reverse-engineered. Their paper publishes its own finishing pipeline in order. After the burst has been aligned and merged (burst photography covers both), the merged raw image passes through thirteen steps, in order of application (Hasinoff et al., 2016, Section 6, read from the rendered page):

1. Black-level subtraction. The offset is measured from optically shielded pixels on the sensor.
2. Lens shading correction. It brightens the corners to undo vignetting and corrects colour that varies with the angle at which light strikes the sensor. It uses a low-resolution RGGB image supplied by the ISP.
3. White balancing. The four RGGB channels are scaled linearly by factors supplied by the ISP.
4. Demosaicking. The output is linear RGB at 12 bits per pixel. It combines techniques from Gunturk et al. (2005): edge-directed interpolation with weighted averaging, constant-hue interpolation, and second-order gradients as correction terms.
5. Chroma denoising. It reduces red and green splotches in dark areas of low-light images, using an approximate bilateral filter (a sparse 3 × 3 non-linear kernel applied in two passes in YUV; compare the bilateral filter).
6. Colour correction. Sensor RGB is converted to linear sRGB by a 3 × 3 matrix supplied by the ISP.
7. Dynamic range compression. Local tone mapping, described below.
8. Dehazing. It reduces veiling glare with a global curve that pushes low values lower while preserving midtones and highlights. At most 0.1% of pixels may be clamped to zero, and only pixels below 7% of the white level are adjusted.
9. Global tone adjustment. An S-shaped contrast curve is concatenated with the standard sRGB transfer function.
10. Chromatic aberration correction. It assumes no lens model. It finds pixels along high-contrast edges and replaces their chroma with that of nearby pixels less likely to be affected.
11. Sharpening. It applies unsharp masking (image sharpening) with a sum of Gaussian kernels built from a three-level convolution pyramid.
12. Hue-specific adjustments. Bluish cyans and purples are shifted towards light blue, and blues and greens are generally saturated, to make skies and vegetation look more appealing.
13. Dithering. Blue noise is added before the reduction from 12 to 8 bits.

The dynamic range compression of step 7 adapts exposure fusion, which normally blends a bracket of captured exposures, to a single high-bit-depth image. The pipeline derives two synthetic exposures, one short and one long, from the merged image by gain and gamma, in greyscale. It fuses them with a fixed weighting function of luma that favours moderately bright pixels, undoes the gamma of the fused greyscale result, and restores colour by copying per-pixel chroma ratios from the linear RGB image (Section 6). The two exposures come from the auto-exposure algorithm, which burst photography describes; high dynamic range imaging covers exposure fusion in general.

The order shows the design choices discussed above:

• White balance before demosaicking, with a demosaicker that uses cross-channel second-order corrections.
• Two kinds of denoising. The merge denoises in raw, temporally, and chroma denoising runs after demosaicking, in YUV.
• Linear colour work, then curves. Colour correction and the dynamic range compression work on linear values, and the global tone curve and gamma come later.
• Correction after the tone curve. Chromatic aberration correction comes after it, and so does sharpening, which Hasinoff et al. (Section 5) note can amplify the ringing of their frequency-domain merge.
• Dithering last. It runs at the bit-depth reduction it protects.

Compression is not one of the thirteen steps. Figure 3’s caption describes the 8-bit output photograph as the input to compression and storage.

Inverse Pipelines

An inverse pipeline runs the stages backwards, from a processed image to an estimate of the raw data that produced it. Brooks et al. and Dewil use it to make training data. Raw images with clean references are scarce, while processed photographs are plentiful, and a network meant for raw input may generalize poorly when trained on data that does not resemble raw data (Brooks et al., 2019, Sections 1–2; Dewil, 2022, Section 1).

Brooks et al. (2019) invert each stage of a simplified pipeline (Sections 3 and 4.1, Fig. 2). To turn an sRGB image into synthetic raw data, they apply, in order:

1. Inverse tone curve. The tone curve is assumed to be the smoothstep $3x^2 - 2x^3$ (Eq. 9), inverted in closed form as $\tfrac12 - \sin\big(\sin^{-1}(1 - 2y)/3\big)$ (Eq. 10). They note that real cameras may use edge-aware local tone mapping, which is difficult to reverse-engineer (Section 3.7).
2. Inverse gamma. The gamma curve $x^{1/2.2}$, with its input clamped at $\epsilon = 10^{-8}$ (Eq. 7), is undone with the exponent 2.2 (Eq. 8).
3. Inverse colour correction. The inverse of a colour correction matrix is applied. The matrix is sampled as a random convex combination of the four fixed matrices of the Darmstadt cameras (Section 3.5).
4. Inverse gains. White balance and digital gain are inverted. The white-balance gains are sampled uniformly from the Darmstadt metadata’s ranges (red 1.9–2.4, blue 1.5–1.9; Section 3.4). The digital gain comes from matching the mean intensities of the synthetic and real data under an exponential model, which gives a ratio of 1.25; inverse gains are sampled from a normal distribution with mean 0.8 and standard deviation 0.1 (Section 3.3).
5. Remosaicking. Two of the three colour values at each pixel are dropped according to the Bayer pattern (Section 3.2).

Shot and read noise are then added with parameters sampled from a model fitted to the Darmstadt data (Section 3.1, Eq. 3).

The inverse gains need care, for the reason given under clipping above. Brooks et al. (2019, Eqs. 5–6) apply a highlight-preserving transformation with threshold $t = 0.9$:

$$\alpha(x) = \left(\frac{\max(x - t, 0)}{1 - t}\right)^2, \qquad f(x, g) = \max\!\left(\frac{x}{g},\; \big(1 - \alpha(x)\big)\frac{x}{g} + \alpha(x)\, x\right).$$

Below $t$ it is plain division. Above $t$, for $g > 1$, it bends up to reach $f(1, g) = 1$, so values that were white stay saturated. The text of Section 3.4 says that $f(1, g) = 1$ holds when $g \le 1$, but Eq. 6 gives $f(1, g) = 1/g > 1$ for $g < 1$, and the paper’s Fig. 4 draws the $g = 0.75$ curve as a straight line that reaches the top of the plot, 1, already at $x = 0.75$. In our reading of the rendered page, the condition should be $g > 1$; the example prints both cases. Dewil (2022, Section 2) writes the same function as $x/g + 100\,x(x - 0.9)^2(1 - 1/g)$ for $x \ge 0.9$ and inverse gain below one, which is Eq. 6 expanded.

Brooks et al. also run the forward pipeline inside training. The network’s output and the clean synthetic raw image both pass through white balance, bilinear demosaicing, the colour matrix and gamma compression before the L1 loss is computed, so that training is aware of the processing that will follow denoising (Section 4.2). On the Darmstadt Noise Dataset, their Table 1 (read from the rendered page) gives 40.35 dB sRGB PSNR for the model trained with this loss in sRGB. The ablations of that model give:

Ablation of the sRGB-loss model	sRGB PSNR (dB)
without unprocessing	39.02
without unprocessing, with four times as many channels per layer	39.35
without the colour matrix, white balance and gain inversion	39.70
without inverting tone mapping and gamma	40.23
noise-blind, trained with additive white Gaussian noise	38.65

The authors conclude that performance is most sensitive to the noise model, that unprocessing contributes substantially, especially on the sRGB metrics, and that a larger network does not make up for omitting it (Section 5).

Dewil (2022) reimplements this unprocessing as an Image Processing On Line demonstration. The article describes each inverse step, gives the algorithm in pseudo-code, adds heteroscedastic Gaussian noise of variance $a u + b$, and exposes the parameters (gamma, gains, a chosen camera’s matrix or a convex combination, black and white levels). Its Section 4 shows how the choice of gains changes the number of saturated areas the highlight-preserving function creates. The article is an MLBriefs paper, whose source code it states has not been reviewed. Two details differ from Brooks et al.:

• Dewil dates Brooks et al.’s proposal to 2018, while the paper appeared at CVPR 2019.
• Dewil gives the global gain’s range as about $[0.5, 1.1]$ in Section 2, but uses $g_{\text{rgb}} = 0.85$ and $2.05$ as the minimum and maximum in Section 4.

An inverse pipeline of this kind is approximate in three ways:

• It inverts a model pipeline, not the camera’s own. The parameters it cannot recover from an image, such as the gains, the matrix and the curve, are sampled from plausible ranges.
• Local tone mapping, sharpening, denoising and compression are not inverted at all. Brooks et al. reduce the effect of noise, quantization, JPEG compression and demosaicing by downsampling their source images by 2 with a Gaussian kernel (Section 4.4).
• Information lost to clipping in the source image cannot be restored. It can only be made statistically plausible, which is the purpose of $f$.

Hasinoff et al. (2016, Section 2) give the converse argument: modelling and reversing an ISP’s proprietary, scene-dependent processing is what starting from raw data avoids.

Examples

Stage order on a simulated pipeline

The script builds a small software pipeline with exact ground truth and measures four interactions between stages. It is an illustration of the properties above under design choices we fixed in a plan before it ran, not a test of a particular camera:

• Gains. The white-balance gains are $(2.15, 1, 1.70)$, inside the ranges that Brooks et al. (2019, Section 3.4) sample.
• Colour matrix. The matrix is the one in Digital Camera’s example. Its rows sum to one, so greys stay grey.
• Curves. Gamma $x^{1/2.2}$ and the smoothstep tone curve follow Brooks et al.’s Eqs. 7 and 9, without the $\epsilon$.
• Noise. Raw noise follows the heteroscedastic Gaussian form of their Eq. 1, with variance $0.004\,x + 10^{-4}$, values we chose.
• Mosaic. The colour filter array is RGGB with red at row 0, column 0, as in Demosaicing, Digital Camera and Image Formation. Arrays are indexed [row, column].
• Border. All image measures drop an 8-pixel border.

There are two scenes, each given as white-balanced linear camera RGB in $[0, 1]$; raw values are the scene divided by the gains:

• Neutral chart. A grey ramp, three grey discs and two sinusoidal gratings, with R = G = B.
• Coffee. scikit-image’s bundled coffee image, unprocessed into linear camera RGB with the inverse of this same pipeline: inverse smoothstep, exponent 2.2, inverse colour matrix, clipping to $[0, 1]$. It is a processed photograph, so it is exact ground truth for the simulation but not a true scene.

The four parts are:

• Part A applies white balance before and after demosaicing, for bilinear interpolation and for Malvar–He–Cutler with Getreuer’s (2011, Fig. 3) filters, taken from the demosaicing article’s code (with blue at green sites assigned as that article does, not as the figure’s labels read). It compares the two orders of each demosaicer, not the demosaicers, which that article does. It prints the largest difference between the orders and each order’s RMSE to the scene, after checking both demosaicers for exactness on a constant and on per-channel linear ramps.
• Part B follows a neutral surface of $k$ times the green saturation level through a sensor that clips each raw channel at 1, then white balance, the colour matrix, a clip to $[0, 1]$, gamma and the tone curve. It does this without and with a second clip to 1 right after white balance, and prints the rendered colour, its spread (largest minus smallest channel) and whether green is the lowest channel.
• Part C adds raw noise for seeds 0, 1 and 2 and compares two places for the same denoiser, scikit-image’s denoise_wavelet (Haar wavelet, soft BayesShrink thresholds, noise level estimated from the finest diagonal wavelet coefficients). In raw, it is applied to each of the four Bayer planes before the pipeline. After the tone curve, it is applied to each rendered channel. The measure is the RMSE to the same pipeline run on the noise-free mosaic, over the image and in dark, mid and bright bands of the noise-free rendering’s mean channel value. The script also prints the noise level the denoiser estimates in each place and, for the rendered channels, the actual standard deviation of the noise.
• Part D unprocesses coffee. It checks the round trip, checks Brooks et al.’s $f$, counts saturated raw pixels after plain inverse gains and after $f$ (inverse gains $1/(1.25\,g_c)$, with the digital gain of their Section 3.3), and processes the $f$-unprocessed image forward again, without and with the clip after white balance.

import os
for v in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread, set before NumPy
import numpy as np
import pywt
from scipy import ndimage, stats
from skimage import data
from skimage.restoration import denoise_wavelet

# Fixed design choices (see the text): white-balance gains inside the ranges that
# Brooks et al. (2019, Sec. 3.4) sample, Digital Camera's example colour matrix
# (rows sum to 1, so greys stay grey), and heteroscedastic Gaussian raw noise
# (Brooks et al., Eq. 1) with our own lambda_shot and lambda_read.
GAINS = np.array([2.15, 1.0, 1.70])           # R, G, B white-balance gains
CCM = np.array([[1.60, -0.45, -0.15], [-0.25, 1.45, -0.20], [0.00, -0.50, 1.50]])
LAM_SHOT, LAM_READ = 0.004, 1e-4              # raw noise variance = 0.004 x + 1e-4
S = (slice(8, -8), slice(8, -8))              # interior: 8-pixel border dropped

# --- Bayer RGGB mosaic and two demosaicers (as in the Demosaicing article) ---
def masks(H, W):
    m = {k: np.zeros((H, W), bool) for k in ("R", "Gr", "Gb", "B")}
    m["R"][0::2, 0::2] = 1; m["Gr"][0::2, 1::2] = 1
    m["Gb"][1::2, 0::2] = 1; m["B"][1::2, 1::2] = 1
    m["G"] = m["Gr"] | m["Gb"]
    return m

def mosaic(rgb, m):
    return rgb[..., 0] * m["R"] + rgb[..., 1] * m["G"] + rgb[..., 2] * m["B"]

def corr(x, k):          # 'mirror' reflects about the edge pixel: keeps CFA parity
    return ndimage.correlate(x, np.asarray(k, float), mode="mirror")

K4 = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) / 4
K8 = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
def bilinear(cfa, m):    # each channel from its own samples
    return np.stack([corr(cfa * m["R"], K8), corr(cfa * m["G"], K4),
                     corr(cfa * m["B"], K8)], -1)

# Malvar-He-Cutler, Getreuer (2011), Fig. 3, coefficients scaled by 8
K_G = np.array([[0, 0, -1, 0, 0], [0, 0, 2, 0, 0], [-1, 2, 4, 2, -1],
                [0, 0, 2, 0, 0], [0, 0, -1, 0, 0]]) / 8
K_ROW = np.array([[0, 0, .5, 0, 0], [0, -1, 0, -1, 0], [-1, 4, 5, 4, -1],
                  [0, -1, 0, -1, 0], [0, 0, .5, 0, 0]]) / 8
K_COL = K_ROW.T
K_DIA = np.array([[0, 0, -1.5, 0, 0], [0, 2, 0, 2, 0], [-1.5, 0, 6, 0, -1.5],
                  [0, 2, 0, 2, 0], [0, 0, -1.5, 0, 0]]) / 8
def mhc(cfa, m):
    g, row, col, dia = (corr(cfa, k) for k in (K_G, K_ROW, K_COL, K_DIA))
    G = np.where(m["G"], cfa, g)
    R = np.where(m["R"], cfa, np.where(m["Gr"], row, np.where(m["Gb"], col, dia)))
    B = np.where(m["B"], cfa, np.where(m["Gb"], row, np.where(m["Gr"], col, dia)))
    return np.stack([R, G, B], -1)

# --- Pipeline stages ---
def wb_mosaic(cfa, m, g):                     # white balance on the mosaic
    return cfa * (g[0] * m["R"] + g[1] * m["G"] + g[2] * m["B"])
gamma = lambda x: x ** (1 / 2.2)              # Brooks et al., Eq. 7 (without eps)
inv_gamma = lambda y: y ** 2.2
smooth = lambda x: 3 * x ** 2 - 2 * x ** 3    # Brooks et al., Eq. 9
inv_smooth = lambda y: 0.5 - np.sin(np.arcsin(1 - 2 * y) / 3)   # Eq. 10
def finish(rgb_wb):                           # colour matrix, clip, gamma, tone curve
    return smooth(gamma(np.clip(rgb_wb @ CCM.T, 0, 1)))
def render(cfa, m):                           # WB, bilinear demosaic, finish
    return finish(bilinear(wb_mosaic(cfa, m, GAINS), m))

def brooks_f(x, g, t=0.9):                    # highlight-preserving inverse gain, Eqs. 5-6
    a = (np.maximum(x - t, 0) / (1 - t)) ** 2
    return np.maximum(x / g, (1 - a) * (x / g) + a * x)

def rmse(a, b):
    return np.sqrt(np.mean((a - b) ** 2))

# --- Scenes: white-balanced linear camera RGB in [0, 1] ---
def neutral_chart(H=256, W=256):
    y, x = np.mgrid[0:H, 0:W].astype(float)
    v = 0.05 + 0.55 * x / (W - 1)                                  # ramp
    for cy, cx, lev in [(64, 48, 0.02), (64, 128, 0.30), (64, 208, 0.85)]:
        v[(y - cy) ** 2 + (x - cx) ** 2 <= 28 ** 2] = lev           # discs
    b1 = (y >= 120) & (y < 170); v[b1] = 0.4 + 0.15 * np.sin(2 * np.pi * x[b1] / 8)
    b2 = (y >= 190) & (y < 240); v[b2] = 0.4 + 0.15 * np.sin(2 * np.pi * y[b2] / 5)
    return np.repeat(v[..., None], 3, -1)

x8 = data.coffee().astype(float) / 255        # 400 x 600 sRGB, bundled
Wraw = inv_gamma(inv_smooth(x8)) @ np.linalg.inv(CCM).T          # unprocess
clipped = ((Wraw < 0) | (Wraw > 1)).any(-1)
Wc = np.clip(Wraw, 0, 1)
SCENES = [("neutral chart", neutral_chart()), ("coffee", Wc)]

# ===== A. White balance and demosaicing: does the order matter? =====
print("A. Checks: max |error| on a constant and on per-channel ramps (64 x 64)")
H = W = 64; m = masks(H, W); yy, xx = np.mgrid[0:H, 0:W].astype(float)
const = np.ones((H, W, 3)) * [0.3, 0.5, 0.7]
ramp = np.stack([0.2 + 0.003 * xx + 0.001 * yy, 0.4 - 0.002 * xx + 0.004 * yy,
                 0.6 + 0.001 * xx - 0.003 * yy], -1)
for nm, f in (("bilinear", bilinear), ("MHC", mhc)):
    e = [np.abs(f(mosaic(t, m), m) - t)[6:-6, 6:-6].max() for t in (const, ramp)]
    print(f"   {nm:9s} constant {e[0]:.1e}  ramps {e[1]:.1e}")
print("A. Order of white balance (WB) and demosaicing, interior")
print(f"   {'scene':14s}{'demosaicer':11s}{'max|diff|':>10s}{'RMSE WB first':>15s}{'RMSE dem. first':>17s}")
for sn, Wsc in SCENES:
    m = masks(*Wsc.shape[:2]); cfa = mosaic(Wsc / GAINS, m)
    for nm, f in (("bilinear", bilinear), ("MHC", mhc)):
        wb_first = f(wb_mosaic(cfa, m, GAINS), m)
        dm_first = f(cfa, m) * GAINS
        d = np.abs(wb_first - dm_first)[S].max()
        print(f"   {sn:14s}{nm:11s}{d:10.1e}{rmse(wb_first[S], Wsc[S]):15.5f}"
              f"{rmse(dm_first[S], Wsc[S]):17.5f}")

# ===== B. A neutral highlight clipped at the sensor =====
print("\nB. Neutral surface of k times the green saturation level")
print(f"   {'k':>5s}  {'after WB':>20s}  {'rendered, no clip after WB':>27s}"
      f"{'spread':>9s}{'G lowest':>9s}  {'spread, clip after WB':>21s}")
for k in [0.5, 0.9, 1.0, 1.2, 1.5, 1.7, 2.0, 2.15, 3.0]:
    wb = np.minimum(k / GAINS, 1.0) * GAINS    # sensor clips raw at 1, then WB
    o_no, o_cl = finish(wb), finish(np.minimum(wb, 1.0))
    low = o_no[1] < min(o_no[0], o_no[2])
    print(f"   {k:5.2f}  {np.array2string(wb, precision=3):>20s}  "
          f"{np.array2string(o_no, precision=3):>27s}{np.ptp(o_no):9.1e}{str(low):>9s}"
          f"  {np.ptp(o_cl):21.1e}")
print(f"   max |after WB - gains| for k >= 2.15: "
      f"{max(np.abs(np.minimum(k / GAINS, 1) * GAINS - GAINS).max() for k in (2.15, 3.0)):.1e}")

# ===== C. Denoising in raw or after the tone curve =====
def sigma_est(img2d):                         # the estimate denoise_wavelet makes
    d = pywt.dwtn(img2d, "db1")["dd"]; d = d[np.nonzero(d)]
    return np.median(np.abs(d)) / stats.norm.ppf(0.75)
def wden(img2d):
    return denoise_wavelet(img2d, wavelet="db1", mode="soft", method="BayesShrink",
                           rescale_sigma=True)
PLANES = {"R": (slice(0, None, 2), slice(0, None, 2)), "Gr": (slice(0, None, 2), slice(1, None, 2)),
          "Gb": (slice(1, None, 2), slice(0, None, 2)), "B": (slice(1, None, 2), slice(1, None, 2))}
def denoise_raw(cfa):                         # each Bayer plane separately
    out = cfa.copy()
    for sl in PLANES.values():
        out[sl] = wden(cfa[sl])
    return out
def denoise_rendered(img):                    # each rendered channel separately
    return np.stack([wden(img[..., c]) for c in range(3)], -1)

print("\nC. Checks")
rng = np.random.default_rng(99)
flat = 0.3 + 0.02 * rng.standard_normal((256, 256))
print(f"   sigma estimate on white noise of 0.02: {sigma_est(flat):.5f}")
xr = 0.2; nz = np.sqrt(LAM_SHOT * xr + LAM_READ) * rng.standard_normal((256, 256))
print(f"   raw noise variance / model at x = 0.2: {nz.var() / (LAM_SHOT * xr + LAM_READ):.4f}")
BANDS = (("dark", 0.0, 0.25), ("mid", 0.25, 0.6), ("bright", 0.6, 1.01))
SEEDS = (0, 1, 2)
for sn, Wsc in SCENES:
    m = masks(*Wsc.shape[:2]); cfa = mosaic(Wsc / GAINS, m)
    ref = render(cfa, m)                      # exact noise-free rendering
    L = ref.mean(-1)[S]
    res = {k: [] for k in ("noisy", "denoised in raw", "denoised after tone curve")}
    sig = {"raw": [], "rendered est": [], "rendered true": []}
    for seed in SEEDS:
        rng = np.random.default_rng(seed)
        noisy_cfa = cfa + np.sqrt(LAM_SHOT * cfa + LAM_READ) * rng.standard_normal(cfa.shape)
        noisy = render(noisy_cfa, m)
        outs = {"noisy": noisy, "denoised in raw": render(denoise_raw(noisy_cfa), m),
                "denoised after tone curve": denoise_rendered(noisy)}
        for k, o in outs.items():
            e = (o - ref)[S]
            res[k].append([rmse(o[S], ref[S])] + [np.sqrt(np.mean(e[(L >= a) & (L < b)] ** 2))
                                                   for _, a, b in BANDS])
        sig["raw"].append([sigma_est(noisy_cfa[sl]) for sl in PLANES.values()])
        sig["rendered est"].append([sigma_est(noisy[..., c]) for c in range(3)])
        sig["rendered true"].append([np.std((noisy - ref)[S][..., c]) for c in range(3)])
    counts = [int(((L >= a) & (L < b)).sum()) for _, a, b in BANDS]
    print(f"C. {sn}: RMSE to the noise-free rendering, mean [min, max] over seeds {SEEDS}")
    print(f"   band pixel counts (dark, mid, bright): {counts}")
    print(f"   {'':26s}" + "".join(f"{h:>26s}" for h in ("all", "dark", "mid", "bright")))
    for k, v in res.items():
        v = np.array(v)
        print(f"   {k:26s}" + "".join(f"{v[:, j].mean():9.5f} [{v[:, j].min():.5f}, {v[:, j].max():.5f}]"
                                       for j in range(4)))
    s = {k: np.mean(v, 0) for k, v in sig.items()}
    print("   sigma estimated per raw plane (R, Gr, Gb, B): " + " ".join(f"{x:.4f}" for x in s["raw"]))
    print("   sigma per rendered channel, estimated / actual: "
          + " ".join(f"{a:.4f}/{b:.4f}" for a, b in zip(s["rendered est"], s["rendered true"])))

# ===== D. The inverse pipeline (unprocessing) =====
print("\nD. Unprocessing coffee with an inverse of this pipeline")
print(f"   round trip max |error| (pixels not clipped by the inverse matrix): "
      f"{np.abs(finish(Wc) - x8)[~clipped].max():.1e}; clipped pixels: {clipped.sum()}")
xs = np.linspace(0, 1, 100001)
GI = GAINS * 1.25                             # WB gains times a digital gain of 1.25
for g in (0.75, 1.25, GI[0], GI[2]):
    lo = xs <= 0.9
    print(f"   g = {g:6.4f}: max |f - x/g| for x <= 0.9: {np.abs(brooks_f(xs[lo], g) - xs[lo] / g).max():.1e};"
          f" f(1, g) = {brooks_f(1.0, g):.4f}")
plain, fb = Wc / GI, brooks_f(Wc, GI)
print(f"   max raw value, plain division: {plain.max():.4f}; with f: {fb.max():.4f}")
print(f"   pixels saturated in raw (any channel at 1): plain {(plain >= 1).any(-1).sum()},"
      f" f {(fb >= 1).any(-1).sum()}; pixels with a channel at 1 before unprocessing: {(Wc >= 1).any(-1).sum()}")
hl = (Wc > 0.9).any(-1)
for nm, v in (("no clip after WB", fb * GI), ("clip after WB", np.minimum(fb * GI, 1))):
    print(f"   forward again, {nm:16s}: RMSE to the original on {hl.sum()} highlight pixels "
          f"{rmse(finish(v)[hl], x8[hl]):.5f}; elsewhere {rmse(finish(v)[~hl], x8[~hl]):.1e}")

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2, PyWavelets with pywt.__version__ 1.8.0):

A. Checks: max |error| on a constant and on per-channel ramps (64 x 64)
   bilinear  constant 0.0e+00  ramps 1.1e-16
   MHC       constant 1.1e-16  ramps 2.2e-16
A. Order of white balance (WB) and demosaicing, interior
   scene         demosaicer  max|diff|  RMSE WB first  RMSE dem. first
   neutral chart bilinear      2.2e-16        0.02614          0.02614
   neutral chart MHC           2.0e-01        0.01006          0.02334
   coffee        bilinear      2.2e-16        0.02594          0.02594
   coffee        MHC           5.7e-01        0.01633          0.02040

B. Neutral surface of k times the green saturation level
       k              after WB   rendered, no clip after WB   spread G lowest  spread, clip after WB
    0.50         [0.5 0.5 0.5]             [0.82 0.82 0.82]  1.1e-16    False                1.1e-16
    0.90         [0.9 0.9 0.9]          [0.994 0.994 0.994]  0.0e+00    False                0.0e+00
    1.00            [1. 1. 1.]                   [1. 1. 1.]  0.0e+00    False                0.0e+00
    1.20         [1.2 1.  1.2]          [1.    0.995 1.   ]  5.1e-03     True                0.0e+00
    1.50         [1.5 1.  1.5]          [1.    0.967 1.   ]  3.3e-02     True                0.0e+00
    1.70         [1.7 1.  1.7]          [1.    0.933 1.   ]  6.7e-02     True                0.0e+00
    2.00         [2.  1.  1.7]          [1.    0.895 1.   ]  1.1e-01     True                0.0e+00
    2.15      [2.15 1.   1.7 ]          [1.    0.872 1.   ]  1.3e-01     True                0.0e+00
    3.00      [2.15 1.   1.7 ]          [1.    0.872 1.   ]  1.3e-01     True                0.0e+00
   max |after WB - gains| for k >= 2.15: 0.0e+00

C. Checks
   sigma estimate on white noise of 0.02: 0.01993
   raw noise variance / model at x = 0.2: 0.9983
C. neutral chart: RMSE to the noise-free rendering, mean [min, max] over seeds (0, 1, 2)
   band pixel counts (dark, mid, bright): [3703, 11901, 41996]
                                                    all                      dark                       mid                    bright
   noisy                       0.06767 [0.06737, 0.06811]  0.06813 [0.06770, 0.06839]  0.07948 [0.07891, 0.07994]  0.06389 [0.06357, 0.06447]
   denoised in raw             0.04248 [0.04209, 0.04311]  0.04434 [0.04385, 0.04521]  0.04686 [0.04648, 0.04748]  0.04098 [0.04043, 0.04170]
   denoised after tone curve   0.06287 [0.06255, 0.06333]  0.06270 [0.06214, 0.06305]  0.07371 [0.07311, 0.07420]  0.05945 [0.05913, 0.06005]
   sigma estimated per raw plane (R, Gr, Gb, B): 0.0269 0.0370 0.0367 0.0290
   sigma per rendered channel, estimated / actual: 0.0268/0.0816 0.0325/0.0506 0.0231/0.0666
C. coffee: RMSE to the noise-free rendering, mean [min, max] over seeds (0, 1, 2)
   band pixel counts (dark, mid, bright): [54803, 139368, 30085]
                                                    all                      dark                       mid                    bright
   noisy                       0.06916 [0.06909, 0.06922]  0.06310 [0.06281, 0.06329]  0.07328 [0.07319, 0.07338]  0.05936 [0.05910, 0.05982]
   denoised in raw             0.05003 [0.04981, 0.05022]  0.04035 [0.04002, 0.04079]  0.05354 [0.05329, 0.05368]  0.04917 [0.04914, 0.04920]
   denoised after tone curve   0.06328 [0.06311, 0.06346]  0.05727 [0.05690, 0.05750]  0.06711 [0.06691, 0.06733]  0.05483 [0.05458, 0.05520]
   sigma estimated per raw plane (R, Gr, Gb, B): 0.0263 0.0291 0.0292 0.0190
   sigma per rendered channel, estimated / actual: 0.0261/0.0717 0.0374/0.0586 0.0320/0.0757

D. Unprocessing coffee with an inverse of this pipeline
   round trip max |error| (pixels not clipped by the inverse matrix): 7.6e-16; clipped pixels: 4
   g = 0.7500: max |f - x/g| for x <= 0.9: 0.0e+00; f(1, g) = 1.3333
   g = 1.2500: max |f - x/g| for x <= 0.9: 0.0e+00; f(1, g) = 1.0000
   g = 2.6875: max |f - x/g| for x <= 0.9: 0.0e+00; f(1, g) = 1.0000
   g = 2.1250: max |f - x/g| for x <= 0.9: 0.0e+00; f(1, g) = 1.0000
   max raw value, plain division: 0.8000; with f: 1.0000
   pixels saturated in raw (any channel at 1): plain 0, f 4; pixels with a channel at 1 before unprocessing: 4
   forward again, no clip after WB: RMSE to the original on 1040 highlight pixels 0.03511; elsewhere 1.0e-16
   forward again, clip after WB   : RMSE to the original on 1040 highlight pixels 0.00183; elsewhere 1.0e-16

Part A. Both demosaicers pass their exactness checks. For bilinear interpolation the two orders agree to rounding error on both scenes, and their RMSEs are identical, as the commutation argument predicts. For Malvar–He–Cutler the orders differ by up to 0.20 on the neutral chart and 0.57 on coffee. White balance first gives the lower error to the scene in both cases. On the neutral chart its RMSE is less than half that of the other order. On coffee the gap is smaller. In our reading the neutral chart is the case the method’s correction is built for, since after white balance its three channels are identical, while before white balance they are scaled copies whose Laplacians the filter mixes with the wrong weights. Coffee’s channels differ after white balance too, so part of the mismatch remains in either order.

Part B. Up to $k = 1$ the surface stays neutral to rounding error, with or without the second clip. Above $k = 1$ the green channel clips first. After white balance the pixel is $(\min(k, 2.15), 1, \min(k, 1.70))$, and from $k = 2.15$ on it is exactly the gain vector. Without a clip after white balance, the colour matrix pulls rendered green below 1 while red and blue clip at 1, so the highlight is magenta: faintly at $k = 1.2$, more strongly as $k$ grows. The tint stops growing once all three raw channels are clipped. With the clip, every case renders neutral white.

Part C. Denoising the raw planes lowers the error to the noise-free rendering much more than denoising the same image after the tone curve, on both scenes, in every band and for every seed; the seed ranges do not overlap. The printed noise levels suggest why. In the rendered channels the denoiser’s estimate is between about a third and two-thirds of the actual noise standard deviation. In our reading, demosaicing and the colour matrix have correlated the noise, so the finest diagonal wavelet coefficients, from which the level is estimated, carry less of it than white noise of the same strength would; the threshold is then too low. In raw, the four planes are independent white noise. That is the denoiser’s model, except that the level varies with the signal while the estimate is one number per plane. The result concerns one denoiser that assumes white noise of one level. A denoiser built for the rendered domain’s noise, or a raw denoiser given the signal-dependent model, could change the margins.

Part D. The inverse pipeline is exact wherever the inverse colour matrix did not push a value outside $[0, 1]$; four pixels were clipped there. Brooks et al.’s $f$ equals plain division below 0.9 for every printed gain. It reaches 1 at $x = 1$ for the gains above 1, and returns $1/0.75$ for $g = 0.75$, which is the discrepancy with the paper’s text noted above. Plain inverse gains cap raw values at $1/1.25 = 0.8$, so no pixel is saturated. With $f$, exactly the pixels that had a channel at 1 before unprocessing are saturated. Processing the $f$-unprocessed image forward again reproduces every non-highlight pixel to rounding error, because $f$ is plain division there. On the 1,040 highlight pixels, the round trip is far closer with a clip after white balance than without it. In our reading, forward gains turn the raised values of $f$ into values above the original ones, above 1 near the top of the range, which the colour matrix then mixes into the other channels; the clip removes the excess above 1 before the mixing.

We wrote the hypotheses and checks in a plan before the script first ran, with each check’s tolerance written beside the value that a planning script reached with the same components on the same scenes; one check, the exact round trip outside the highlights, rested on a derivation instead. Every check held: exactness of both demosaicers, the bilinear commutation, neutral rendering for $k \le 1$ and with the clip after white balance, the clipped neutral equal to the gains, the noise-level estimate within 5% of a known white-noise level, the simulated raw noise variance within 5% of its model, the round trip, the properties of $f$, the saturation counts, and the exact round trip outside the highlights. The hypotheses’ outcomes were as follows:

• The two orders of bilinear demosaicing agree: held.
• The two orders of Malvar–He–Cutler differ by more than $10^{-3}$: held on both scenes.
• White balance first gives the lower Malvar–He–Cutler RMSE: held on the neutral chart and held on coffee.
• Without the clip after white balance, the highlight is non-neutral (spread above $10^{-3}$) with green lowest for every $k > 1$: held at every printed $k$. The spread part of this hypothesis is a consequence of the construction, not a test. The planning script had already printed these spreads before the bar of $10^{-3}$ was written. The green-lowest part had not been computed.
• With the clip, the highlight is neutral: held.
• Raw-domain denoising beats denoising after the tone curve over the whole image, for every seed: held on both scenes.
• Both denoisers beat the noisy rendering: held on both scenes.
• The same ranking holds in each brightness band: held in all three bands on both scenes.
• The forward round trip of the $f$-unprocessed highlights is closer with the clip than without: held.

The pilot run was approved without changes, and nothing changed after the first run. The plan said that a neutral highlight clipped at the sensor “takes on the inverse white-balance gains”; by construction it takes on the gains themselves, as Part B prints. The script runs in about a second, and two runs gave identical output.

Common Misconceptions

• “There is a standard order of stages.” The sources disagree, and the table above lists four orders. What can be said is which swaps are harmless (per-channel gains and per-channel linear filters) and which are not (gains and clipping, gains and cross-channel filters, curves and spatial filters).
• “White balance is just a multiplication, so it can go anywhere.” It commutes with per-channel linear stages, but not with clipping or with demosaicers that borrow detail across channels, as Parts A and B show.
• “Undoing the gamma recovers raw data.” The gamma is one of several stages. The tone curve, the colour matrix, the gains, local tone mapping and clipped values also stand between a processed image and raw data. Inverse pipelines invert a model of the pipeline with sampled parameters (digital camera shows that undoing a standard gamma alone leaves the response non-linear).
• “A hardware ISP and a software pipeline produce the same image.” In HDR+ the hardware viewfinder and the software photograph differ visibly, and the authors say so (Hasinoff et al., 2016, Section 2).
• “Denoising can always be left to the end.” The noise a late denoiser sees has been correlated and reshaped by the stages before it. In the example, a white-noise denoiser underestimated it after the tone curve.

Where It Is Used

• Cameras and phones. Every digital camera runs a pipeline, in hardware, in software, or split between them (digital camera, computational photography). Burst systems such as HDR+ insert align-and-merge before the finishing stages (burst photography).
• Raw converters and stage algorithms. Each stage is its own field: demosaicing, color constancy for the white-balance gains, image denoising, tone mapping and high dynamic range imaging, gamma correction, image sharpening and image compression. The pipeline decides what input each of them receives.
• Training data for raw-domain networks. Inverse pipelines turn processed photographs into synthetic raw training data for denoising (Brooks et al., 2019; Dewil, 2022). Raw-to-sRGB networks for extreme low light, which replace much of the pipeline with one learned mapping, are described in low-light image enhancement.
• Vision systems. Methods that measure light, compare intensities across images or assume a noise model depend on which stages ran before them. Image formation follows the physical chain, and Digital Camera lists where processing breaks the assumptions of vision methods.