Lesson 25: Image Formation¶

Up to now we have been processing images without thinking about where they come from. This lesson unpacks that mystery by considering: geometry (how the projection matrix is derived from a pinhole camera model), optics (what a real lens does that a pinhole camera fails to model — focus and blur), and sensing/color (how a real sensor captures only one color per pixel and needs to reconstruct the rest, along with the two classic pitfalls — gamma and white balance — that trip up naive image processing).

In [1]:

import numpy as np
import cv2
import matplotlib.pyplot as plt

The pinhole projection matrix, derived¶

A pinhole camera lets through only the single ray of light passing through one point (the center of projection). By similar triangles, a 3D point $(X, Y, Z)$ (camera-centered coordinates) lands on the image plane at

$$x = f\frac{X}{Z}, \qquad y = f\frac{Y}{Z}$$

where $f$ is the distance from the pinhole to the image plane. Converting to pixel units (accounting for pixel size and where the origin is) folds in the intrinsics $K$; accounting for the camera's position and orientation relative to the world folds in the extrinsics $[R | t]$. Put together, in homogeneous coordinates:

$$\underbrace{\begin{bmatrix}u\\v\\1\end{bmatrix}}_{\text{pixel}} \;\propto\; \underbrace{\begin{bmatrix}f_x&0&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}}_{K}\underbrace{\begin{bmatrix}R & t\end{bmatrix}}_{\text{extrinsics}}\begin{bmatrix}X\\Y\\Z\\1\end{bmatrix}_{\text{world}}$$

This is exactly the $P = K [R | t]$ that Lessons 26-28 (Epipolar Geometry, Camera Calibration, and Structure from Motion) will use freely. Every symbol can be traced back to a physical cause: $f_x, f_y$ to focal length and pixel size, $(c_x, c_y)$ to where the optical axis hits the sensor, $R, t$ to the camera's pose.

In [2]:

K = np.array([[500, 0, 320], [0, 500, 240], [0, 0, 1]], dtype=np.float64)
R, t = np.eye(3), np.array([0, 0, 5.0])
P = K @ np.hstack([R, t.reshape(3, 1)])

cube_corners = np.array([[x, y, z] for x in (-1, 1) for y in (-1, 1) for z in (-1, 1)], dtype=np.float64)
homogeneous = np.hstack([cube_corners, np.ones((8, 1))])
projected = (P @ homogeneous.T).T
pixels = projected[:, :2] / projected[:, 2:3]

edges = [(0, 1), (0, 2), (0, 4), (1, 3), (1, 5), (2, 3), (2, 6),
         (3, 7), (4, 5), (4, 6), (5, 7), (6, 7)]

plt.figure(figsize=(5, 4))
for i, j in edges:
    plt.plot([pixels[i, 0], pixels[j, 0]], [pixels[i, 1], pixels[j, 1]], color='tab:blue')
plt.scatter(*pixels.T, color='red', zorder=5)
plt.gca().invert_yaxis()
plt.title('A 3D cube, projected through P = K[R|t]')
plt.show()

No description has been provided for this image

Real lenses: focus and defocus blur¶

A pinhole is an idealization — it lets through so little light that a real pinhole camera needs absurdly long exposures. Real cameras use a lens with a wide aperture instead, which gathers much more light but focuses sharply only for objects at one particular distance (the focal plane). Points nearer or farther than that spread their light over a small disk on the sensor (the circle of confusion) instead of a point — exactly the blurring convolution kernels from Lessons 9-10, just applied with a kernel size that depends on depth.

In [3]:

scene = np.zeros((200, 300, 3), dtype=np.uint8)
cv2.rectangle(scene, (20, 20), (90, 180), (255, 120, 30), -1)    # near object
cv2.rectangle(scene, (110, 20), (190, 180), (30, 200, 255), -1)  # object at the focal plane
cv2.rectangle(scene, (210, 20), (280, 180), (120, 255, 60), -1)  # far object

depths_m = {'near': 2.0, 'mid': 5.0, 'far': 9.0}
focus_depth_m = 5.0

def defocus_kernel_size(depth, focus, strength=15.0):
    radius = strength * abs(depth - focus) / focus
    ksize = int(radius) * 2 + 1
    return max(1, ksize)

def blur_with_context(image, y0, y1, x0, x1, ksize, pad):
    """Blur a region using a padded crop, so the blur can pull in real neighboring
    pixels (including the background) instead of only the region's own solid color."""
    if ksize <= 1:
        return image[y0:y1, x0:x1].copy()
    yy0, yy1 = max(0, y0 - pad), min(image.shape[0], y1 + pad)
    xx0, xx1 = max(0, x0 - pad), min(image.shape[1], x1 + pad)
    blurred = cv2.GaussianBlur(image[yy0:yy1, xx0:xx1], (ksize, ksize), 0)
    return blurred[y0 - yy0:y0 - yy0 + (y1 - y0), x0 - xx0:x0 - xx0 + (x1 - x0)]

