Rotate a vector by a randomly oriented angle

Question

We start off with a unit vector $\mathbf{v}$ randomly oriented in 3D space and we want to generate another unit vector $\mathbf{w}$ so that

$$ \mathbf{w}\cdot \mathbf{v} = \cos \beta $$

where $\beta$ is a given fixed angle. The new vector $\mathbf{w}$ must be randomly oriented, so that on average the cross product

$$ \langle \mathbf{w}\times \mathbf{v} \rangle = 0 $$

Normalization is $|\mathbf{v}|^2 = |\mathbf{w}|^2 = 1 $

How would you generate a particular vector $\mathbf{w}$ which satisfies the above conditions, using only 1 random number, and requiring as little computations as possible?

I have a "dirty" solution using polar coordinates, where I use an orthonormal set of vectors

$$ \mathbf{v} = \begin{pmatrix} \sin \theta \cos \phi\\ \sin \theta \sin \phi\\ \cos \theta \end{pmatrix} \quad \mathbf{v}^* = \begin{pmatrix} \cos \theta \cos \phi\\ \cos \theta \sin \phi\\ -\sin \theta \end{pmatrix} \quad \mathbf{v}^{\dagger} = \begin{pmatrix} -\sin \phi\\ \cos \phi\\ 0 \end{pmatrix} $$

or, in Cartesian components

$$ \mathbf{v} = \begin{pmatrix} x\\ y\\ z \end{pmatrix} \quad \mathbf{v}^* = \frac{1}{\sqrt{1-z^2}}\begin{pmatrix} xz\\ yz\\ z^2-1 \end{pmatrix} \quad \mathbf{v}^{\dagger} = \frac{1}{\sqrt{1-z^2}} \begin{pmatrix} -y\\ x\\ 0 \end{pmatrix} $$

The random orientation is then selected by a linear superposition \begin{equation} \mathbf{a} = \mathbf{v}^* \cos(2\pi \alpha) + \mathbf{v}^{\dagger} \sin(2\pi \alpha) \end{equation} where $\alpha \in (0,1)$ is a random number drawn from a uniform distribution. The rotation is given simply by

$$ \mathbf{w} = \mathbf{v}\cos \beta + \mathbf{a}\sin \beta $$

Mathematically this seems perfect, except in the situation where $z=1$ or close to it. Also, if I just keep rotating my vector, the normalization drifts away from 1 due to roundoff errors.

Do you have some ideas on a more general formula which would be numerically robust and maintain normalization more rigorously, preferably with not too many computations?

Answer 1

First, set up an orthonormal coordinate system $\{ \hat{\bf x}, \hat{\bf y}, {\bf v} \}$ with ${\bf v}$ as one axis. Pick any random vector ${\bf a}$ which is not collinear with ${\bf v}$. Then either choice of unit vector $$\hat{\bf x} := \frac{{\bf a} \times {\bf v}}{|{\bf a} \times {\bf v}|} \quad \text{or} \quad \hat{\bf x} := \frac{{\bf a} - ({\bf a} \cdot {\bf v}) {\bf v}}{|{\bf a} - ({\bf a} \cdot {\bf v}) {\bf v}|}$$ will be orthogonal to ${\bf v}$. $\hat{\bf y} := {\bf v} \times \hat{\bf x}$ completes the orthonormal coordinate frame. Now you can simply use spherical coordinates with ${\bf v}$ as the axis that is usually denoted as the $z$-axis: just let $${\bf w} = \sin \beta \cos \phi\, \hat{\bf x} + \sin \beta \sin \phi\, \hat{\bf y} + \cos \beta\, {\bf v},$$ where $\alpha$ is drawn uniformly from $[0,2 \pi)$.

Answer 2

If you want to use a single random number, you first need to find "a" solution (one vector w that meets the criterion), then rotate that vector about $\mathbf{v}$ by a randomly generated angle.

