Scipy ndimage形态运算符会使我的计算机内存RAM饱和(8GB) [英] Scipy ndimage morphology operators saturate my computer memory RAM (8GB)

问题描述

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

I need to compute morphological opening for 3D array of shape (400,401,401), size 64320400 bytes using a 3D structure element with a radius of 17 or greater. The size of structure element ndarray is 42875 bytes. Using scipy.ndimage.morphology.binary_opening, the whole process consumes 8GB RAM.

我已经在GitHub上阅读了scipy/ndimage/morphology.py，据我所知，形态侵蚀运算符是用纯C实现的.我很难理解ni_morphology.c的源代码，所以我没有找到该代码的任何部分都会导致如此巨大的内存利用率.添加更多的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.

重现该问题:

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.

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语言实现的，因此修改它们是一个困难的前景.

I too do openings of increasing radius for granulometry, and I ran into this same problem. In fact, the memory usage increases as roughly R^6 where R is the radius of the spherical kernel. That's quite a rate of increase! I did some memory profiling, including splitting the opening into an erosion and then a dilation (the definition of opening), and found that the large memory usage comes from SciPy's binaries and is cleared as soon as the result is returned to the calling Python script. SciPy's morphology codes are mostly implemented in C, so modifying them is a difficult prospect.

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

Anyway the OP's last comment: "After some researche I turned to Opening implementation using convolution -> multiplication of Fourier transforms - O(n log n), and no so big memory overhead." helped me figure out the solution, so thanks for that. The implementation however was not obvious at first. For anyone else who happens upon this problem I am going to post the implementation here.

我将开始讨论膨胀，因为二进制腐蚀只是二进制图像的补码(逆)的膨胀，然后将结果求逆.

I will start talking about dilation, because binary erosion is just the dilation of the complement (inverse) of a binary image, and then the result is inverted.

简而言之:根据由Kosheleva等人撰写的白皮书，膨胀可以看作是数据集A与结构化元素(球形核)B的卷积，阈值高于一定的价值.卷积也可以在频率空间中进行(通常快得多)，因为频率空间中的乘法与实际空间中的卷积相同.因此，通过先对A和B进行傅立叶变换，然后将它们相乘，然后对结果进行逆变换，然后对那个进行阈值大于0.5的阈值运算，即可得出A与B的关系.我链接的白皮书表示阈值大于0，但是大量测试表明该错误结果带有许多伪影；

In short: according to this white paper by Kosheleva et al, dilation can be viewed as a convolution of the dataset A with the structuring element (spherical kernel) B, thresholded above a certain value. Convolutions can also be done (often much faster) in frequency space, since a multiplication in frequency space is the same as convolution in real space. So by taking the Fourier transform of A and B first, multiplying them, and then inverse-transforming the result, and then thresholding that for values above 0.5, you get the dilation of A with B. (Note that the white paper I linked says to threshold above 0, but much testing showed that that gave wrong results with many artifacts; another white paper by Kukal et al. gives the threshold value as >0.5, and that gave identical results as scipy.ndimage.binary_dilation for me. I'm not sure why the discrepancy, and I wonder if I missed some detail of ref 1's nomenclature)

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

Proper implementation of that involves padding for size, but luckily for us, it's already been done in scipy.signal.fftconvolve(A,B,'same') - this function does what I just described and takes care of padding for you. Giving the third option as 'same' will return a result the same size as A, which is what we want (otherwise it will be padded out by the size of 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. (请注意，代码可以用更少的行完成，但是我试图在此处进行说明.)

Erosion in principal is this: you invert A, dilate it by B as above, and then re-invert the result. But it requires a slight trick to match exactly the results from scipy.ndimage.binary_erosion - you must pad the inversion with 1s out to at least the radius R of the spherical kernel B. So erosion can be implemented thusly to get identical results to scipy.ndimage.binary_erosion. (Note that the code could be done in fewer lines but I'm trying to be illustrative here.)

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])

