如何在dask阵列上执行“加窗”操作

Question

我有一个具有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函数，但我不认为我完全理解如何使用，特别是因为输出与输入的维度不同；i、 e.输出仅为整个NxN面片的中心像素

Answer 1

通过一点实验，我想出了一种执行此操作的方法

我重新编写了距离度量，使其与dask兼容（即，在直方图中添加了range参数，因为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参数以重叠的方式进行实验

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