defocused = scene.copy()
regions = {'near': (20, 180, 20, 90), 'mid': (20, 180, 110, 190), 'far': (20, 180, 210, 280)}
for name, (y0, y1, x0, x1) in regions.items():
    k = defocus_kernel_size(depths_m[name], focus_depth_m)
    print(f'{name:>5} object: depth={depths_m[name]}m, blur kernel size={k}')
    defocused[y0:y1, x0:x1] = blur_with_context(scene, y0, y1, x0, x1, k, pad=k)

fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
axes[0].imshow(cv2.cvtColor(scene, cv2.COLOR_BGR2RGB))
axes[0].set_title('All in focus (pinhole idealization)')
axes[1].imshow(cv2.cvtColor(defocused, cv2.COLOR_BGR2RGB))
axes[1].set_title(f'Focused at {focus_depth_m}m\n(near/far objects blurred)')
for ax in axes:
    ax.axis('off')
plt.tight_layout()
plt.show()

 near object: depth=2.0m, blur kernel size=19
  mid object: depth=5.0m, blur kernel size=1
  far object: depth=9.0m, blur kernel size=25

This is why a wider aperture (more light, but a shallower depth of field) trades off against a narrower one (less light, but more of the scene in focus) — the classic photographic aperture/depth-of-field tradeoff, and the reason portrait photos often blur the background while keeping the subject sharp.

The Bayer color filter array: one color per pixel¶

Most camera sensors don't measure red, green, and blue at every pixel. A single monochrome photosensor array sits under a mosaic of tiny color filters — the Bayer filter (Bayer, 1976) — so each individual pixel physically records only one of the three colors. The most common arrangement (RGGB) repeats a 2x2 tile: one red, one blue, and two green filters, since the human eye is most sensitive to green and luminance detail. Every pixel in the raw sensor output is missing two-thirds of its color information by construction.

In [4]:

rng = np.random.default_rng(0)
bayer_scene = np.zeros((120, 160, 3), dtype=np.uint8)
cv2.rectangle(bayer_scene, (10, 10), (70, 100), (40, 180, 220), -1)
cv2.circle(bayer_scene, (110, 60), 40, (200, 90, 40), -1)
bayer_scene = np.clip(bayer_scene.astype(np.float64) + rng.normal(0, 4, bayer_scene.shape), 0, 255).astype(np.uint8)

def make_bayer_rggb(img):
    # (row, col) parity: (even,even)=R, (even,odd)=G, (odd,even)=G, (odd,odd)=B
    H, W = img.shape[:2]
    mosaic = np.zeros((H, W), dtype=np.uint8)
    mosaic[0::2, 0::2] = img[0::2, 0::2, 2]  # R  (img is BGR, so channel 2 = R)
    mosaic[0::2, 1::2] = img[0::2, 1::2, 1]  # G
    mosaic[1::2, 0::2] = img[1::2, 0::2, 1]  # G
    mosaic[1::2, 1::2] = img[1::2, 1::2, 0]  # B
    return mosaic

def colorize_bayer_rggb(mosaic):
    """Put each pixel's raw reading into its OWN filter's color channel, zeroing the
    other two -- this is what actually lets you see the RGGB pattern, since a plain
    grayscale view of the mosaic just shows intensities with no color information at all."""
    H, W = mosaic.shape
    colored = np.zeros((H, W, 3), dtype=np.uint8)
    colored[0::2, 0::2, 0] = mosaic[0::2, 0::2]  # R filter
    colored[0::2, 1::2, 1] = mosaic[0::2, 1::2]  # G filter
    colored[1::2, 0::2, 1] = mosaic[1::2, 0::2]  # G filter
    colored[1::2, 1::2, 2] = mosaic[1::2, 1::2]  # B filter
    return colored

bayer = make_bayer_rggb(bayer_scene)
bayer_colored = colorize_bayer_rggb(bayer)

y0, y1, x0, x1 = 40, 56, 20, 36  # a 16x16 crop, fully inside a colored shape

