# Scaling and Resampling

> Nearest, bilinear and box filtering, the half-pixel bug that shifts your whole image, and why downscaling without a filter turns detail into a solid bar.

Every filter so far kept the grid it was given: one output pixel for one input
pixel, in the same place. Changing the size breaks that. Output pixels no longer
line up with input pixels, and the question becomes what colour sits at a
position that is not a pixel.

```zig
const std = @import("std");
const canvas = @import("_canvas.zig");

var pixels: [canvas.width * canvas.height]u32 = undefined;

/// Coordinates are in 16.16 fixed point: 16 bits of pixel index and 16 bits of
/// position inside the pixel. Integers throughout keeps this exactly
/// reproducible, which matters more here than anywhere else in the section
/// because a resampler that rounds differently on two machines produces
/// different files.
const frac_bits = 16;
const one = 1 << frac_bits;

const Image = struct {
    pixels: []const u32,
    width: usize,
    height: usize,

    fn at(self: Image, x: i64, y: i64) u32 {
        const cx: usize = @intCast(std.math.clamp(x, 0, @as(i64, @intCast(self.width)) - 1));
        const cy: usize = @intCast(std.math.clamp(y, 0, @as(i64, @intCast(self.height)) - 1));
        return self.pixels[cy * self.width + cx];
    }
};

/// Nearest neighbour: round the sample position to the closest pixel centre.
fn nearest(src: Image, fx: i64, fy: i64) u32 {
    return src.at((fx + one / 2) >> frac_bits, (fy + one / 2) >> frac_bits);
}

/// Bilinear: the four surrounding pixels, weighted by how close the sample
/// point is to each.
fn bilinear(src: Image, fx: i64, fy: i64) u32 {
    // Sample positions are relative to pixel *centres*, so shift back by half a
    // pixel before splitting into index and fraction.
    const px = fx - one / 2;
    const py = fy - one / 2;
    const x0 = px >> frac_bits;
    const y0 = py >> frac_bits;
    const tx = px & (one - 1);
    const ty = py & (one - 1);

    const c00 = canvas.channels(src.at(x0, y0));
    const c10 = canvas.channels(src.at(x0 + 1, y0));
    const c01 = canvas.channels(src.at(x0, y0 + 1));
    const c11 = canvas.channels(src.at(x0 + 1, y0 + 1));

    return canvas.rgb(
        mix4(c00.r, c10.r, c01.r, c11.r, tx, ty),
        mix4(c00.g, c10.g, c01.g, c11.g, tx, ty),
        mix4(c00.b, c10.b, c01.b, c11.b, tx, ty),
    );
}

fn mix4(a: u8, b: u8, c: u8, d: u8, tx: i64, ty: i64) u8 {
    const top = @as(i64, a) * (one - tx) + @as(i64, b) * tx;
    const bottom = @as(i64, c) * (one - tx) + @as(i64, d) * tx;
    const value = (top * (one - ty) + bottom * ty) >> (frac_bits * 2);
    return @intCast(std.math.clamp(value, 0, 255));
}

/// Box filter: average every source pixel that falls inside the output pixel.
/// For downscaling this is the one that is actually correct, because it looks
/// at all the input rather than a sample of it.
fn box(src: Image, dst_x: usize, dst_y: usize, dst_w: usize, dst_h: usize) u32 {
    const x_start = (dst_x * src.width) / dst_w;
    const x_end = @max(x_start + 1, ((dst_x + 1) * src.width) / dst_w);
    const y_start = (dst_y * src.height) / dst_h;
    const y_end = @max(y_start + 1, ((dst_y + 1) * src.height) / dst_h);

    var r: u32 = 0;
    var g: u32 = 0;
    var b: u32 = 0;
    var n: u32 = 0;
    for (y_start..y_end) |y| {
        for (x_start..x_end) |x| {
            const c = canvas.channels(src.at(@intCast(x), @intCast(y)));
            r += c.r;
            g += c.g;
            b += c.b;
            n += 1;
        }
    }
    return canvas.rgb(@intCast(r / n), @intCast(g / n), @intCast(b / n));
}

const Sampler = enum { nearest, bilinear, box };

/// Resample `src` into a `dst_w` x `dst_h` buffer.
///
/// The mapping is the part that goes wrong. An output pixel covers a range of
/// the input, and the point to sample is the *centre* of that range:
///
///     src = (dst + 0.5) * scale - 0.5
///
/// Dropping either half-pixel gives `src = dst * scale`, which lines the two
/// images up by their top-left corners instead of their centres and shifts the
/// result by half a source pixel. That is the single most common resampling
/// bug, and on a 2x upscale it is a visible offset.
fn resample(
    dst: []u32,
    dst_w: usize,
    dst_h: usize,
    src: Image,
    sampler: Sampler,
    correct_centres: bool,
) void {
    const x_scale = @divTrunc(@as(i64, @intCast(src.width)) << frac_bits, @as(i64, @intCast(dst_w)));
    const y_scale = @divTrunc(@as(i64, @intCast(src.height)) << frac_bits, @as(i64, @intCast(dst_h)));

    for (0..dst_h) |y| {
        for (0..dst_w) |x| {
            if (sampler == .box) {
                dst[y * dst_w + x] = box(src, x, y, dst_w, dst_h);
                continue;
            }
            const fx = if (correct_centres)
                (((@as(i64, @intCast(x)) << frac_bits) + one / 2) * x_scale >> frac_bits) - one / 2
            else
                @as(i64, @intCast(x)) * x_scale;
            const fy = if (correct_centres)
                (((@as(i64, @intCast(y)) << frac_bits) + one / 2) * y_scale >> frac_bits) - one / 2
            else
                @as(i64, @intCast(y)) * y_scale;

            dst[y * dst_w + x] = switch (sampler) {
                .nearest => nearest(src, fx, fy),
                .bilinear => bilinear(src, fx, fy),
                .box => unreachable,
            };
        }
    }
}

/// Print a small buffer as luma values, since at these sizes the ASCII dump has
/// nothing to show.
fn printGrid(out: *std.Io.Writer, buf: []const u32, w: usize, h: usize) !void {
    for (0..h) |y| {
        try out.writeAll("   ");
        for (0..w) |x| try out.print("{d:>4}", .{canvas.luma(buf[y * w + x])});
        try out.writeByte('\n');
    }
}

pub fn main(init: std.process.Init) !void {
    var buf: [32768]u8 = undefined;
    var file_writer = std.Io.File.stdout().writerStreaming(init.io, &buf);
    const out = &file_writer.interface;

    canvas.scene(&pixels);
    const src = Image{ .pixels = &pixels, .width = canvas.width, .height = canvas.height };

    // A 4x1 strip upscaled to 8x1, so every number can be checked by hand.
    // Source luma: two dark, two bright.
    var strip = [_]u32{
        canvas.rgb(0, 0, 0),
        canvas.rgb(0, 0, 0),
        canvas.rgb(255, 255, 255),
        canvas.rgb(255, 255, 255),
    };
    const strip_src = Image{ .pixels = &strip, .width = 4, .height = 1 };
    var up: [8]u32 = undefined;

    try out.writeAll("4 pixels -> 8, source luma 0 0 255 255\n");
    resample(&up, 8, 1, strip_src, .nearest, true);
    try out.writeAll("  nearest, centres correct:\n");
    try printGrid(out, &up, 8, 1);
    resample(&up, 8, 1, strip_src, .bilinear, true);
    try out.writeAll("  bilinear, centres correct:\n");
    try printGrid(out, &up, 8, 1);
    resample(&up, 8, 1, strip_src, .bilinear, false);
    try out.writeAll("  bilinear, half-pixel bug (note the shift):\n");
    try printGrid(out, &up, 8, 1);

    // Downscale the scene to a quarter of its size, three ways.
    const small_w = canvas.width / 4;
    const small_h = canvas.height / 4;
    var small: [(canvas.width / 4) * (canvas.height / 4)]u32 = undefined;

    resample(&small, small_w, small_h, src, .nearest, true);
    try out.writeAll("\n64x32 -> 16x8, nearest:\n");
    try printGrid(out, &small, small_w, small_h);

    resample(&small, small_w, small_h, src, .box, true);
    try out.writeAll("\n64x32 -> 16x8, box filtered:\n");
    try printGrid(out, &small, small_w, small_h);

    // Aliasing, with nowhere to hide. A row of alternating black and white
    // columns has a real average of 127 everywhere. Nearest neighbour samples
    // every fourth pixel, and every fourth pixel of an alternating pattern is
    // the *same* colour, so the answer is a solid bar of whichever phase the
    // grid happened to land on. The detail is not merely lost; it has become a
    // large, confident, wrong feature.
    var stripes: [64]u32 = undefined;
    for (&stripes, 0..) |*p, i| {
        p.* = if (i % 2 == 0) canvas.rgb(0, 0, 0) else canvas.rgb(255, 255, 255);
    }
    const stripe_src = Image{ .pixels = &stripes, .width = 64, .height = 1 };
    var thin: [16]u32 = undefined;

    try out.writeAll("\nalternating black/white columns, 64 -> 16 (true average is 127)\n");
    resample(&thin, 16, 1, stripe_src, .nearest, true);
    try out.writeAll("  nearest:\n");
    try printGrid(out, &thin, 16, 1);
    resample(&thin, 16, 1, stripe_src, .box, true);
    try out.writeAll("  box filtered:\n");
    try printGrid(out, &thin, 16, 1);

    try out.flush();
}
```

