使用 dask 进行 3D 体积处理

Question

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

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

现在是最重要的部分：

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

你觉得-可能吗？

为了用图像来解释它——请将其视为 3D 域分解（这是来自 DWT——但有助于说明 from here）：

illistration of domain decomposition

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

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

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

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

Answer 1

我听到几个问题。我会一一回复

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

是的。 Lazy Dask 操作产生一个依赖图。在 dask.arrays 的情况下，该图是按块计算的。如果你的输出只依赖于图的一个子集，那么 Dask 将删除不需要的任务。这方面的深入文档特别是 here and the 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 代码那样高效。不过，我会把性能问题留给你。

Answer 2

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

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

至于在 Dask 中使用 Sobel 运算符，可以使用此 sobel function from dask-ndfilters（许可许可）在 Dask 阵列上执行此操作。它处理引擎盖下块上的适当光环，返回一个新的 Dask 阵列。

由于 SciPy 的 sobel 函数（以及 dask-ndfilters 的 sobel 也是）在一维上运行，因此需要映射每个轴并堆栈以获得完整的 Sobel 算子结果。也就是说，这很简单。下面是一个简短的片段，展示了如何在随机 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()

使用 dask 进行 3D 体积处理

3D volume processing using dask

python

filtering

convolution

multidimensional-array

dask