Apply 3-argument function to 3D numpy array

Question

I have a 3D numpy array A of shape (2133, 3, 3). Basically this is a list of 2133 lists with three 3D points. Furthermore I have a function which takes three 3D points and returns one 3D point, x = f(a, b, c), with a, b, c, x numpy arrays of length 3. Now I want to apply f to A, so that the output is an array of shape (2133, 3). So something like numpy.array([f(*A[0]),...,f(*A[2132])).

I tried numpy.apply_along_axis and numpy.vectorize without success.

To be more precise the function f I consider is given by:

def f(a, b, c, r1, r2=None, r3=None):
    a = np.asarray(a)
    b = np.asarray(b)
    c = np.asarray(c)

    if np.linalg.matrix_rank(np.matrix([a, b, c])) != 3:
        # raise ValueError('The points are not collinear.')
        return None

    a, b, c, = sort_triple(a, b, c)

    if any(r is None for r in (r2, r3)):
        r2, r3 = (r1, r1)

    ex = (b - a) / (np.linalg.norm(b - a))
    i = np.dot(ex, c - a)
    ey = (c - a - i*ex) / (np.linalg.norm(c - a - i*ex))
    ez = np.cross(ex, ey)
    d = np.linalg.norm(b - a)
    j = np.dot(ey, c - a)

    x = (pow(r1, 2) - pow(r2, 2) + pow(d, 2)) / (2 * d)
    y = ((pow(r1, 2) - pow(r3, 2) + pow(i, 2) + pow(j, 2)) / (2*j)) - ((i/j)*x)
    z_square = pow(r1, 2) - pow(x, 2) - pow(y, 2)
    if z_square >= 0:
        z = np.sqrt(z_square)
        intersection = a + x * ex + y*ey + z*ez
        return intersection

A = np.array([[[131.83, 25.2, 0.52], [131.51, 22.54, 0.52],[133.65, 23.65, 0.52]], [[13.02, 86.98, 0.52], [61.02, 87.12, 0.52],[129.05, 87.32, 0.52]]])

r1 = 1.7115

Seeing what you tried with np.vectorize and np.apply_along_axis would be helpful. Also, what exactly does the function f do? You might be able to replace it with a vectorized version. — user2699
– user2699, Commented Aug 24, 2017 at 14:17
@user2699 I provided the function f. I think the problem with np.apply_along_axis is that it is applied to a 1D-slice along the given axis. np.vectorize does not work because the unpacking is not done in the correct way. Even if I rewrite f, so that it takes an array of shape (3,3), how do I tell numpy to iterate over the first axis and take the 3x3 subarrays? — patrik
– patrik, Commented Aug 24, 2017 at 14:36
There is no automatic way to do that, you need to write f in a way that it works with arrays of points instead of single points. I think it is almost okay, although you'd need to change the first if, add an axis parameter to the np.linalg.norm calls and maybe a few things more (sort_triple, which I assume is another if your functions, might need to be adapted too). Btw, consider replacing every pow(x, 2) call with np.square(x). — javidcf
– javidcf, Commented Aug 24, 2017 at 14:37
Hi @jdehesa what exactly do you mean with "you need to write f in a way that it works with arrays of points instead of single points." Replacing a, b, c as arguments through p where p is an array of shape (3,3), and then unpacking p does not seem to work. Do you mean the shape of the input argument of f has to be of the same shape as A, at least hast three dimensions? Thanks for the hint withnp.square(x). — patrik
– patrik, Commented Aug 24, 2017 at 17:04
@patrik No, I meant the function must be designed to take three matrices with shape (2133, 3) (or (X, 3) for whatever X) and return directly the (2133, 3) result. For example, you'd need to change ex = (b - a) / (np.linalg.norm(b - a)) with ex = (b - a) / (np.linalg.norm(b - a, axis=0)), etc. Each statement must work with all the points at the same time. There is no general way to automate it (efficiently, that is; you can use loops or comprehension for convenience, like in the answer, if performance is not a big issue). — javidcf
– javidcf, Commented Aug 24, 2017 at 17:14

patrik · Accepted Answer · 2017-09-04 08:05:58Z

Thanks to the great help of @jdehesa I was able to produce an alternative solution to the one given by @hpaulj. I am not sure if this solution is the most elegant one but it worked so far. Comments are appreciated.

def sort_triple(a, b, c):
    pts = np.stack((a, b, c), axis=1)
    xSorted = pts[np.arange(pts.shape[0])[:, None], np.argsort(pts[:, :, 0])]
    orientation = np.cross(xSorted[:, 1] - xSorted[:, 0], xSorted[:, 2] -
                           xSorted[:, 0])[:, 2] >= 0
    xSorted_flipped = np.stack((xSorted[:, 0], xSorted[:, 2], xSorted[:, 1]),
                               axis=1)
    xSorted = np.where(orientation[:, np.newaxis, np.newaxis], xSorted,
                       xSorted_flipped)
    return map(np.squeeze, np.split(xSorted, 3, axis=1))


def f(A, r1, r2=None, r3=None):

    a, b, c = map(np.squeeze, np.split(A, 3, axis=1))
    a, b, c = sort_triple(a, b, c)

    if any(r is None for r in (r2, r3)):
        r2, r3 = (r1, r1)

    ex = (b - a) / (np.linalg.norm(b - a, axis=1))[:, np.newaxis]
    i = inner1d(ex, (c - a))
    ey = ((c - a - i[:, np.newaxis]*ex) /
          (np.linalg.norm(c - a - i[:, np.newaxis]*ex, axis=1))[:, np.newaxis])
    ez = np.cross(ex, ey)
    d = np.linalg.norm(b - a, axis=1)
    j = inner1d(ey, c - a)

    x = (np.square(r1) - np.square(r2) + np.square(d)) / (2 * d)
    y = ((np.square(r1) - np.square(r3) + np.square(i) + np.square(j)) / (2*j) -
         i/j*x)
    z_square = np.square(r1) - np.square(x) - np.square(y)
    mask = z_square < 0
    z_square[mask] *= 0
    z = np.sqrt(z_square)
    z[mask] = np.nan
    intersection = (a + x[:, np.newaxis] * ex + y[:, np.newaxis] * ey +
                    z[:, np.newaxis] * ez)
    return intersection

Probably the map parts in each function could be done better. Maybe also the excessive use of np.newaxis.

hpaulj · Accepted Answer · 2017-08-24 16:52:48Z

0

This works fine (after commenting out sort_triple):

res = [f(*row,r1) for row in A]
print(res)

producing:

[array([ 132.21182324,   23.80481826,    1.43482849]), None]

That looks like one row produced a (3,) array, the other had some sort of problem and produced None. I don't know if that None was due to removing the sort or not. But in any case, turning a mix of arrays and None back into an array would be a problem. If all items of res were matching arrays, we could stack them back into a 2d array.

There are ways of getting modest speed improvements (compared to this list comprehension). But with a complex function like this, the time spent in the function (called 2000 times) dominates the time spent by the iteration mechanism.

And since you are iterating on the 1st dimension, and passing the other 2 (as 3 arrays), this explicit loop is a lot easier to use than vectorize, frompyfunc or apply_along/over...

To get significant time savings you have to write f() to work with the 3d array directly.

edited Aug 24, 2017 at 16:52

answered Aug 24, 2017 at 16:46

hpaulj

233k14 gold badges260 silver badges392 bronze badges

3 Comments

patrik Over a year ago

Thanks @hpaulj. It can indeed happen that f produces a None. f computes the trilateration point of three points, and sometimes it just does not exist for a given value of r1. What is the proper way of dealing with such a case?

hpaulj Over a year ago

For a list it is fine to return a None. The question is - what do you want in a 2d array? Skip that case? Fill it with some sort of dummy values? This list comprehension approach gives you a lot of flexibility.

patrik Over a year ago

Thanks I like the list approach, because of its simplicity. Maybe I'll rewrite f as practice and then I think I will set np.nan for all three coordinates if no trilateration point exists.

Collectives™ on Stack Overflow

Apply 3-argument function to 3D numpy array

2 Answers 2

Comments

3 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

3 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related