scipy.ndimage.interpolation.zoom使用类最近邻算法进行缩容

Question

在测试scipy的缩放功能时，我发现向下搜索数组的结果类似于最近邻算法，而不是平均。这极大地增加了噪声，并且对于许多应用程序来说通常是次优的。

是否有一种替代方案不使用类似最近邻的算法，并且可以在缩减规模时适当地对阵列进行平均？虽然粗粒化适用于整数比例因子，但我也需要非整数比例因子。

测试用例:创建一个随机的100*M x 100*M数组，对于M= 2..20，按M的三种方式缩小数组：

1)通过取MxM块中的平均值2)通过使用具有缩放因子1/M的scipy缩放3)通过在

得到的数组具有相同的均值和形状，但scipy的数组的方差与最近邻的数组的方差一样高。对scipy.zoom采取不同的顺序并没有真正的帮助。

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()

​

​

编辑: scipy.ndimage.zoom在真实噪声图像上的性能也很差

​

​

原图在这里http://wiz.mit.edu/lena_noisy.png

生成它的代码：

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

没有人发布一个有效的答案，所以我将发布一个我目前使用的解决方案。不是最优雅的，但很管用。

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, e.g. 100x100 ->
    25x25, you just need to do block-averaging, that's 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我发现order=1至少比default或order=3更好地保留了平均值(真的如预期的)