Image resampling filters

Deep dive

How it works

In: a source image (rows of f32 or u8 pixels), the output size, a kernel (box, triangle, Catmull-Rom, Lanczos3) · Out: the resized image · State: per-axis weight tables, built once per size pair · Cost: two 1-D passes; order of 1 ns per source pixel at 4x shrink (measured below)

Output pixel o of out from in pixels, scale s = in/out:

1. Pixel centres sit at +0.5: center = (o + 0.5) * s.
2. Kernel width: fscale = max(s, 1), support = kernel.support * fscale. Only shrinking widens the kernel.
3. Taps are source pixels i in [center - support + 0.5, center + support + 0.5), clipped to the image.
4. Weight w_i = kernel((i + 0.5 - center) / fscale); divide by the sum so weights add to 1 (this also fixes the borders, by renormalising over the clipped taps).
5. Horizontal pass into a temporary out_w x in_h, then vertical pass. Order is mathematically irrelevant; do first the axis that shrinks more.

Kernels (x in source pixels):

box       1            |x| <= 0.5
triangle  1 - |x|      |x| < 1
cubic     (a+2)|x|^3 - (a+3)|x|^2 + 1                 |x| < 1
          a|x|^3 - 5a|x|^2 + 8a|x| - 4a               1 <= |x| < 2    (a = -0.5: Catmull-Rom/Keys)
lanczos3  sinc(x) * sinc(x/3)   |x| < 3,  sinc(x) = sin(pi x)/(pi x)

Mitchell-Netravali is the same cubic family with B = C = 1/3: blurrier, almost no ringing.

Code

use std::f64::consts::PI;

#[derive(Clone, Copy)]
pub enum Filter { Box, Triangle, CatmullRom, Lanczos3 }

impl Filter {
    pub fn support(self) -> f64 {
        match self { Filter::Box => 0.5, Filter::Triangle => 1.0, Filter::CatmullRom => 2.0, Filter::Lanczos3 => 3.0 }
    }
    pub fn weight(self, x: f64) -> f64 {
        let x = x.abs();
        match self {
            Filter::Box => if x <= 0.5 { 1.0 } else { 0.0 },
            Filter::Triangle => (1.0 - x).max(0.0),
            Filter::CatmullRom => { // Keys cubic, a = -0.5
                const A: f64 = -0.5;
                if x < 1.0 { ((A + 2.0) * x - (A + 3.0)) * x * x + 1.0 }
                else if x < 2.0 { (((x - 5.0) * x + 8.0) * x - 4.0) * A }
                else { 0.0 }
            }
            Filter::Lanczos3 => if x < 3.0 { sinc(x) * sinc(x / 3.0) } else { 0.0 },
        }
    }
}

fn sinc(x: f64) -> f64 {
    if x == 0.0 { 1.0 } else { let p = PI * x; p.sin() / p }
}

/// For each output index: first source index and normalized weights (flat, fixed stride).
pub struct Coeffs { pub start: Vec<usize>, pub len: Vec<usize>, pub w: Vec<f32>, pub stride: usize }

pub fn coeffs(in_size: usize, out_size: usize, f: Filter) -> Coeffs {
    let scale = in_size as f64 / out_size as f64;
    let fscale = scale.max(1.0);              // widen the kernel only when shrinking
    let support = f.support() * fscale;
    let stride = support.ceil() as usize * 2 + 1;
    let mut c = Coeffs { start: vec![], len: vec![], w: vec![0.0; out_size * stride], stride };
    for o in 0..out_size {
        let center = (o as f64 + 0.5) * scale; // pixel centers at +0.5
        let lo = ((center - support + 0.5).max(0.0)) as usize;
        let hi = (((center + support + 0.5) as usize)).min(in_size);
        let row = &mut c.w[o * stride..(o + 1) * stride];
        let mut sum = 0.0;
        for (k, i) in (lo..hi).enumerate() {
            let w = f.weight((i as f64 + 0.5 - center) / fscale);
            row[k] = w as f32;
            sum += w;
        }
        for k in 0..hi - lo { row[k] = (row[k] as f64 / sum) as f32; } // partition of unity, also at borders
        c.start.push(lo);
        c.len.push(hi - lo);
    }
    c
}

/// One channel, row-major f32 (linear light, premultiplied if it has alpha).
pub fn resize(src: &[f32], w: usize, h: usize, ow: usize, oh: usize, f: Filter) -> Vec<f32> {
    let (cx, cy) = (coeffs(w, ow, f), coeffs(h, oh, f));
    let mut tmp = vec![0.0f32; ow * h];       // horizontal pass: dot products along rows
    for y in 0..h {
        let row = &src[y * w..(y + 1) * w];
        for x in 0..ow {
            let k = &cx.w[x * cx.stride..x * cx.stride + cx.len[x]];
            let s = &row[cx.start[x]..cx.start[x] + cx.len[x]];
            tmp[y * ow + x] = s.iter().zip(k).map(|(a, b)| a * b).sum();
        }
    }
    let mut out = vec![0.0f32; ow * oh];      // vertical pass: weighted sum of whole rows (cache friendly)
    for y in 0..oh {
        let orow = &mut out[y * ow..(y + 1) * ow];
        for (j, &wt) in cy.w[y * cy.stride..y * cy.stride + cy.len[y]].iter().enumerate() {
            let r = &tmp[(cy.start[y] + j) * ow..(cy.start[y] + j + 1) * ow];
            for (o, &v) in orow.iter_mut().zip(r) { *o += wt * v; }
        }
    }
    out
}

