Scipy ndimage形态学运算符占满我的8GB内存

5 投票

2 回答

1828 浏览

提问于 2025-04-18 15:19

我需要对一个形状为 (400, 401, 401) 的三维数组进行形态学开运算，这个数组的大小是 64320400 字节。我使用的三维结构元素的半径需要在 17 或更大。这个结构元素的大小是 42875 字节。使用 scipy.ndimage.morphology.binary_opening 进行处理时，整个过程消耗了 8GB 的内存。

我在 GitHub 上查看了 scipy/ndimage/morphology.py 的代码，按照我的理解，形态学的腐蚀操作是用纯 C 语言实现的。对我来说，理解 ni_morphology.c 的源代码有点困难，所以我没有找到导致如此巨大的内存使用的代码部分。增加内存并不是一个可行的解决方案，因为随着结构元素半径的增大，内存使用可能会呈指数级增长。

要重现这个问题：

import numpy as np
from scipy import ndimage

arr_3D = np.ones((400,401,401),dtype="bool")

str_3D = ndimage.morphology.generate_binary_structure(3,1)
big_str_3D = ndimage.morphology.iterate_structure(str_3D,20)

arr_out_3D = ndimage.morphology.binary_opening(arr_3D, big_str_3D)

这个过程大约需要 7GB 的内存。

有没有人能给我一些建议，如何在上述例子中计算形态学运算？

内存管理图像处理计算效率三维数组 C语言实现结构元素形态学运算腐蚀操作

2 个回答

我猜这个代码可能是在尝试分解一个结构元素，然后进行几次并行计算。每次计算都用的是整个原始数据的一个副本。其实400x400x400的数据量并不算大……

据我所知，因为你只进行了一次开运算或闭运算，所以它最多只会使用原始数据的三倍内存：原始数据 + 膨胀/腐蚀 + 最终结果……

你可以试着自己手动实现一下……虽然可能会慢一些，但这个代码其实很简单，应该能让你更好地理解这个问题……

回答于 2025-04-18 由 Python大师

分享举报

我也在做粒度分析时遇到了同样的问题，使用的开口半径逐渐增大。实际上，内存使用量大约是半径 R 的六次方增长的，这个增长速度可真不小！我做了一些内存分析，包括把开口操作分成侵蚀和膨胀（开口的定义），发现大内存使用主要来自 SciPy 的二进制文件，并且一旦结果返回给调用的 Python 脚本，内存就会被释放。SciPy 的形态学代码大多是用 C 语言实现的，所以修改起来比较困难。

不过，原帖最后的评论：“经过一些研究，我转向了使用卷积实现开口操作——即傅里叶变换的乘法——O(n log n)，并且内存开销不大。”这让我找到了问题的解决方案，感谢这个提示。不过，最开始实现起来并不明显。对于其他遇到这个问题的人，我会在这里分享实现的过程。

我会先讲讲膨胀，因为二进制侵蚀实际上就是对二进制图像的补集（反转）进行膨胀，然后再反转结果。

简单来说：根据Kosheleva 等人的这篇白皮书，膨胀可以看作是数据集 A 与结构元素（球形核）B 的卷积，结果会在某个值以上进行阈值处理。卷积也可以在频域中进行（通常速度更快），因为频域中的乘法等同于实域中的卷积。因此，先对 A 和 B 进行傅里叶变换，然后相乘，最后进行逆变换，再对结果进行阈值处理（大于 0.5），就得到了 A 和 B 的膨胀。（注意，我链接的白皮书说要阈值处理大于 0，但经过多次测试发现这样会产生很多伪影；Kukal 等人的另一篇白皮书给出的阈值是 >0.5，这样的结果与 scipy.ndimage.binary_dilation 完全一致。我不太清楚为什么会有这个差异，可能是我没理解第一篇文献中的某些细节)

