我也为粒度测量做了增加半径的开口,我遇到了同样的问题。事实上,内存使用量大约增加了 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 的临界大小并在那个时候切换方法。
参考文献,尽管我很抱歉它们不容易引用,因为它们都没有提供年份:
- O. Kosheleva、SD Cabrera、GA Gibson、M. Koshelev使用快速傅里叶变换快速实现形态学运算。http://www.cs.utep.edu/vladik/misha5.pdf
- J. Kukal, D. Majerova, A. Prochazka用球形面具对灰色图像进行膨胀和侵蚀。http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf