Demosaicing

Demosaicing, also spelled demosaicking (the spelling of most of the papers cited here), reconstructs a full-colour image from the output of a single-chip colour image sensor. A colour filter array (CFA) over the sensor lets each photosite measure only one of red, green and blue, so the raw image is a mosaic with one value per pixel. Demosaicing estimates the two missing values at every pixel. It runs on every photograph from a single-sensor camera, and its errors appear as colour fringes, “zipper” patterns along edges and false colour in fine texture. Because two thirds of the values in a demosaiced image are estimated rather than measured, the method also decides how the noise of neighbouring pixels and channels becomes correlated.

The simplest method interpolates each channel from its own samples, which image formation shows on bars. This article covers the methods that use the other channels: gradient-corrected linear filters, directional interpolation, colour-difference and residual interpolation, non-local methods, and convolutional networks trained for demosaicing alone or jointly with denoising. A worked example compares four of them on images with known ground truth.

Problem Definition

On the Bayer pattern, green photosites form a quincunx (checkerboard) lattice covering half of the array, and red and blue photosites each cover a quarter, spaced every two pixels. Rows alternate between “red rows” (R, G, R, G, …) and “blue rows” (G, B, G, B, …) (Getreuer’s IPOL article on Malvar–He–Cutler, 2011a, Section 1). Jin, Guo, Morel and Facciolo (2021, eq. 1) write the mosaic as a masked sum of the full-colour image,

$$Q = M_R \cdot R + M_G \cdot G + M_B \cdot B,$$

where $M_R$, $M_G$ and $M_B$ are binary masks that are 1 where the CFA passes that colour, and the products are taken pixel by pixel. Demosaicing estimates $(R, G, B)$ at every pixel from $Q$ and the masks. Ehret and Facciolo (2019, Section 1) describe it as an ill-posed inverse problem: the mosaic has a third as many values as the image it must produce.

Sampling explains why the problem is hard. Red and blue are sampled every second row and column, so along the axes their Nyquist limit is a quarter of a cycle per pixel, a period of 4 pixels; green, on its quincunx lattice, resolves twice that frequency along the axes (it follows from the lattice spacing). Detail between the two limits is still recorded in green but aliases in red and blue (sampling and aliasing). A method that interpolates each channel separately can only produce the aliased red and blue, and the mismatch with green appears as colour. The methods beyond bilinear interpolation therefore borrow detail from the densely sampled green, or from all channels at once, under an assumption about how the channels are related. Jin et al. (2021, Section 1) state it as strong correlation between the channels, which amounts to the differences between R, G and B being smooth.

Input and Output

• Input. A single-channel mosaic of linear raw values, normally after black-level subtraction (the image signal processing pipeline describes where demosaicing sits among the other stages and how sources disagree on its order relative to white balance), together with the CFA layout and its phase: which of the four positions of the $2 \times 2$ Bayer tile holds red. The four phases are named here by the top-left tile read row by row, RGGB, GRBG, GBRG and BGGR; Jin et al.’s demonstration lets the user select the Bayer pattern (Section 8). A demosaicer run with the wrong phase interpolates every channel from the wrong samples.
• Convention in this article. RGGB with red at row 0, column 0: even rows are red rows and odd rows blue rows. This is the convention of the code in image formation and digital camera. Arrays are indexed [row, column]. Jin et al. (2021, Section 3.1) index pixels as $(i, j)$ with $i$ the horizontal direction, and Getreuer’s (2011a) formula for red at a green site in a red row implies the same, so in their formulas $i$ is the column.
• Output. A full-resolution RGB image in the camera’s own colour space, which later stages white-balance, colour-correct and tone-map. Jin et al.’s Hamilton–Adams algorithms and the four methods of the example below return the CFA value itself at the pixel where it was measured.

Variants

• Other filter arrays. Fujifilm’s X-Trans pattern repeats on a $6 \times 6$ tile instead of $2 \times 2$ (Ehret & Facciolo, 2019, Fig. 2). Methods written for the Bayer tile do not apply to it directly; Gharbi et al.’s network, for instance, has a separate X-Trans variant without its initial downsampling, which its authors judged too aggressive for a $6 \times 6$ pattern (as Ehret & Facciolo, 2019, Section 2.1, describe).
• Joint demosaicking and denoising. Raw data is noisy, and demosaicing changes the noise (image noise): image denoising notes that raw data has a simple noise model but needs denoising before demosaicing. Joint methods solve both problems at once; the noisy variant of Gharbi et al.’s network takes the noise standard deviation as an extra input channel (Ehret & Facciolo, 2019, Section 2.1), and Zhao et al.’s study of losses for joint demosaicking and denoising networks is summarized in loss function and SSIM.
• Multi-frame reconstruction. A burst of slightly shifted raw frames samples each colour at more positions than one frame. Wronski et al. (2019) replace demosaicing with a merge of raw burst frames directly onto an RGB grid, as computational photography and burst photography describe.

Challenges

• Aliasing. Detail above the red and blue Nyquist limit cannot be recovered from those channels alone. Getreuer’s Zhang–Wu article shows a fence whose slats red and blue sample well below the Nyquist rate, on which bilinear interpolation reaches a PSNR of 20.42 dB against 38.29 dB for the Zhang–Wu method (Getreuer, 2011b, Section 5).
• Artifacts that PSNR averages away. Getreuer (2011b, Section 5) names the typical artifacts: zipper patterns (sequences of adjacent pixels of alternating brightness), false colours and desaturation. Gharbi et al. argued that L2 and PSNR are not much affected by demosaicking artifacts, and that because difficult patches are rare a network would hardly learn to avoid them without a training set biased towards them (as Ehret & Facciolo, 2019, Section 2.2, report).
• Channel correlation varies. Duran and Buades (2015, Section 1) contrast the Kodak images, which have few saturated colours and high inter-channel correlation, with the IMAX collection, which has many more saturated colours and edges between coloured regions, and note that methods designed for one set give state-of-the-art results on it but usually poor results on the other. Getreuer (2011b, Section 5) reports the same for the Zhang–Wu method: it works very well on Kodak, which he calls somewhat grey and blurred, and less well on IMAX and McMaster; on the peppers image, whose colour edges violate the correlation assumption, bilinear interpolation scores higher than Hamilton–Adams, Gunturk et al. and Zhang–Wu.
• Noise. Correction terms borrowed from other channels are differences of samples, which add noise. The example below measures how much.

