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

如何解决如何对 dask 数组执行“窗口化”操作

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

基本上取一对，说前两个“时间”图像，然后在它们的 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 块的中心像素。

解决方法

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

我重新编写了我的距离度量以使其与 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),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),dtype='float64').compute()

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

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