BookNotesSlidesP-sets Syllabusdraft · v0.2.49.1
◀▲⌂▶
💬Comments welcome. To leave a note, select any text and click the note / highlight button that pops up — or open the panel with the tab at the top-right (‹). Notes are visible only inside our private review group.

15.2 Advanced traditional demosaicking⧉

The methods below, published between the late 1980s and the mid-2000s, combine the two ideas of Chapter 15.1, Basic traditional demosaicking, edge direction and color differences, in different proportions, and hand-designed methods of this family are still what most cameras and phones run in their image signal processors, where a demosaicker must run on every preview and video frame, fast, small and predictable. Learned demosaicking has shipped mainly in desktop raw converters, such as Adobe's Enhance Details and DxO's DeepPRIME, which also denoises (Chapter 15.6, Learned demosaicking). Some phone pipelines avoid single-frame demosaicking altogether by merging a burst of raw frames directly onto a full-color grid (Chapter 14.2, Multi-frame super-resolution). Figure 15.2.1 runs ten of them side by side on the same mosaic, including the residual interpolation and ARI of this chapter, so each method's assumption can be seen where it fails.

fig-algo-lab-demosaic
Figure 15.2.1. Interactive: the demosaicking lab. One Bayer mosaic reconstructed by up to four methods at once, chosen from ten classical methods: nearest neighbor, bilinear, smooth hue transition, Freeman's median of color differences (on bilinear, or on the chapter's edge-based demosaicking), Malvar–He–Cutler, Hamilton–Adams, variable-number-of-gradients (VNG), homogeneity-directed (AHD), residual interpolation (RI) and adaptive residual interpolation (ARI), in a mouse-driven four-way split with a magnifier on the pointer, per-pane parameter sliders and live PSNR and SSIM against the original. The last six can be followed by Freeman's median, 3 × 3 or 5 × 5. A pane can instead show the ground truth, the raw mosaic, the difference of two panes, or a pane's amplified error; a sensor-noise slider corrupts the mosaic before reconstruction, a resolution slider gives the sensor fewer pixels over the same scene, and the zone-plate and Siemens-star inputs make false color and zipper unmissable; the zone plate and the star are band-limited: their finest rings and spokes fade to gray at the sensor's sampling limit, so the ground truth carries no aliasing of its own.

Nearest neighbor copies each missing value from a nearby sample of the right color, with no arithmetic. The name is loose, because most sites have several equally near samples of a given color: a red site has four green neighbors at the same distance, and four blue ones on the diagonals. An implementation breaks the tie with a fixed rule, usually by working tile by tile. In the lab's version, every pixel of a 2 × 2 tile takes that tile's one red and one blue, a green site keeps its own green, and the other two sites take the same one of the tile's greens. That makes it a close relative of the downsampling of Downsampling and bilinear interpolation, and how they fail: downsampling turns each tile into one pixel, and nearest neighbor is almost that same small image blown back up to full size by copying each value into the four pixels of its tile. The only differences are that the measured greens stay where they were and nothing is averaged. Nearest neighbor therefore keeps the pixel count but adds almost no detail over downsampling, and it keeps downsampling's flaw: each color is attributed to a place up to a pixel away from where it was measured, so every edge becomes a colored staircase, the strongest zipper of any method here. It costs nothing, which is why it survives in live previews and thumbnails. Bilinear interpolation (Downsampling and bilinear interpolation, and how they fail) averages instead of copying and fixes the misregistration in smooth regions, not at edges.

Figure 15.2.2 works through every method of this section on a small mosaic, one step per stage, starting with nearest neighbor.

fig-demosaick-methods-steps
Figure 15.2.2, stage 2. The methods of this section, one step at a time. A small scene is area-sampled onto a grid and kept through an RGGB mosaic; the middle grid shows the numbers of the current step, bold where measured and italic where interpolated, and each method's result joins the column on the right once its last step is reached, with its mean error against the scene. This stage is nearest neighbor: each 2 × 2 tile's red and blue fill the whole tile, a green site keeps its green, and the tile's top-right green fills its red and blue sites. The input is a corner, a sine grating, a Siemens star or a zone plate, in three pairs of colors, on a grid of 8 × 8 with the values or of 16 × 16 or 32 × 32 without them, which the star and the zone plate need. The zone plate and the star are band-limited: their finest rings and spokes fade to gray at the sensor's sampling limit, so the ground truth carries no aliasing of its own. Interactive: choose the input, its colors and the grid, set the corner's or the grating's angle, the grating's period, the star's spokes or the zone plate's highest frequency, drag the corner or the pattern's center, and step through each method with the stage buttons.

Smooth hue transition Cok 1987 (smooth hue transition), one of the first methods to use one channel to reconstruct another, interpolates green first and then interpolates hue rather than red and blue: the ratio $R/G$ at the red sites, spread to its neighbors and multiplied back by green. It assumes that hue changes slowly across an object even where brightness does not. The choice of hue coordinate matters, and the lab lets you switch it: a ratio becomes unstable where green is near zero, a difference $R-G$ does not, and a log ratio sits in between. The difference form is the recipe of Interpolate differences, not channels.

