The Magic Kernel

Happy to answer questions or address criticisms.

There is an article [0] that points this out, although later it goes into recommending odd filters, which insidiously break many applications where you would want to use downscaling by shifting the image by a fraction of a pixel.

[0] https://cbloomrants.blogspot.com/2011/03/03-24-11-image-filt...

Fourier properties

of the Magic Kernel Sharp algorithm, and discovered that it possesses high-order zeros at multiples of the sampling frequency, which explains why it produces such amazingly clear results: these high-order zeros greatly suppress any potential aliasing artifacts that might otherwise appear at low frequencies..."

[...]

>"Following that work, I analytically derived the

Fourier transform

of the Magic Kernel in closed form, and found, incredulously, that

it is simply the cube of the sinc function.

This implies that the Magic Kernel is just the

rectangular window function convolved with itself twice

-which, in retrospect, is completely obvious. This observation, together with a precise definition of the requirement of the Sharp kernel, allowed me to obtain an analytical expression for the exact Sharp kernel, and hence also for the exact Magic Kernel Sharp kernel, which I recognized is just the third in a sequence of fundamental resizing kernels. These findings allowed me to explicitly show why Magic Kernel Sharp is superior to any of the Lanczos kernels. It also allowed me to derive further members of this fundamental sequence of kernels, in particular the sixth member, which has the same computational efficiency as Lanczos-3, but has far superior properties."

Paper: http://www.johncostella.com/magic/mks.pdf

PDS: My thoughts: Absolutely fascinating! (Also, that the the "Magic Kernel is just the rectangular window function convolved with itself twice" -- was not obvious to me, so it's good you included that!)

> Computationally, however, nothing short of Fast Fourier Transforms or hardware acceleration can make a large blurring job efficient: for large values of blur, the support of any blur kernel is large, and not practical for some at-scale applications.

The standard way to do Gaussian blur is to approximate it as a B-spline, specifically a cascade of simple moving average filters, also known as a Hogenauer filter or CIC filter in the time domain. Typically a cascade of three filters is used, precisely as Costella describes. This is both simpler and faster than doing the task with a Fourier transform, and moreover it is an online algorithm — it produces output samples at some processing lag behind the input samples.

As you can easily verify, the following program uses this approach to convolve its input buffer with a Gaussian kernel with support of 766 samples, constructed in just the way Costella describes (as a B-spline), in a single pass over the input buffer, with only six additions and subtractions per sample, not counting the additions and subtractions used for addressing.

    #include <stdio.h>

    enum {lag = 256, image_size = 4096};

    unsigned long buf[image_size] = {0};

    int
    main(int argc, char **argv)
    {
      buf[image_size/2] = 1;

      for (size_t i = 0; i < image_size - 3 * lag - 3; i++) {
        buf[i + 3 * lag + 2] += buf[i + 3 * lag + 1];
        buf[i + 3 * lag + 1] += buf[i + 3 * lag];
        buf[i + 3 * lag] += buf[i + 3 * lag - 1];
        buf[i + 2 * lag] -= buf[i + 3 * lag];
        buf[i + lag] -= buf[i + 2 * lag];
        buf[i] = buf[i + lag] - buf[i];
      }

      for (size_t i = 0; i < image_size - 3 * lag - 3; i++) {
        if (i) printf((i%8) ? ", " : ",\n");
        printf("%6ld", buf[i]);
      }
      printf("\n");

      return 0;
    }

Adding an outer loop reveals that on my nine-year-old laptop it takes 6.5 μs to run over the whole buffer, or 1.6 nanoseconds per sample. In theory this performance does not depend on the lag and thus the blur radius, but I've only tested a single lag.

The output samples are in some sense scaled by 2²⁴ = 16777216, or by 2¹⁶ = 65536 if you're interpolating a signal that's being upsampled by zero insertion by a factor of 256. This scale factor does depend on the lag (it is the cube of the lag).

This version computes the results in place starting from output that is already in memory, but you can also use a ring buffer with a capacity of at least 3lag + 2 = 770 samples. If you're downsampling, an alternative is to decimate the samples as soon as they've passed through the cascade of three integrators at the top of the loop, thus reducing both the amount of memory needed and the number of operations to do.

Overflow of intermediate results can happen, but it is not a problem as long as the final output is within range, because exact integer addition mod k is a group. If you were doing this in floating point, you might prefer to interleave the integrators with the comb filters to prevent intermediate overflow and consequent loss of precision leading to numerical instability.

Only the parts of the result buffer over which all six pipeline stages have passed contain valid results. The simplest solution to this is padding, but you can also use special logic to run the integrators but not the combs over the beginning of the input.

I would argue that 17 instructions or 1.6 nanoseconds per sample is always "efficient" and "practical" for nearly all "at-scale applications". If you were blurring a 3840×2160 display at 120 fps, one gigapixel per second, it would take 13 milliseconds of my laptop CPU cores per second for each of the R, the G, and the B; you could do it on about 6 cores of a modern CPU, without a GPU.