Approaches

Bilinear interpolation

Each missing value is the average of the nearest samples of the same colour: green at a red or blue site from its four axial neighbours, red and blue at a green site from their two neighbours, and red at a blue site (or blue at red) from the four diagonal neighbours (Getreuer, 2011a, Section 2, gives the axial and diagonal cases). It is the linear interpolation of each channel on its own lattice, it is exact on constant and linear images, and it ignores the other channels. Getreuer (2011a, Section 2) calls its quality generally quite poor: because the channels are interpolated independently, misalignments near edges give strong colour distortions and zipper artifacts.

Gradient-corrected linear interpolation: Malvar–He–Cutler

Malvar, He and Cutler (2004) kept the method linear but added cross-channel Laplacian corrections, following Pei and Tam (as Getreuer, 2011a, Section 2, describes; the original paper was not read). In Getreuer’s notation, green at a red site is

$$\hat G(i,j) = \hat G_{\text{bl}}(i,j) + \alpha\, \Delta_R(i,j), \qquad \Delta_R(i,j) = R(i,j) - \tfrac14\big[R(i-2,j) + R(i+2,j) + R(i,j-2) + R(i,j+2)\big],$$

where $\hat G_{\text{bl}}$ is the bilinear estimate and $\Delta_R$ a five-point Laplacian of the red samples. Red at a green site adds $\beta$ times a nine-point Laplacian of the green samples to the bilinear red, and red at a blue site adds $\gamma$ times the five-point Laplacian of blue; blue is symmetric. The gains were set to the values that minimized the mean squared error over the Kodak suite and rounded to dyadic rationals, $\alpha = \tfrac12$, $\beta = \tfrac58$, $\gamma = \tfrac34$, so that the filters can be computed with integers and bit shifts; Getreuer reports that they approximate the optimal $5 \times 5$ Wiener filters within 5% in mean squared error (Getreuer, 2011a, Section 2). The result is eight $5 \times 5$ filters, of four shapes (Getreuer, 2011a, Fig. 3, coefficients scaled by 8). Green at a red or blue site, and red at a green site in a red row, are

$$\frac18\begin{bmatrix} 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\end{bmatrix}, \qquad \frac18\begin{bmatrix} 0&0&\tfrac12&0&0\\ 0&-1&0&-1&0\\ -1&4&5&4&-1\\ 0&-1&0&-1&0\\ 0&0&\tfrac12&0&0\end{bmatrix}.$$

The second filter’s two 4s fall on the red neighbours to the left and right; its transpose serves where the wanted colour lies above and below, which for red is a green site in a blue row. Blue at a green site uses the same two shapes with the rows exchanged: in a blue row the blue neighbours are to the left and right. Getreuer also writes the second filter in integer form, a sum with weights 1, −2, 8 and 10 divided by 16 (Getreuer, 2011a, Section 2).

Two details of Getreuer’s figure matter to an implementer. Its bottom row labels the left-and-right shape “B at green locations in red rows” and the above-and-below shape “B at green locations in blue rows”, read from the rendered page. In a red row (R, G, R, G) the horizontal neighbours of a green site are red, not blue, so these two labels appear to be exchanged; the example below runs the filters as labelled and shows that they are then not exact on linear images. Second, Getreuer compares the method with Hamilton–Adams on one image and finds it visually sharper than bilinear interpolation but not as sharp as Hamilton–Adams, which he attributes to Hamilton–Adams’s nonlinear strategy (Getreuer, 2011a, Section 3, Fig. 5).

Directional interpolation: Hamilton–Adams

Laroche and Prescott chose between horizontal and vertical interpolation of green by comparing gradients, and Hamilton and Adams (1997) added second-order corrections from the red and blue channels (as Jin et al., 2021, Section 1, describe; the patent was not read). Jin et al. (2021, Section 3) derive the method from a second-order Taylor expansion and the assumption that, along a red row, the second horizontal derivatives of red and green are nearly equal. With $i$ horizontal, the horizontal estimate of green at a red site is

$$\tilde G_H(i,j) = \frac{G(i-1,j) + G(i+1,j)}{2} - \frac{R(i-2,j) - 2R(i,j) + R(i+2,j)}{4},$$

and $\tilde G_V$ is the vertical analogue. A classifier measures how flat each direction is,

$$CL_H(i,j) = \big|G(i+1,j) - G(i-1,j)\big| + \big|R(i-2,j) - 2R(i,j) + R(i+2,j)\big|,$$

and the green estimate is $\tilde G_H$ where $CL_H < CL_V$, $\tilde G_V$ where $CL_H > CL_V$, and their mean on ties (Jin et al., 2021, Section 3.1, Algorithm 1). Red and blue then follow with the interpolated green as a guide. In Jin et al.’s Algorithm 2, red at a green site in a red row is the average of its two red neighbours minus a quarter of the three-point horizontal Laplacian of the interpolated green (the vertical analogue in a blue row). Red at a blue site is the same along the diagonal or the anti-diagonal, whichever has the smaller classifier; blue is symmetric (Section 3.2).

Read from the rendered pages, Jin et al.’s derivation does not reach its own formula exactly. Their eq. (8) gives green as the average of its neighbours minus half the second derivative, and eq. (9) estimates the second derivative of red as $(R(i-2,j) - 2R(i,j) + R(i+2,j))/4$. Combining the two gives a coefficient of $\tfrac18$ on the red difference, while the stated formula and Algorithm 1 use $\tfrac14$. In eq. (10), the second line writes the same estimate with the kernel $\tilde\Delta_H = [-1\ 0\ 2\ 0\ -1]$, which reverses the sign of the correction relative to the first line; Algorithm 1 uses $[1, 0, -2, 0, 1]$ and agrees with the first line. The example below implements Algorithm 1 as printed. In our reading the factor $\tfrac14$, like Malvar–He–Cutler’s fitted gains, is larger than the Taylor value, so neither method is exact on a quadratic image, which the example confirms.