fig, axes = plt.subplots(1, 3, figsize=(11, 3.5))
axes[0].imshow(cv2.cvtColor(bayer_scene, cv2.COLOR_BGR2RGB))
axes[0].add_patch(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, edgecolor='white', facecolor='none', linewidth=1.5))
axes[0].set_title('true scene (each pixel: R, G, and B)')
axes[1].imshow(bayer[y0:y1, x0:x1], cmap='gray')
axes[1].set_title('raw mosaic\n(grayscale readout -- no color info visible)')
axes[2].imshow(bayer_colored[y0:y1, x0:x1])
axes[2].set_title('same patch, colorized by filter\n(the actual RGGB pattern)')
for ax in axes:
    ax.axis('off')
plt.tight_layout()
plt.show()

Demosaicking: reconstructing RGB from the mosaic¶

Demosaicking (or "debayering") fills in each pixel's two missing colors by interpolating from same-color neighbors — the simplest version is bilinear: average the nearby known samples of a color, weighted by distance, exactly Lesson 9's bilinear interpolation, just interpolating across a sparse, patterned grid of known samples instead of a dense downsampled image.

In [5]:

def demosaick_bilinear_rggb(mosaic):
    H, W = mosaic.shape
    R = np.zeros((H, W)); G = np.zeros((H, W)); B = np.zeros((H, W))
    R[0::2, 0::2] = mosaic[0::2, 0::2]
    G[0::2, 1::2] = mosaic[0::2, 1::2]; G[1::2, 0::2] = mosaic[1::2, 0::2]
    B[1::2, 1::2] = mosaic[1::2, 1::2]

    Rmask = np.zeros((H, W)); Rmask[0::2, 0::2] = 1
    Gmask = np.zeros((H, W)); Gmask[0::2, 1::2] = 1; Gmask[1::2, 0::2] = 1
    Bmask = np.zeros((H, W)); Bmask[1::2, 1::2] = 1

    def fill(channel, mask):
        # weighted average of known same-color neighbors within a 3x3 window
        kernel = np.array([[0.25, 0.5, 0.25], [0.5, 1.0, 0.5], [0.25, 0.5, 0.25]])
        weight_sum = cv2.filter2D(mask, -1, kernel, borderType=cv2.BORDER_REFLECT)
        value_sum = cv2.filter2D(channel * mask, -1, kernel, borderType=cv2.BORDER_REFLECT)
        filled = channel.copy()
        empty = mask == 0
        filled[empty] = value_sum[empty] / np.clip(weight_sum[empty], 1e-6, None)
        return filled

    R_full, G_full, B_full = fill(R, Rmask), fill(G, Gmask), fill(B, Bmask)
    return np.clip(np.stack([B_full, G_full, R_full], axis=-1), 0, 255).astype(np.uint8)

manual_demosaick = demosaick_bilinear_rggb(bayer)
cv2_demosaick = cv2.cvtColor(bayer, cv2.COLOR_BayerBG2BGR)

manual_err = np.abs(manual_demosaick.astype(np.float64) - bayer_scene.astype(np.float64)).mean()
cv2_err = np.abs(cv2_demosaick.astype(np.float64) - bayer_scene.astype(np.float64)).mean()
print(f'manual bilinear demosaick MAE vs ground truth: {manual_err:.2f}')
print(f'cv2 demosaick MAE vs ground truth:              {cv2_err:.2f}')

fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))
for ax, im, title in zip(axes, [bayer_scene, manual_demosaick, cv2_demosaick],
                          ['ground truth', 'manual bilinear\ndemosaick', 'cv2 demosaick']):
    ax.imshow(cv2.cvtColor(im, cv2.COLOR_BGR2RGB))
    ax.set_title(title, fontsize=9)
    ax.axis('off')
plt.tight_layout()
plt.show()

manual bilinear demosaick MAE vs ground truth: 3.37
cv2 demosaick MAE vs ground truth:              3.43

Both methods reconstruct the ground truth closely (the residual error is mostly just the sensor noise already present in the mosaic, not a demosaicking artifact) — a from-scratch bilinear fill and OpenCV's built-in demosaicking land within a fraction of a pixel value of each other. Real demosaicking algorithms are more sophisticated than pure bilinear (edge-aware interpolation that avoids blurring across object boundaries, since naive bilinear demosaicking is exactly what causes the color-fringing "zipper" artifacts visible along sharp edges in cheap cameras), but the core idea — interpolate each missing color from its nearest same-color neighbors — is the same.

The pitfall: getting the mosaic pattern wrong¶

The R/G/B arrangement matters, and there's more than one convention (RGGB, BGGR, GRBG, GBRG, depending on the sensor). Demosaicking with the wrong assumed pattern doesn't fail gracefully — it silently produces a plausible-looking but badly wrong-colored image, since every pixel gets a value, just the wrong one.

In [6]:

