Mar*_*urz 5 python numpy image-processing scipy mathematical-morphology
我需要使用半径为 17 或更大的 3D 结构元素来计算 3D 形状数组 (400,401,401)、大小为 64320400 字节的形态学开度。结构元素 ndarray 的大小为 42875 字节。使用scipy.ndimage.morphology.binary_opening，整个过程消耗 8GB RAM。
I have read scipy/ndimage/morphology.py on GitHub, and as far as I can tell, the morphology erosion operator is implemented in pure C. It is to difficult for me to understand the ni_morphology.c source, so I haven't found any part of this code which leads to such enormous memory utilization. Adding more RAM is not a workable solution, since memory usage may increase exponentially with the structure element radius.
To reproduce the problem:
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)
This takes approximately 7GB RAM.
Does anyone have some suggestions for how compute morphology in example described above?
我也做粒度增加半径的开口，我遇到了同样的问题。事实上，内存使用量大约增加了 R^6，其中 R 是球形内核的半径。这是一个相当大的增长率！我做了一些内存分析，包括将开孔拆分为腐蚀，然后是膨胀（开孔的定义），发现大量内存使用来自 SciPy 的二进制文件，并在结果返回到调用 Python 脚本后立即清除. SciPy 的形态代码大部分是用 C 实现的，因此修改它们是一个困难的前景。
无论如何，OP 的最后一条评论是：“经过一些研究，我转向使用卷积的 Opening 实现 -> 傅立叶变换的乘法 - O(n log n)，并且没有那么大的内存开销。” 帮我找出解决方案，所以谢谢你。然而，实施起初并不明显。对于遇到此问题的任何其他人，我将在此处发布实现。
我将开始谈论膨胀，因为二值腐蚀只是二值图像的补码（逆）的膨胀，然后将结果反转。
简而言之：根据Kosheleva 等人的这篇白皮书，扩张可以被视为数据集 A 与结构元素（球核）B 的卷积，阈值高于某个值。卷积也可以在频率空间中完成（通常快得多），因为频率空间中的乘法与实际空间中的卷积相同。因此，通过采取傅立叶变换A和B的第一，乘以它们，然后逆变换的结果，然后进行阈值是对大于0.5的值，你得到的扩张与B（注白纸我联系说阈值高于 0，但大量测试表明，这会产生错误的结果和许多伪影；Kukal 等人的另一篇白皮书. 给出阈值 > 0.5，这给了我与 scipy.ndimage.binary_dilation 相同的结果。我不确定为什么会出现这种差异，我想知道我是否遗漏了 ref 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
原则上的侵蚀是这样的：你把 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
我希望这可以帮助其他遇到此问题的人。关于我得到的结果的说明：
另外两个快速说明：
首先：考虑我在关于 erode_v1 的中间部分讨论的填充。用 1s 填充倒数基本上允许从数据集的边缘以及数据集中的任何界面发生侵蚀。根据您的系统和您尝试执行的操作，您可能需要考虑这是否真正代表了您想要的处理方式。如果没有，您可以考虑使用“反射”边界条件进行填充，这将模拟边缘附近任何特征的延续。我建议尝试使用不同的边界条件（膨胀和侵蚀）并对结果进行可视化和量化，以确定最适合您的系统和目标的方法。
第二：这种基于频率的方法不仅在内存方面更好，而且在速度方面也更好——在大多数情况下。对于小内核 B，原始方法更快。然而，无论如何，小内核运行得非常快，所以为了我自己的目的，我不在乎。如果您这样做（就像您多次执行小内核一样），您可能希望找到 B 的临界大小并在此时切换方法。
参考文献，尽管我很抱歉它们不容易引用，因为它们都没有提供年份：
一个大胆的猜测是，代码试图以某种方式分解结构元素并进行多个并行计算。每个计算都有其自己的整个原始数据的副本。400x400x400 并不是那么大...
AFAIK，由于您正在执行一次打开/关闭操作，因此它最多应该使用原始数据内存的 3 倍：原始 + 膨胀/腐蚀 + 最终结果......
您可以尝试自己手动实现...它可能会慢一些，但是代码很简单并且应该可以深入了解问题...