We have some flexibility when finding that first vector $w_0$ - we just need to find "a" vector at right angles to $v$ and rotate about that axis by $\beta$. We can choose any vector that is not $v$ and take the cross product with $v$ to get a vector that is at right angles. The general equation for rotating about a vector by an angle is given by the Rodriguez formula - if we are rotating about one of the unit axes, that simplifies considerably and allows you to get the initial vector.

Rotating that initial vector about $\mathbf{v}$ is going to take a bit of work. It might be convenient to use the formulation given here which uses quaternions to simplify the calculation of the rotation matrix. It means that you have to first express the rotation about the axis as a quaternion (where $\phi$ is a randomly generated angle between 0 and $2\pi$):

$$q = \begin{bmatrix}{ q_r\\ q_i\\ q_j\\ q_k} \end{bmatrix}= \begin{bmatrix}{ \cos\frac{\phi}{2}\\ v_1 \sin{\frac{\phi}{2}}\\ v_2 \sin{\frac{\phi}{2}}\\ v_3 \sin{\frac{\phi}{2}}} \end{bmatrix}$$

After which you define the rotation as

$$\mathbf{R} = \begin{bmatrix} 1-2s(q_j^2+q_k^2) & 2s(q_i q_j - q_k q_r) & 2s(q_i q_k + q_j q_r)\\ 2s(q_i q_j + q_k q_r) & 1 - 2s(q_i^2 + q_k^2) & 2s(q_j q_k - q_i q_r)\\ 2s(q_i q_k - q_j q_r) & 2s(q_j q_k + q_i q_r) & 1 - 2s(q_i^2 + q_j^2) \end{bmatrix}$$

where $s=1$ if you constructed your quaternion to have unit length (more generally, $s = ||\mathbf{q}||^{-2}$).

Computationally, a lot of terms appear in pairs so you only need to do the multiplication once.

It's possible there is a clever shortcut - but I don't know of one.

Update

I decided to see if I could implement the above in a piece of C code, and the result is here. No warranties that this is bug free, but some simple testing suggested it is "mostly correct". Might be a good starting point (note - in latest edit I made some changes to the calculation of w0 as the original was wrong...)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define TWOPI (4*acos(0))
#define NUM 1000000

float normalize(float *a) {
// normalize the vector given
  float temp = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
  for(int ii=0; ii<3; ii++) a[ii]/=temp;
  return temp;
}

void cross(float *a, float *b, float *c) {
// return cross product of a and b in c
  c[0] = a[1]*b[2] - b[1]*a[2];
  c[1] = a[2]*b[0] - a[0]*b[2];
  c[2] = a[0]*b[1] - a[1]*b[0];
  return;
}

float dot(float *a, float *b) {
  float sum=0;
  for(int ii=0; ii<3; ii++) sum+=a[ii]*b[ii];
  return sum;
}

void printmat(float r[][3]) {
  for(int ii=0; ii<3; ii++) {
    for(int jj=0; jj<3; jj++) {
      printf("%.4f\t", r[ii][jj]);
    }
    printf("\n");
  }
  printf("\n");
  return;
}

void printquat(float q[]) {
  for(int ii=0; ii<4; ii++) printf("%.4f\n", q[ii]);
  printf("\n");
  return;
}

void rotvec(float *axis, float *vec, float angle, float* result) {
// rotate vec about axis by angle
// return result
  float q[4]; // quaternion representation of vector
  float R[3][3]; // rotation matrix
  float s = sin(angle/2), c = cos(angle/2);
  // construct quaternion:
  q[0] = c;
  for(int ii=0; ii<3; ii++) q[ii+1]=axis[ii]*s;

  #ifdef VERBOSE
    printf("quaternion:\n");
    printquat(q);
  #endif
  // construct rotation matrix:
  float qii, qij, qik, qir, qjj, qjk, qjr, qkk, qkr;
  qii=q[1]*q[1]; qij=q[1]*q[2]; qik=q[1]*q[3]; qir=q[1]*q[0];
  qjj=q[2]*q[2]; qjk=q[2]*q[3]; qjr=q[2]*q[0];
  qkk=q[3]*q[3]; qkr=q[3]*q[0];
  R[0][0]=1-2*(qjj+qkk);
  R[0][1]=2*(qij-qkr);
  R[0][2]=2*(qik+qjr);
  R[1][0]=2*(qij+qkr);
  R[1][1]=1-2*(qii+qkk);
  R[1][2]=2*(qjk-qir);
  R[2][0]=2*(qij-qjr);
  R[2][1]=2*(qjk+qir);
  R[2][2]=1-2*(qii+qjj);

  #ifdef VERBOSE
    printf("rotation matrix: \n");
    printmat(R);
  #endif

  for(int ii=0; ii<3; ii++) {
    result[ii]=0;
    for(int jj=0; jj<3; jj++) {
      result[ii]+=R[ii][jj]*vec[jj];
    }
  }
  return;
}