wrong_demosaick = cv2.cvtColor(bayer, cv2.COLOR_BayerRG2BGR)  # wrong pattern assumption
wrong_err = np.abs(wrong_demosaick.astype(np.float64) - bayer_scene.astype(np.float64)).mean()
print(f'correct pattern (BayerBG2BGR) MAE: {cv2_err:.2f}')
print(f'WRONG pattern (BayerRG2BGR) MAE:   {wrong_err:.2f}')

fig, axes = plt.subplots(1, 2, figsize=(6, 3.5))
axes[0].imshow(cv2.cvtColor(cv2_demosaick, cv2.COLOR_BGR2RGB)); axes[0].set_title('correct pattern', fontsize=9)
axes[1].imshow(cv2.cvtColor(wrong_demosaick, cv2.COLOR_BGR2RGB)); axes[1].set_title('wrong pattern', fontsize=9)
for ax in axes:
    ax.axis('off')
plt.tight_layout()
plt.show()

correct pattern (BayerBG2BGR) MAE: 3.43
WRONG pattern (BayerRG2BGR) MAE:   64.69

The error jumps by more than an order of magnitude with the wrong pattern — every red sample gets treated as blue and vice versa, producing a strongly color-shifted (not just slightly-off) image. This is the practical reason RAW image pipelines always need to know a camera's exact sensor layout: demosaicking is one of the few image-processing steps where a metadata mistake corrupts every single pixel at once, and the result still looks like a photograph, just the wrong one.

Gamma correction: pixel values are not linear light¶

Cameras and displays don't store brightness proportionally to physical light intensity (linear radiance). Instead, values are gamma-encoded: encoded = linear ** (1/gamma) (typically gamma ≈ 2.2, close to the sRGB standard), which allocates more encoded levels to darker tones — matching human vision's greater sensitivity to shadows than highlights, and historically matching how CRT displays responded to voltage. Decoding reverses it: linear = encoded ** gamma.

In [7]:

def encode_gamma(linear, gamma=1 / 2.2):
    return np.clip(linear, 0, 1) ** gamma

def decode_gamma(encoded, gamma=2.2):
    return np.clip(encoded, 0, 1) ** gamma

linear_values = np.linspace(0, 1, 100)
plt.plot(linear_values, linear_values, '--', color='gray', label='identity (no gamma)')
plt.plot(linear_values, encode_gamma(linear_values), label='gamma-encoded (stored pixel value)')
plt.xlabel('linear light intensity')
plt.ylabel('encoded value')
plt.legend(fontsize=8)
plt.title('Gamma encoding curve')
plt.show()

The pitfall: averaging encoded pixels is not averaging light¶

Blending, resizing, and blurring (Lessons 8-10) all implicitly assume a linear quantity is being averaged. Pixel values usually aren't linear — they're gamma-encoded — so naively averaging two encoded pixel values gives the wrong physical brightness, sometimes badly so.

In [8]:

black_encoded, white_encoded = encode_gamma(0.0), encode_gamma(1.0)

# WRONG: average the encoded (stored) pixel values directly
naive_blend_encoded = (black_encoded + white_encoded) / 2
naive_blend_as_linear_light = decode_gamma(naive_blend_encoded)

# RIGHT: decode to linear light first, average there, then re-encode for storage/display
correct_blend_linear = (decode_gamma(black_encoded) + decode_gamma(white_encoded)) / 2
correct_blend_encoded = encode_gamma(correct_blend_linear)

print(f'true 50% linear-light gray, correctly encoded: {correct_blend_encoded:.3f}')
print(f'naive average of encoded black/white:           {naive_blend_encoded:.3f}')
print()
print(f'naive result represents only {naive_blend_as_linear_light:.1%} of the true linear brightness '
      f'(should be 50%)')

swatch = np.zeros((80, 240), dtype=np.float64)
swatch[:, :80] = naive_blend_encoded
swatch[:, 80:160] = correct_blend_encoded
swatch[:, 160:] = 0.5
plt.imshow(swatch, cmap='gray', vmin=0, vmax=1)
plt.title('naive blend  |  correct blend  |  encoded 0.5 for reference')
plt.axis('off')
plt.show()

true 50% linear-light gray, correctly encoded: 0.730
naive average of encoded black/white:           0.500

naive result represents only 21.8% of the true linear brightness (should be 50%)

The naive blend is visibly, dramatically darker than the physically correct 50%-gray result — a real and common bug when image-processing code (blurring, mipmap generation, alpha blending) operates directly on gamma-encoded pixel values instead of linearizing first.

White balance: color depends on the light source¶

