python - Scipy ndimage 形态学运算符使我的计算机内存 RAM (8GB) 饱和

Question

我需要使用半径为 17 或更大的 3D 结构元素来计算形状为 (400,401,401)、大小为 64320400 字节的 3D 数组的形态学开度。结构元素ndarray的大小为
42875 字节。使用 scipy.ndimage.morphology.binary_opening ，整个过程消耗8GB RAM。

我已阅读 scipy/ndimage/morphology.py在 GitHub 上，据我所知，形态侵 eclipse 算子是用纯 C 语言实现的。我很难理解 ni_morphology.c源，所以我没有发现这段代码的任何部分会导致如此巨大的内存利用率。添加更多 RAM 不是一个可行的解决方案，因为内存使用可能会随着结构元素半径呈指数增长。

要重现问题:

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 RAM。

有没有人对上述示例中如何计算形态学有一些建议？

Answer 1

我也做粒度增加半径的开口，我遇到了同样的问题。事实上，内存使用量大约增加了 R^6，其中 R 是球形内核的半径。这是一个相当大的增长率!我做了一些内存分析，包括将开孔拆分为腐 eclipse ，然后是膨胀(开孔的定义)，发现大量内存使用来自 SciPy 的二进制文件，并在结果返回到调用 Python 脚本后立即清除. SciPy 的形态代码大部分是用 C 实现的，因此修改它们是一个困难的前景。

无论如何，OP 的最后一条评论是:“经过一些研究，我转向使用卷积的 Opening 实现 -> 傅立叶变换的乘法 - O(n log n)，并且没有那么大的内存开销。”帮我找出解决方案，所以谢谢你。然而，实现起初并不明显。对于遇到此问题的任何其他人，我将在此处发布实现。

我将开始谈论膨胀，因为二值腐 eclipse 只是二值图像的补码(逆)的膨胀，然后将结果反转。

简而言之:根据this white paper by Kosheleva et al ，膨胀可以被视为数据集 A 与结构元素(球核)B 的卷积，阈值高于某个值。卷积也可以在频率空间中完成(通常快得多)，因为频率空间中的乘法与实际空间中的卷积相同。因此，通过首先对 A 和 B 进行傅立叶变换，将它们相乘，然后对结果进行逆变换，然后对大于 0.5 的值进行阈值处理，就可以得到 A 与 B 的膨胀。(请注意，我链接的白皮书说阈值高于 0，但大量测试表明，这给出了许多工件的错误结果；another white paper by Kukal et al . 将阈值设为 >0.5，这对我来说给出了与 scipy.ndimage.binary_dilation 相同的结果。我不知道为什么差异，我想知道我是否错过了引用文献 1 命名法的一些细节)

正确的实现涉及填充大小，但对我们来说幸运的是，它已经在 scipy.signal.fftconvolve(A,B,'same') 中完成了。 - 此函数执行我刚刚描述的操作并为您处理填充。将第三个选项设为“相同”将返回与 A 大小相同的结果，这正是我们想要的(否则它将被 B 的大小填充)。

所以膨胀是:

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

原则上的侵 eclipse 是这样的:你把 A 取反，用 B 扩大它，然后重新取反结果。但是它需要一个小技巧来精确匹配来自 scipy.ndimage.binary_erosion 的结果 - 你必须用 1s 填充反演到至少球形内核 B 的半径 R。因此可以实现侵 eclipse 以获得与 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])

您可以通过另一种方式获得相同的侵 eclipse 结果，如 the white paper by Kukal et al 所示。 - 他们指出 A 和 B 的卷积可以通过 > m-0.5 的阈值化来进行侵 eclipse ，其中 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 的中间部分讨论的填充。用 1s 填充倒数基本上允许从数据集的边缘以及数据集中的任何界面发生侵 eclipse 。根据您的系统和您尝试执行的操作，您可能需要考虑这是否真正代表了您想要的处理方式。如果没有，您可以考虑使用“反射”边界条件进行填充，这将模拟边缘附近任何要素的延续。我建议尝试使用不同的边界条件(膨胀和侵 eclipse )并对结果进行可视化和量化，以确定最适合您的系统和目标的方法。

第二:这种基于频率的方法不仅在内存方面更好，而且在速度方面也更好——在大多数情况下。对于小内核 B，原始方法更快。然而，无论如何，小内核运行得非常快，所以为了我自己的目的，我不在乎。如果您这样做(就像您多次执行小内核一样)，您可能希望找到 B 的临界大小并在此时切换方法。

引用文献，尽管我很抱歉它们不容易引用，因为它们都没有提供年份:

使用快速傅立叶变换快速实现形态学操作 作者:O. Kosheleva、S. D. Cabrera、G. A. Gibson、M. Koshelev。 http://www.cs.utep.edu/vladik/misha5.pdf

J. Kukal、D. Majerova、A. Prochazka 对带有球形蒙版的灰色图像的膨胀和侵 eclipse 。 http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf