我需要计算形状(400,401,401)的3D数组的形态开口，大小为64320400字节，使用半径为17或更大的3D结构元素。结构元素narray的大小为42875字节。使用scipy.ndimage.morphology.binary_opening，整个进程消耗8GB内存。
我在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的内存。
有没有人对上面描述的例子中如何计算形态学有一些建议?
我也在做粒度测定时增加半径的开口，我遇到了同样的问题。实际上，内存使用大约以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. nimage相同。我的Binary_dilation。我不确定为什么会有差异，我想知道我是否错过了ref 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. nimage的结果。binary_erosion -你必须用1s填充反转到至少球面核b的半径R，因此可以实现侵蚀，从而获得与scipy. nimage .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获得的相同答案相同。nimage形态学操作(以及skimage操作，后者在后端只调用nimage操作)。
- 对于我最大的内核(R=21)，内存使用减少了30倍!速度也快了20倍。
- 我只对二进制图像进行了测试-我只是不知道灰度，但是在下面的第二个参考文献中有一些讨论。
还有两个提示:
首先:考虑我在中间部分讨论的关于erode_v1的填充。用1s填充逆基本上允许从数据集的边缘以及数据集中的任何接口发生侵蚀。根据您的系统和您想要做的事情，您可能需要考虑这是否真正代表了您想要的处理方式。如果没有，您可以考虑用"反射"边界条件填充，这将模拟边缘附近任何特征的延续。我建议尝试不同的边界条件(膨胀和侵蚀)，并将结果可视化和量化，以确定最适合您的系统和目标。
第二:这种基于频率的方法不仅在内存上更好，而且在速度上也更好——在大多数情况下。对于小核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。
- 球面掩模灰度图像的膨胀和侵蚀[/em]http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf
一种大胆的猜测是代码试图以某种方式分解结构元素并进行若干并行计算。每个计算都有一个完整原始数据的副本。400x400x400没那么大tbh…
我知道，由于您正在进行单个打开/关闭，它应该最多使用原始数据内存的3倍:原始数据+膨胀/侵蚀+最终结果…
你可以试着自己动手实现…它可能会慢一些，但是代码很简单，应该能让你对问题有一些了解……