Colour differences and directional LMMSE: Zhang–Wu

Zhang and Wu (2005) work with the primary difference signals $G - R$ and $G - B$, which they assume to be smooth (as Getreuer, 2011b, Section 4, describes). Along each row and column, second-order Laplacian interpolation (the formula of Hamilton–Adams’s directional estimate) gives a first estimate of the difference; each directional estimate is denoised by linear minimum mean-square-error (LMMSE) estimation with local statistics, and the horizontal and vertical estimates are fused with weights that minimize the estimated error, which yields green. Red and blue then come from the differences: $G - R$, known at red sites, is averaged from four diagonal neighbours at blue sites and then from four axial neighbours at green sites, and red is green minus the interpolated difference (Getreuer, 2011b, Sections 4.1–4.3). The example below uses this last step with Hamilton–Adams green.

Threshold-free and residual interpolation

Jin et al. (2021) trace a line of methods that all interpolate green first:

• GBTF (gradient-based threshold-free interpolation, Pekkucuksen and Altunbasak) computes Hamilton–Adams estimates of the colour differences in four directions, north, south, east and west, and combines them with weights instead of a hard choice, so that a pixel near an edge can draw on the one side where interpolation is reliable (Section 4).
• RI (residual interpolation, Kiku et al., 2013) keeps GBTF’s framework but replaces its Hamilton–Adams step with guided-filter interpolation improved by a residual technique (Sections 1 and 5).
• MLRI and MLRI+wei compute the tentative estimates by minimizing the Laplacian of the image with a guided filter, and the latter takes a weighted average of the guided filter’s linear coefficients (Section 1).
• ARI (adaptive residual interpolation) combines RI and MLRI adaptively and chooses the number of iterations at each pixel (Section 1).

On their own implementations, Jin et al. (2021, Section 7.2, Table 3) report that GBTF, RI, MLRI and MLRI+wei all improve on Hamilton–Adams and that their average difference in colour PSNR over the four datasets is less than 0.25 dB, and that ARI is the best of the interpolation methods. The table spreads wider than that sentence suggests: without noise the four averages run from 38.44 to 38.89 dB, against 35.03 dB for Hamilton–Adams, and on the rows labelled Imax GBTF falls slightly below Hamilton–Adams at $\sigma = 3$ and 10. ARI is also the slowest: in their Table 4, Hamilton–Adams takes 0.30 s for a $500 \times 500$ image and ARI 17.17 s.

Non-local and adaptive methods

Self-similarity can resolve what local neighbourhoods cannot. Duran and Buades (2015, Sections 1–3) first interpolate green along four directions, reconstruct red and blue from each by bilinear interpolation of the colour differences, and choose among the results by their chromatic smoothness; a parameter sets how much inter-channel correlation the interpolation assumes, chosen automatically from the image’s chromatic gradient. A second step filters the colour differences, rather than the channel values, with non-local means adapted in the same way. LSSC (Mairal et al.) is a dictionary-learning method that learns its parameters for every image; Jin et al. (2021, Section 7.2) find its PSNR almost the same as ARI’s, at 734 s for a $500 \times 500$ image.

Learned demosaicing

Convolutional networks trained on full-colour images and their simulated mosaics score higher than the hand-designed methods on noise-free mosaics in both IPOL comparisons cited here. Two are reviewed in detail by Ehret and Facciolo (2019), who reran both with the authors’ weights:

• Gharbi et al. (2016) downsample the mosaic into a four-channel image, apply 15 convolutional layers, upsample to an RGB image of full size, concatenate it with the mosaic split into its channels, and finish with two more convolutions; the network is trained with an L2 loss. After pretraining on 1.3 million ImageNet images and 1 million MirFlickr images, it was fine-tuned on patches mined for difficulty: luminance artifacts found with the HDR-VDP-2 metric and moiré found as excess Fourier energy in the demosaicked patch (Ehret & Facciolo, 2019, Sections 2.1–2.2).
• Tan et al. (2017) start from a bilinear interpolation and add two residual stages, the first estimating green and the second the full image, with a loss that adds a green-only term to the full-image term (Ehret & Facciolo, 2019, Section 3).

In Ehret and Facciolo’s Table 1 (PSNR in dB, read from the rendered page), Tan et al. reach 41.99 on Kodak and Gharbi et al. 39.08 on McMaster, against 41.44 and 37.48 for the best hand-designed entries in those columns (LSSC on Kodak, ARI on McMaster). The authors conclude that the CNNs beat the best human-crafted methods by almost 2 dB while being an order of magnitude faster; the pooled “All” column gives 38.56 dB for Tan et al. against 36.86 dB for LSSC, the best hand-designed entry. Their Table 2 qualifies the second claim: on one Kodak image, on a CPU, both networks take about 14 s, against 61 s for ARI and 736 s for LSSC but under 2 s for RI and MLRI. Jin et al. (2021, Table 3) reran the two networks through Ehret and Facciolo’s implementations, as RCNN (Tan et al.) and JCNN (Gharbi et al.). Without noise, JCNN has the best average colour PSNR over their four datasets, 41.53 dB against 39.82 dB for ARI. At noise $\sigma = 10$, JCNN, trained for noisy input, leads with 33.66 dB, and ARI (29.20 dB) is ahead of RCNN (28.64 dB). Their text says that at $\sigma = 10$ and above ARI performs better than the deep-learning methods, but their table shows that this holds against RCNN and not against JCNN.

Examples

Four demosaicers on images with known ground truth

The script compares four methods on the same RGGB mosaics:

