像素化位置和通量的算法

Algorithm to pixelate positions and fluxes

如果你在一个平面上有 n 个粒子(位置为 (x_n,y_n)),具有一定的通量流感x_n,你会怎么做?像素化这些粒子,所以你必须从 (x,y) 到 (pixel_i, pixel_j) space 并且你必须总结落入的 m 个粒子的通量每个像素?有什么建议么?谢谢!

有几种方法可以解决您的问题。

假设:您的位置已存储为两个 numpy array 形状 (N, ),即 n 的位置 x_n(或 y_n)在 [0, N),我们称它们为 xy。通量存储在一个numpy array中,形状相同,fluxes.

1 - 强化病例 创建看起来像网格的东西:

#get minimums and maximums position
mins = int(x.min()), int(y.min())
maxs = int(x.max()), int(y.max())
#actually you can also add and subtract 1 or more unit 
#in order to have a grid larger than the x, y extremes
#something like mins-=epsilon and maxs += epsilon

#create the grid
xx = np.arange(mins[0], maxs[0])
yy = np.arange(mins[1], maxs[1])

现在你可以执行双for loop,每次添加xxyy的两个连续元素,要做到这一点,你可以简单地采取:

x1 = xx[:-1] #excluding the last element
x2 = xx[1:]  #excluding the first element

#the same for y:
y1 = yy[:-1] #excluding the last element
y2 = yy[1:]  #excluding the first element

fluxes_grid = np.zeros((xx.shape[0], yy.shape[0]))
for i, (x1_i, x2_i) in enumerate(zip(x1, x2)):
    for j, (y1_j, y2_j) in enumerate(zip(y1, y2)):
        idx = np.where((x>=x1_i) & (x<x2_i) & (y>=y1_j) & (y<y2_j))[0]
        fluxes_grid[i,j] = np.sum(fluxes[idx])

在此循环结束时,您将得到一个网格,其元素是表示通量总和的像素。

2 - 使用像 K-NN 这样的量化算法

如果您有很多 o 点,以至于循环需要数小时,会发生什么情况? 更快的解决方案是使用量化方法,例如 K 最近邻,刚性网格上的 KNN。 运行 KNN 有很多方法(包括已经实现的版本,例如 sklearn KNN)。但如果您可以利用 GPU,则效率会有所不同。例如,这是我的 tensorflow (vs 2.1) 实现。定义方形网格后:

_min, maxs = min(mins), max(maxs)
xx = np.arange(_min, _max)
yy = np.arange(_min, _max)

您可以构建矩阵 grid 和位置矩阵 X: 网格 = np.column_stack([xx, yy]) X = np.column_stack([x, y])

那么你必须定义一个矩阵欧氏成对距离函数:

@tf.function
def pairwise_dist(A, B):  
    # squared norms of each row in A and B
    na = tf.reduce_sum(tf.square(A), 1)
    nb = tf.reduce_sum(tf.square(B), 1)
    # na as a row and nb as a co"lumn vectors
    na = tf.reshape(na, [-1, 1])
    nb = tf.reshape(nb, [1, -1])
    # return pairwise euclidead difference matrix
    D = tf.sqrt(tf.maximum(na - 2*tf.matmul(A, B, False, True) + nb, 0.0))
    return D

因此:

#compute the pairwise distances:
D = pairwise_dist(grid, X)
D = D.numpy() #get a numpy matrix from a tf tensor
#D has shape M, N, where M is the number of points in the grid and N the number of positions.
#now take a rank and from this the best K (e.g. 10)
ranks = np.argsort(D, axis=1)[:, :10]
#for each point in the grid you have the nearest ten.

现在你必须取这 10 个位置对应的通量并求和。

我已经避免进一步指定第二种方法,我不知道你目录的维度,如果你有或没有 GPU 或者你是否想使用这种优化。 如果你愿意,我可以改进这个解释,前提是你有兴趣。