Its three steps are the smooth-hue stages of Figure 15.2.2.

fig-demosaick-methods-steps
Figure 15.2.2, stage 4. Smooth hue transition: $R/G$ and $B/G$ at the red and blue sites, against the bilinear green; the ratios interpolated; and multiplied back by green. Beside a black region the green is near zero and the ratio is unstable; the lab's guard, one level added to both, keeps it finite. Interactive: the same controls; the stages step through the three steps.

Freeman's method Freeman 1988 (median demosaicking) demosaicks with bilinear interpolation, then cleans up the hue: it applies a 3 × 3 median to the difference channels $R-G$ and $B-G$ and adds green back. Within an object the color difference is nearly constant, so a false-color pixel is an isolated outlier in a flat channel, which is exactly what a median removes. A 5 × 5 median often does a better job, because false color on fine texture tends to come in clusters wider than a 3 × 3 window can outvote, though a larger median also starts to carry color across object boundaries. The median does nothing for zippering in the green channel, so although the method is described on top of bilinear, it is mostly used as a clean-up step on top of a better demosaicker. Figure 15.2.1 offers both bases: on bilinear, the median removes most of the fringes and leaves the zippering; on the chapter's edge-based demosaicking, the green is already good and the median only has the remaining false color to remove. The same pass survives as the final clean-up of many later methods, AHD among them. (Freeman's patent differs in its details: it describes two colors alternating along a line, interpolates each linearly, and applies a one-dimensional median to their difference, five pixels wide in its example, with each pixel keeping its own measured value. The two-dimensional Bayer version above is the one usually implemented.)

Figure 15.2.2 shows the three steps on the grid.

fig-demosaick-methods-steps
Figure 15.2.2, stage 7. Freeman's method in the Bayer form the text describes: $R-G$ and $B-G$ from the bilinear planes at every pixel, drawn as the color they add to gray; their 3 × 3 median; and green added back. A fourth stage puts the median on the chapter's edge-based demosaicking instead of bilinear and adds that green back. There the false color comes in bands along slanted edges and in patches on fine patterns rather than in lone pixels, so the median there is 5 × 5, which removes most of what a 3 × 3 median keeps. Interactive: the same controls; the stages step through the three steps, then the median on the edge-based result.

Gradient-corrected linear interpolation Malvar et al. 2004, by Malvar, He and Cutler, keeps a fixed linear filter but adds the Laplacian of the channel measured at the center pixel, scaled by a gain, to the bilinear estimate of each missing channel. A red pixel's curvature, for example, corrects the green estimated there, on the assumption that the channels vary together. The result is a set of fixed 5 × 5 filters with no decisions, so it is fast, easy to run on a GPU, and free of the unstable direction choices that make mazes; it is, for example, what MATLAB's demosaic function implements. Because it does not adapt, it still blurs across edges where an adaptive method would not.

fig-malvar-kernels
Figure 15.2.3. The four Malvar–He–Cutler filters, in eighths: green at a red or blue site, red at a green site in a red row and in a blue row, and red at a blue site (blue at a red site is the same by symmetry). Each cell shows its weight and the mosaic color it applies to. The bilinear part is shaded blue and the Laplacian correction from the center channel orange, with the published gains α = 1/2, β = 5/8 and γ = 3/4.

The same correction, pixel by pixel, is in Figure 15.2.2.

fig-demosaick-methods-steps
Figure 15.2.2, stage 11. Malvar–He–Cutler: the correction each missing channel receives, the Laplacian of the measured channel scaled by the published gains; then the bilinear estimate plus that correction. Near an edge the correction overshoots, and the values leave the range 0 to 1. Interactive: the same controls, and a checkbox that follows the result with Freeman's 5 × 5 median; the stages step through the two steps.

Hamilton–Adams Hamilton & Adams 1997, met in Interpolate along the edge, not across it, combines both ideas of the basic chapter, and it does so in three passes over the mosaic, each using the channels already filled in.

  1. Green at the red and blue sites. Measure a horizontal and a vertical activity at the site, each the sum of a green gradient across the site and the second difference of the center color along the same line. Interpolate green along the less active direction, as the average of its two green neighbors plus a quarter of the center color's second difference, and use all four neighbors when the two activities are equal. This is the pseudocode of Interpolate along the edge, not across it.
  2. Red and blue at the green sites. Now that green exists everywhere, a missing red at a green site is the average of its two red neighbors, horizontal or vertical depending on the row, corrected by half of green's second difference along that line. It is color-difference interpolation in disguise: the result equals the site's green plus the average of the two neighbors' $R-G$, up to the correction.
  3. Red at the blue sites, and blue at the red sites. The four nearest samples of the missing color sit on the diagonals, so the method chooses between the two diagonals with the same kind of test, a gradient of the missing color along the diagonal plus green's second difference along it, and interpolates along the calmer one with the same half-weight correction.
Algorithm 15.2.1
# Hamilton-Adams, steps 2 and 3 (step 1, green, is in the edge-directed section)
# step 2: red at a green site in a red row (blue rows: use the vertical pair)
R[x,y] = (R[x-1,y] + R[x+1,y])/2 + (2*G[x,y] - G[x-1,y] - G[x+1,y])/2

# step 3: red at a blue site, choosing a diagonal
dN = abs(R[x-1,y-1] - R[x+1,y+1]) + abs(2*G[x,y] - G[x-1,y-1] - G[x+1,y+1])   # NW-SE
dP = abs(R[x+1,y-1] - R[x-1,y+1]) + abs(2*G[x,y] - G[x+1,y-1] - G[x-1,y+1])   # NE-SW
if   dN < dP: R[x,y] = (R[x-1,y-1] + R[x+1,y+1])/2 + (2*G[x,y] - G[x-1,y-1] - G[x+1,y+1])/2
elif dP < dN: R[x,y] = (R[x+1,y-1] + R[x-1,y+1])/2 + (2*G[x,y] - G[x+1,y-1] - G[x-1,y+1])/2
else:         R[x,y] = (R[x-1,y-1] + R[x+1,y+1] + R[x+1,y-1] + R[x-1,y+1])/4 \
                       + (4*G[x,y] - G[x-1,y-1] - G[x+1,y+1] - G[x+1,y-1] - G[x-1,y+1])/4
# blue at a red site: the same, with R and B exchanged

The corrections are what separate it from bilinear. Each one assumes that the missing channel curves the same way as the channel that was measured or already reconstructed at the site, so a sharp edge in green is copied into red and blue at the same place instead of being smeared into a fringe. The direction tests are what separate it from Malvar–He–Cutler, which applies the same kind of correction everywhere, in every direction. Hamilton–Adams costs a small, fixed amount of work per pixel, and its green step is the starting point of many later methods, including the residual interpolation of the next chapter; its failures are those of any two-way choice: diagonal texture, and fine patterns where both directions look equally active.

Figure 15.2.2 runs the three passes on the grid; the grating near 45 degrees shows the failure.

fig-demosaick-methods-steps
Figure 15.2.2, stage 13. Hamilton–Adams in its three passes: green along the calmer axis with the curvature correction, the bar marking the axis averaged; red and blue at the green sites from the pair on that line, corrected by half of green's second difference; and red at the blue sites and blue at the red ones along the calmer diagonal, the bar marking it. Interactive: the same controls, and a checkbox that follows the result with Freeman's 5 × 5 median; the stages step through the three passes.

Variable number of gradients Chang, Cheung & Pan 1999 (VNG), VNG, by Chang, Cheung and Pan, replaces the two-way choice with a vote among eight directions, and lets the image decide how many of them to trust. At each pixel it works in the 5 × 5 neighborhood around the site, in four steps.

  1. Eight gradients. For each compass direction, north, northeast, east and so on, sum the absolute differences between pairs of same-color samples that lie along that direction on that side of the pixel. Pairing same-color samples keeps the gradient free of the differences between channels. A direction that points into a smooth region gets a small gradient; one that crosses an edge gets a large one.
  2. A threshold. From the smallest gradient $g_{\min}$ and the largest $g_{\max}$, set $T = 1.5\,g_{\min} + 0.5\,(g_{\max}-g_{\min})$. The threshold is relative: it asks which directions are nearly as smooth as the smoothest one here, not which are below a fixed number, so the same rule works in shadows and highlights.
  3. A vote. Keep every direction whose gradient is below $T$. For each kept direction, average the samples of each color that lie on that side of the pixel, and sum these averages over the kept directions, one sum per color.
  4. A color difference. Estimate each missing color from the center sample and the difference of the sums: at a red site, for instance, $\hat G = R + (S_G - S_R)/n$, where $S_G$ and $S_R$ are the sums for green and red over the $n$ kept directions. This is the color-difference idea of Interpolate differences, not channels, averaged over the directions the vote allowed.
Algorithm 15.2.2
# VNG at one pixel, in outline (per-site sample weights omitted)
g = {d: gradient(d) for d in [N, NE, E, SE, S, SW, W, NW]}       # same-color pairs along d
T = 1.5*min(g.values()) + 0.5*(max(g.values()) - min(g.values()))
keep = [d for d in g if g[d] <= T]
S = {c: sum(average_on_side(c, d) for d in keep) for c in [R, G, B]}  # one sum per color
for c in missing_colors:
    value[c] = center_value + (S[c] - S[center_color]) / len(keep)

The number of directions that pass is what gives the method its name and its behavior. In a flat region all eight pass and it averages widely, which suppresses noise; along a straight edge only the two directions parallel to it pass, and it averages along the edge like Hamilton–Adams; at a corner or a line end it can keep three or four, which a two-way choice cannot express (Figure 15.2.4). The cost is about eight gradients per pixel over a 5 × 5 window, several times Hamilton–Adams, and the relative threshold can pass too many directions in fine, low-contrast texture, where everything looks equally smooth. VNG was long one of the modes of the dcraw raw decoder.

fig-vng-directions
Figure 15.2.4. Variable number of gradients at one red site of a real mosaic, the NO PARKING sign of problem set 3. The eight compass gradients are drawn as arrows scaled by their size, with the threshold $T = 1.5\,\min + 0.5\,(\max-\min)$ marked. In a flat patch all eight pass and the method averages in every direction; on a letter stroke only north and south pass, the directions along the stroke. The gradients follow the published idea, not its exact per-site weights.

Figure 15.2.2 runs VNG on every pixel of the grid.

fig-demosaick-methods-steps
Figure 15.2.2, stage 16. Variable number of gradients: at each pixel, how many of the eight directions fall below the threshold and which, drawn as spokes; then the missing colors as the center sample plus the color differences averaged over them. The gradients follow the published idea, not its exact per-site weights. Interactive: the same controls, and a checkbox that follows the result with Freeman's 5 × 5 median; the stages step through the two steps.

Two further refinements are common. Iteration: once green and the color differences are in, re-estimate the differences against the interpolated green and repeat. Alternating projections Gunturk, Altunbasak & Mersereau 2002 (alternating projections) formalize this as projecting in turn onto two sets, the images that agree with the measured samples and the images whose red and blue have the same fine detail as green. And directional estimation with explicit noise models: Zhang and Wu's directional linear minimum mean-square-error method Zhang & Wu 2005 (DLMMSE) fuses the horizontal and vertical estimates with weights set by how noisy each one is, rather than choosing one.

Across the eight panes, three things change from method to method. The assumption about images: smooth channels, smooth hue, correlated channels, or edges along one of a few directions. The decision made at each pixel: none for bilinear and Malvar–He–Cutler, a direction for Hamilton–Adams, VNG and AHD (#Adaptive homogeneity-directed demosaicking (AHD), next). And the cost, from nothing for nearest neighbor to two full reconstructions and a color-space conversion for AHD. The residual-interpolation methods (#Residual interpolation and ARI) take the next step by replacing the color difference with the residual of a guided fit.

For reproducible comparisons, state the implementation and demosaicking mode. LibRaw and dcraw expose several algorithms, but their defaults are not interchangeable with the Malvar–He–Cutler method. Compare against a strong classical baseline such as adaptive residual interpolation Monno et al. 2015 (adaptive residual interpolation) as well as a simple linear method.

The sections that follow take up three families that go further: AHD, which builds whole reconstructions in two directions and keeps the more homogeneous one; residual interpolation and ARI, which interpolate what a guided prediction leaves over; and frequency-domain demultiplexing, which treats the mosaic as a modulated signal. The chapter ends with the situations that defeat all of them.

15.2.1 Adaptive homogeneity-directed demosaicking (AHD)⧉

If you want to implement one classical demosaicker that is simple and works well, this is probably your best bet, followed by Freeman's median of the color differences (earlier in this chapter) to remove the false color that remains. Every method so far decides the interpolation direction before interpolating, from gradients in the mosaic. Adaptive homogeneity-directed demosaicking Hirakawa & Parks 2005 (AHD), by Hirakawa and Parks, decides afterward. It reconstructs the whole image twice, once interpolating only horizontally and once only vertically, and at each pixel keeps whichever reconstruction looks more plausible there. The criterion is the key idea: interpolating in the wrong direction averages across an edge and creates colors that exist on neither side of it, so the wrong candidate is less homogeneous, locally, than the right one.

  1. Two candidates. Build a full image with every interpolation along rows: green with a Hamilton–Adams-style correction, then red and blue by interpolating color differences. Build a second with every interpolation along columns.
  2. A perceptual space. Convert both candidates to CIELAB, so that lightness and color can be compared separately and a distance means roughly what a viewer would notice (Chapter 4.5, Advanced color: opponents, perceptually-uniform, and preferred).
  3. Homogeneity. For each pixel of each candidate, count the neighbors in a small window that are close to it in both lightness and color, within thresholds $\varepsilon_L$ and $\varepsilon_C$. The thresholds are set at every pixel from the differences to its nearest neighbors along each candidate's own direction, so the test adapts to local contrast rather than using one number for shadows and highlights.
  4. The choice. Sum each candidate's counts over a 3 × 3 window and keep, at each pixel, the candidate with the larger sum; average the two on a tie.
  5. Clean-up. Finish with a few passes of a median over the color differences $R-G$ and $B-G$, Freeman's clean-up.
Algorithm 15.2.3
# AHD, in outline
Ih, Iv = demosaic_along_rows(mosaic), demosaic_along_columns(mosaic)
Lh, Lv = to_cielab(Ih), to_cielab(Iv)
eps_L, eps_C = adaptive_thresholds(Lh, Lv)             # per pixel, from nearest-neighbor differences
Hh = count_close_neighbors(Lh, eps_L, eps_C)           # lightness and color both within threshold
Hv = count_close_neighbors(Lv, eps_L, eps_C)
Hh, Hv = box_sum_3x3(Hh), box_sum_3x3(Hv)
out = where(Hh > Hv, Ih, where(Hv > Hh, Iv, (Ih + Iv) / 2))
out = median_color_differences(out, passes=3)

Figure 15.2.5 shows the homogeneity test at one pixel.

fig-ahd-homogeneity-criterion
Figure 15.2.5. AHD's homogeneity criterion at one pixel. Left, the mosaic. Then the horizontal and the vertical candidate, the ringed pixel's four neighbors outlined green when they are within both thresholds, $\varepsilon_L$ in lightness and $\varepsilon_C$ in color (CIELAB), and red when they are not; beneath each, the differences, the count and its sum over 3 × 3. The thresholds are the same for both candidates, set from the pixel's own neighbors along each candidate's direction. Right, the result: each pixel framed orange where the horizontal candidate won, blue where the vertical one did, gray on a tie. Next to an edge, the candidate that interpolated across it has made colors that neither side has, and its neighbors stop matching. The inputs are the corner, a sine grating, a Siemens star and a zone plate (both band-limited), the maze's colored stripes, and crops of four photographs, on a grid of 10 × 10 or 16 × 16. Interactive: choose the input and the grid, set the pattern's angle, period, spokes or frequency, drag on the mosaic to move the pattern or pan the photograph, and click a pixel to test it.

The figure below runs it on a real facade, where the cornice's horizontal lines choose the horizontal candidate and the window columns the vertical one (Figure 15.2.6), and Figure 15.2.2 makes the same choice cell by cell on the chapter's small inputs.

fig-ahd-homogeneity
Figure 15.2.6. AHD on a crop of the Prudential facade, from a mosaic sampled from the photograph as in the lab: the horizontal candidate, the vertical candidate, their homogeneity maps, the per-pixel choice and the result. The cornice's horizontal lines choose the horizontal candidate and the window columns the vertical one. Where the candidates differ, the choice picks the more accurate one at about two pixels in three, and the result's mean error is below either candidate's. Each picture enlarges on its own.

Figure 15.2.2 steps through the same algorithm on the small inputs.

fig-demosaick-methods-steps
Figure 15.2.2, stage 18. AHD: the horizontal candidate, the vertical candidate, the choice at each pixel from their homogeneity counts summed over 3 × 3, and the result. The published method finishes with a median of the color differences; the checkbox adds one, Freeman's 5 × 5. Interactive: the same controls, and a checkbox that follows the result with Freeman's 5 × 5 median; the stages step through the four steps.

AHD costs two full reconstructions, two color-space conversions and the counting, several times Hamilton–Adams. It was the default of the dcraw raw decoder and of many converters built on it for years, and it is still a strong classical baseline. Its limits are the other methods' limits. It still chooses between two directions, so diagonal texture gives both candidates errors and the choice cannot fix them. In fine texture near the sampling limit both candidates are inhomogeneous and the choice becomes erratic, which is where mazes come from. And noise lowers the homogeneity of both candidates everywhere, so the choice degrades as the noise rises, one reason to denoise first or jointly (Chapter 15.8, Joint demosaicking and denoising).

15.2.2 Residual interpolation and ARI⧉

The best-performing classical demosaickers, and the ones to compare against today, come from the residual-interpolation family, and adaptive residual interpolation (ARI) is its strongest member. It appeared a decade after AHD, outperforms it clearly in published comparisons on the standard test sets, and needs no training. It is also built directly on the guided filter of Chapter 10.6, Guided image filtering, so it is a payoff of that chapter rather than a new tool. It takes more work to implement than AHD, which remains the simpler choice when that matters.

The idea changes what gets interpolated. Every classical method interpolates something, and each step of the classical line chose a flatter quantity. Bilinear interpolates the channel itself and fails at edges. Hamilton–Adams and its successors interpolate the color difference $G-R$, which is smooth wherever the two channels differ by an offset, and fails where that assumption does: at a saturated color boundary, at a specular highlight, wherever the ratio between the channels changes rather than their offset. Residual interpolation Kiku et al. 2013 (residual interpolation), by Kiku, Monno, Tanaka and Okutomi, goes one step further. Rather than hope that some quantity is flat, it constructs one: it fits a local prediction of the sparse channel from the dense one and interpolates only what the prediction gets wrong.

$$ \underbrace{R}_{\text{full of structure}} \quad\to\quad \underbrace{G - R}_{\text{smoother, if the offset model holds}} \quad\to\quad \underbrace{R - \tilde R}_{\text{potentially smoother}} $$

Take red, sampled on a quarter of the pixels, and green, already reconstructed everywhere because it is sampled twice as densely. In a small window, model red as a linear function of green,

$$ \tilde R(\mathbf p) = a\, G(\mathbf p) + b, $$

with $a$ and $b$ estimated from the red samples that exist in that window. This is a guided filter with green as the guide He, Sun & Tang 2013 (guided filter), restricted to the measured sites. In a window $\Omega$, with $\mathcal R$ the set of sites where red was measured, solve

$$\min_{a,b}\sum_{q\in\Omega\cap\mathcal R}(R_q-aG_q-b)^2+\epsilon a^2,$$

where $\epsilon>0$ stabilizes the slope. Unmeasured sites must not enter the fit as zeros; a window with no red sample has to grow or fall back to a simpler estimate; and since every pixel lies in many windows, the overlapping predictions are averaged, which is the guided filter's own rule. The tentative estimate $\tilde R$ is a full-resolution red built entirely out of green's structure, and it follows a change of contrast between the channels that a difference cannot. The residual is the part of red that green cannot explain,

$$ \Delta R = R - \tilde R, \qquad \text{known only where red was measured,} $$

and the reconstruction interpolates it, by any simple method, and adds it back:

$$ \hat R = \tilde R + \operatorname{interp}(\Delta R). $$

The prediction carries the shared detail, and the residual corrects its errors where measurements constrain them. When the local model explains the structure the channels share, the residual is smaller and flatter than the color difference, and the same cheap interpolator does much less damage to it. That flatness is an assumption to test, not a property the fit enforces: a small residual at the measured sites does not prove that the residual between them is small or smooth, and a misleading guide, a sharp color boundary or a poorly conditioned fit can transfer the wrong structure into the result.

fig-residual-interpolation-scanline
Figure 15.2.7. Interactive. One scanline across a hard color boundary, with four candidate quantities plotted on the same axes: the red channel, the color difference $G-R$, the ratio $R/G$, and the residual $R - \tilde R$ against a fitted guide. The channel swings, the difference is flat until the boundary and then jumps, the ratio is worse where green is small, and the residual stays flat across the boundary. Under the plot, the reconstruction each choice produces on the full image, with its error map, so the flatness of a curve and the quality of a result are on screen together. A scene selector moves the scanline onto a saturated flower, a specular highlight and a neutral edge, which are the three places the older assumptions fail. Interactive: choose the scene and the scanline, toggle the four quantities, set the guide window size, and compare reconstructions and error maps.

The full algorithm applies the same idea twice, first to green and then to red and blue.

  1. Green first. Estimate green along rows and along columns separately, each time guiding by the measured red or blue along that line and interpolating the residual along it, then combine the two directional estimates with weights that favor the direction whose color differences vary less locally.
  2. Tentative red. With the full green as the guide, fit a local linear model $\tilde R = aG + b$ in each window, using only the sites where red was measured, and evaluate it everywhere.
  3. Residual. At the measured red sites, form $R-\tilde R$.
  4. Interpolate and add back. Interpolate the residual bilinearly to every site and add it to $\tilde R$. Blue is the same with blue's samples.

On a red flower against green leaves, the residual is close to zero where the color difference is large, and residual interpolation scores 3.3 dB above AHD (Figure 15.2.8).

fig-ari-residual
Figure 15.2.8. Residual interpolation on saturated color edges: red bracts against green leaves, a crop of a photograph whose RGGB mosaic is sampled from it, so every result is scored against the original. For red, the tentative estimate $\tilde R = aG + b$ comes from a guided filter with the reconstructed green as guide, fitted on the red samples only. Beside it, on the same color scale, the color difference $G-R$ that earlier methods interpolate and the residual $R-\tilde R$ that residual interpolation interpolates: the difference is large wherever the flower is, while the residual is close to zero except at a few edge pixels, and its mean absolute Laplacian, the roughness a bilinear interpolator suffers from, is 6.0 against 9.2. The results: AHD 32.8 dB, residual interpolation 36.2 dB, MLRI 36.5 dB and ARI 37.7 dB, with error maps on one scale; the errors that remain sit on the petal edges. Residual interpolation and MLRI are our implementation from the papers; ARI is the authors' code, ported.

Minimized-Laplacian residual interpolation, MLRI Kiku et al. 2014 (minimized-Laplacian residual interpolation), changes only how $a$ and $b$ are fitted: instead of least squares on the values, it minimizes the Laplacian energy of the residual, which is what the final bilinear interpolation needs to be small. It changes how the regression coefficients are estimated, not which channel serves as the guide. ARI Monno et al. 2015 (adaptive residual interpolation), by Monno, Kiku, Tanaka and Okutomi, iterates both, using each output as the next guide, and at each pixel selects the variant and the number of iterations whose residual is smoothest, judged by the local fit. That is the direction-selection idea of Hamilton–Adams operating one level up: not which direction to interpolate in, but which method to trust here. A later version extends it to multispectral filter arrays Monno et al. 2017 (ARI, multispectral). The authors' reference code has been used widely enough in later work that comparisons against it can be trusted.

Algorithm 15.2.4
# Residual interpolation for red, green already reconstructed (blue is the same)
def residual_interpolation(G, R_mosaic, red_mask, radius, eps):
    a, b = masked_guided_fit(guide=G, target=R_mosaic, mask=red_mask, radius=radius, eps=eps)
    R_tent = box_mean(a, radius) * G + box_mean(b, radius)     # guided filter output
    residual = where(red_mask, R_mosaic - R_tent, 0)
    return R_tent + bilinear_from_samples(residual, red_mask)

# MLRI: the same, with a and b fitted to minimize the Laplacian of the residual
# ARI, in outline: at each pixel keep the (method, iteration) whose residual is smoothest
estimates = []
guide = G
for k in range(K):
    for method in [RI, MLRI]:
        R_k = method(guide, R_mosaic, red_mask)
        estimates.append((R_k, laplacian_energy(R_mosaic - R_k, red_mask)))
    guide = R_k
R = per_pixel_argmin(estimates)          # smoothest residual wins, in a local window

Figure 15.2.2 steps through residual interpolation cell by cell on the chapter's small inputs, from the directional green to ARI's per-pixel choice, and the demosaicking lab at the start of this chapter (Figure 15.2.1) runs the same code on photographs.

fig-demosaick-methods-steps
Figure 15.2.2, stage 22. Residual interpolation on the small inputs. Green is estimated along rows and along columns from residuals, and the two are weighted by which direction's color difference is smoother; a bar marks the heavier direction. The following stages show the tentative red $\tilde R = aG + b$ with the fitted $a$ and $b$ per cell, the residual beside the color difference $G-R$ at the red and blue sites, the result, and a map of ARI's choice at each pixel between RI and MLRI and between the first and second pass. On a colored edge the residual is smaller than the color difference. ARI here is the published idea in outline, the second pass guided by the first's tentative estimate and the smoothest residual chosen, and on inputs this small it does not always beat plain residual interpolation. Interactive: the same controls as the other stages; step through the five residual-interpolation stages, and the checkbox adds Freeman's median to the result.

The limits are the guided filter's. A local linear model cannot follow channels that are related nonlinearly, at clipping or under strong chromatic aberration, and fitting $a$ and $b$ from a few noisy samples is unstable, so noise should be removed first, accounted for within the fit, or handled by joint reconstruction (Chapter 15.8, Joint demosaicking and denoising).

This family, not bilinear interpolation, is the baseline a learned demosaicker should be measured against. It embodies three separate decisions: predict one channel from a guide, interpolate what the prediction missed, and choose locally between useful variants. A learned model can change any of the three, and a comparison with bilinear alone leaves those established improvements unaccounted for, so it cannot tell which part of a network's gain required learning. Compare actual accuracy, noise assumptions and runtime; no ranking of classical methods, and no fixed advantage in decibels, carries over across datasets and settings. A network can represent similar cross-channel predictions, but it chooses what to interpolate implicitly, which is one reason it is hard to say what it learned: a benchmark gain does not show that it learned this algorithm, and attributing the gain requires controlled experiments rather than an analogy between architectures (Chapter 15.6, Learned demosaicking).

15.2.3 Demosaicking as demultiplexing: frequency-domain methods⧉

One family of methods treats the mosaic the way a radio engineer treats a signal. As Chapter 15.3, Mosaics and aliasing shows, a Bayer mosaic is exactly a luma component $L$ plus two chroma components modulated onto carriers at the highest frequencies the grid supports: $C_1$ on the horizontal and vertical carriers, $C_2$ on the diagonal one,

$$y_{ij}=L_{ij}+C_{1,ij}\big[(-1)^i+(-1)^j\big]+C_{2,ij}(-1)^{i+j}.$$

Demosaicking is then demultiplexing: isolate each modulated chroma band with a bandpass filter, shift it back to baseband by multiplying by its carrier, and take the luma as what remains. Alleysson, Süsstrunk and Hérault Alleysson, Süsstrunk & Hérault 2005 (frequency-domain demosaicing) arrived at this from the side of human vision, where luma and chroma are carried separately; Dubois Dubois 2005 (frequency-domain demosaicking) made it adaptive. Nowhere else in this chapter is the algorithm itself the frequency argument, rather than a heuristic that happens to be analyzable in frequency terms.

  1. Diagonal chroma. Filter the mosaic with a bandpass filter centered on the diagonal carrier and multiply by $(-1)^{i+j}$ to bring $C_2$ back to baseband.
  2. Axial chroma, twice. $C_1$ is carried twice, on the horizontal and on the vertical carrier. Extract it from each with its own bandpass filter and demodulate.
  3. Choose the cleaner copy. Luma detail leaks into the two copies differently: a vertical texture puts its energy near the horizontal carrier and spoils that copy, while the vertical copy stays clean. Dubois weights the two estimates at each pixel by the local energy of the crosstalk in each, which is the edge-direction idea of this chapter moved into the frequency domain.
  4. Luma. Remodulate the chroma estimates and subtract them from the mosaic; what remains is $L$.
  5. Back to RGB. $G=L-C_2$, $R=L+C_2+2C_1$, $B=L+C_2-2C_1$.
Algorithm 15.2.5
# frequency-domain demosaicking (Dubois style), for an RGGB mosaic y
C2  = lowpass(bandpass(y, carrier=(pi, pi)) * alternate(i + j))
C1h = lowpass(bandpass(y, carrier=(pi, 0))  * alternate(i))
C1v = lowpass(bandpass(y, carrier=(0, pi))  * alternate(j))
eh, ev = local_energy(crosstalk_h), local_energy(crosstalk_v)   # luma leaking into each copy
C1  = (ev * C1h + eh * C1v) / (eh + ev)                         # trust the cleaner copy
L   = y - C1 * (alternate(i) + alternate(j)) - C2 * alternate(i + j)
G, R, B = L - C2, L + C2 + 2*C1, L + C2 - 2*C1

Figure 15.2.2 runs the same demultiplexing on the small inputs.

fig-demosaick-methods-steps
Figure 15.2.2, stage 27. Frequency-domain demosaicking on the small inputs, with the convention of the text: the magnitude of the mosaic's Fourier transform with the three carriers ringed, then the demodulated $C_2$, the two copies of $C_1$ with the weight given to the horizontal one, the luma left after subtracting the remodulated chroma, and the result. A vertical grating shows the point of the per-pixel choice: its luma leaks into the horizontal copy of $C_1$, and the weight moves to the vertical copy. The weights come from each copy's local energy of detail beyond a wider smoothing, a stand-in for Dubois's crosstalk estimate. Interactive: the same controls; step through the five frequency-domain stages, and try the zone plate at 32 × 32 to see the carriers clearly.

The filters are where the quality lives. Leung, Jeon and Dubois Leung, Jeon & Dubois 2011 (least-squares luma-chroma demultiplexing) design them by least squares on training images rather than by hand, the least-squares luma-chroma demultiplexing that is the version to cite for performance: a small learned element inside an entirely classical structure. The limit is the one Chapter 15.3, Mosaics and aliasing names: where luma and chroma overlap in frequency, as on fine texture near the carriers, no linear filter can separate them, and the leftover luma shows up as false color.

On the upper floors of a tower facade, where the window grid approaches the carriers, choosing the cleaner copy of $C_1$ per pixel removes much of the false color that averaging the two copies leaves, and the method scores above AHD (Figure 15.2.9). The same implementation falls below bilinear interpolation on the flower of Figure 15.2.8, whose saturated edges carry chroma at frequencies the low-pass filters remove.

fig-frequency-demosaic
Figure 15.2.9. Demosaicking as demultiplexing, on a facade whose window grid is fine enough to put luma detail near the carriers; the RGGB mosaic is sampled from the photograph. The mosaic's spectrum shows luma around the origin and the chroma carriers ringed: $C_1$ on the axes at $\pi$, $C_2$ in the corners. Demodulating each carrier and low-pass filtering gives $C_1$ and $C_2$; luma is what remains once they are remodulated and subtracted. Averaging the two copies of $C_1$ (Alleysson and colleagues' version) scores 37.6 dB, the same as AHD; choosing the copy with less luma leakage at each pixel (after Dubois) scores 39.8 dB. The false-color maps show where luma near a carrier was read as color, in vertical streaks when the copies are averaged, and the enlarged box shows it in the windows. Our implementation, with our own low-pass filters rather than the least-squares filters of Leung, Jeon and Dubois.

15.2.4 Where the classical methods stop⧉

Four situations defeat this family, and they are worth listing because the learned methods of Chapter 15.6, Learned demosaicking are judged on exactly them.

Diagonal edges can challenge a horizontal/vertical decision. More directions or a broader model can help; failure is not inevitable on every diagonal.

Repetitive texture near the sampling limit can alias into apparently consistent false structure. Wider context helps some cases, but a single mosaic may leave genuine ambiguity.

Clipped samples no longer report the incident signal linearly. Mark them and avoid treating their color differences as ordinary measurements. The amount of recoverable highlight color depends on which channels remain useful.

Noise perturbs gradients and can switch interpolation decisions. The error depends on signal-to-noise ratio, not ISO alone. Chapter 15.8, Joint demosaicking and denoising incorporates this uncertainty in reconstruction.