5
$\begingroup$

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?

$\endgroup$

migrated from physics.stackexchange.com Oct 2 '17 at 10:37

This question came from our site for active researchers, academics and students of physics.

3
$\begingroup$

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)$.

$\endgroup$
  • $\begingroup$ Very nice. When you have to generate lots of these vectors, this method greatly simplifies the calculation as you can precompute most of it. I may have a go comparing my original implementation against this one for speed comparison... $\endgroup$ – Floris Sep 30 '17 at 22:20
  • $\begingroup$ One improvement I would like to suggest: don't pick a "random" vector $\mathbf{a}$ - instead, choose one of [1 0 0], [0 1 0], [0 0 1] by looking at which of the coefficients of $\mathbf{v}$ is smallest: that gives you a "mostly orthogonal" first guess, and ensures that the cross product $\mathbf{a}\times\mathbf{v}$ is well conditioned. $\endgroup$ – Floris Oct 2 '17 at 13:14
2
$\begingroup$

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;
}
$\endgroup$
1
$\begingroup$

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$.

$\endgroup$
  • 2
    $\begingroup$ Your second paragraph is needlessly complicated. Once you've found $(\hat{x}, \hat{y}, \hat{v})$, just let ${\bf w} = \sin \beta \cos \alpha\, \hat{x} + \sin \beta \sin \alpha\, \hat{y} + \cos \beta\, \hat{v}$, where $\alpha$ is drawn uniformly from $[0, 2 \pi)$. $\endgroup$ – tparker Sep 30 '17 at 18:04
  • $\begingroup$ @tparker your solution is very elegant. Why don't you write it as an answer? $\endgroup$ – Floris Sep 30 '17 at 20:01
  • $\begingroup$ @Floris Okay, done. $\endgroup$ – tparker Sep 30 '17 at 22:07
0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ okay, but what is the distribution from which you would draw the new random component? It must be demonstrated that on average the cross product $\langle \mathbf{w} \times \mathbf{v} \rangle = 0$, i.e. the new vector is equally likely to fall anywhere on the cone at an angle $\beta$ to the old vector. $\endgroup$ – Airidas Korolkovas Sep 30 '17 at 16:10
  • $\begingroup$ Why not [0,1]? This is a bijection. It should give you what you want. $\endgroup$ – occamsrazor Sep 30 '17 at 16:20
  • 2
    $\begingroup$ The length of an individual component of $w$ (which is the randomly generated vector) does not follow a uniform distribution. This method won't work. $\endgroup$ – Floris Sep 30 '17 at 16:42
  • $\begingroup$ @Floris Yeah, you're right! This is so cool! Should it be sinusoidal? $\endgroup$ – occamsrazor Sep 30 '17 at 16:44
  • $\begingroup$ I think $\cos(rand(-\pi/2,\pi/2))$ might get you part of the way; but you still need to decide on the sign of the other components or you will end up in just some quadrant - and that means generating more than one random number. I suppose that generating a random number between $-\pi$ and $\pi$ might suffice - you could then take the sine and the cosine to give you the two pieces you need. Maybe. $\endgroup$ – Floris Sep 30 '17 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.