【发布时间】：2022-12-05 10:53:01
【问题描述】：
首先，我以前从未问过关于 stackoverflow 的问题，我会尽力遵守网站指南，但如果我应该更改我的帖子，请告诉我。
我正在尝试编写一个可以从二进制 3D 图像中快速提取孔径分布的函数。我通过计算图像的局部厚度来做到这一点，其方式与在 ImageJ 的局部厚度插件中实现的方式类似。我理想地需要此函数在 1 秒内运行，因为我在模拟退火过程中调用它约 200000 次。它部分在 CPU（第 12 代 Intel(R) Core(TM) i7-12700KF、20 核、16GB RAM）上执行，部分在 GPU（RTX GeForce 3050、8GB）上执行。
该功能有效，但发生了一些事情我认为在后端，这是人为地减慢它的速度。这可能与线程、GPU 到 CPU 开销或某种“冷却”期有关。
函数分为三部分：
-
欧氏距离变换 - 在 CPU 上执行，使用 edt 包并行执行。目前在 250^3 二进制图像上需要大约 0.25 秒
-
3d 骨架化 - 使用 skimage.morphology.skeletonize_3d 在 CPU 上执行，但使用 dask 将图像分割成块。此实现由 porespy.filters.chunked_func 提供。将骨架乘以距离变换以获得骨架，其值等于到最近背景体素的最小距离。这个过程需要 0.45 到 0.5 秒。
-
使用半径等于骨架体素值的球形结构元素扩大骨架上的每个体素。这是在 for 循环中完成的，从最大结构元素大小开始，并按降序排列。较大的球体不会被较小的球体覆盖。使用 cupyx.scipy.signal.signaltools.convolve 在 GPU 上使用 fft 卷积完成膨胀，这需要大约 0.005 秒。
然而，重现我所看到的效果需要更少的代码。关键部分是按顺序执行许多 fft 卷积。
一个最小的可重现示例如下：
import skimage import time import cupy as cp from cupyx.scipy.signal.signaltools import convolve # Generate a binary image im = cp.random.random((250,250,250)) > 0.4 # Generate spherical structuring kernels for input to convolution structuring_kernels = {} for r in range(1,21): structuring_kernels.update({r: cp.array(skimage.morphology.ball(r))}) # run dilation process in loop for i in range(10): s = time.perf_counter() for j in range(20,0,-1): convolve(im, structuring_kernels[j], mode='same', method='fft') e = time.perf_counter() # time.sleep(2) print(e-s)按原样运行时，在前几个循环之后，每个扩张循环在我的计算机上需要大约 1.8 秒。如果我取消注释
time.sleep(2)行（即在每个循环之间暂停 2 秒），那么循环函数调用只需要 0.05 秒。我怀疑这与线程或 GPU 使用有关，因为它需要几个循环才能达到 1.8 秒，然后它会保持稳定在该值。当我监控我的 GPU 使用情况时，3D 监视器迅速达到 100% 并保持接近该水平。如果我只是受到 GPU 容量的限制，为什么前几个循环运行得更快？会发生内存泄漏吗？有谁知道为什么会这样，如果有办法防止它，可能在 cupy 中使用后端控件？
我不确定这是否有必要，但我的整体局部厚度函数如下：
import porespy as ps from skimage.morphology import skeletonize_3d import time import numpy as np import cupy as cp from edt import edt from cupyx.scipy.signal.signaltools import convolve def local_thickness_cp(im, masks=None, method='fft'): """ Parameters ---------- im: 3D voxelized image for which the local thickness map is desired masks: (optional) A dictionary of the structuring elements to be used method: 'fft' or 'direct' Returns ------- The local thickness map """ s = time.perf_counter() # Calculate the euclidean distance transform using edt package dt = cp.array(edt(im, parallel=15)) e = time.perf_counter() # print(f'EDT took {e - s}') s = time.perf_counter() # Calculate the skeleton of the image and multiply by dt skel = cp.array(ps.filters.chunked_func(skeletonize_3d, overlap=17, divs=[2, 3, 3], cores=20, image=im).astype(bool)) * dt e = time.perf_counter() # print(f'skeletonization took {e - s} seconds') r_max = int(cp.max(skel)) s = time.perf_counter() if not masks: masks = {} for r in range(int(r_max), 0, -1): masks.update({r: cp.array(ps.tools.ps_ball(r))}) e = time.perf_counter() # print(f'mask creation took {e - s} seconds') # Initialize the local thickness image final = cp.zeros(cp.shape(skel)) time_in_loop = 0 s = time.perf_counter() for r in range(r_max, 0, -1): # Get a mask of where the skeleton has values between r-1 and r skel_selected = ((skel > r - 1) * (skel <= r)).astype(int) # Perform dilation on the mask using fft convolve method, and multiply by radius of pore size dilation = (convolve(skel_selected, masks[r], mode='same', method=method) > 0.1) * r # Add dilation to local thickness image, where it is still zero (ie don't overwrite previous inserted values) final = final + (final == 0) * dilation e = time.perf_counter() # print(f'Dilation loop took {e - s} seconds') return final现在，理论上，该函数应该需要大约 0.80 秒的时间来计算。但是，当在单独的图像上循环调用时，大约需要 1.5 秒。但是，如果我在每次函数调用后添加一个
time.sleep(1)，那么该函数确实需要大约 0.8 秒。
【问题讨论】：
标签： python fft convolution cupy slowdown