*Runnable: compiled to WebAssembly and executed by CI against Zig master. (`09-graphics.resample`)*

## The mapping, and the bug in it

An output pixel covers a range of the input. The point to sample is the centre
of that range:

```
src = (dst + 0.5) * scale - 0.5
```

Both halves matter. `dst + 0.5` moves from the output pixel's index to its
centre; `- 0.5` moves from the input's centre back to its index. Drop them and
you get `src = dst * scale`, which aligns the two images by their top-left
corners rather than their centres and shifts the result by half a source pixel.

The snippet upscales four pixels to eight, with and without the correction:

```
4 pixels -> 8, source luma 0 0 255 255
  bilinear, centres correct:
      0   0   0   0  63 191 255 255
  bilinear, half-pixel bug (note the shift):
      0   0   0   0 127 255 255 255
```

The correct version puts a symmetric two-sample ramp across the edge. The buggy
one puts a one-sample ramp starting later, and the whole image sits half a
source pixel to the left of where it should. On a 2x upscale that is a visible
offset, and it is the single most common resampling bug: easy to write, hard to
see in a single image, obvious the moment you blend the result against the
original.

Fixed point rather than floats here, because a resampler that rounds differently
on two machines produces different files, and a difference of one in the last
bit is enough to fail a checksum.

