5

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

我在 GitHub 上读过scipy/ndimage/morphology.py,据我所知,形态侵蚀算子是用纯 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。

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

4

2 回答 2

5

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

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

我将开始谈论膨胀,因为二进制侵蚀只是二进制图像的补码(逆)的膨胀,然后结果被反转。

简而言之:根据Kosheleva 等人的这篇白皮书,膨胀可以看作是数据集 A 与结构元素(球形核)B 的卷积,阈值高于某个值。卷积也可以在频率空间中完成(通常更快),因为频率空间中的乘法与实际空间中的卷积相同。因此,通过首先对 A 和 B 进行傅里叶变换,将它们相乘,然后对结果进行逆变换,然后对大于 0.5 的值进行阈值处理,就可以得到 A 与 B 的膨胀。(请注意,我链接的白皮书说阈值高于 0,但大量测试表明,这给出了许多伪影的错误结果;Kukal 等人的另一篇白皮书. 将阈值设为 >0.5,这对我来说给出了与 scipy.ndimage.binary_dilation 相同的结果。我不确定为什么会出现这种差异,我想知道我是否错过了参考 1 命名法的一些细节)

正确的实现涉及到大小的填充,但对我们来说幸运的是,它已经完成了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 的结果 - 你必须用 1s 填充反演到至少球形核 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 的临界大小并在那个时候切换方法。

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

  1. O. Kosheleva、SD Cabrera、GA Gibson、M. Koshelev使用快速傅里叶变换快速实现形态学运算。http://www.cs.utep.edu/vladik/misha5.pdf
  2. J. Kukal, D. Majerova, A. Prochazka用球形面具对灰色图像进行膨胀和侵蚀。http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf
于 2018-07-03T16:42:10.563 回答
1

一个疯狂的猜测是代码试图以某种方式分解结构元素并进行多个并行计算。每个计算都有自己的整个原始数据的副本。400x400x400 没那么大...

AFAIK,由于您正在进行单次打开/关闭,它最多应该使用原始数据的 3 倍内存:原始 + 膨胀/腐蚀 + 最终结果......

如果您自己手动实现,您可以尝试实现......它可能会更慢,但代码很简单,应该对问题有所了解......

于 2015-09-18T13:01:42.397 回答