您可以通过另一种方式获得相同的侵蚀结果，如

You can get identical erosion results another way, as shown in the white paper by Kukal et al - they point out that the convolution of A and B can be made into an erosion by thresholding by > m-0.5 , where m is the "size" of B (which turns out to be the volume of the sphere, not the volume of the array). I showed erode_v1 first because it's slightly easier to understand, but the results are the same here:

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

我希望这对其他有此问题的人有所帮助.关于我得到的结果的注释:

I hope this helps anyone else having this problem. Notes about the results I got:

• 我在2D和3D中都对此进行了测试，所有结果均与scipy.ndimage形态学运算(以及skimage运算，在后端仅称为ndimage)得到的相同答案相同.
• 对于我最大的内核(R = 21)，内存使用量减少了30倍！速度也快了20倍.
• 尽管我只在二进制图像上进行了测试-我只是不了解灰度，但是下面的第二个参考文献对此进行了一些讨论.

另外两个快速注释:

首先:考虑我在中间部分讨论的关于erode_v1的填充.用1填充反数基本上可以使腐蚀从数据集的边缘以及数据集中的任何界面发生.根据您的系统和要执行的操作，您可能要考虑这是否真正代表了您希望其处理的方式.如果不是，则可以考虑使用反射"边界条件进行填充，该条件将模拟边缘附近所有特征的延续.我建议您在不同的边界条件(包括膨胀和侵蚀)上进行试验，并对结果进行可视化和量化，以确定哪种最适合您的系统和目标.

First: Consider the padding I discuss in the middle section about erode_v1. Padding the inverse out with 1s basically allows erosion to occur from the edges of the dataset as well as from any interface in the dataset. Depending on your system and what you are trying to do, you may want to consider whether or not this truly represents the way you want it handled. If not, you might consider padding out with the 'reflect' boundary condition, which would simulate a continuation of any features near the edge. I recommend playing around with different boundary conditions (on both dilation and erosion) and visualizing and quantifying the results to determine what suits your system and goals the best.

第二:大多数情况下，这种基于频率的方法不仅在内存方面而且在速度方面都更好.对于小内核B，原始方法更快.但是，无论如何，小内核都可以非常快速地运行，因此出于我自己的目的，我不在乎.如果这样做(就像您多次执行小内核一样)，则可能需要找到B的临界大小并在此时切换方法.

Second: This frequency-based method is not only better in memory but also in speed - for the most part. For small kernels B, the original method is faster. However, small kernels run very quickly anyway, so for my own purposes I don't care. If you do (like if you are doing a small kernel many times), you may want to find the critical size of B and switch methods at that point.

参考文献，尽管我很抱歉，因为它们都不提供年份，所以不容易被引用:

References, though I apologize that they are not easy to cite as they provide neither year:

O. Kosheleva，S.D. Cabrera，G.A. Gibson，M.Koshelev使用
1. 使用快速傅立叶变换快速实现形态学运算. http://www.cs.utep.edu/vladik/misha5.pdf
J. Kukal，D.Majerova和A. Prochazka的
2. 带球形罩的灰度图像的扩张和腐蚀. http://http％3A％2F％2Fwww2. humusoft.cz％2Fwww％2Fpapers％2Ftcp07％2F001_kukal.pdf

1. Fast Implementation of Morphological Operations Using Fast Fourier Transform by O. Kosheleva, S. D. Cabrera, G. A. Gibson, M. Koshelev. http://www.cs.utep.edu/vladik/misha5.pdf
2. Dilation and Erosion of Gray Images with Spherical Masks by J. Kukal, D. Majerova, A. Prochazka. http://http%3A%2F%2Fwww2.humusoft.cz%2Fwww%2Fpapers%2Ftcp07%2F001_kukal.pdf

这篇关于Scipy ndimage形态运算符会使我的计算机内存RAM饱和(8GB)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！