## Nearest and bilinear

**Nearest** rounds the sample position to the closest pixel centre. One read per
output pixel, no arithmetic on colours at all, and it preserves exact values,
which is why it is the right choice for pixel art, indexed colour, and any image
whose values are labels rather than measurements. A segmentation mask
interpolated bilinearly produces class 7.5, which is not a class.

**Bilinear** reads the four surrounding pixels and weights them by how close the
sample point is to each. It is a genuine interpolation, and for upscaling
photographs it is the cheapest thing that does not look blocky. Bicubic reads 16
and gives a slightly sharper result; Lanczos reads more still. The pattern
continues, with each step buying less.

## Downscaling is a different problem

Upscaling asks for values between the ones you have. Downscaling asks you to
throw values away, and the way you throw them away decides whether the result is
a smaller version of the picture or a different picture.

The snippet makes this unambiguous with a row of alternating black and white
columns downscaled to a quarter width. The true average of that row is 127
everywhere:

```
alternating black/white columns, 64 -> 16 (true average is 127)
  nearest:
      0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  box filtered:
    127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127
```

Nearest neighbour samples every fourth pixel, and every fourth pixel of an
alternating pattern is the same colour. The answer is a solid black bar. Not
approximately right, not noisy: confidently, uniformly wrong, and it would have
been a solid white bar had the grid landed one pixel over.

That is aliasing. The input contained detail finer than the output grid can
represent, and point-sampling turned it into a large false feature rather than
losing it gracefully. The same effect is why a striped shirt strobes on video
and why a brick wall in the distance turns into moiré.

The fix is to remove the detail before discarding the samples: average every
input pixel that falls inside the output pixel, which is the box filter above,
or blur first and then point-sample. Either way, the filter has to look at all
the input. Any downscaler that reads one pixel per output pixel is guessing.

## Where this connects

[Antialiasing](https://www.ziglang.in/learn/graphics/antialiasing/) took sixteen samples inside one pixel to
turn a hard inside/outside test into a coverage fraction. This is the same idea
arriving from the other direction: a pixel is not a point, it is an area, and
treating it as a point is what produces both the staircase there and the solid
bar here.

Two practical notes:

**Filter in linear light.** Averaging sRGB bytes makes a downscaled image
darker than it should be, and the error is largest exactly where contrast is
highest. This is the same caveat [alpha blending](https://www.ziglang.in/learn/graphics/blending/) sets out,
and downscaling is where it is most visible, because a checkerboard of black and
white averages to 127 in sRGB and to about 187 done correctly.

**Use premultiplied alpha.** Interpolating colours that carry their own alpha
without premultiplying pulls the colour of fully transparent pixels into the
result, which puts a halo of whatever the transparent regions happened to
contain around every edge.