A camera records the product of a surface's reflectance and the illuminant's color, not the surface's "true" color alone — a white sheet of paper looks orange under warm tungsten light and blue-ish under overcast sky, even though our visual system usually (mostly) compensates for this automatically (color constancy). Cameras have to do the same correction deliberately: white balancing.

In [9]:

rng = np.random.default_rng(0)
true_scene = np.zeros((150, 150, 3), dtype=np.uint8)
cv2.rectangle(true_scene, (20, 20), (130, 130), (150, 150, 150), -1)  # a neutral gray card
cv2.circle(true_scene, (75, 75), 30, (60, 120, 200), -1)
true_scene = np.clip(true_scene.astype(np.float64) + rng.normal(0, 5, true_scene.shape), 0, 255).astype(np.uint8)

tungsten_gains = np.array([1.3, 1.0, 0.6])  # BGR: boosts red, cuts blue -- a warm cast
cast_scene = np.clip(true_scene.astype(np.float64) * tungsten_gains, 0, 255).astype(np.uint8)

fig, axes = plt.subplots(1, 2, figsize=(6, 3.5))
axes[0].imshow(cv2.cvtColor(true_scene, cv2.COLOR_BGR2RGB))
axes[0].set_title('True colors')
axes[1].imshow(cv2.cvtColor(cast_scene, cv2.COLOR_BGR2RGB))
axes[1].set_title('Under warm tungsten light')
for ax in axes:
    ax.axis('off')
plt.tight_layout()
plt.show()

The gray-world assumption¶

A simple, classic white-balancing algorithm assumes that, averaged over an entire real-world scene, colors roughly cancel out to neutral gray. Under that assumption, any systematic difference between a photo's per-channel averages reveals the illuminant's color cast — so rescaling each channel to equalize the averages should approximately undo it.

In [10]:

def gray_world_white_balance(image):
    image_f = image.astype(np.float64)
    channel_means = image_f.reshape(-1, 3).mean(axis=0)
    target_gray = channel_means.mean()
    gains = target_gray / channel_means
    return np.clip(image_f * gains, 0, 255).astype(np.uint8)

corrected_scene = gray_world_white_balance(cast_scene)

error_before = np.abs(cast_scene.astype(int) - true_scene.astype(int)).mean()
error_after = np.abs(corrected_scene.astype(int) - true_scene.astype(int)).mean()
print(f'mean abs pixel error before correction: {error_before:.2f}')
print(f'mean abs pixel error after correction:  {error_after:.2f}')

fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))
for ax, im, title in zip(axes, [true_scene, cast_scene, corrected_scene],
                          ['True colors', 'Color cast', 'Gray-world corrected']):
    ax.imshow(cv2.cvtColor(im, cv2.COLOR_BGR2RGB))
    ax.set_title(title, fontsize=9)
    ax.axis('off')
plt.tight_layout()
plt.show()

mean abs pixel error before correction: 18.98
mean abs pixel error after correction:  7.43

The correction substantially reduces the error, but not to zero — the gray-world assumption is only approximate here (this small scene isn't perfectly neutral on average, since it's dominated by the colored circle), which is exactly why real cameras often combine it with other cues (an actual detected gray/white reference patch, learned scene statistics, or user-specified presets) rather than relying on the gray-world assumption alone.

Exercise¶

Change strength in defocus_kernel_size to make the depth-of-field effect much shallower (larger strength) or much deeper (smaller strength). At strength=0, what should happen, and does the code produce that?
Modify make_bayer_rggb and the Rmask/Gmask/Bmask assignments in demosaick_bilinear_rggb to build and decode a BGGR mosaic instead of RGGB (swap which corner is red vs. blue). Find the matching OpenCV code by trying each of cv2.COLOR_BayerBG2BGR, cv2.COLOR_BayerGB2BGR, cv2.COLOR_BayerGR2BGR, cv2.COLOR_BayerRG2BGR and checking which one drops the MAE back down to the noise floor — OpenCV's naming convention for these codes is notoriously easy to get backwards, which is itself a small demonstration of this lesson's point about silent pattern-mismatch errors.
Repeat the gamma-blending experiment, but blend 0.2 and 0.8 (instead of pure black and white) in both the naive and correct ways. Is the discrepancy between them larger or smaller than the black/white case, and why might that make sense given the shape of the gamma curve?
Modify the white-balance demo so the scene is dominated by a large saturated red region instead of a neutral gray card (i.e. make gray-world's core assumption clearly false). Does gray_world_white_balance still improve the result, make no difference, or actively make it worse?

← Lesson 24: Image Stitching and Mosaicking

Lesson 26: Epipolar Geometry →