title: Compute rotation matrices from arrays of "aim" and "up" vectors
body: <p>I want to compute an array of rotation matrices from given arrays of <code>aim</code> and <code>up</code> vectors.</p>
  
  <p>For simplicity I'll assume the <code>aim</code> axes will correspond to the matrices's <code>x</code> components, and <code>up</code> axes to the matrices's <code>y</code> components.</p>
  
  <p>The only way that i know of is by doing a series of cross products:</p>
  
  <pre><code>import cProfile
  import numpy as np
  from numpy.core.umath_tests import inner1d
  
  normalize = lambda V: V/(inner1d(V,V)**0.5)[:,np.newaxis] # inner1d is faster than np.linalg.norm on large arrays
  
  def vectorToMatrix(X,Y):
      X = normalize(X) # make sure X is normalized
      Z = normalize(np.cross(X,Y)) # Z is the normal to the XY plane
  
      # Re-adjust Y to keep matrices orthogonal
      Y = np.cross(Z,X)
      return np.dstack((X,Y,Z)).swapaxes(2,1)
  </code></pre>
  
  <p>Running this on a million random items</p>
  
  <pre><code>np.random.seed(30)
  n = 10**6
  X = np.random.random((n,3))
  Y = np.random.random((n,3))
  M = vectorToMatrix(X,Y)
  
  cProfile.run('vectorToMatrix(X,Y)') # 61 function calls in 0.255 seconds
  </code></pre>
  
  <p>I'm looking for other methods, preferably leveraging <code>numpy/scipy</code>, that could help calculate the same result given by <code>vectorToMatrix</code> with a performance boost. </p>