Fast Interpolation / Resample of Numpy Array - Python

TLDR; Check out Skimage's pyramid_gaussian. On a single image of (512, 512) it shows an order of 0.3M times speed-up. (163 ms/471 ns = 346072). Pillow-SIMD does super-fast resamples/resize but requires you to uninstall PIL, Pillow before you install it. It uses parallel processing (single instruction, multiple data - SIMD) and better algorithms such as replacing a convolution-based Gaussian blur with a sequential-box one. Advice to use this for production environments after setting up a separate venv.

There are multiple ways of upsampling and downsampling your image. I will add some benchmarks of the methods that I have worked with. I will keep updating this answer as I come across more methods so that this can act as a reference for others on SO.

#Utility function for plotting original, upsampled, and downsampled image

def plotit(img, up, down):
    fig, axes = plt.subplots(1,3, figsize=(10,15))
    axes[0].imshow(img)
    axes[1].imshow(up)
    axes[2].imshow(down)
    axes[0].title.set_text('Original')
    axes[1].title.set_text('Upsample')
    axes[2].title.set_text('Downsample')

IIUC, this is somewhat your pipeline -

from scipy.ndimage import zoom
from skimage.data import camera

img = camera() #(512,512)
up = zoom(img,2) #upsample image

#some code
...

down = zoom(up,0.5) #downsample the upsampled image

plotit(img, up, down)

Methods & Benchmarks (in no specific order)

1. Scipy zoom (order=3)

The array is zoomed using spline interpolation of a given order, in this case, the default is order = 3.

%%timeit
#from scipy.ndimage import zoom
up = zoom(img,2)
down = zoom(up,0.5)

#163 ms ± 12.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

2. Scipy zoom (order=0)

The array is zoomed using spline interpolation of a given order, in this case, order = 0.

%%timeit
#from scipy.ndimage import zoom
up = zoom(img,2, order=0)
down = zoom(up,0.5, order=0)

#18.7 ms ± 950 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

3. Skimage pyramid_gaussian

In a Gaussian pyramid, subsequent images are weighted down using a Gaussian average (Gaussian blur) and scaled-down. Applying a Gaussian blur to an image is the same as convolving the image with a Gaussian function. The amount of blur depends on the standard deviation size (sigma).

%%timeit
#from skimage.transform import import pyramid_gaussian
up = pyramid_gaussian(img,2)
down = pyramid_gaussian(up,0.5)

#471 ns ± 30.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

4. Skimage pyramid_expand and pyramid_reduce

An image pyramid is a type of multi-scale image representation in which an image is subject to repeated smoothing and subsampling. The first function smooths and then upsamples the image, while the second does the same but instead downsamples, both of these use spline order=1 by default.

%%timeit
#from skimage.transform import import pyramid_expand, pyramid_reduce
up = pyramid_expand(img,2)
down = pyramid_reduce(up,2)

#120 ms ± 3.08 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

5. Skimage rescale

Scales an image by a certain factor. Performs spline interpolation (order=1 by default) to up-scale or down-scale N-dimensional images. Note that anti-aliasing should be enabled when down-sizing images to avoid aliasing artifacts.

%%timeit
#from skimage.transform import import rescale
up = rescale(img,2)
down = rescale(up,0.5)

#83 ms ± 3.69 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

6. PIL resize with Nearest pixel filter for resampling

Returns a resized copy of this image. Picks one nearest pixel from the input image. Ignores all other input pixels.

%%timeit
#from PIL import Image
im = Image.fromarray(img)
up = im.resize((im.width*2, im.height*2),resample=Image.NEAREST)
down = up.resize((up.width//2, up.height//2),resample=Image.NEAREST)

#704 µs ± 29.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

7. PIL resize with BILINEAR filter for resampling

Returns a resized copy of this image. For resize calculates the output pixel value using linear interpolation on all pixels that may contribute to the output value. For other transformations, linear interpolation over a 2x2 environment in the input image is used.

%%timeit
#from PIL import Image
im = Image.fromarray(img)
up = im.resize((im.width*2, im.height*2),resample=Image.BILINEAR)
down = up.resize((up.width//2, up.height//2),resample=Image.BILINEAR)

#10.2 ms ± 877 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

8. PIL resize with BICUBIC filter for resampling

Returns a resized copy of this image. For resize calculates the output pixel value using cubic interpolation on all pixels that may contribute to the output value. For other transformations cubic interpolation over a 4x4 environment in the input image is used.

%%timeit
#from PIL import Image
im = Image.fromarray(img)
up = im.resize((im.width*2, im.height*2),resample=Image.BICUBIC)
down = up.resize((up.width//2, up.height//2),resample=Image.BICUBIC)

#12.3 ms ± 326 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

9. PIL resize with Lanczos filter for resampling

Returns a resized copy of this image. Calculates the output pixel value using a high-quality Lanczos filter (a truncated sinc) on all pixels that may contribute to the output value.

%%timeit
#from PIL import Image
im = Image.fromarray(img)
up = im.resize((im.width*2, im.height*2),resample=Image.LANCZOS)
down = up.resize((up.width//2, up.height//2),resample=Image.LANCZOS)

#15.7 ms ± 184 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Check the benchmarks in https://python-pillow.org/pillow-perf/

TL;DR: for most scenarios the fastest way in Python is using pillow-simd.