将浮点坐标添加到 numpy 数组

Question

我想通过将基于坐标质心的强度拆分到相邻像素来将浮点坐标添加到 numpy 数组。

以整数为例：

import numpy as np

arr = np.zeros((5, 5), dtype=float)

coord = [2, 2]
arr[coord[0], coord[1]] = 1

arr
>>> array([[0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0.],
           [0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0.]])

但是，当coord 是浮点数据时，我想在相邻像素之间分配强度，例如。 coord = [2.2, 1.7].

我考虑过使用高斯，例如：

grid = np.meshgrid(*[np.arange(i) for i in arr.shape], indexing='ij')

out = np.exp(-np.dstack([(grid[i]-c)**2 for i, c in enumerate(coord)]).sum(axis=-1) / 0.5**2)

效果很好，但对于 3d 数据和数千个点来说会变得很慢。

任何建议或想法将不胜感激，谢谢。

根据@rpoleski 的建议，取一个局部区域并按距离应用权重。这是一个好主意，虽然我的实现没有保持坐标的原始质心，例如：

from scipy.ndimage import center_of_mass

coord = [2.2, 1.7]

# get region coords
grid = np.meshgrid(*[range(2) for i in coord], indexing='ij')
# difference Euclidean distance between coords and coord
delta = np.linalg.norm(np.dstack([g-(c%1) for g, c, in zip(grid, coord)]), axis=-1)

value = 3 # pixel value of original coord
# create final array by 1/delta, ie. closer is weighted more
# normalise by sum of 1/delta
out = value * (1/delta) / (1/delta).sum()

out.sum()
>>> 3.0 # as expected

# but
center_of_mass(out)
>>> (0.34, 0.63) # should be (0.2, 0.7) in this case, ie. from coord

有什么想法吗？

Answer 1

这是一个简单（因此很可能足够快）的解决方案，它保持质心并且总和 = 1：

arr = np.zeros((5, 5), dtype=float)

coord = [2.2, 0.7]

indexes = np.array([[x, y] for x in [int(coord[0]), int(coord[0])+1] for y in [int(coord[1]), int(coord[1])+1]])
values = [1. / (abs(coord[0]-index[0]) * abs(coord[1]-index[1])) for index in indexes]
sum_values = sum(values)
for (value, index) in zip(values, indexes):
    arr[index[0], index[1]] = value / sum_values
print(arr)
print(center_of_mass(arr))

导致：

[[0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.  ]
 [0.   0.24 0.56 0.   0.  ]
 [0.   0.06 0.14 0.   0.  ]
 [0.   0.   0.   0.   0.  ]]
(2.2, 1.7)

注意：我使用的是出租车距离 - 它们适用于质心计算。

Answer 2

对于任何需要此功能的人，感谢@rpoleski，我想出了这个，它使用Numba 来加快计算速度。

@numba.njit
def _add_floats_to_array_2d(coords, arr, values):
    """

    Distribute float values around neighbouring pixels in array whilst maintinating center of mass.
    Uses Manhattan (taxicab) distances for center of mass calculation.

    This function uses numba to speed up the calculation but is limited to exactly 2D.

    Parameters
    ----------
    coords: (N, ndim) ndarray
        Floats to distribute into array.
    arr: ndim ndarray
        Floats will be distributed into this array.
        Array is modified in place.
    values: (N,) arraylike
        The total value of each coord to distribute into arr.

    """
    indices_local = np.array([[[0, 0], [1, 0]], [[0, 1], [1, 1]]])

    for i, c in enumerate(coords):
        temp_abs = np.abs(indices_local - np.remainder(c, 1))
        temp = 1.0 / (temp_abs[..., 0] * temp_abs[..., 1])
        # handle perfect integers
        for j in range(temp.shape[0]):
            for k in range(temp.shape[1]):
                if np.isinf(temp[j, k]):
                    temp[j, k] = 0
        arr[int(c[0]) : int(c[0]) + 2, int(c[1]) : int(c[1]) + 2] += (
            values[i] * temp / temp.sum()
        )

一些测试：

arr = np.zeros((256, 256))

coords = np.random.rand(10000, 2) * arr.shape[0]
values = np.ones(len(coords))

%timeit arr[tuple(coords.astype(int).T)] = values
>>> 106 µs ± 4.08 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit _add_floats_to_array_2d(coords, arr, values)
>>> 13.5 ms ± 546 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

事实上，最好将其与缓冲函数进行比较，因为第一个测试将覆盖任何先前的值而不是累加：

%timeit np.add.at(arr, tuple(coords.astype(int).T), values)
>>> 1.23 ms ± 178 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)