使用 dask 进行 3D 体积处理答案

【问题标题】：3D volume processing using dask使用 dask 进行 3D 体积处理
【发布时间】：2017-08-13 20:37:18
【问题描述】：

我现在正在使用 dask 使用一些简单的模板来探索 3D 交互式体积卷积。让我解释一下我的意思：

假设您有一个 3D 数据，您希望通过 Sobel 变换对其进行处理（例如获得 L1 或 L2 梯度）。
然后您将输入的 3D 图像划分为子体积（具有一些重叠边界 - 对于 3x3x3 模板 Sobel，它将需要 +2 个样本重叠/填充）
现在让我们假设您在整个 3D 体积上创建了 Sobel 3D 变换的延迟计算，但尚未执行它。

现在是最重要的部分：

我想编写一个函数，该函数将从虚拟转换的数据中提取某些特定的 2D 部分。
最后让 dask 计算所有内容：
- 但是我需要做的不是为我计算整个变换然后提供一个部分。
  - 我只需要它来执行计算特定 2D 转换图像切片所需的那些任务。

你认为这可能吗？

为了用图像来解释它——请认为这是一个 3D 域分解（这是来自 DWT——但很好的说明from here）：

illistration of domain decomposition

并假设有一个函数使用 dask 计算整个体积的 3D 变换。但我想得到——例如——转换后的 3D 数据的 2D 图像，它由 LLL1、LLH1、HLH1、HLL1 平面组成（基本上是一个切片） .

重要的部分不是计算整个子立方体——而是让 dask 以某种方式自动跟踪计算图中的依赖关系并仅评估那些。

请不要担心计算 vs.复制时间。假设它有完美的比例。

如果需要更多说明，请告诉我！感谢您的帮助！

【问题讨论】：

标签： python multidimensional-array filtering convolution dask

【解决方案1】：

不要减损@MRocklin 精心布置的答案，而是要添加更多内容。

我还经常发现自己需要对大规模阵列数据进行边缘检测和其他图像处理技术。由于 Dask 是一个非常好的库，用于在大型阵列数据上构建和探索此类计算工作流，因此在一个名为 dask-image 的 GitHub 组织中为一些常见的图像处理技术组合了一些实用程序库。它们的设计主要是为了模仿 SciPy 的 ndimage API。

关于在 Dask 中使用 Sobel 运算符，可以使用来自 dask-ndfilters（许可许可）的 sobel 函数在 Dask 阵列上执行此操作。它处理引擎盖下方的块上的正确光环，返回一个新的 Dask Array。

由于 SciPy 的 sobel 函数（以及 dask-ndfilters 的 sobel）在一维上运行，因此需要映射每个轴和堆栈以获得完整的 Sobel 运算符结果。也就是说，这很简单。下面是一个简短的 sn-p，展示了如何在随机 Dask Array 上执行此操作。还包括沿 XZ 平面进行切片。尽管人们可以轻松地获取任何其他切片或对数据执行其他操作。

希望这会有所帮助。 :)

import dask.array as da
import dask_ndfilters as da_ndfilt

d = da.random.random((100, 120, 140), chunks=(25, 30, 35))
ds = da.stack([da_ndfilt.sobel(d, axis=i) for i in range(d.ndim)])
dsp = ds[:, :, 0, :]
asp = dsp.compute()

【讨论】：

【解决方案2】：

我听到几个问题。我会逐一回答

Dask 能否跟踪输出子集所需的任务并仅计算这些任务？

是的。 Lazy Dask 操作生成依赖图。在 dask.arrays 的情况下，这个图是每块的。如果您的输出仅取决于图表的子集，则 Dask 将删除不必要的任务。对此的深入文档是 here 和 cull 特别优化。

以这个 100,000 x 100,000 数组为例

>>> x = da.random.random((100000, 100000), chunks=(1000, 1000))

假设我从中添加了几个 1d 切片

>>> y = x[5000, :] + x[:, 5000].T

生成的优化图仅足以计算输出

>>> graph = y._optimize(y.dask, y._keys())  # you don't need to do this
>>> len(graph)                              # it happens automatically
301

而且我们可以很快地计算出结果：

In [8]: %time y.compute()
CPU times: user 3.18 s, sys: 120 ms, total: 3.3 s
Wall time: 936 ms
Out[8]: 
array([ 1.59069994,  0.84731881,  1.86923216, ...,  0.45040813,
        0.86290539,  0.91143427])

现在，这并不完美。它确实必须生成我们两个切片接触的所有 1000x1000 块。但是你可以在那里控制粒度。

简答：Dask 将自动检查图表并仅运行评估输出所需的那些任务。你不需要做任何特别的事情来做到这一点。

使用 dask.array 进行重叠数组计算是个好主意吗？

也许吧。相关文档页面位于Overlapping Blocks with Ghost Cells。 Dask.array 具有方便的功能，可以轻松写下来。但是，它将创建内存中的副本。许多像你这样的人发现 memcopy 太慢了。 Dask 通常不支持就地计算，因此我们不能像正确的 MPI 代码那样高效。不过，我会把性能问题留给你。

【讨论】：