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

Posted 2023-02-23

【中文标题】Python 中值滤波器应用于 3D 数组以产生 2D 结果

【英文标题】：Python median filter applied to 3D array to produce 2D result

【发布时间】：2018-09-19 07:10:50

【问题描述】：

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

我有一个尺寸为 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 个是非零值。在处理过程中必须丢弃零，因此，我没有使用稀疏数组表示。

【问题讨论】：

数据的真实形状/数组顺序是什么？在文本中它是 (10000x12000x750)，但在 memmap 示例中它的顺序不同（750、12000、10000）？为什么在这里使用非分块数组，有几个缺点（只能有效访问变化最快的维度）？

对不起，我写的解释有误。我的实际数据数组有维度（750、12000、10000）。我提交了不同维度的分配代码，因为我正在测试使用不同的步幅是否有任何收益，但我会对其进行编辑以避免混淆。是的，我正在使用分块访问，但只是为了读取来自 memmap 的块。我一次将整个数组复制到 memmap，因为在我的测试中，对原始数据结构 (gdal.Dataset) 的分块访问比 memmap 慢一点，但如果你提出更好的方法，我可以更改它。

如果你在数据数组中的稀疏度是 98%，那么你真的应该换个角度思考这个问题。实际上，您对 x,y 中的大多数位置进行了多次测量（但可能对于某些没有数据的位置），并且您正在尝试估计局部平均值。一定要转换为稀疏格式来做到这一点。

如果您的支持大致一致（数据点的数量不会因地区而异），那么几乎任何基于网格的插值都可以。如果它确实有所不同，我建议使用 k-最近邻平均值（我有一个加权 k-最近邻平均值实现 here，但是在 x,y 中获得 k-最近数据点的中位数也很简单） )。

感谢您对 Paul 的评论和建议。是的，我的数据大约 98% 是稀疏的，但仅在一个轴上。也许更容易想到维度为 (12000, 10000, 750) 的数据，因为这样一来，所有 (x,y) 都有值，但是在 z 轴上它并没有所有 750 个值，而只有 20 个非零值(2%)。

【参考方案1】：

这对评论来说太长了：

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

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

【讨论】：

感谢您的建议。你是正确的应用平均值。我的算法的第一个版本是使用均值实现的，结果还不错。但我注意到中位数改善了结果，因为我有许多被噪声污染的值。我肯定会得到不同的结果，因为应该丢弃许多零（我编辑了原始帖子以包含这些重要信息），但轴 (z) 的中位数上的中位数 (x, y) 可能具有中间行为在全中位数 (x, y, z) 和均值 (x, y, mean (z)) 的结果之间。