NumPy: Execute function over each ndarray element

arrays, multidimensional-array, numpy, python

Solution

List comprehensions are a very inefficient way to deal with numpy arrays. They're an especially poor choice for the distance calculation.

To find the difference between your data and a point, you'd just do `data - point`. You can then calculate the distance using `np.hypot`, or if you'd prefer, square it, sum it, and take the square root.

It's a bit easier if you make it an Nx2 array for the purposes of the calculation though.

Basically, you want something like this:

import numpy as np

data = np.array([[[1704, 1240],
                  [1745, 1244],
                  [1972, 1290],
                  [2129, 1395],
                  [1989, 1332]],

                 [[1712, 1246],
                  [1750, 1246],
                  [1964, 1286],
                  [2138, 1399],
                  [1989, 1333]],

                 [[1721, 1249],
                  [1756, 1249],
                  [1955, 1283],
                  [2145, 1399],
                  [1990, 1333]]])

point = [1989, 1332]

#-- Calculate distance ------------
# The reshape is to make it a single, Nx2 array to make calling `hypot` easier
dist = data.reshape((-1,2)) - point
dist = np.hypot(*dist.T)

# We can then reshape it back to AxBx1 array, similar to the original shape
dist = dist.reshape(data.shape[0], data.shape[1], 1)
print dist

This yields:

array([[[ 299.48121811],
        [ 259.38388539],
        [  45.31004304],
        [ 153.5219854 ],
        [   0.        ]],

       [[ 290.04310025],
        [ 254.0019685 ],
        [  52.35456045],
        [ 163.37074401],
        [   1.        ]],

       [[ 280.55837182],
        [ 247.34186868],
        [  59.6405902 ],
        [ 169.77926846],
        [   1.41421356]]])

Now, removing the closest element is a bit harder than simply getting the closest element.

With numpy, you can use boolean indexing to do this fairly easily.

However, you'll need to worry a bit about the alignment of your axes.

The key is to understand that numpy "broadcasts" operations along the last axis. In this case, we want to brodcast along the middle axis.

Also, `-1` can be used as a placeholder for the size of an axis. Numpy will calculate the permissible size when `-1` is put in as the size of an axis.

What we'd need to do would look a bit like this:

#-- Remove closest point ---------------------
mask = np.squeeze(dist) != dist.min(axis=1)
filtered = data[mask]

# Once again, let's reshape things back to the original shape...
filtered = filtered.reshape(data.shape[0], -1, data.shape[2])

You could make that a single line, I'm just breaking it down for readability. The key is that `dist != something` yields a boolean array which you can then use to index the original array.

So, Putting it all together:

import numpy as np

data = np.array([[[1704, 1240],
                  [1745, 1244],
                  [1972, 1290],
                  [2129, 1395],
                  [1989, 1332]],

                 [[1712, 1246],
                  [1750, 1246],
                  [1964, 1286],
                  [2138, 1399],
                  [1989, 1333]],

                 [[1721, 1249],
                  [1756, 1249],
                  [1955, 1283],
                  [2145, 1399],
                  [1990, 1333]]])

point = [1989, 1332]

#-- Calculate distance ------------
# The reshape is to make it a single, Nx2 array to make calling `hypot` easier
dist = data.reshape((-1,2)) - point
dist = np.hypot(*dist.T)

# We can then reshape it back to AxBx1 array, similar to the original shape
dist = dist.reshape(data.shape[0], data.shape[1], 1)

#-- Remove closest point ---------------------
mask = np.squeeze(dist) != dist.min(axis=1)
filtered = data[mask]

# Once again, let's reshape things back to the original shape...
filtered = filtered.reshape(data.shape[0], -1, data.shape[2])

print filtered

Yields:

array([[[1704, 1240],
        [1745, 1244],
        [1972, 1290],
        [2129, 1395]],

       [[1712, 1246],
        [1750, 1246],
        [1964, 1286],
        [2138, 1399]],

       [[1721, 1249],
        [1756, 1249],
        [1955, 1283],
        [2145, 1399]]])

On a side note, if more than one point is equally close, this won't work. Numpy arrays have to have the same number of elements along each dimension, so you'll need to re-do your grouping in that case.

Problem

Original source