Rotation of a vector distribution to align with a normal vector

I generate a distribution of random points on a unit hemisphere whose pole is on the positive z-axis (the base lies in the x-y plane). Each point represents a directional vector $v$ in which a ray will be fired.

I wish to rotate this distribution of points to align the 'hemisphere' with a vector which is the surface normal vector ($V_n$) of a given surface.

Currently I complete this by converting the cartesian coordinates of $N$ to spherical coordinates as follows:

$\theta = \arctan(N_y/N_x)$

$\phi = \arccos(N_z)$

I then form a rotation matrix for rotating around the z and x axis and use it to transform the original vector $v$ to get my desired vector $v'$. $v'$ will now lie on a hemisphere which is aligned with the surface normal $N$.

Is there a better way to complete this operation without having to calculate the spherical coordinates of the normal vector then compose a rotation matrix and preform a matrix multiplication. This function is used in a tight loop of some computationally expensive code so it would be ideal to preform some fundamental optimisation.

Perhaps some properties of vectors could be used it would be most beneficial if I could avoid having to calculate $\arccos$ and $\arctan$.

-
1. It might be better to use the two-argument arctangent for your application. 2. You do know that $\sin(\arctan(x,y))=\frac{y}{\sqrt{x^2+y^2}}$ and $\cos(\arctan(x,y))=\frac{x}{\sqrt{x^2+y^2}}$, don't you? –  Ｊ. Ｍ. Sep 3 '11 at 8:58
I do use the two argument arctan however I could use $sin(arctan(x,y))=\cfrac{y}{\sqrt{x^2+y^2}}$ to calculate the rotation matrix directly rather than calculating the angle from arctan then the matrix coefficient using sin. –  cubiclewar Sep 3 '11 at 10:04
Note that the convention of your computing environment might be backwards from the convention I use. The atan2() in your system might have the y first and the x last. –  Ｊ. Ｍ. Sep 3 '11 at 10:07

Yes, you can do this without any trigonometric functions, and even without any square roots. Take a look at Rodrigues' formula for a rotation matrix given an angle $\theta$ and an axis along a unit vector $k$:

$$R=I\cos\theta+\sin\theta[k]_\times+(1-\cos\theta)kk^{\text T}\;,$$

where $I$ is the identity matrix and

$$[k]_\times=\begin{pmatrix}0&-k_z&k_y\\k_z&0&-k_x\\-k_y&k_x&0\end{pmatrix}\;.$$

You have not an axis and an angle, but a vector $N$ that you want to rotate the $z$ axis into. The axis should be perpendicular to this vector and the $z$ axis, and the angle should be the angle between this vector and the $z$ axis. Since the $z$ component of your vector is the cosine of the angle it forms with the $z$ axis, you already have $\cos\theta=N_z$. Now consider

$$N_\perp=\begin{pmatrix}0\\0\\1\end{pmatrix}\times N=\begin{pmatrix}-N_y\\N_x\\0\end{pmatrix}\;.$$

This vector is perpendicular to both $e_z$ and $N$, so its direction is along the desired rotation axis. Its magnitude is $\sin\theta$. So

$$\begin{eqnarray} N_\perp&=&\sin\theta k\;,\\ [N_\perp]_\times&=&\sin\theta[k]_\times\;,\\ \end{eqnarray}$$

where $k$ is a unit vector along the rotation axis. Now we have most of the ingredients of Rodrigues' formula. The only problem left is that we don't have $kk^{\text T}$ but $\sin^2\theta kk^{\text T}$. But $\sin^2\theta=1-\cos^2\theta$, and we have $\cos\theta$, so we can correct for that using only elementary arithmetic. Putting it all together, we have

$$\begin{eqnarray} R &=& I\cos\theta+[N_\perp]_\times+(1-\cos\theta)N_\perp N_\perp^{\text T}/(1-\cos^2\theta) \\ &=& I\cos\theta+[N_\perp]_\times+N_\perp N_\perp^{\text T}/(1+\cos\theta) \\ &=& IN_z+ \begin{pmatrix}0&0&N_x\\0&0&N_y\\-N_x&-N_y&0\end{pmatrix} +\frac1{1+N_z} \begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix} \\ &=& \begin{pmatrix}N_z&0&N_x\\0&N_z&N_y\\-N_x&-N_y&N_z\end{pmatrix} +\frac1{1+N_z} \begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix} \;. \end{eqnarray}$$

You can check that this is an orthogonal matrix that rotates a unit vector along the $z$ axis into $N$. If you organize the operations efficiently, you only need two negations, three additions, one division and five multiplications to calculate the rotation matrix.

Note that the result becomes undefined at $N_z=-1$, corresponding to the fact that there's no preferred choice of axis in this case. You can either arbitrarily choose an axis specifically for that case, e.g. the $x$ axis, or, since this might also cause numerical problems near $N_z=-1$, you can instead solve the problem for $-N_z$ whenever $N_z<0$, and then invert the resulting vectors. This will produce a mirror image of your distribution, but I'd assume that this is rotationally invariant around the $z$ axis anyway, and in that case the inversion wouldn't make a difference.

-
Wow, thank you very much, from omitting the trigonometric functions I reduced my run time by around 40% . After a few days of searching for a method to simplify this I didn't really come up with anything. I good mathematician will save you every time! –  cubiclewar Sep 3 '11 at 10:28
+1. This is a great answer. –  Geoff Mar 18 '13 at 19:21