title: scipy.ndimage.interpolation.zoom uses nearest-neighbor-like algorithm for scaling-down
body: <p>While testing scipy's zoom function, I found that the results of scailng-down an array are similar to the nearest-neighbour algorithm, rather than averaging. This increases noise drastically, and is generally suboptimal for many application. </p>
  
  <p>Is there an alternative that does not use nearest-neighbor-like algorithm and will properly average the array when downsizing? While coarsegraining works for integer scaling factors, I would need non-integer scaling factors as well. </p>
  
  <p>Test case: create a random 100*M x 100*M array, for M = 2..20 
  Downscale the array by the factor of M three ways: </p>
  
  <p>1) by taking the mean in MxM blocks 
  2) by using scipy's zoom with a scaling factor 1/M 
  3) by taking a first point within a</p>
  
  <p>Resulting arrays have the same mean, the same shape, but scipy's array has the variance as high as the nearest-neighbor. Taking a different order for scipy.zoom does not really help.  </p>
  
  <pre><code>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()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/eaaIv.png" rel="noreferrer"><img src="https://i.stack.imgur.com/eaaIv.png" alt="enter image description here"></a></p>
  
  <p>EDIT: Performance of scipy.ndimage.zoom on a real noisy image is also very poor </p>
  
  <p><a href="https://i.stack.imgur.com/NuC6Z.png" rel="noreferrer"><img src="https://i.stack.imgur.com/NuC6Z.png" alt="enter image description here"></a></p>
  
  <p>The original image is here <a href="http://wiz.mit.edu/lena_noisy.png" rel="noreferrer">http://wiz.mit.edu/lena_noisy.png</a></p>
  
  <p>The code that produced it: </p>
  
  <pre><code>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()
  </code></pre>