2

I have a basis for a plane in 3 dimensions: (u, v).

I would like to obtain all linear combinations of this basis to basically go through my whole plane:

for i in [0, 512[ and j in [0, 512[, get all (i * u + j * v).

I need this to be fast so for loops are not really an option. How can I do that with numpy broadcasting?

Had a look at https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html and my impression is that it's not possible to do...

Tried:

# This is an orthonormal basis but there is no guarantee it is
u = np.array([1, 0, 0])
v = np.array([0, 1, 0])
tmp = np.arange(512)
factors = itertools.combinations(tmp, 2)
pixels = factors[0] * u + factors[1] + v

But obviously it does not work.

Is there a solution to this problem? And if yes then how?

Improve this question

2 Answers 2

Reset to default
2

Multiply (u, v) with a 2D index grid:

ind = np.indices((512, 512))
pixels = ind[0, ..., np.newaxis] * u + ind[1, ..., np.newaxis] * v

>>> %timeit ind = np.indices((512, 512)); pixels = ind[0, ..., np.newaxis] * u + ind[1, ..., np.newaxis] * v
8.06 ms ± 69.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Multiply u with a 1D index range, multiply v with a 1D index range, broadcast and combine to 2D:

i512 = np.arange(512)[:, np.newaxis]
pixels = (i512 * u)[:, np.newaxis, :] + (i512 * v)[np.newaxis, :, :]

>>> %timeit i512 = np.arange(512)[:, np.newaxis]; pixels = (i512 * u)[:, np.newaxis, :] + (i512 * v)[np.newaxis, :, :]
4.06 ms ± 58.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Improve this answer
Sign up to request clarification or add additional context in comments.

5 Comments

I'm not sure I understand your example but it seems that this work in the particular case of my example. What if the basis is ([1, 1, 1], [1, -1, 1]) ?
The first example is a direct translation of the itertools combinations code you tried into only using numpy's np.indices. The second example uses the linear independence of the the first and second terms of i * u + j * v to speed up the computations. It should work for any values of u and v. Where exactly are you struggling to understand the solution?
Well the fact that you are using a "newaxis" on the x axis on u and "newaxis" on the y axis on v. I'm still relatively new to numpy so I need to review those concept.
docs.scipy.org/doc/numpy-1.13.0/reference/arrays.indexing.html is a good resource for that, worth studying as well as the broadcasting page you mentioned.
Your solution works. In order to have something completely answering the question, like @blue_note's answer does, you just need to add a pixels.reshape((512*512), 3) at the end. Your solutions is the fatest.
1

A few edits to your code

factors = itertools.product(tmp, tmp)
factors = list(zip(*factors))
factors = np.array(factors)
pixels = factors[0][:, None] * u + factors[1][:, None] + v

First, a math error: you want the Cartesian product, not the combinations.

Now, the actual syntax error: itertools produces a list of [(i1, j1), (i2, j2)..], not a numpy array. So, for your final line to work you have to

  1. zip(*) to convert the list to (i1, i2), (j1, j2) format
  2. make it a numpy array
  3. in the factors[0], factors[1] vectors, which are 1D, add [:, None] to convert the to columns, so that broadcasting works.
Improve this answer

2 Comments

Your solution works and is easy to understand. Unfortunately it is significantly slower than @w-m solution.
@LukeSkywalker: I know, I didn't write it as a proposed method, in terms of efficiency. Just wanted to explain why you get an error.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.