Scipy.ndimage.interpolation.zoom uses a close neighbor algorithm to scale

Question

When testing the scipy zoom function, I found that the results of the scailng-down array are similar to the nearest neighbor algorithm, not averaging. This significantly increases noise and, as a rule, is not optimal for many applications.

Is there an alternative that does not use an algorithm close to the neighbor, and will correctly average the array during reduction? Although coarsegraining works for integer scaling factors, I also need non-integer scaling factors.

Test case: create a random array of 100 * M x 100 * M, for M = 2..20 Scaling an array with coefficient M in three ways:

1) by taking the average value in MxM blocks 2) using scipy zoom with a scaling factor of 1 / M 3) by taking the first point inside

The resulting arrays have the same average value, the same shape, but the scipy array has a variance as high as its nearest neighbor. Taking a different order for scipy.zoom really doesn't help.

import scipy.ndimage.interpolation import numpy as np import matplotlib.pyplot as plt mean1, mean2, var1, var2, var3 = [],[],[],[],[] values = range(1,20) # down-scaling factors for M in values: N = 100 # size of an array a = np.random.random((N*M,N*M)) # large array b = np.reshape(a, (N, M, N, M)) b = np.mean(np.mean(b, axis=3), axis=1) assert b.shape == (N,N) #coarsegrained array c = scipy.ndimage.interpolation.zoom(a, 1./M, order=3, prefilter = True) assert c.shape == b.shape d = a[::M, ::M] # picking one random point within MxM block assert b.shape == d.shape mean1.append(b.mean()) mean2.append(c.mean()) var1.append(b.var()) var2.append(c.var()) var3.append(d.var()) plt.plot(values, mean1, label = "Mean coarsegraining") plt.plot(values, mean2, label = "mean scipy.zoom") plt.plot(values, var1, label = "Variance coarsegraining") plt.plot(values, var2, label = "Variance zoom") plt.plot(values, var3, label = "Variance Neareset neighbor") plt.xscale("log") plt.yscale("log") plt.legend(loc=0) plt.show() 

EDIT: scipy.ndimage.zoom performance in real noise image is also very poor

Original image here http://wiz.mit.edu/lena_noisy.png

The code that created it:

 from PIL import Image import numpy as np import matplotlib.pyplot as plt from scipy.ndimage.interpolation import zoom im = Image.open("/home/magus/Downloads/lena_noisy.png") im = np.array(im) plt.subplot(131) plt.title("Original") plt.imshow(im, cmap="Greys_r") plt.subplot(132) im2 = zoom(im, 1 / 8.) plt.title("Scipy zoom 8x") plt.imshow(im2, cmap="Greys_r", interpolation="none") im.shape = (64, 8, 64, 8) im3 = np.mean(im, axis=3) im3 = np.mean(im3, axis=1) plt.subplot(133) plt.imshow(im3, cmap="Greys_r", interpolation="none") plt.title("averaging over 8x8 blocks") plt.show() 

Answer 1

No one posted a working answer, so I will post the solution I am currently using. Not the most elegant, but it works.

 import numpy as np import scipy.ndimage def zoomArray(inArray, finalShape, sameSum=False, zoomFunction=scipy.ndimage.zoom, **zoomKwargs): """ Normally, one can use scipy.ndimage.zoom to do array/image rescaling. However, scipy.ndimage.zoom does not coarsegrain images well. It basically takes nearest neighbor, rather than averaging all the pixels, when coarsegraining arrays. This increases noise. Photoshop doesn't do that, and performs some smart interpolation-averaging instead. If you were to coarsegrain an array by an integer factor, eg 100x100 -> 25x25, you just need to do block-averaging, that easy, and it reduces noise. But what if you want to coarsegrain 100x100 -> 30x30? Then my friend you are in trouble. But this function will help you. This function will blow up your 100x100 array to a 120x120 array using scipy.ndimage zoom Then it will coarsegrain a 120x120 array by block-averaging in 4x4 chunks. It will do it independently for each dimension, so if you want a 100x100 array to become a 60x120 array, it will blow up the first and the second dimension to 120, and then block-average only the first dimension. Parameters ---------- inArray: n-dimensional numpy array (1D also works) finalShape: resulting shape of an array sameSum: bool, preserve a sum of the array, rather than values. by default, values are preserved zoomFunction: by default, scipy.ndimage.zoom. You can plug your own. zoomKwargs: a dict of options to pass to zoomFunction. """ inArray = np.asarray(inArray, dtype=np.double) inShape = inArray.shape assert len(inShape) == len(finalShape) mults = [] # multipliers for the final coarsegraining for i in range(len(inShape)): if finalShape[i] < inShape[i]: mults.append(int(np.ceil(inShape[i] / finalShape[i]))) else: mults.append(1) # shape to which to blow up tempShape = tuple([i * j for i, j in zip(finalShape, mults)]) # stupid zoom doesn't accept the final shape. Carefully crafting the # multipliers to make sure that it will work. zoomMultipliers = np.array(tempShape) / np.array(inShape) + 0.0000001 assert zoomMultipliers.min() >= 1 # applying scipy.ndimage.zoom rescaled = zoomFunction(inArray, zoomMultipliers, **zoomKwargs) for ind, mult in enumerate(mults): if mult != 1: sh = list(rescaled.shape) assert sh[ind] % mult == 0 newshape = sh[:ind] + [sh[ind] // mult, mult] + sh[ind + 1:] rescaled.shape = newshape rescaled = np.mean(rescaled, axis=ind + 1) assert rescaled.shape == finalShape if sameSum: extraSize = np.prod(finalShape) / np.prod(inShape) rescaled /= extraSize return rescaled myar = np.arange(16).reshape((4,4)) rescaled = zoomArray(myar, finalShape=(3, 5)) print(myar) print(rescaled) 

Answer 2

FWIW I found that order = 1 at least keeps the average value much better than the default value or order = 3 (as expected)