Lanczos resampling of pixel maps in C

• 27 November 2024
• C/C++,
• Graphics

Lanczos resampling is an anti-aliasing technique originating from signal processing that is used in filtering and improving image quality. At a high-level, it filters out high-frequency artifacts for smooth interpolation. We'll look at how this works shortly.

We use PPM (Portable Pixel Map) image files because it's very basic and easier to read or write to than JPEG or many other major formats. Let's start with the reconstruction kernel which is defined in terms of $sinc$ for windowing and normalization.

$$\text{sinc}(x)=\{\begin{array}{ll}1& \text{if }x=0,\\ \frac{\sin (\pi x)}{\pi x}& \text{otherwise.}\end{array}$$

$$L(x,a)=\{\begin{array}{ll}\text{sinc}(x)\cdot \text{sinc}\left(\frac{x}{a}\right)& \text{if }|x|<a,\\ 0& \text{otherwise.}\end{array}$$

The kernel limits its output to a small neighborhood by clamping $sinc$ to a window size (lossless).

static double sinc(double x) {
  return x == 0.0 ? 1.0 : sin(M_PI * x) / (M_PI * x);
}

static double lanczos_kernel(double x, int a) {
  return fabs(x) < a ? sinc(x) * sinc(x / a) : 0.0;
}

At a given channel $c$, we have an output pixel $I^{\prime }(ox,oy)$ which is formed from the weighted sum of nearby pixels under interpolation from our kernel. We have:

$$I^{\prime }(ox,oy,c)=\frac{\sum _{x=x_s}^{x_e}\sum _{y=y_s}^{y_e}I(x,y,c)\cdot L(s_x-x,a)\cdot L(s_y-y,a)}{\sum _{x=x_s}^{x_e}\sum _{y=y_s}^{y_e}L(s_x-x,a)\cdot L(s_y-y,a)},$$

such that:

$$s_x=ox\cdot \frac{iw}{ow},s_y=oy\cdot \frac{ih}{oh}$$

$$x_s=\lfloor s_x\rfloor -a+1,x_e=\lfloor s_x\rfloor +a$$

$$y_s=\lfloor s_y\rfloor -a+1,y_e=\lfloor s_y\rfloor +a$$

and $a$ is the Lanczos kernel window size.

We start with $I(x,y,c)$, the input channels. Lanczos kernel provides weights based on the distance between the original pixel location $(x,y)$ and the fractional source coordinates $(s_x,s_y)$.

void lanczos_resample(const uint8_t* input, int iw, int ih, int ic,
                      uint8_t* output, int ow, int oh, int a) {
  double x_ratio = (double)iw / ow;
  double y_ratio = (double)ih / oh;

  for (int oy = 0; oy < oh; oy++) {
    for (int ox = 0; ox < ow; ox++) {
      for (int c = 0; c < ic; c++) {
        double acc = 0.0, w_s = 0.0;

        // fractional source coords

        double s_x = ox * x_ratio;  // <--
        double s_y = oy * y_ratio;  // <--
      }
    }
  }
}

*_ratio is provided for scaling the input image itself which we pass around as a buffer. The looping is row-major since PPM files store the buffer sequentially.

$$⟶(row)\downarrow (col)$$

We find the neighborhood by determining the range of pixels $(xs,xe,ys,ye)$ in the input image that influence the current output pixel (i.e. per-pixel, so it runs in the innermost loop). This range depends on the kernel size $a$.

int x_s = (int)floor(s_x) - a + 1;
int x_e = (int)floor(s_x) + a;
int y_s = (int)floor(s_y) - a + 1;
int y_e = (int)floor(s_y) + a;

The Lanczos kernel covers a symmetric range of input pixels around the central position which is why we consistently add $-a+1$ and $a$ for start and end indices respectively.

Why?

When computing the starting index $(xs,ys)$, subtracting $a$ would include $a$ pixels before the central position. However, pixel indices start at 0, and the inclusion of $+1$ forces the starting index to stay within the inclusive range of pixels that the kernel is applied to.

Consider $s_x=5.7$ (fractional coordinate):

• $\lfloor s_x\rfloor =5$ (nearest integer below $s_x$).
• ${\ast }_s$​ starts from $\lfloor s_x\rfloor -a+1=5-3+1=3$.
• ${\ast }_e$​ ends at $\lfloor s_x\rfloor +a=5+3=8$.

This specific range $[3,4,5,6,7]$ includes $2a-1$ pixels symmetrically around $s_x$ and perfectly aligns with the kernel’s range. Wikipedia has more resources on the implications of this.

• For each pixel $(x,y)$ in the neighborhood:

  • Compute its weight as the product of two Lanczos kernels: $$w=L(sr{c}_x-x,a)\cdot L(sr{c}_y-y,a)$$
  • Accumulate the weighted contribution to the output pixel value: $$\text{acc}\text{+=}I(x,y,c)\cdot w$$
  • Accumulate the total weight: $$w_s\text{+=}w$$

$$I^{\prime }(ox,oy,c)=\frac{\text{acc}}{w_s}$$

for (int y = y_s; y <= y_e; y++) {
  for (int x = x_s; x <= x_e; x++) {
    if (x >= 0 && x < iw && y >= 0 && y < ih) {
      double dx = s_x - x;
      double dy = s_y - y;
      double w = lanczos_kernel(dx, a) * lanczos_kernel(dy, a);

      acc += input[(y * iw + x) * ic + c] * w;
      w_s += w;
    }
  }
}

We essentially write to the output channel using the input buffer. It's extremely important to not select any pixels outside the range for visual clarity. The resulting image will appear with some level of blur and the pixelated parts get aliased.

output[(oy * ow + ox) * ic + c] = w_s > 0.0 ? (uint8_t)(acc / w_s) : 0;

I tried this with a sample image. The original image is at 42.9 kB, resampled output was at 171.7 kB (roughly 75% gain in size and this is important since PPM does not support compression). There's really not much else besides reading and writing to a PPM and selecting a scale factor (2 or 4 work well, find the final source below).

The PNG equivalents above do not carry all the details but still show a discrepancy in file size, resampled output image converted to PNG has a larger file size.

Acknowledgements #

Cornell CS 664 - Computer Vision Prof. Dan Huttenlocher, Fall 2003, sample house image.

Resources #

• Lancsoz resampling - Wikipedia

Cf. final source.

• Previous: Endianness in the JVM bytecode layout