如何在dask数组中进行径向求和？

Question

我正在尝试对 dask 数组的值进行径向求和，在该数组中我保留了分块数据，并对每个半径求和。将总和标准化为“像素”总和可能很有用。

这是dask数组的结构： dask.array 其中，shape=(264, 256, 1500)，dtype=int16，chunksize=(264, 256, 992)，chunktype=numpy.ndarray 数据包含：强度（计数）和维度： x: 尺寸 (264,) 的距离 (nm) y: 尺寸 (256,) 的距离 (nm) energy_scale: 能量 (eV) 大小 (1500,)

我尝试过用这段代码来掩盖一个圆并创建一个新数组，它效果很好，但需要径向求和。

import numpy as np

def cartesian_to_polar(x, y, center):

    x_rel = x - center[1]
    y_rel = y - center[0]
    radius = np.sqrt(x_rel**2 + y_rel**2)
    theta = np.arctan2(y_rel, x_rel)
    return radius, theta

def get_polar_coordinates_grid(center, max_radius, dimensions):
    x = np.arange(dimensions[1])
    y = np.arange(dimensions[0])
    xv, yv = np.meshgrid(x, y)
    rv, thetav = cartesian_to_polar(xv, yv, center)
    
    # Mask to select points within the specified radius
    mask = rv <= max_radius
    
    mask_3D = np.repeat(mask[:, :, np.newaxis], dimensions[2], axis=2)
    
    print('mask_3D shape',mask_3D.shape)
    
    return mask_3D, rv, thetav

def polar_profile(dataset, center, max_radius):
    print('center', center[0],center[1])

    mask_3D, rv, thetav = get_polar_coordinates_grid(center, max_radius, dataset.shape[:3])
    radial_dataset_array = np.repeat(rv[:, :, np.newaxis], dataset.shape[2], axis=2)

    masked_dataset = np.where(mask_3D, dataset, 0)

    masked_r_dataset = np.where(mask_3D, radial_dataset_array, 0)

从这里我不确定索引需要求和的像素的最佳路线是什么。一个简单的 numpy 数组的例子也会很有帮助。

Answer 1

实际上比这简单得多。一种方法是这样做（这里我创建了一个示例数据，将其替换为您自己的数据）：

import dask.array as da
import numpy as np


dask_data = da.random.random((264, 256, 1500), chunks=(264, 256, 992))

center = (dask_data.shape[0] // 2, dask_data.shape[1] // 2)
max_radius = min(center)

x = da.arange(dask_data.shape[0])
y = da.arange(dask_data.shape[1])
xv, yv = da.meshgrid(x - center[0], y - center[1], indexing='ij')
radii = da.sqrt(xv**2 + yv**2).astype(np.int32)

def sum_radial(image, radii, max_radius):
    results = []
    for r in range(max_radius+1):
        mask = radii == r
        masked_data = da.where(mask[..., None], image, 0)
        sum_data = masked_data.sum(axis=(0, 1))
        count = mask.sum()
        results.append((sum_data, count))
    return results

results = sum_radial(dask_data, radii, max_radius)
final_results = [(res[0] / res[1]).compute() for res in results if res[1] > 0]

返回

[array([0.56819324, 0.48555518, 0.07560507, ..., 0.57646385, 0.54231102,
       0.24472355]), array([0.69557989, 0.55792832, 0.37600816, ..., 0.59541941, 0.47944983,
       0.44399378]), array([0.54944245, 0.37145509, 0.57202096, ..., 0.5346811 , 0.46998448,
       0.48209734]), array([0.45247991, 0.48744803, 0.56118335, ..., 0.40804105, 0.57898455,
       0.48248958]), array([0.45528472, 0.43234064, 0.50973776, ..., 0.43357677, 0.54758867,
       0.48335612]), array([0.50817523, 0.54346835, 0.48714012, ..., 0.46681618, 0.58126478,
       0.41308635]), array([0.4513112 , 0.58265886, 0.51799577, ..., 0.50165861, 0.52393633,
       0.52736444]), array([0.47323004, 0.48509912, 0.50154742, ..., 0.44656244, 0.45017714,
       0.45944342]), array([0.51737857, 0.43069107, 0.50277041, ..., 0.50754802, 0.45124549,
       0.48241854]), array([0.47561275, 0.47823851, 0.45425869, ..., 0.48548346, 0.46774252,
       0.53712027]), array([0.45032079, 0.46940781, 0.54524277, ..., 0.54348768, 0.47486265,
       0.57205173]), array([0.41538367, 0.57499891, 0.53469115, ..., 0.49469812, 0.52803723,
       0.45382559]), array([0.46701844, 0.52337592, 0.52955901, ..., 0.50932645, 0.5619902 ,
       0.5217232 ]), array([0.51389081, 0.45033721, 0.52085882, ..., 0.46968577, 0.44170561,
       0.48609703]), array([0.52183791, 0.48536374, 0.46433651, ..., 0.47977694, 0.51830913,
...
       0.50629353]), array([0.50827807, 0.50630572, 0.48872231, ..., 0.50363937, 0.49413712,
       0.5133742 ]), array([0.51127133, 0.51176922, 0.48647271, ..., 0.49598711, 0.50950352,
       0.51000234]), array([0.51796942, 0.46370582, 0.48640914, ..., 0.51021386, 0.5104863 ,
       0.49745391])]