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

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

要删除双重 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数组）吗？

Edit1 真实数据数组有许多零值。当考虑单个轴时，在 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)。