正确实现这一点需要进行大小填充，但幸运的是，这已经在scipy.signal.fftconvolve(A,B,'same')中完成了——这个函数做的就是我刚才描述的，并且会为你处理填充。将第三个参数设为 'same' 会返回与 A 相同大小的结果，这正是我们想要的（否则结果会被 B 的大小填充）。

所以，膨胀的实现是：

from scipy.signal import fftconvolve
def dilate(A,B):
    return fftconvolve(A,B,'same')>0.5

侵蚀的原理是这样的：你先反转 A，然后像上面那样用 B 进行膨胀，最后再反转结果。但为了与 scipy.ndimage.binary_erosion 的结果完全一致，需要稍微做个小技巧——你必须将反转后的图像用 1 填充到至少球形核 B 的半径 R。这样就可以实现与 scipy.ndimage.binary_erosion 完全相同的结果。（注意，代码可以写得更简洁，但我这里想要更清楚地说明。）

from scipy.signal import fftconvolve
import numpy as np
def erode_v1(A,B,R):
    #R should be the radius of the spherical kernel, i.e. half the width of B
    A_inv = np.logical_not(A)
    A_inv = np.pad(A_inv, R, 'constant', constant_values=1)
    tmp = fftconvolve(A_inv, B, 'same') > 0.5
    #now we must un-pad the result, and invert it again
    return np.logical_not(tmp[R:-R, R:-R, R:-R])

你也可以通过另一种方式获得相同的侵蚀结果，正如Kukal 等人的白皮书所指出的那样——他们提到 A 和 B 的卷积可以通过阈值处理 > m-0.5 来实现侵蚀，其中 m 是 B 的“大小”（实际上是球体的体积，而不是数组的体积）。我先展示了 erode_v1，因为它稍微容易理解，但结果是一样的：

from scipy.signal import fftconvolve
import numpy as np
def erode_v2(A,B):
    thresh = np.count_nonzero(B)-0.5
    return fftconvolve(A,B,'same') > thresh

希望这能帮助到其他遇到这个问题的人。关于我得到的结果的一些说明：

我在 2D 和 3D 中都进行了测试，所有结果都与 scipy.ndimage 的形态学操作（以及 skimage 操作，后者实际上只是调用了 ndimage 的操作）完全一致。
对于我最大的核（R=21），内存使用量减少了 30 倍！速度也快了 20 倍。
不过我只在二进制图像上进行了测试——对于灰度图像我就不太清楚了，但在下面的第二个参考文献中有一些讨论。

再补充两点：

首先：考虑我在中间部分关于 erode_v1 讨论的填充。用 1 填充反转的边界，基本上允许侵蚀从数据集的边缘以及数据集中的任何接口进行。根据你的系统和你想要做的事情，你可能需要考虑这是否真正代表你想要的处理方式。如果不是，你可以考虑用“反射”边界条件进行填充，这样可以模拟边缘附近特征的延续。我建议尝试不同的边界条件（在膨胀和侵蚀上都试试），并可视化和量化结果，以确定哪种最适合你的系统和目标。

第二：这种基于频率的方法不仅在内存上更好，在速度上也是如此——大部分情况下。对于小的核 B，原始方法更快。然而，小核的运行速度本来就很快，所以对我来说并不重要。如果你在多次使用小核时关心速度，可能需要找到 B 的临界大小，并在那个点切换方法。

参考文献，虽然我很抱歉它们不容易引用，因为没有提供年份：

Fast Implementation of Morphological Operations Using Fast Fourier Transform by O. Kosheleva, S. D. Cabrera, G. A. Gibson, M. Koshelev. http://www.cs.utep.edu/vladik/misha5.pdf
Dilation and Erosion of Gray Images with Spherical Masks by J. Kukal, D. Majerova, A. Prochazka. http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf

回答于 2025-04-18 由 Python大师

分享举报

Scipy ndimage形态学运算符占满我的8GB内存

2 个回答

撰写回答