Pixel interpolation(binning?)

Question

I am trying to write up a pixel interpolation(binning?) algorithm (I want to, for example, take four pixels and take their average and produce that average as a new pixel). I've had success with stride tricks to speed up the "partitioning" process, but the actual calculation is really slow. For a 256x512 16-bit grayscale image I get the averaging code to take 7s on my machine. I have to process from 2k to 20k images depending on the data set. The purpose is to make the image less noisy (I am aware my proposed method decreases resolution, but this might not be a bad thing for my purposes).

import numpy as np
from numpy.lib.stride_tricks import as_strided
from scipy.misc import imread
import matplotlib.pyplot as pl
import time

def sliding_window(arr, footprint):
    """ Construct a sliding window view of the array"""
    t0 = time.time()
    arr = np.asarray(arr)
    footprint = int(footprint)
    if arr.ndim != 2:
        raise ValueError("need 2-D input")
    if not (footprint > 0):
        raise ValueError("need a positive window size")
    shape = (arr.shape[0] - footprint + 1,
             arr.shape[1] - footprint + 1, footprint, footprint)
    if shape[0] <= 0:
        shape = (1, shape[1], arr.shape[0], shape[3])
    if shape[1] <= 0:
        shape = (shape[0], 1, shape[2], arr.shape[1])
    strides = (arr.shape[1]*arr.itemsize, arr.itemsize,
               arr.shape[1]*arr.itemsize, arr.itemsize)

    t1 = time.time()
    total = t1-t0
    print "strides"
    print total

    return as_strided(arr, shape=shape, strides=strides)




def binning(w,footprint):

    #the averaging block
    #prelocate memory
    binned = np.zeros(w.shape[0]*w.shape[1]).reshape(w.shape[0],w.shape[1])
    #print w
    t2 = time.time()
    for i in xrange(w.shape[0]):
        for j in xrange(w.shape[1]):
            binned[i,j] = w[i,j].sum()/(footprint*footprint + 0.0)

    t3 = time.time()
    tot = t3-t2
    print tot

    return binned

Output:
5.60283660889e-05
7.00565886497

Is there some built-in/optimized function that would to the same thing I want, or should I just try and make a C extension (or even something else) ?

Below is the additional part of the code just for completeness, since I think the functions are the most important here. Image plotting is slow, but I think there is a way to improve it for example here

for i in range(2000):
    arr = imread("./png/frame_" + str("%05d" % (i + 1) ) + ".png").astype(np.float64)
    w = sliding_window(arr,footprint)
    binned = binning(w,footprint)
    pl.imshow(binned,interpolation = "nearest")
    pl.gray()
    pl.savefig("binned_" + str(i) + ".png")
    enter code here

What I am looking for could be called interpolation. I just used the term the person who advised me to do this used. Probably that is the reason why I was finding histogram related stuff !

Apart from the median_filter I tried generic_filter from scipy.ndimage but those did not give me the results I wanted (they had no "valid" mode like in convolution i.e. these relied on going out of bounds of the array when moving the kernel arround). I asked in code review and it seems that stackoverflow would be a more suitable place for this question.

Answer 1

Without diving into your code, I think what you what is just to resize the image with interpolation. You should use an image library for this operation, as it will have heavily optimized code.

Since you are using SciPy, you might want to start with PIL, the Python Imaging Library. Use the resize method, were you can pass the desired interpolation parameter, probably Image.BILINEAR in your case.

It should look something like this:

import Image
im = Image.fromarray(your_numpy)
im.resize((w/2, h/2), Image.BILINEAR)

Edit: I just noticed, you can do it even with scipy itself, look at the documentation for

scipy.misc.imresize

a = np.array([[0,1,2],[3,4,5],[6,7,8],[9,10,11]], dtype=np.uint8)
res = scipy.misc.imresize(a, (3,2), interp="bilinear")

Answer 2

For the average look at scipy.ndimage.filters.uniform_filter, in scipy.ndimage.filters, you have lots kernels for convolution that are much faster than a direct convolution with scipy.convolve.

Answer 3

For completeness, there's an interesting solution on scipy blog. The idea is to reshape the array to higher dimensions, then apply mean over one dimension and again, reshape into a smaller array. Supposedly it's fast as well.

https://scipython.com/blog/binning-a-2d-array-in-numpy/

Answer 4

For the ones looking for true binning, rather than interpolation or decimation: this is also provided in the Pillow module with the function Image.reduce. The output of Image.reduce is equal to the rebin method from scipython.com linked by @Tilen K.

image = np.arange(16).astype(float).reshape(4,4)

array([[ 0., 1., 2., 3.], [ 4., 5., 6., 7.], [ 8., 9., 10., 11.], [12., 13., 14., 15.]])

np.asarray(Image.fromarray(image).reduce(2))

array([[ 2.5, 4.5], [10.5, 12.5]], dtype=float32)