void printvec(float *v) {
  for(int ii=0; ii<3; ii++) printf("%f\n", v[ii]);
  printf("\n");
  return;
}


void findw0(float *v, float b, float *w0) {
  // find any vector that has an angle b with vector v
  int mi=0;
  float minax=fabs(v[0]);
  for(int ii=1; ii<3; ii++) {
    if(fabs(v[ii])<minax) {
      mi=ii;
      minax = fabs(v[ii]);
    }
  }
  float X[]={0,0,0};
  X[mi]=1;
  float A[3]; // the axis we want to rotate about at right angles to v
  cross(v, X, A); // this finds "a" vector at right angles
  normalize(A); // make it unit length
  rotvec(A, v, b, w0); // rotate v about A by angle b
  return;
}

int main(void) {
// given an input vector v
// and an angle b
// find a random vector w that gives dot(v,w)=cos(b)

  float v[]={0,0,1};
  float b = acos(0)/2; // pi/4
  normalize(v);
  // find "a vector" w0 that has the right relationship:
  float w0[3];
  findw0(v, b, w0);
  printf("vector w0 that is pi/4 from Z axis:\n");
  printvec(w0);
  float w[3];
  float crossSum[]={0,0,0};
  float crossTemp[3];
  float bsum = 0;
  for(int ii=0; ii<NUM; ii++) {
    // generate random angle between -pi and pi:
    b = TWOPI*(rand()*1.0/RAND_MAX-0.5);
    bsum+=b;
    rotvec(v, w0, b, w);
    cross(v, w, crossTemp);
    for(int jj=0; jj<3; jj++) crossSum[jj]+=crossTemp[jj]/NUM;
  }
  printf("after %d iterations, mean cross product is\n", NUM);
  printvec(crossSum);
  printf("mean angle of rotation is %.3f\n", bsum/NUM);
  return 1;
}

Answer 3

Find a vector $\vec a$ such that $\vec v \times \vec a \ne 0$. Define $\hat x = (\vec v\,\times \vec a)/||\vec v \times \vec a||$. By construction, $\hat x$ is a unit vector and is orthogonal to $\vec v$. Define $\hat y = \vec v \times \vec x$. The triple $(\hat v, \hat x, \hat y)$ form a right handed coordinate system.

Form a second coordinate system $(\hat v, \hat x', \hat y')$ by rotating about $\hat v$ by an angle $\alpha$, where $\alpha$ is drawn from $U(0,2\pi)$. Finally, form a third coordinate system $(\hat w, \hat x', \hat y'')$ by rotating about $\hat x'$ by the desired angle $\beta$. By construction, $\hat w$ is a unit vector and $\hat w \cdot \hat v = \cos \beta$. Thanks to the original random rotation about $\hat v$, $\langle \hat w \times \hat v \rangle = 0$.

Answer 4

You have $$\mathbf{v}\cdot\mathbf{w} = v_1w_1 + v_2w_2 + v_3w_3 = \cos\beta$$ and $$ v_1^2 + v_2^2 + v_3^2 = 1\,.$$ Randomly generate one component of $\mathbf{v}$ and use the two equations above to solve for the other components. Does this work?