1. Bilinear interpolation of each channel.
2. MHC: the four Malvar–He–Cutler filter shapes of Getreuer (2011a, Fig. 3), with blue at a green site in a blue row using the left-and-right shape.
3. HA: Hamilton–Adams green (Jin et al., 2021, Algorithm 1) and Hamilton–Adams red and blue with the Laplacian of the interpolated green (Algorithm 2).
4. HA+CD: the same Hamilton–Adams green, then red and blue by bilinear interpolation of the colour differences $G - R$ and $G - B$ (Getreuer, 2011b, Section 4.3).

How to complete red and blue after Hamilton–Adams green is a choice the sources leave open: Jin et al.’s Algorithm 2 uses Laplacian corrections, while Zhang–Wu, GBTF and Duran and Buades use colour differences. The script therefore prints both. No library in SciPy, scikit-image 0.25.2 or Pillow provides a demosaicer, so the implementations are checked against inputs with known answers. The script has five parts:

• Part A checks exactness on a constant colour, on three independent linear ramps (one per channel) and on a quadratic bowl shared by the channels. It also runs Getreuer’s Fig. 3 with its bottom-row labels taken literally.
• Part B checks Malvar–He–Cutler’s integer form against its kernel form. It compares sums of integers divided by 16, Getreuer’s written-out formula for red at a green site in a red row, and the decomposition of each filter into bilinear interpolation plus a gain times a Laplacian.
• Part C measures colour PSNR (CPSNR) on the four bundled colour photographs of scikit-image.
• Part D measures false colour on a grey zone plate, whose local frequency rises linearly from the centre to half a cycle per pixel at the edge of a disc. Red and blue reach their Nyquist limit along the axes at a quarter of a cycle per pixel.
• Part E adds white Gaussian noise to the mosaic before demosaicing. It measures each method’s variance gain on flat grey, compares it for the linear methods with the gain predicted from the filter taps, and repeats the CPSNR measurement at two noise levels with three seeds.

CPSNR follows Jin et al. (2021, Section 7.1): ten times the base-10 logarithm of the squared peak value over the mean of the three channels’ mean squared errors, which is the common colour definition of PSNR. Jin et al. state no border crop. Dropping 8 pixels at each border is our choice, made because every method’s estimates near the border depend on how the image is extended; the script prints the mean without the crop as well. Outputs are clipped to $[0, 1]$ before they are scored. The boundary rule is SciPy’s mode="mirror", a reflection about the edge pixel, which maps a pixel at $-1$ to $+1$ and so keeps the CFA phase of the reflected samples.

The bundled images (astronaut, coffee, chelsea and rocket, with odd sizes trimmed to even) are themselves processed photographs; rocket is stored as a JPEG, and scikit-image’s documentation does not say how any of them was captured or processed. They are exact ground truth for the simulated mosaic, not true full-colour scenes, and their numbers are not comparable with Kodak or McMaster results. The script uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 and runs in about 2 seconds.

We wrote the checks and hypotheses in a plan before the script was first run, computing each criterion’s reachable value with an earlier version of the same script. The pilot run printed one value, the sum of a Laplacian’s taps, in full precision; that line was changed to print it in scientific notation with two significant digits before the first run, a print-only change the editor approved.

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
from scipy import ndimage
from skimage import data

# Bayer RGGB, as in Image Formation and Digital Camera: rows 0, 2, ... are
# "red rows" (R G R G ...), rows 1, 3, ... are "blue rows" (G B G B ...).
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")

# 1. Bilinear: each channel from its own samples
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):
    return np.stack([corr(cfa * m["R"], K8), corr(cfa * m["G"], K4),
                     corr(cfa * m["B"], K8)], -1)

# 2. 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      # G at R or B
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  # wanted colour left/right
K_COL = K_ROW.T                                               # wanted colour above/below
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  # R at B, B at R
def mhc(cfa, m, as_labelled=False):
    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)))
    if as_labelled:      # Fig. 3's bottom row taken literally: row form in red rows
        B = np.where(m["B"], cfa, np.where(m["Gr"], row, np.where(m["Gb"], col, dia)))
    else:                # in a blue row the blue neighbours are left and right
        B = np.where(m["B"], cfa, np.where(m["Gb"], row, np.where(m["Gr"], col, dia)))
    return np.stack([R, G, B], -1)

# 3. Hamilton-Adams, Jin et al. (2021), Section 3, Algorithms 1 and 2
def ha_green(cfa, m):
    h = lambda k: np.asarray(k, float)[None, :]
    v = lambda k: np.asarray(k, float)[:, None]
    lap5, kk, d = [1, 0, -2, 0, 1], [0, 1, 0, 1, 0], [-1, 0, 1]
    gi = [0.5 * a - 0.25 * b for a, b in zip(kk, lap5)]
    GH, GV = corr(cfa, h(gi)), corr(cfa, v(gi))
    CH = np.abs(corr(cfa, h(d))) + np.abs(corr(cfa, h(lap5)))
    CV = np.abs(corr(cfa, v(d))) + np.abs(corr(cfa, v(lap5)))
    G = np.where(CV > CH, GH, np.where(CV < CH, GV, (GH + GV) / 2))
    return np.where(m["G"], cfa, G)

def ha_rb(cfa, m, G, c):              # c = "R" or "B": Algorithm 2, and its blue analogue
    X = cfa * m[c]
    near_h = m["Gr"] if c == "R" else m["Gb"]   # green sites with X left and right
    near_v = m["Gb"] if c == "R" else m["Gr"]
    other = m["B"] if c == "R" else m["R"]
    KH, LH = np.array([[1, 0, 1]]) / 2, np.array([[1, -2, 1]]) / 4
    KP, LP = np.diag([1, 0, 1]) / 2, np.diag([1, -2, 1]) / 4
    KN, LN = np.fliplr(KP), np.fliplr(LP)
    DP, DN = np.diag([-1, 0, 1]), np.fliplr(np.diag([-1, 0, 1]))
    XH = corr(X, KH) - corr(G, LH)
    XV = corr(X, KH.T) - corr(G, LH.T)
    XP = corr(X, KP) - corr(G, LP)
    XN = corr(X, KN) - corr(G, LN)
    CP = np.abs(corr(X, DP)) + np.abs(corr(G, 4 * LP))
    CN = np.abs(corr(X, DN)) + np.abs(corr(G, 4 * LN))
    XD = np.where(CP > CN, XN, np.where(CP < CN, XP, (XN + XP) / 2))
    return np.where(m[c], cfa, np.where(near_h, XH, np.where(near_v, XV,
                    np.where(other, XD, 0))))

