0

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
patrik
  • 13
  • 3
  • 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 Aug 24 '17 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 Aug 24 '17 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)`. – jdehesa Aug 24 '17 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 with`np.square(x)`. – patrik Aug 24 '17 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). – jdehesa Aug 24 '17 at 17:14
  • Thank you @jdehesa probably you meant `axis=1`. I had to use `np.newaxis` very often. Furthermore I ran into problems for the last step. It could happen that z_square is negative in that case intersection should be zero. In the end I want an array shape `(n, 3)` where some rows can be None, and not `[None, None, None]`. Is this possible? – patrik Aug 28 '17 at 14:10
  • @patrik Yes, sorry, you're right, it was `axis=1`. Umh, I don't think you can have rows with a single element, but you can put `[nan, nan, nan]` in those rows and then do something like `np.logical_and.reduce(~np.isnan(result), axis=1)` (hope I got it right this time) to get a mask for the "good" values only. – jdehesa Aug 29 '17 at 08:33
  • @jdehesa Thanks so much, you helped me a lot and I learned a lot. Before I post my final 'vectorized' version of `f` one more question. What is the proper pythonic way of applying `np.sqrt` to an array with possibly negative values for which I expect `np.nan`. Can I just apply it? In ipython I got a warning `RuntimeWarning: invalid value encountered in sqrt`. Can I just ignore it? Thanks again. – patrik Aug 30 '17 at 08:41
  • @patrik You can, it is a [documented feature](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.sqrt.html), you can [hide the warnings](https://stackoverflow.com/questions/9031783/hide-all-warnings-in-ipython) if you want. If you don't like it, though, an prefer to be more explicit about it, you can do something like `m = arr < 0; arr[m] *= 0; res = np.sqrt(arr); res[m] = np.nan`. – jdehesa Aug 30 '17 at 09:09
  • Thanks @jdehesa for this great `numpy` lecture. I learned a lot. I postet a solution based on your input. I guess it's not the most elegant, but maybe not that bad. ;-) Probably the `map` part is not that pythonic. Thanks again. – patrik Sep 04 '17 at 08:08

2 Answers2

1

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.

patrik
  • 13
  • 3
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.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • 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? – patrik Aug 28 '17 at 14:00
  • 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. – hpaulj Aug 29 '17 at 02:16
  • 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. – patrik Aug 30 '17 at 08:47