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 存储的真实数据，我的输入数组被细分为几个块，但是为了提高我的测试速度，我将在这篇文章中使用一个简化的数组 (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

为了去掉double-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) 中应用中值滤波器（甚至可能是均值滤波器），与精确计算中值相比，会导致不可接受的错误同时过滤 (x,y,z)。