Python 中值滤波器应用于 3D 数组以产生 2D 结果

Question

我在这个论坛上看到过几个关于在移动窗口中应用中值滤波器的讨论，但我的应用程序有一个特殊之处。

我有一个维度为 750x12000x10000 的 3D 数组，我需要应用中值滤波器来生成 2D 数组 (12000x10000)。为此，每个中值计算都应考虑固定的邻域窗口(通常为 100x100)和所有 z 轴 值。矩阵中有一些零值，在计算中位数时不应考虑它们。为了处理真实数据，我使用了 numpy.memmap:

fp = np.memmap(filename, dtype='float32', mode='w+', shape=(750, 12000, 10000))

为了处理用 memmap 存储的真实数据，我的输入数组被分割为几个 block ，但为了提高我的测试速度，我将在这篇文章中使用一个简化的数组 (11, 200, 300) 和一个较小的窗口 (11, 5, 5) 或 (11, 50, 50)，我希望得到一个结果矩阵 (200, 300):

import numpy as np
from timeit import default_timer as timer

zsize, ysize, xsize = (11, 200, 300)
w_size = 5 #to generate a 3D window (all_z, w_size, w_size)
#w_size = 50 #to generate a 3D window (all_z, w_size, w_size)

m_in=np.arange(zsize*ysize*xsize).reshape(zsize, ysize, xsize)
m_out = np.zeros((ysize, xsize))

首先，我尝试了蛮力法，但它如预期的那样非常慢(即使对于小数组也是如此):

start = timer()
for l in range(0, ysize):
    i_l = max(0, l - w_size/2)
    o_l = min(ysize, i_l+w_size/2)
    for c in range(0, xsize):
        i_c = max(0, c - w_size/2)
        o_c = min(xsize, i_c+w_size/2)
        values = m_in[:, i_l:o_l, i_c:o_c]
        values = values[np.nonzero(values)]
        value = np.median(values)
        m_out[l, c] = value
end = timer()
print("Time elapsed: %f seconds"%(end-start))
#11.7 seconds with 50 in z, 7.9 seconds with 5 in z

要删除双 for，我尝试使用 itertools.product，但它仍然很慢:

from itertools import product
for l, c in product(range(0, ysize), range(0, xsize)):
    i_l = max(0, l - w_size/2)
    o_l = min(ysize, i_l+w_size/2)
    i_c = max(0, c - w_size/2)
    o_c = min(xsize, i_c+w_size/2)
    values = m_in[:, i_l:o_l, i_c:o_c]
    values = values[np.nonzero(values)]
    value = np.median(values)
    m_out[l, c] = value
#11.7 seconds with 50 in z, 2.3 seconds with 5

所以我尝试利用numpy的矩阵运算的性能，所以我尝试用scipy.ndimage:

from scipy import ndimage
m_all = ndimage.median_filter(m_in, size=(zsize, w_size, w_size))
m_out[:] = m_all[0] #only first layer of 11, considering all the same
#a lot of seconds with 50 in z, 7.9 seconds with 5

和scipy.signal也是:

m_all = signal.medfilt(m_in, kernel_size=(zsize, w_size, w_size))
m_out[:] = m_all[0] #only first layer of 11, considering all the same
#a lot of seconds with 50 in z, 7.8 seconds with 5 in z

但在这两种 scipy 情况下，都存在处理浪费，因为该函数应用于输入矩阵的所有 3D 位置，但是，它只能应用于第一层，使用维度为 (all_z, w_size, w_size)。

在我所有的测试中，即使我使用了缩减矩阵和窗口 ((11, 200, 300) 和 (11, 50, 50))，我的执行时间也并不快。使用我的真实数据(750x12000x10000 的数组和 750x100x100 的窗口)，性能将更加关键。

拜托，谁能帮助我以更好的 pythonic 方式应用中值滤波器(3D 数组到 2D 数组)？

编辑1 真实数据数组有很多零值。当考虑单个轴时，在 750 个值中，大约有 15 个是非零值。在处理过程中必须丢弃零，因此，我没有使用稀疏数组表示。

Answer 1

这对于评论来说太长了:

如果您应用均值滤波器，这个问题将是微不足道的:您将在 z 轴上取均值，然后在 2D 中应用均值滤波器；这将完全等同于一次计算整个 (x,y,z) 邻域的平均值，因为平均值操作是关联的(如果这是术语；我的意思是:f(f(a,b), c) = f(a, b, c))。

原则上，中位数并非如此。但是，由于您在 (x,y) 和 z 中的邻域都相当大，我认为关联性仍然大致成立(除非您的数据是从一个古怪的分布中提取的，这可能不是因为这看起来像某种成像数据).如果我是你，我会测试一些测试数据，如果先在 z 中应用中值，然后在 (x,y) 中应用中值滤波器(或者甚至可能是均值滤波器)，与准确计算中值相比，会导致 Not Acceptable 错误同时过滤 (x,y,z)。