Add float coordinates to numpy array

Question

I would like to add float coordinates to a numpy array by splitting the intensity based on the centre of mass of the coordinate to neighbouring pixels.

As an example with integers:

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.]])

However I would like to distribute the intensity across neighbouring pixels when coord is float data, eg. coord = [2.2, 1.7].

I have considered using a gaussian, eg:

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)

which gives good results, but becomes slow for 3d data and thousands of points.

Any advice or ideas would be appreciated, thanks.

Based on @rpoleski suggestion, take a local region and apply weighting by distance. This is a good idea although the implementation I have does not maintain the original centre of mass of the coordinates, for example:

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

Any ideas?

Yes the idea I had does work but it's too slow for practical use. It may be quicker to slice the array and only evaluate the gaussian for a small slice, maybe 3 or 4 sigma, around the point. I was wondering whether anyone else had come across this problem, or whether it has a good solution, as I imagine it has been run into a many times but I couldn't find anything in my search. — Paddy Harrison
– Paddy Harrison, Commented Apr 27, 2020 at 14:36
You don't say that gaussian is required. If it's not, then you can distribute intensity only to the 4 neighbouring: arr[int(c[0]), int(c[1])], arr[int(c[0])+1, int(c[1])] etc. with values being proportional to distance from c. In your solution, I think your code is slow because you distribute the signal over whole array and np.exp() calculation is slow. Take only nearby points and calculate -distance_from_coords**2/2 and use that as argument for np.exp(). — rpoleski
– rpoleski, Commented Apr 30, 2020 at 6:48
@gnodab yes the Euclidean distance doesn't maintain CoM of the input coordinates, but, as @rpoleski showed, using the taxicab or Manhattan distance does. There is also the cityblock function in the scipy library with this functionality for anyone interested. I didn't know about this before and it's always good to pass along new information: docs.scipy.org/doc/scipy/reference/generated/… — Paddy Harrison
– Paddy Harrison, Commented May 18, 2020 at 11:18

rpoleski · Accepted Answer · 2020-04-30 15:54:46Z

2

+50

Here is a simple (and hence most probably fast enough) solution that keeps the center of mass and has sum = 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))

which results in:

[[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)

Note: I'm using taxicab distances - they're good for center of mass calculations.

edited Apr 30, 2020 at 15:54

answered Apr 30, 2020 at 15:36

rpoleski

9986 silver badges12 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Paddy Harrison Over a year ago

Thanks for the answer- your solution works great. It appear that the problem I was having was because of using the Euclidean distance. Your taxicab distances work much better, I was unfamiliar with them before. Is this a well known solution?

rpoleski Over a year ago

I'm glad I could help. Yes, cab drivers know it :) Seriously, it fulfils all requirements for being mathematical distance (e.g. triangle inequality), so you can use it for many applications. It has some funny properties, e.g., taxicab circle looks like a square. I'm guessing extending my answer to gaussian should be easy.

Paddy Harrison · Accepted Answer · 2021-03-05 17:57:29Z

For anyone needing this functionality, and thanks to @rpoleski, I came up with this which uses Numba to speed up the calculation.

@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()
        )

Some testing:

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)

In fact it is better to compare to compare to a buffered function as the first test will overwrite any previous values instead of accumulating:

%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)

Collectives™ on Stack Overflow

Add float coordinates to numpy array

2 Answers 2

2 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related