Tested: constants survive any resize (partition of unity), identity size reproduces the input, triangle upscaling reproduces a linear ramp, box 2x downscale equals 2x2 block averages, Catmull-Rom rings on a step.

Complexity and constants

• Weights: O(out x taps) kernel evaluations (sin for Lanczos), once per axis and size pair, not per pixel or channel.
• Pixels: horizontal pass out_w * in_h * taps_x, vertical out_w * out_h * taps_y multiply-adds, where taps = 2*ceil(support*fscale)+1 slots (about 2*support*fscale non-zero).
• Memory traffic: the source is read once, the temporary (out_w * in_h * 4 bytes for f32) written and read once. A 2-D kernel would instead re-read overlapping source windows for every output pixel.
• Measured here (one channel f32, one thread, plain scalar Rust, no SIMD, this Mac): 4096x4096 to 1024x1024 took order of 1 ns per source pixel for triangle, Catmull-Rom and Lanczos3 (11-18 ms); 1024 to 2048 took order of 1-1.5 ns per output pixel.

Efficiency

• Cost: measured here (one channel f32, one thread, plain scalar Rust, no SIMD, this Mac): 4096x4096 to 1024x1024 took order of 1 ns per source pixel for triangle, Catmull-Rom and Lanczos3 (11-18 ms); 1024 to 2048 took order of 1-1.5 ns per output pixel. Bound: memory bandwidth for large images, multiply-adds for wide kernels.
• Levers:
  • separable passes — Lanczos3 upscaling needs at most 6 + 6 taps per output pixel instead of 36; shrinking by s costs s instead of s^2 — always
  • precomputed weights — no sin or branches in the inner loop; the vertical pass is a streaming out += w * row that auto-vectorises — any size
  • pre-shrink by an integer box average, or decode a JPEG at 1/2, 1/4 or 1/8 scale — fewer pixels reach the kernel (Pillow reducing_gap: larger is closer to exact, smaller is faster) — shrink factors above about 2
  • fixed-point weights — Pillow's 8-bit path uses 22 fractional bits and i32 accumulation, so SIMD handles several pixels per instruction — 8-bit images
• Naive → standard → state of the art:
  • 2-D kernel per output pixel · 36 taps for Lanczos3, rereads overlapping windows
  • separable Lanczos3 with precomputed weights · about 1 ns per source pixel here
  • SIMD with pre-shrink and fixed point (fast_image_resize, Pillow) · bounded by memory bandwidth; bilinear or box/area for previews and exact integer shrinks, Lanczos for final output

Implementation notes

• Layout: planar (one channel at a time) or interleaved RGBA with the channel loop innermost; interleaved with u8 and SIMD is what fast libraries do. Vertical pass should walk whole rows, never columns (stride access misses cache on every tap).
• Tile the work so the temporary fits in L2 (a band of rows per vertical output row), or process in strips for huge images (tiled pipeline).
• Colour: decode sRGB to linear, resample, encode back (colour spaces). Many libraries skip this by default and average in gamma space, which darkens fine detail.
• Alpha: premultiply before, divide after, then clamp (alpha compositing).
• Clamp the result to the valid range; cubic and Lanczos overshoot at edges. Do not clamp intermediate rows when the data is float and later passes need the true values.
• Cache coefficient tables keyed by (in, out, filter); thumbnailers resize thousands of images to the same size.
• Cropping/subpixel offset: shift center by the crop origin; do not crop first.

Pitfalls

• Using a fixed-width kernel when shrinking (OpenCV/GL linear/cubic, PyTorch without antialias=True): aliasing and moire. OpenCV's docs recommend INTER_AREA for shrinking.
• Forgetting the +0.5 centre convention: images shift by half a pixel per resize and drift over repeated resizes.
• Not renormalising weights at borders: edge pixels darken.
• Resampling in gamma-encoded sRGB: bright thin lines lose energy; alpha non-premultiplied: coloured fringes around transparent pixels.
• Negative weights with u8 accumulators: wrap-around instead of clamping gives black/white speckles.
• Lanczos on upscaled pixel art or text: visible halos; use a cubic with less ringing or nearest.
• Resizing repeatedly (e.g. every pipeline stage): each pass blurs; resample once from the original.

Variants and alternatives

In production

• Web and photo: Pillow (Image.resize, thumbnail), libvips (shrink then kernel reduce), the Rust image crate (imageops::resize) and fast_image_resize (SIMD), ImageMagick.
• Remote sensing and GIS: gdalwarp -r selects near, bilinear, cubic, cubicspline, lanczos, average and more when reprojecting rasters (docs).
• Video: FFmpeg's scaler offers bilinear, bicubic and lanczos with a configurable tap count (docs).
• Machine learning: torchvision's resize has an antialias option for downscaling (docs); OpenCV cv::resize is the common alternative.
• Games and GPUs: texture samplers (bilinear/trilinear hardware, with custom shaders for bicubic and Lanczos).

Implementations

Sources

• Robert G. Keys, Cubic Convolution Interpolation for Digital Image Processing, IEEE Trans. ASSP 29(6), 1981. paper
• Don P. Mitchell, Arun N. Netravali, Reconstruction Filters in Computer Graphics, SIGGRAPH 1988. paper
• Ken Turkowski, Filters for Common Resampling Tasks, Graphics Gems I, 1990.
• Pillow, Resample.c (coefficient computation used above) and the Image.resize docs: handbook