我在这个论坛上看过几个关于将中值滤波器应用于移动窗口的讨论,但我的应用程序有一个特殊的特性.

我有一个尺寸为750x12000x10000的3D阵列,我需要应用中值滤波器来生成2D阵列(12000×10000).为此,每个中值计算应考虑固定的邻域窗口(通常为100×100)和所有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阵列)吗？

EDIT1
真实数据阵列有许多零值.当考虑单个轴时,在750个值中,大约15个是非零值.必须在处理中丢弃零,因此,我没有使用稀疏数组表示.

解决方法:

这最终导致评论太长：

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

原则上,中位数不是这样.然而,由于(x,y)和z中的邻域都相当大,我认为关联性仍然大致保持(除非你的数据来自一个可能不是因为这看起来像某种成像数据的笨拙的分布).如果我是你,我会测试一些测试数据,如果首先在z中应用中位数然后在(x,y)中的中值滤波器(或者甚至是均值滤波器)导致不可接受的误差与计算中值完全相比同时过滤(x,y,z).