如何对 dask 数组执行“窗口化”操作答案

【问题标题】：How to perform a "windowed" operation on a dask array如何对 dask 数组执行“窗口化”操作
【发布时间】：2021-01-08 04:34:09
【问题描述】：

我有一个具有 3 个维度（时间、x 和 y）的 Xarray - 它们基本上是一堆图像。我想以“窗口”方式对图像对进行操作。

基本上取一对，比如前两个“时间”图像，然后在它们上方的5x5 窗口上调用一个函数。在 numpy 中，我会通过简单地切片并通过首先形成“补丁”来计算我的指标来做到这一点：

params = []
for i in range(0,patch1.shape[0],1):
    for j in range(0,patch1.shape[1],1):
        window1 = np.copy(imga[i:i+N,j:j+N]).flatten()
        window2 = np.copy(imgb[i:i+N,j:j+N]).flatten()
        params.append((window1, window2))

这里N 是窗口大小，例如5。然后计算我的指标：

def f(param):
    return metric_function(*param)

with Pool(4) as p:
    r = list(tqdm.tqdm(p.imap(f,params), total=len(params)))

但是，我很难将它翻译成 dask，我需要一些帮助。我的第一个直觉是使用map_overlap 函数，但我认为我并不完全理解如何使用，特别是因为输出与输入的维度不同；即输出只会是整个 NxN 补丁的中心像素。

【问题讨论】：

听起来你想做的是一个模板，stackoverflow.com/questions/40117237/… 可能有帮助？
接近了！并且绝对可以帮助我取得一些进展。我想通过实验，我设法想出了一个解决方案。让我在这里发布作为答案，等待更多回复！
也可以考虑map_overlapdocs.dask.org/en/latest/…
@mdurant 结合使用 map_overlap 和 skimage.util.shape.view_as_windows（并循环遍历块），我能够以相当快的方式实现我的目标。我将在下面更新我的答案以纠正这一点！

标签： python dask python-xarray dask-distributed scientific-computing

【解决方案1】：

通过一些实验，我想出了一种方法来执行此操作。

我将我的距离度量重新编写为与dask 兼容（即，将范围参数添加到直方图，因为 dask.array.histogram 期望这一点，而 numpy 没有）。这是距离度量：

def cauchy_schwartz(chunk):
    
    imga = chunk[0]
    imgb = chunk[1]
    
    p, _ = np.histogram(np.ravel(imga), bins=20, range=[imga.min(),imgb.max()])
    p = p/np.sum(p)
    q, _ = np.histogram(np.ravel(imgb), bins=20, range=[imgb.min(),imgb.max()])
    q = q/np.sum(q)

    n_d = np.array(da.sum(p * q)) 
    d_d = np.sqrt(np.sum(np.power(p, 2)) * np.sum(np.power(q, 2)))
    return np.array([-1.0 * np.log10( n_d/ d_d)])[None, None, None]

这里的关键是返回一个形状为(1,1,1) 的数组，后面会解释。

然后我简单地将我的数据分块到所需的窗口大小（在本例中为 `9）。

dcube2 = dcube.chunk((2,9,9)).persist()

我能够通过 9x9 窗口进行处理：

output = da.map_blocks(cauchy_schwartz, dcube2.data, chunks=(1,1,1), dtype='float64').compute()

我正在尝试使用map_blocks 和depth 参数以重叠的方式执行此操作。

这可能不是最有效的解决方案，但它是我能想到的最好的解决方案。

【讨论】：