def ha(cfa, m):
    G = ha_green(cfa, m)
    return np.stack([ha_rb(cfa, m, G, "R"), G, ha_rb(cfa, m, G, "B")], -1)

# 4. Hamilton-Adams green, then colour differences (Getreuer's Zhang-Wu, Section 4.3)
def ha_cd(cfa, m):
    G = ha_green(cfa, m)
    out = [None, G, None]
    for i, c, o in ((0, "R", "B"), (2, "B", "R")):
        D = (G - cfa) * m[c]                       # G - X, known at X sites
        D = D + corr(D, [[1, 0, 1], [0, 0, 0], [1, 0, 1]]) / 4 * m[o]
        D = D + corr(D, [[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4 * m["G"]
        out[i] = G - D
    return np.stack(out, -1)

METHODS = {"bilinear": bilinear, "MHC": mhc, "HA": ha, "HA+CD": ha_cd}
def demosaic(name, cfa, m):
    return np.clip(METHODS[name](cfa, m), 0, 1)

def cpsnr(x, ref, crop=8):
    s = (slice(crop, -crop),) * 2 if crop else (slice(None),) * 2
    return 10 * np.log10(1 / np.mean((x[s] - ref[s]) ** 2))

# A. Checks: exactness on constants and independent ramps, quadratic, MHC forms
H = W = 64; m = masks(H, W); y, x = 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 * x + 0.001 * y, 0.4 - 0.002 * x + 0.004 * y,
                 0.6 + 0.001 * x - 0.003 * y], -1)
quad = np.repeat((0.3 + 2e-4 * ((x - 32) ** 2 + (y - 32) ** 2))[..., None], 3, -1)
inner = (slice(6, -6), slice(6, -6))
print("A. max |error| in the interior (6-px border dropped)")
print(f"{'method':<24}{'constant':>10}{'ramps':>10}{'quadratic':>11}")
tests = list(METHODS.items()) + [("MHC, Fig. 3 as labelled", lambda c, mm: mhc(c, mm, True))]
for name, f in tests:
    e = [np.abs(f(mosaic(t, m), m) - t)[inner].max() for t in (const, ramp, quad)]
    print(f"{name:<24}{e[0]:>10.1e}{e[1]:>10.1e}{e[2]:>11.1e}")

img = data.astronaut()[:512, :512]
mi = masks(512, 512)
F = mosaic(img.astype(np.int64), mi)                  # integer CFA, 0..255
kern = mhc(F.astype(float) / 1, mi)                     # kernel form, floats
ints = {k: np.rint(16 * K).astype(np.int64) for k, K in
        (("G", K_G), ("row", K_ROW), ("col", K_COL), ("dia", K_DIA))}
def icorr(k): return ndimage.correlate(F, ints[k], mode="mirror")  # integer sums
iR = np.where(mi["R"], 16 * F, np.where(mi["Gr"], icorr("row"),
              np.where(mi["Gb"], icorr("col"), icorr("dia"))))
iG = np.where(mi["G"], 16 * F, icorr("G"))
iB = np.where(mi["B"], 16 * F, np.where(mi["Gb"], icorr("row"),
              np.where(mi["Gr"], icorr("col"), icorr("dia"))))
integer = np.stack([iR, iG, iB], -1)
print("B. MHC integer form (sums of integers, then /16) against kernel form:",
      f"max diff {np.abs(integer / 16 - kern).max():.1e};",
      f"all 16x sums integers: {all(np.array_equal(v, np.rint(16 * v / 16)) for v in ints.values())}")
# Getreuer's written-out formula for R at green in a red row, (i, j) = (column, row)
P = np.pad(F, 2, mode="reflect")                        # numpy 'reflect' = scipy 'mirror'
def Fp(di, dj): return P[2 + dj:2 + dj + 512, 2 + di:2 + di + 512]
num = (Fp(0, -2) + Fp(0, 2)
       - 2 * (Fp(-1, -1) + Fp(1, -1) + Fp(-2, 0) + Fp(2, 0) + Fp(-1, 1) + Fp(1, 1))
       + 8 * (Fp(-1, 0) + Fp(1, 0)) + 10 * Fp(0, 0))
print("   written formula (R at G, red row) vs kernel form:",
      f"max diff {np.abs(num / 16 - kern[..., 0])[mi['Gr']].max():.1e};",
      f"bit-shift >>4 vs exact /16: max diff {np.abs((num >> 4) - num / 16)[mi['Gr']].max():.4f}")
# Laplacian form: bilinear plus a gain times the Laplacian of the measured channel
lap5 = np.array([[0, 0, -.25, 0, 0], [0] * 5, [-.25, 0, 1, 0, -.25], [0] * 5, [0, 0, -.25, 0, 0]])
Gbl = np.pad(K4, 1); Gbl[2, 2] = 0                     # four axial neighbours / 4
Dbl = np.pad(K8, 1) * 0; Dbl[1:4:2, 1:4:2] = 0.25      # four diagonal neighbours / 4
Hbl = np.zeros((5, 5)); Hbl[2, 1] = Hbl[2, 3] = 0.5    # left and right neighbours / 2
print("   K_G  = bilinear + 1/2 Laplacian:", np.allclose(K_G, Gbl + 0.5 * lap5),
      "| K_DIA = bilinear + 3/4 Laplacian:", np.allclose(K_DIA, Dbl + 0.75 * lap5))
L9 = (K_ROW - Hbl) / (5 / 8)
print("   (K_ROW - bilinear) / (5/8): centre", L9[2, 2], f"sum {L9.sum():.1e}",
      "non-zero taps", int((L9 != 0).sum()))

# C. CPSNR on bundled images (8-px border dropped)
names = ["astronaut", "coffee", "chelsea", "rocket"]
imgs = {}
for n in names:
    a = getattr(data, n)().astype(float) / 255
    imgs[n] = a[:a.shape[0] // 2 * 2, :a.shape[1] // 2 * 2]
print("C. CPSNR (dB), 8-px border dropped")
print(f"{'image':<12}" + "".join(f"{k:>10}" for k in METHODS))
res = {}
for n, a in imgs.items():
    mm = masks(*a.shape[:2]); c = mosaic(a, mm)
    res[n] = [cpsnr(demosaic(k, c, mm), a) for k in METHODS]
    print(f"{n:<12}" + "".join(f"{v:>10.2f}" for v in res[n]))
print(f"{'mean':<12}" + "".join(f"{np.mean([res[n][i] for n in names]):>10.2f}" for i in range(4)))
print(f"{'mean, crop 0':<12}" + "".join(
    f"{np.mean([cpsnr(demosaic(k, mosaic(a, masks(*a.shape[:2])), masks(*a.shape[:2])), a, 0) for a in imgs.values()]):>10.2f}"
    for k in METHODS))

# D. False colour on a grey zone plate; local frequency f = k r cycles/pixel
N = 512; yy, xx = np.mgrid[0:N, 0:N] - N // 2
r = np.hypot(xx, yy); k = 0.5 / (N // 2)
zp = np.repeat((0.5 + 0.4 * np.cos(np.pi * k * r ** 2))[..., None], 3, -1)
mz = masks(N, N); cz = mosaic(zp, mz); f = k * r
bands = [(0, .125), (.125, .25), (.25, .375), (.375, .5)]
print("D. data check: grey", np.array_equal(zp[..., 0], zp[..., 2]) and np.array_equal(zp[..., 0], zp[..., 1]),
      f"| range {zp.min():.3f}-{zp.max():.3f} | pixels per band",
      [int(((f >= a) & (f < b)).sum()) for a, b in bands])
print("   zone plate: mean false colour (max - min over channels), by local frequency")
print(f"{'method':<11}" + "".join(f"{f'{a}-{b}':>12}" for a, b in bands))
for name in METHODS:
    out = demosaic(name, cz, mz); fc = out.max(-1) - out.min(-1)
    print(f"{name:<11}" + "".join(f"{fc[(f >= a) & (f < b)].mean():>12.4f}" for a, b in bands))

# E. Noise: variance gain on flat grey, and CPSNR with noise added to the mosaic
taps = {"G at R/B": K_G, "row form": K_ROW, "column form": K_COL, "R at B / B at R": K_DIA}
print("E. MHC filters, sum of squared taps:",
      ", ".join(f"{n} {(K ** 2).sum():.4f}" for n, K in taps.items()))
print(f"   bilinear: axial 4-avg {(K4[K4 != 1] ** 2).sum():.4f}, 2-avg 0.5000, diagonal 4-avg 0.2500")
pred = {"bilinear": np.mean([1, .5, .5, .25, 1, 1, .25, .25, 1, .5, .5, .25]),
        "MHC": np.mean([1, (K_ROW ** 2).sum(), (K_COL ** 2).sum(), (K_DIA ** 2).sum(),
                        1, 1, (K_G ** 2).sum(), (K_G ** 2).sum(),
                        1, (K_ROW ** 2).sum(), (K_COL ** 2).sum(), (K_DIA ** 2).sum()])}
print("   predicted mean variance gain from taps (3 channels x 4 sites):",
      ", ".join(f"{k} {v:.4f}" for k, v in pred.items()))
S = 256; ms = masks(S, S); flat = np.full((S, S, 3), 0.5); sig = 0.02
print("   output noise variance / input variance, flat grey, sigma 0.02, seeds 0-2")
for name in METHODS:
    g = []
    for seed in range(3):
        n = np.random.default_rng(seed).normal(0, sig, (S, S))
        out = METHODS[name](mosaic(flat, ms) + n, ms)
        g.append(((out - 0.5)[8:-8, 8:-8] ** 2).mean() / sig ** 2)
    print(f"   {name:<10} {min(g):.4f} - {max(g):.4f}   (input noise std, seed 2: {n.std():.4f})")
print("   mean CPSNR over the four images (dB), min - max over seeds 0-2")
for s in (0.01, 0.03):
    row = []
    for name in METHODS:
        vals = []
        for seed in range(3):
            rng = np.random.default_rng(seed); v = []
            for a in imgs.values():
                mm = masks(*a.shape[:2]); c = mosaic(a, mm) + rng.normal(0, s, a.shape[:2])
                v.append(cpsnr(demosaic(name, c, mm), a))
            vals.append(np.mean(v))
        row.append(f"{name} {min(vals):.3f}-{max(vals):.3f}")
    print(f"   sigma {s}: " + "; ".join(row))

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2):

A. max |error| in the interior (6-px border dropped)
method                    constant     ramps  quadratic
bilinear                   0.0e+00   1.1e-16    4.0e-04
MHC                        1.1e-16   2.2e-16    2.0e-04
HA                         5.6e-17   1.1e-16    2.0e-04
HA+CD                      5.6e-17   2.2e-16    2.0e-04
MHC, Fig. 3 as labelled    4.0e-01   3.6e-01    2.0e-04
B. MHC integer form (sums of integers, then /16) against kernel form: max diff 0.0e+00; all 16x sums integers: True
   written formula (R at G, red row) vs kernel form: max diff 0.0e+00; bit-shift >>4 vs exact /16: max diff 0.9375
   K_G  = bilinear + 1/2 Laplacian: True | K_DIA = bilinear + 3/4 Laplacian: True
   (K_ROW - bilinear) / (5/8): centre 1.0 sum -1.1e-16 non-zero taps 9
C. CPSNR (dB), 8-px border dropped
image         bilinear       MHC        HA     HA+CD
astronaut        30.48     34.53     34.73     35.65
coffee           29.44     33.17     33.08     33.87
chelsea          33.96     38.43     37.67     39.22
rocket           29.99     31.02     31.84     32.26
mean             30.97     34.29     34.33     35.25
mean, crop 0     31.00     34.34     34.37     35.29
D. data check: grey True | range 0.100-0.900 | pixels per band [12849, 38580, 64348, 90080]
   zone plate: mean false colour (max - min over channels), by local frequency
method          0-0.125  0.125-0.25  0.25-0.375   0.375-0.5
bilinear         0.0344      0.1500      0.3036      0.4182
MHC              0.0147      0.0387      0.1225      0.2991
HA               0.0165      0.0669      0.1316      0.3649
HA+CD            0.0065      0.0192      0.0312      0.2761
E. MHC filters, sum of squared taps: G at R/B 0.5625, row form 0.9922, column form 0.9922, R at B / B at R 0.9531
   bilinear: axial 4-avg 0.2500, 2-avg 0.5000, diagonal 4-avg 0.2500
   predicted mean variance gain from taps (3 channels x 4 sites): bilinear 0.5833, MHC 0.9167
   output noise variance / input variance, flat grey, sigma 0.02, seeds 0-2
   bilinear   0.5757 - 0.5868   (input noise std, seed 2: 0.0199)
   MHC        0.9067 - 0.9178   (input noise std, seed 2: 0.0199)
   HA         0.7287 - 0.7399   (input noise std, seed 2: 0.0199)
   HA+CD      0.9093 - 0.9188   (input noise std, seed 2: 0.0199)
   mean CPSNR over the four images (dB), min - max over seeds 0-2
   sigma 0.01: bilinear 30.638-30.639; MHC 33.213-33.221; HA 33.442-33.449; HA+CD 33.916-33.925
   sigma 0.03: bilinear 28.735-28.739; MHC 29.163-29.180; HA 29.760-29.773; HA+CD 29.378-29.391

What the checks show

All four methods reproduce constants and independent linear ramps to rounding error. This follows from their construction: each correction is a discrete Laplacian of one channel, which vanishes on a linear function, and each method uses only samples of the channel it corrects with. None is exact on the quadratic bowl. Read as labelled, Getreuer’s Fig. 3 gives an error of 0.36 on the ramps and 0.40 on the constant colour, the size of the difference between channel values. In that reading, blue at a green site is interpolated from red neighbours, which is consistent with the two labels being exchanged.

The integer form of Malvar–He–Cutler and its kernel form give the same output at every pixel, and so does Getreuer’s written-out formula. A right shift by 4 bits rounds down and therefore differs from exact division by up to 15/16 of a level. The green and red-at-blue filters are exactly bilinear interpolation plus $\tfrac12$ and $\tfrac34$ of the five-point Laplacian. Subtracting bilinear interpolation from the red-at-green filter and dividing by $\tfrac58$ leaves nine taps that sum to zero with a centre of 1, the “nine-point Laplacian” Getreuer names without printing its coefficients.

What the comparison shows

Every method that uses the other channels beats bilinear interpolation on every image, by 1.0 to 2.3 dB on rocket and by 3.6 to 5.3 dB on the other three. Malvar–He–Cutler and Hamilton–Adams with Laplacian-corrected red and blue split the four images between them, with nearly equal means. Hamilton–Adams green followed by colour differences is best on all four, by 0.4 to 1.6 dB over the same green with Jin et al.’s red and blue. Rocket gains least from any method. Dropping the 8-pixel border changes the printed means by 0.03 to 0.05 dB and leaves the order unchanged.

What the zone plate shows

On the grey zone plate every method makes little false colour below a quarter of a cycle per pixel and much more above it, where red and blue alias. Bilinear interpolation makes the most in every band; between 0.25 and 0.375 cycles per pixel it makes about 10 times as much as HA+CD. HA+CD makes the least in every band. On a grey image $G - R$ and $G - B$ are exactly zero, so the colour-difference step adds no colour of its own, and any false colour comes from errors in the green estimate. In the top band, from 0.375 to 0.5 cycles per pixel, all methods produce strong false colour. In our reading this is because green aliases there too: its quincunx lattice resolves frequencies with $|f_x| + |f_y| \le \tfrac12$, which along the diagonals is only about 0.35 cycles per pixel.

What the noise shows

Linear filters change the variance of white noise by the sum of their squared taps, as image sharpening explains. Here the gain differs from one CFA site to another: measured samples keep gain 1, while Malvar–He–Cutler’s interpolated values have gains from 0.56 to 0.99, against 0.25 to 0.5 for bilinear interpolation. Averaged over the three channels and four site types, the measured gain on flat grey matches the prediction for both linear methods. Malvar–He–Cutler multiplies the noise variance by about 0.92, against 0.58 for bilinear interpolation. Hamilton–Adams, which is nonlinear and has no taps to sum, measures about 0.73, and HA+CD about 0.91.

The CPSNR results follow. Malvar–He–Cutler’s 3.3 dB lead over bilinear interpolation without noise falls to 2.6 dB at $\sigma = 0.01$ (about 2.55 levels of 255) and to under 0.5 dB at $\sigma = 0.03$. The two methods with the lowest variance gain lose the least as the noise grows. At $\sigma = 0.03$ the order of the two Hamilton–Adams variants is reversed: Laplacian-corrected red and blue now beats colour differences. The ranges over three seeds are under two hundredths of a dB, far smaller than the differences between methods. A method that wins on clean mosaics therefore does not necessarily win on noisy ones; Jin et al. (2021, Section 7.2) likewise add noise at several levels to test the robustness of the methods they compare.

Outcomes against the plan

The plan fixed five checks, data checks and six hypotheses before the first run. Its criteria came from values that the earlier version of the script reached on the same set-up, so a “held” confirms the derivations and the implementations rather than an independent prediction. All held:

• Checks.
  • Exactness on constants and independent ramps within $10^{-12}$ for all four methods: held (largest error $2.2 \times 10^{-16}$).
  • Errors above $10^{-6}$ on the quadratic for every method: held.
  • Getreuer’s Fig. 3 read as labelled has an error above 0.01 on the ramps: held (0.36).
  • The integer, written-out and kernel forms of Malvar–He–Cutler agree within $10^{-12}$, its filters equal bilinear interpolation plus $\tfrac12$ and $\tfrac34$ Laplacians, and the implied nine-point Laplacian has centre 1, zero sum and nine taps: held.
  • The measured variance gains lie within 0.02 of the tap prediction for bilinear interpolation and Malvar–He–Cutler: held.
  • The data checks (grey zone plate, value range, populated bands, input noise level) held.
• Hypotheses.
  • MHC, HA and HA+CD each beat bilinear interpolation on all four images: held.
  • Neither MHC nor HA wins on all four images: held.
  • With the same green, colour-difference red and blue beat Laplacian-corrected red and blue on all four images: held.
  • On the zone plate, every method has more false colour in each band above a quarter cycle per pixel than in each band below: held. Bilinear interpolation has the most false colour in every band: held. HA+CD has the least in every band: held.
  • Malvar–He–Cutler’s variance gain exceeds bilinear interpolation’s: held. Its CPSNR lead over bilinear interpolation shrinks as $\sigma$ rises from 0 to 0.01 to 0.03: held.
  • Bilinear interpolation and Hamilton–Adams, the two methods with the lowest measured variance gain, lose the least CPSNR from $\sigma = 0$ to $\sigma = 0.03$: held.

Datasets and Benchmarks

• Kodak. 24 noiseless full-colour images of $768 \times 512$ or $512 \times 768$ pixels, downsampled from $3072 \times 2048$ originals (Jin et al., 2021, Section 7.1). Malvar–He–Cutler’s gains were fitted on it (Getreuer, 2011a), so in our reading results on Kodak partly reflect that fit.
• McMaster and IMAX. McMaster has 18 images of $500 \times 500$ pixels cropped from $2310 \times 1814$ originals (Jin et al., 2021, Section 7.1). Jin et al.’s tables label these rows “Imax”, while their text names the McMaster dataset. Duran and Buades (2015, Section 1) describe an IMAX collection with many more saturated colours and edges between coloured regions than Kodak, and Getreuer (2011b, Section 5) names IMAX and McMaster as two datasets, so the sources do not settle whether the two names denote the same set.
• Microsoft Demosaicing Dataset. 57 Canon and 500 Panasonic images of $210 \times 318$ pixels, built from raw images taken with a Canon EOS 550D and a Panasonic Lumix DMC-LX3 (Jin et al., 2021, Section 7.1).
• Flickr500 subset. 14 images from Syu et al.’s Flickr500, used by Ehret and Facciolo (2019, Section 4).

Protocols differ in ways that change the numbers. Getreuer (2011b, Section 5) removes a 5-pixel margin, Ehret and Facciolo (2019, Section 4) pad each image by 48 pixels with symmetry before mosaicking and remove the padding before scoring, and Jin et al. state no crop. The same method can therefore score differently in different papers: ARI has 39.90 dB colour PSNR on Kodak in Jin et al.’s Table 3 and 39.24 dB PSNR on Kodak in Ehret and Facciolo’s Table 1, from different implementations and protocols. Compare numbers only within one table.

Evaluation Metrics

• Colour PSNR (CPSNR). The PSNR of the three channels’ mean squared error, averaged over channels before the logarithm (Jin et al., 2021, Section 7.1); see PSNR and mean squared error. Per-channel PSNR is also reported. Jin et al. (2021, Tables 1–2) give green and red separately and state that blue behaves like red.
• Visual artifacts. Zipper patterns, false colour and moiré are judged on crops of difficult images. Getreuer’s IPOL articles and Ehret and Facciolo (2019, Figs. 6–9) show such crops. Gharbi et al. mined training patches for these artifacts with a perceptual metric and a Fourier criterion (as Ehret & Facciolo, 2019, Section 2.2, describe).
• Noise robustness. Results at several noise levels added to the mosaic; Jin et al. (2021, Section 7.2) use $\sigma = 1$, 3 and 10, with CPSNR computed for a peak of 255.
• Run time. Reported per image size and machine (Jin et al., 2021, Table 4; Ehret & Facciolo, 2019, Table 2).

Applications

• Camera pipelines. Every single-sensor camera demosaics, in hardware or in software; digital camera places the stage in the processing chain and notes that camera demosaicers are mostly proprietary, and image sharpening that sharpening is often combined with demosaicing.
• Raw processing and computational photography. Raw converters demosaic on a computer, and burst pipelines can merge raw frames before or instead of demosaicing (computational photography).
• Vision on processed images. Methods that assume independent pixel noise (image noise) or exact colour at edges inherit the demosaicer’s choices.
• Training data for restoration. Simulated mosaics of full-colour images train and test demosaicers and joint demosaicking–denoising networks: Gharbi et al. and Tan et al. made their training mosaics by masking (Ehret & Facciolo, 2019, Sections 2.2 and 3.2), with the caveat that the reference images were themselves processed.

Open Problems

• Difficult content. Ehret and Facciolo (2019, Figs. 7 and 9) show that every compared method, the CNNs included, still leaves some moiré on shutters, and that only Gharbi et al.’s network does well on their hardest example.
• Noise and real raw data. In Jin et al.’s Table 3 the advantage of RCNN, a learned method without a denoising function, over ARI shrinks as noise grows and reverses at $\sigma = 10$, and the authors observe that the gap between deep-learning and interpolation methods narrows at higher noise. Joint methods trained for noise lead there. Evaluations that add Gaussian noise to mosaics of finished photographs do not model real sensor noise (image noise).
• Ground truth. The test images are processed photographs, Kodak’s downsampled from larger originals and McMaster’s cropped from them (Jin et al., 2021, Section 7.1), not the scenes themselves. In our reading, scores on simulated mosaics measure how well a method reproduces these processed references as much as how well it reconstructs a scene.
• Other arrays and multi-frame capture. Non-Bayer filter arrays such as X-Trans need their own methods or networks, and burst methods that bypass demosaicing change what a single-frame demosaicer is for.